]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split up evaluation_kernels.h into cell and face parts 15983/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Tue, 12 Sep 2023 12:26:06 +0000 (14:26 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Wed, 13 Sep 2023 16:35:15 +0000 (18:35 +0200)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_kernels_face.h [new file with mode: 0644]
include/deal.II/matrix_free/evaluation_template_face_factory.templates.h
include/deal.II/matrix_free/evaluation_template_factory_internal.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/fe_point_evaluation.h
include/deal.II/matrix_free/tensor_product_kernels.h

index e63dc157b9d61cc03295d0467754caedd3508006..afdf62854a2b16c3ca7881ea20a0fa09fd89d8d4 100644 (file)
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/base/geometry_info.h>
-#include <deal.II/base/ndarray.h>
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/vectorization.h>
 
-#include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/fe_evaluation_data.h>
 #include <deal.II/matrix_free/shape_info.h>
@@ -1825,164 +1822,6 @@ namespace internal
       Number                                *values_dofs_actual,
       FEEvaluationData<dim, Number, false>  &fe_eval,
       const bool                             add_into_values_array = false);
-
-    template <int  direction,
-              bool contract_over_rows,
-              bool symmetric_evaluate = true,
-              int  extra_stride       = 1>
-    static void
-    work_normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
-                const Number                                            *in,
-                Number                                                  *out,
-                const bool add_into_result  = false,
-                const int  subface_index_1d = 0)
-    {
-      AssertIndexRange(direction, dim);
-      AssertDimension(fe_degree, data.fe_degree);
-      AssertDimension(n_q_points_1d, data.n_q_points_1d);
-      constexpr int  n_rows    = fe_degree + 1;
-      constexpr int  n_columns = n_q_points_1d;
-      constexpr int  mm        = contract_over_rows ? n_rows : n_columns;
-      constexpr int  nn        = contract_over_rows ? n_columns : n_rows;
-      const Number2 *shape_data =
-        symmetric_evaluate ?
-          data.shape_values_eo.data() :
-          data.values_within_subface[subface_index_1d].data();
-      Assert(shape_data != nullptr, ExcNotInitialized());
-      Assert(contract_over_rows == false || !add_into_result,
-             ExcMessage("Cannot add into result if contract_over_rows = true"));
-
-      constexpr int n_blocks1  = Utilities::pow(fe_degree, direction);
-      constexpr int n_blocks2  = Utilities::pow(fe_degree, dim - direction - 1);
-      constexpr int stride_in  = contract_over_rows ? 1 : extra_stride;
-      constexpr int stride_out = contract_over_rows ? extra_stride : 1;
-      constexpr EvaluatorVariant variant =
-        symmetric_evaluate ? evaluate_evenodd : evaluate_general;
-
-      for (int i2 = 0; i2 < n_blocks2; ++i2)
-        {
-          for (int i1 = 0; i1 < n_blocks1; ++i1)
-            {
-              if (contract_over_rows == false && add_into_result)
-                apply_matrix_vector_product<variant,
-                                            EvaluatorQuantity::value,
-                                            n_rows,
-                                            n_columns,
-                                            n_blocks1 * stride_in,
-                                            n_blocks1 * stride_out,
-                                            contract_over_rows,
-                                            true>(shape_data, in, out);
-              else
-                apply_matrix_vector_product<variant,
-                                            EvaluatorQuantity::value,
-                                            n_rows,
-                                            n_columns,
-                                            n_blocks1 * stride_in,
-                                            n_blocks1 * stride_out,
-                                            contract_over_rows,
-                                            false>(shape_data, in, out);
-
-              in += stride_in;
-              out += stride_out;
-            }
-          in += n_blocks1 * (mm - 1) * stride_in;
-          out += n_blocks1 * (nn - 1) * stride_out;
-        }
-    }
-
-    template <int  direction,
-              int  normal_direction,
-              bool contract_over_rows,
-              bool symmetric_evaluate = true,
-              int  extra_stride       = 1>
-    static void
-    work_tangential(
-      const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
-      const Number                                            *in,
-      Number                                                  *out,
-      const int subface_index_1d = 0)
-    {
-      AssertIndexRange(direction, dim);
-      AssertDimension(fe_degree - 1, data.fe_degree);
-      AssertDimension(n_q_points_1d, data.n_q_points_1d);
-      static_assert(direction != normal_direction,
-                    "Cannot interpolate tangentially in normal direction");
-
-      constexpr int  n_rows    = fe_degree;
-      constexpr int  n_columns = n_q_points_1d;
-      const Number2 *shape_data =
-        symmetric_evaluate ?
-          data.shape_values_eo.data() :
-          data.values_within_subface[subface_index_1d].data();
-      Assert(shape_data != nullptr, ExcNotInitialized());
-
-      constexpr int n_blocks1 =
-        (direction > normal_direction) ?
-          Utilities::pow(n_q_points_1d, direction) :
-          (direction > 0 ?
-             (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
-             1);
-      constexpr int n_blocks2 =
-        (direction > normal_direction) ?
-          Utilities::pow(fe_degree, dim - 1 - direction) :
-          ((direction + 1 < dim) ?
-             (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
-             1);
-
-      constexpr EvaluatorVariant variant =
-        symmetric_evaluate ? evaluate_evenodd : evaluate_general;
-
-      // Since we may perform an in-place interpolation, we must run the step
-      // expanding the size of the basis backward ('contract_over_rows' aka
-      // 'evaluate' case), so shift the pointers and decrement during the loop
-      if (contract_over_rows)
-        {
-          in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
-          out += extra_stride *
-                 ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
-          for (int i2 = 0; i2 < n_blocks2; ++i2)
-            {
-              for (int i1 = 0; i1 < n_blocks1; ++i1)
-                {
-                  apply_matrix_vector_product<variant,
-                                              EvaluatorQuantity::value,
-                                              n_rows,
-                                              n_columns,
-                                              n_blocks1,
-                                              n_blocks1 * extra_stride,
-                                              true,
-                                              false>(shape_data, in, out);
-
-                  --in;
-                  out -= extra_stride;
-                }
-              in -= n_blocks1 * (n_rows - 1);
-              out -= n_blocks1 * (n_columns - 1) * extra_stride;
-            }
-        }
-      else
-        {
-          for (int i2 = 0; i2 < n_blocks2; ++i2)
-            {
-              for (int i1 = 0; i1 < n_blocks1; ++i1)
-                {
-                  apply_matrix_vector_product<variant,
-                                              EvaluatorQuantity::value,
-                                              n_rows,
-                                              n_columns,
-                                              n_blocks1 * extra_stride,
-                                              n_blocks1,
-                                              false,
-                                              false>(shape_data, in, out);
-
-                  in += extra_stride;
-                  ++out;
-                }
-              in += n_blocks1 * (n_columns - 1) * extra_stride;
-              out += n_blocks1 * (n_rows - 1);
-            }
-        }
-    }
   };
 
 
@@ -1999,7 +1838,7 @@ namespace internal
       const EvaluationFlags::EvaluationFlags evaluation_flag,
       Number                                *values_dofs,
       FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array)
+      const bool                             add)
   {
     Assert(dim == 2 || dim == 3,
            ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas "
@@ -2025,6 +1864,9 @@ namespace internal
 
     if (integrate)
       {
+        EvaluatorTensorProductAnisotropic<dim, fe_degree, n_q_points_1d, false>
+          eval;
+
         const bool do_values = evaluation_flag & EvaluationFlags::values;
         if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
           integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
@@ -2032,12 +1874,9 @@ namespace internal
                                                               gradients,
                                                               do_values);
         if constexpr (dim > 2)
-          work_tangential<2, 0, false>(shape_data[1], values, values);
-        work_tangential<1, 0, false>(shape_data[1], values, values);
-        work_normal<0, false>(shape_data[0],
-                              values,
-                              values_dofs,
-                              add_into_values_array);
+          eval.template tangential<2, 0>(shape_data[1], values, values);
+        eval.template tangential<1, 0>(shape_data[1], values, values);
+        eval.template normal<0>(shape_data[0], values, values_dofs, add);
 
         values += n_points;
         gradients += n_points * dim;
@@ -2049,12 +1888,9 @@ namespace internal
                                                               gradients,
                                                               do_values);
         if constexpr (dim > 2)
-          work_tangential<2, 1, false>(shape_data[1], values, values);
-        work_tangential<0, 1, false>(shape_data[1], values, values);
-        work_normal<1, false>(shape_data[0],
-                              values,
-                              values_dofs,
-                              add_into_values_array);
+          eval.template tangential<2, 1>(shape_data[1], values, values);
+        eval.template tangential<0, 1>(shape_data[1], values, values);
+        eval.template normal<1>(shape_data[0], values, values_dofs, add);
 
         if constexpr (dim > 2)
           {
@@ -2067,20 +1903,19 @@ namespace internal
                                                                   values,
                                                                   gradients,
                                                                   do_values);
-            work_tangential<1, 2, false>(shape_data[1], values, values);
-            work_tangential<0, 2, false>(shape_data[1], values, values);
-            work_normal<2, false>(shape_data[0],
-                                  values,
-                                  values_dofs,
-                                  add_into_values_array);
+            eval.template tangential<1, 2>(shape_data[1], values, values);
+            eval.template tangential<0, 2>(shape_data[1], values, values);
+            eval.template normal<0>(shape_data[0], values, values_dofs, add);
           }
       }
     else
       {
-        work_normal<0, true>(shape_data[0], values_dofs, values);
-        work_tangential<1, 0, true>(shape_data[1], values, values);
+        EvaluatorTensorProductAnisotropic<dim, fe_degree, n_q_points_1d, true>
+          eval;
+        eval.template normal<0>(shape_data[0], values_dofs, values);
+        eval.template tangential<1, 0>(shape_data[1], values, values);
         if constexpr (dim > 2)
-          work_tangential<2, 0, true>(shape_data[1], values, values);
+          eval.template tangential<2, 0>(shape_data[1], values, values);
         if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
           evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
                                                              values,
@@ -2090,10 +1925,10 @@ namespace internal
         gradients += n_points * dim;
         values_dofs += dofs_per_component;
 
-        work_normal<1, true>(shape_data[0], values_dofs, values);
-        work_tangential<0, 1, true>(shape_data[1], values, values);
+        eval.template normal<1>(shape_data[0], values_dofs, values);
+        eval.template tangential<0, 1>(shape_data[1], values, values);
         if constexpr (dim > 2)
-          work_tangential<2, 1, true>(shape_data[1], values, values);
+          eval.template tangential<2, 1>(shape_data[1], values, values);
         if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
           evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
                                                              values,
@@ -2105,9 +1940,9 @@ namespace internal
             gradients += n_points * dim;
             values_dofs += dofs_per_component;
 
-            work_normal<2, true>(shape_data[0], values_dofs, values);
-            work_tangential<0, 2, true>(shape_data[1], values, values);
-            work_tangential<1, 2, true>(shape_data[1], values, values);
+            eval.template normal<2>(shape_data[0], values_dofs, values);
+            eval.template tangential<0, 2>(shape_data[1], values, values);
+            eval.template tangential<1, 2>(shape_data[1], values, values);
             if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
               evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
                                                                  values,
@@ -2118,22 +1953,6 @@ namespace internal
 
 
 
-  /**
-   * Helper function to specify whether transformation to collocation should
-   * be used: It should give correct results (first condition), we need to be
-   * able to initialize the fields in shape_info.templates.h from the
-   * polynomials (second condition), and it should be the most efficient
-   * choice in terms of operation counts (third condition).
-   */
-  constexpr bool
-  use_collocation_evaluation(const unsigned int fe_degree,
-                             const unsigned int n_q_points_1d)
-  {
-    return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
-           (n_q_points_1d <= 3 * fe_degree / 2 + 1);
-  }
-
-
   /**
    * This class chooses an appropriate evaluation/integration strategy based on
    * the template parameters and the shape_info variable which contains runtime
@@ -2384,3055 +2203,187 @@ namespace internal
 
 
 
-  template <bool symmetric_evaluate,
-            int  dim,
-            int  fe_degree,
-            int  n_q_points_1d,
-            typename Number>
-  struct FEFaceEvaluationImpl
+  /**
+   * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation,
+   * using an FEEvaluationData argument.
+   */
+  template <int dim, typename Number>
+  struct CellwiseInverseMassMatrixImplBasic
   {
-    // We enable a transformation to collocation for derivatives if it gives
-    // correct results (first two conditions), if it is the most efficient
-    // choice in terms of operation counts (third condition) and if we were
-    // able to initialize the fields in shape_info.templates.h from the
-    // polynomials (fourth condition).
     using Number2 =
-      typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
-
-    using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
-                                                             evaluate_general,
-                                        dim - 1,
-                                        fe_degree + 1,
-                                        n_q_points_1d,
-                                        Number,
-                                        Number2>;
-
-    static Eval
-    create_evaluator_tensor_product(
-      const MatrixFreeFunctions::UnivariateShapeData<Number2> &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);
-        }
-    }
+      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
 
-    static void
-    evaluate_in_face(
-      const unsigned int                                       n_components,
-      const EvaluationFlags::EvaluationFlags                   evaluation_flag,
-      const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
-      Number                                                  *values_dofs,
-      Number                                                  *values_quad,
-      Number                                                  *gradients_quad,
-      Number                                                  *hessians_quad,
-      Number                                                  *scratch_data,
-      const unsigned int                                       subface_index)
+    template <int fe_degree, int = 0>
+    static bool
+    run(const unsigned int                          n_components,
+        const FEEvaluationData<dim, Number, false> &fe_eval,
+        const Number                               *in_array,
+        Number                                     *out_array)
     {
-      Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
-      Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
+      const unsigned int given_degree =
+        (fe_degree > -1) ? fe_degree :
+                           fe_eval.get_shape_info().data.front().fe_degree;
 
-      const std::size_t n_dofs = 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);
+      const unsigned int dofs_per_component =
+        Utilities::pow(given_degree + 1, dim);
 
-      // keep a copy of the original pointer for the case of the Hessians
-      Number *values_dofs_ptr = values_dofs;
+      Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
+      Assert(fe_eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
+             ExcNotImplemented());
 
-      if ((evaluation_flag & EvaluationFlags::values) != 0u &&
-          ((evaluation_flag & EvaluationFlags::gradients) == 0u))
-        for (unsigned int c = 0; c < n_components; ++c)
-          {
-            switch (dim)
-              {
-                case 3:
-                  eval0.template values<0, true, false>(values_dofs,
-                                                        values_quad);
-                  eval1.template values<1, true, false>(values_quad,
-                                                        values_quad);
-                  break;
-                case 2:
-                  eval0.template values<0, true, false>(values_dofs,
-                                                        values_quad);
-                  break;
-                case 1:
-                  values_quad[0] = values_dofs[0];
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
-            // Note: we always keep storage of values, 1st and 2nd derivatives
-            // in an array
-            values_dofs += 3 * n_dofs;
-            values_quad += n_q_points;
-          }
-      else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
-        for (unsigned int c = 0; c < n_components; ++c)
-          {
-            switch (dim)
-              {
-                case 3:
-                  if (symmetric_evaluate &&
-                      use_collocation_evaluation(fe_degree, n_q_points_1d))
-                    {
-                      eval0.template values<0, true, false>(values_dofs,
-                                                            values_quad);
-                      eval0.template values<1, true, false>(values_quad,
-                                                            values_quad);
-                      EvaluatorTensorProduct<evaluate_evenodd,
-                                             dim - 1,
-                                             n_q_points_1d,
-                                             n_q_points_1d,
-                                             Number,
-                                             Number2>
-                        eval_grad({}, data.shape_gradients_collocation_eo, {});
-                      eval_grad.template gradients<0, true, false, 3>(
-                        values_quad, gradients_quad);
-                      eval_grad.template gradients<1, true, false, 3>(
-                        values_quad, gradients_quad + 1);
-                    }
-                  else
-                    {
-                      // grad x
-                      eval0.template gradients<0, true, false>(values_dofs,
-                                                               scratch_data);
-                      eval1.template values<1, true, false, 3>(scratch_data,
-                                                               gradients_quad);
-
-                      // grad y
-                      eval0.template values<0, true, false>(values_dofs,
-                                                            scratch_data);
-                      eval1.template gradients<1, true, false, 3>(
-                        scratch_data, gradients_quad + 1);
-
-                      if ((evaluation_flag & EvaluationFlags::values) != 0u)
-                        eval1.template values<1, true, false>(scratch_data,
-                                                              values_quad);
-                    }
-                  // grad z
-                  eval0.template values<0, true, false>(values_dofs + n_dofs,
-                                                        scratch_data);
-                  eval1.template values<1, true, false, 3>(scratch_data,
-                                                           gradients_quad + 2);
-
-                  break;
-                case 2:
-                  eval0.template values<0, true, false, 2>(values_dofs + n_dofs,
-                                                           gradients_quad + 1);
-                  eval0.template gradients<0, true, false, 2>(values_dofs,
-                                                              gradients_quad);
-                  if ((evaluation_flag & EvaluationFlags::values) != 0u)
-                    eval0.template values<0, true, false>(values_dofs,
-                                                          values_quad);
-                  break;
-                case 1:
-                  values_quad[0]    = values_dofs[0];
-                  gradients_quad[0] = values_dofs[1];
-                  break;
-                default:
-                  AssertThrow(false, ExcNotImplemented());
-              }
-            values_dofs += 3 * n_dofs;
-            values_quad += n_q_points;
-            gradients_quad += dim * n_q_points;
-          }
+      EvaluatorTensorProduct<evaluate_evenodd,
+                             dim,
+                             fe_degree + 1,
+                             fe_degree + 1,
+                             Number,
+                             Number2>
+        evaluator({},
+                  {},
+                  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
+                  given_degree + 1,
+                  given_degree + 1);
 
-      if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
+      for (unsigned int d = 0; d < n_components; ++d)
         {
-          values_dofs = values_dofs_ptr;
-          for (unsigned int c = 0; c < n_components; ++c)
-            {
-              switch (dim)
-                {
-                  case 3:
-                    // grad xx
-                    eval0.template hessians<0, true, false>(values_dofs,
-                                                            scratch_data);
-                    eval1.template values<1, true, false>(scratch_data,
-                                                          hessians_quad);
-
-                    // grad yy
-                    eval0.template values<0, true, false>(values_dofs,
-                                                          scratch_data);
-                    eval1.template hessians<1, true, false>(scratch_data,
-                                                            hessians_quad +
-                                                              n_q_points);
-
-                    // grad zz
-                    eval0.template values<0, true, false>(values_dofs +
-                                                            2 * n_dofs,
-                                                          scratch_data);
-                    eval1.template values<1, true, false>(scratch_data,
-                                                          hessians_quad +
-                                                            2 * n_q_points);
-
-                    // grad xy
-                    eval0.template gradients<0, true, false>(values_dofs,
-                                                             scratch_data);
-                    eval1.template gradients<1, true, false>(scratch_data,
-                                                             hessians_quad +
-                                                               3 * n_q_points);
-
-                    // grad xz
-                    eval0.template gradients<0, true, false>(values_dofs +
-                                                               n_dofs,
-                                                             scratch_data);
-                    eval1.template values<1, true, false>(scratch_data,
-                                                          hessians_quad +
-                                                            4 * n_q_points);
-
-                    // grad yz
-                    eval0.template values<0, true, false>(values_dofs + n_dofs,
-                                                          scratch_data);
-                    eval1.template gradients<1, true, false>(scratch_data,
-                                                             hessians_quad +
-                                                               5 * n_q_points);
-
-                    break;
-                  case 2:
-                    // grad xx
-                    eval0.template hessians<0, true, false>(values_dofs,
-                                                            hessians_quad);
-                    // grad yy
-                    eval0.template values<0, true, false>(
-                      values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
-                    // grad xy
-                    eval0.template gradients<0, true, false>(
-                      values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
-                    break;
-                  case 1:
-                    hessians_quad[0] = values_dofs[2];
-                    break;
-                  default:
-                    AssertThrow(false, ExcNotImplemented());
-                }
-              values_dofs += 3 * n_dofs;
-              hessians_quad += dim * (dim + 1) / 2 * n_q_points;
-            }
+          const Number *in  = in_array + d * dofs_per_component;
+          Number       *out = out_array + d * dofs_per_component;
+          // Need to select 'apply' method with hessian slot because values
+          // assume symmetries that do not exist in the inverse shapes
+          evaluator.template hessians<0, true, false>(in, out);
+          if (dim > 1)
+            evaluator.template hessians<1, true, false>(out, out);
+          if (dim > 2)
+            evaluator.template hessians<2, true, false>(out, out);
         }
-    }
-
-    static void
-    integrate_in_face(
-      const unsigned int                                       n_components,
-      const EvaluationFlags::EvaluationFlags                   integration_flag,
-      const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
-      Number                                                  *values_dofs,
-      Number                                                  *values_quad,
-      Number                                                  *gradients_quad,
-      Number                                                  *hessians_quad,
-      Number                                                  *scratch_data,
-      const unsigned int                                       subface_index)
-    {
-      Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
-      Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
-
-      const std::size_t n_dofs =
-        fe_degree > -1 ?
-          Utilities::pow(fe_degree + 1, dim - 1) :
-          (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 ((integration_flag & EvaluationFlags::values) != 0u &&
-          (integration_flag & EvaluationFlags::gradients) == 0u)
-        for (unsigned int c = 0; c < n_components; ++c)
-          {
-            switch (dim)
-              {
-                case 3:
-                  eval1.template values<1, false, false>(values_quad,
-                                                         values_quad);
-                  eval0.template values<0, false, false>(values_quad,
-                                                         values_dofs);
-                  break;
-                case 2:
-                  eval0.template values<0, false, false>(values_quad,
-                                                         values_dofs);
-                  break;
-                case 1:
-                  values_dofs[0] = values_quad[0];
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
-            values_dofs += 3 * n_dofs;
-            values_quad += n_q_points;
-          }
-      else if ((integration_flag & EvaluationFlags::gradients) != 0u)
-        for (unsigned int c = 0; c < n_components; ++c)
-          {
-            switch (dim)
-              {
-                case 3:
-                  // grad z
-                  eval1.template values<1, false, false, 3>(gradients_quad + 2,
-                                                            scratch_data);
-                  eval0.template values<0, false, false>(scratch_data,
-                                                         values_dofs + n_dofs);
-                  if (symmetric_evaluate &&
-                      use_collocation_evaluation(fe_degree, n_q_points_1d))
-                    {
-                      EvaluatorTensorProduct<evaluate_evenodd,
-                                             dim - 1,
-                                             n_q_points_1d,
-                                             n_q_points_1d,
-                                             Number,
-                                             Number2>
-                        eval_grad({}, data.shape_gradients_collocation_eo, {});
-                      if ((integration_flag & EvaluationFlags::values) != 0u)
-                        eval_grad.template gradients<1, false, true, 3>(
-                          gradients_quad + 1, values_quad);
-                      else
-                        eval_grad.template gradients<1, false, false, 3>(
-                          gradients_quad + 1, values_quad);
-                      eval_grad.template gradients<0, false, true, 3>(
-                        gradients_quad, values_quad);
-                      eval0.template values<1, false, false>(values_quad,
-                                                             values_quad);
-                      eval0.template values<0, false, false>(values_quad,
-                                                             values_dofs);
-                    }
-                  else
-                    {
-                      if ((integration_flag & EvaluationFlags::values) != 0u)
-                        {
-                          eval1.template values<1, false, false>(values_quad,
-                                                                 scratch_data);
-                          eval1.template gradients<1, false, true, 3>(
-                            gradients_quad + 1, scratch_data);
-                        }
-                      else
-                        eval1.template gradients<1, false, false, 3>(
-                          gradients_quad + 1, scratch_data);
-
-                      // grad y
-                      eval0.template values<0, false, false>(scratch_data,
-                                                             values_dofs);
-
-                      // grad x
-                      eval1.template values<1, false, false, 3>(gradients_quad,
-                                                                scratch_data);
-                      eval0.template gradients<0, false, true>(scratch_data,
-                                                               values_dofs);
-                    }
-                  break;
-                case 2:
-                  eval0.template values<0, false, false, 2>(gradients_quad + 1,
-                                                            values_dofs +
-                                                              n_dofs);
-                  eval0.template gradients<0, false, false, 2>(gradients_quad,
-                                                               values_dofs);
-                  if ((integration_flag & EvaluationFlags::values) != 0u)
-                    eval0.template values<0, false, true>(values_quad,
-                                                          values_dofs);
-                  break;
-                case 1:
-                  values_dofs[0] = values_quad[0];
-                  values_dofs[1] = gradients_quad[0];
-                  break;
-                default:
-                  AssertThrow(false, ExcNotImplemented());
-              }
-            values_dofs += 3 * n_dofs;
-            values_quad += n_q_points;
-            gradients_quad += dim * n_q_points;
-          }
-
-      if ((integration_flag & EvaluationFlags::hessians) != 0u)
+      for (unsigned int q = 0; q < dofs_per_component; ++q)
         {
-          values_dofs = values_dofs_ptr;
-          for (unsigned int c = 0; c < n_components; ++c)
-            {
-              switch (dim)
-                {
-                  case 3:
-                    // grad xx
-                    eval1.template values<1, false, false>(hessians_quad,
-                                                           scratch_data);
-                    if ((integration_flag & (EvaluationFlags::values |
-                                             EvaluationFlags::gradients)) != 0u)
-                      eval0.template hessians<0, false, true>(scratch_data,
-                                                              values_dofs);
-                    else
-                      eval0.template hessians<0, false, false>(scratch_data,
-                                                               values_dofs);
-
-                    // grad yy
-                    eval1.template hessians<1, false, false>(hessians_quad +
-                                                               n_q_points,
-                                                             scratch_data);
-                    eval0.template values<0, false, true>(scratch_data,
-                                                          values_dofs);
-
-                    // grad zz
-                    eval1.template values<1, false, false>(hessians_quad +
-                                                             2 * n_q_points,
-                                                           scratch_data);
-                    eval0.template values<0, false, false>(scratch_data,
-                                                           values_dofs +
-                                                             2 * n_dofs);
-
-                    // grad xy
-                    eval1.template gradients<1, false, false>(hessians_quad +
-                                                                3 * n_q_points,
-                                                              scratch_data);
-                    eval0.template gradients<0, false, true>(scratch_data,
-                                                             values_dofs);
-
-                    // grad xz
-                    eval1.template values<1, false, false>(hessians_quad +
-                                                             4 * n_q_points,
-                                                           scratch_data);
-                    if ((integration_flag & EvaluationFlags::gradients) != 0u)
-                      eval0.template gradients<0, false, true>(scratch_data,
-                                                               values_dofs +
-                                                                 n_dofs);
-                    else
-                      eval0.template gradients<0, false, false>(scratch_data,
-                                                                values_dofs +
-                                                                  n_dofs);
-
-                    // grad yz
-                    eval1.template gradients<1, false, false>(hessians_quad +
-                                                                5 * n_q_points,
-                                                              scratch_data);
-                    eval0.template values<0, false, true>(scratch_data,
-                                                          values_dofs + n_dofs);
-
-                    break;
-                  case 2:
-                    // grad xx
-                    if ((integration_flag & (EvaluationFlags::values |
-                                             EvaluationFlags::gradients)) != 0u)
-                      eval0.template hessians<0, false, true>(hessians_quad,
-                                                              values_dofs);
-                    else
-                      eval0.template hessians<0, false, false>(hessians_quad,
-                                                               values_dofs);
-
-                    // grad yy
-                    eval0.template values<0, false, false>(
-                      hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
-                    // grad xy
-                    if ((integration_flag & EvaluationFlags::gradients) != 0u)
-                      eval0.template gradients<0, false, true>(
-                        hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
-                    else
-                      eval0.template gradients<0, false, false>(
-                        hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
-                    break;
-                  case 1:
-                    values_dofs[2] = hessians_quad[0];
-                    if ((integration_flag & EvaluationFlags::values) == 0u)
-                      values_dofs[0] = 0;
-                    if ((integration_flag & EvaluationFlags::gradients) == 0u)
-                      values_dofs[1] = 0;
-                    break;
-                  default:
-                    AssertThrow(false, ExcNotImplemented());
-                }
-              values_dofs += 3 * n_dofs;
-              hessians_quad += dim * (dim + 1) / 2 * n_q_points;
-            }
+          const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
+          for (unsigned int d = 0; d < n_components; ++d)
+            out_array[q + d * dofs_per_component] *= inverse_JxW_q;
         }
+      for (unsigned int d = 0; d < n_components; ++d)
+        {
+          Number *out = out_array + d * dofs_per_component;
+          if (dim > 2)
+            evaluator.template hessians<2, false, false>(out, out);
+          if (dim > 1)
+            evaluator.template hessians<1, false, false>(out, out);
+          evaluator.template hessians<0, false, false>(out, out);
+        }
+      return false;
     }
   };
 
 
 
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  struct FEFaceEvaluationImplRaviartThomas
+  /**
+   * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation
+   * with user-provided coefficients at quadrature points (in contrast to
+   * CellwiseInverseMassMatrixImplBasic, which implicitly uses `1/(|J|xW)' as
+   * coefficient).
+   */
+  template <int dim, typename Number>
+  struct CellwiseInverseMassMatrixImplFlexible
   {
     using Number2 =
-      typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
+      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
 
-    /**
-     * Apply the sum factorization kernels within the face for Raviart-Thomas
-     * elements for either evaluation or integration
-     */
-    template <bool do_integrate>
-    static inline void
-    evaluate_or_integrate_in_face(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      const std::vector<MatrixFreeFunctions::UnivariateShapeData<Number2>>
-                        &shape_data,
-      Number            *values_dofs_in,
-      Number            *values_quad,
-      Number            *gradients_quad,
-      Number            *scratch_data,
-      const unsigned int subface_index,
-      const unsigned int face_direction)
+    template <int fe_degree, int = 0>
+    static bool
+    run(const unsigned int                          n_desired_components,
+        const FEEvaluationData<dim, Number, false> &fe_eval,
+        const ArrayView<const Number>              &inverse_coefficients,
+        const bool                                  dyadic_coefficients,
+        const Number                               *in_array,
+        Number                                     *out_array)
     {
-      AssertDimension(shape_data.size(), 2);
-
-      const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
-      const int n_rows_n = degree + 1;
-      const int n_rows_t = degree;
-      const dealii::ndarray<int, 3, 3> dofs_per_direction{
-        {{{n_rows_n, n_rows_t, n_rows_t}},
-         {{n_rows_t, n_rows_n, n_rows_t}},
-         {{n_rows_t, n_rows_t, n_rows_n}}}};
-
-      (void)scratch_data;
-      (void)subface_index;
-      // TODO: This is currently not implemented, but the test
-      // matrix_vector_rt_face_03 apparently works without it -> check
-      // if (subface_index < GeometryInfo<dim - 1>::max_children_per_cell)
-      //  Assert(false, ExcNotImplemented());
-
-      using EvalAniso =
-        FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
-                         dim - 1,
-                         fe_degree,
-                         n_q_points_1d,
-                         Number>;
-      using Eval = EvaluatorTensorProduct<evaluate_evenodd,
-                                          dim - 1,
-                                          fe_degree,
-                                          n_q_points_1d,
-                                          Number,
-                                          Number2>;
-
-      std::array<int, dim> values_dofs_offsets = {};
-      for (unsigned int comp = 0; comp < dim - 1; ++comp)
+      const unsigned int given_degree =
+        (fe_degree > -1) ? fe_degree :
+                           fe_eval.get_shape_info().data.front().fe_degree;
+
+      const unsigned int dofs_per_component =
+        Utilities::pow(given_degree + 1, dim);
+
+      Assert(inverse_coefficients.size() > 0 &&
+               inverse_coefficients.size() % dofs_per_component == 0,
+             ExcMessage(
+               "Expected diagonal to be a multiple of scalar dof per cells"));
+
+      if (!dyadic_coefficients)
         {
-          if (dim == 2)
-            values_dofs_offsets[comp + 1] =
-              values_dofs_offsets[comp] +
-              3 * dofs_per_direction[comp][(face_direction + 1) % dim];
-          else
-            values_dofs_offsets[comp + 1] =
-              values_dofs_offsets[comp] +
-              3 * dofs_per_direction[comp][(face_direction + 1) % dim] *
-                dofs_per_direction[comp][(face_direction + 2) % dim];
+          if (inverse_coefficients.size() != dofs_per_component)
+            AssertDimension(n_desired_components * dofs_per_component,
+                            inverse_coefficients.size());
+        }
+      else
+        {
+          AssertDimension(n_desired_components * n_desired_components *
+                            dofs_per_component,
+                          inverse_coefficients.size());
         }
 
-      // Jacobians on faces are reordered to enable simple access with the
-      // regular evaluators; to get the RT Piola transform right, we need to
-      // pass through the values_dofs array in a permuted right order
-      std::array<unsigned int, dim> components;
-      for (unsigned int comp = 0; comp < dim; ++comp)
-        components[comp] = (face_direction + comp + 1) % dim;
+      Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
+      Assert(fe_eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
+             ExcNotImplemented());
 
-      for (const unsigned int comp : components)
-        {
-          Number *values_dofs = values_dofs_in + values_dofs_offsets[comp];
+      EvaluatorTensorProduct<evaluate_evenodd,
+                             dim,
+                             fe_degree + 1,
+                             fe_degree + 1,
+                             Number,
+                             Number2>
+        evaluator({},
+                  {},
+                  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
+                  given_degree + 1,
+                  given_degree + 1);
+
+      const Number *in  = in_array;
+      Number       *out = out_array;
+
+      const Number *inv_coefficient = inverse_coefficients.data();
+
+      const unsigned int shift_coefficient =
+        inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
+                                                           0;
 
-          std::array<int, 2> n_blocks{
-            {dofs_per_direction[comp][(face_direction + 1) % dim],
-             (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] :
-                        1)}};
+      const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
+      const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
 
-          if constexpr (dim == 3)
+      for (unsigned int d = 0; d < n_comp_outer; ++d)
+        {
+          for (unsigned int di = 0; di < n_comp_inner; ++di)
+            {
+              const Number *in_  = in + di * dofs_per_component;
+              Number       *out_ = out + di * dofs_per_component;
+              evaluator.template hessians<0, true, false>(in_, out_);
+              if (dim > 1)
+                evaluator.template hessians<1, true, false>(out_, out_);
+              if (dim > 2)
+                evaluator.template hessians<2, true, false>(out_, out_);
+            }
+          if (dyadic_coefficients)
             {
-              EvaluatorTensorProduct<evaluate_evenodd,
-                                     dim - 1,
-                                     n_q_points_1d,
-                                     n_q_points_1d,
-                                     Number,
-                                     Number2>
-                eval_g({},
-                       shape_data[0].shape_gradients_collocation_eo.data(),
-                       {});
-              if (!do_integrate)
+              const auto n_coeff_components =
+                n_desired_components * n_desired_components;
+              if (n_desired_components == dim)
                 {
-                  // Evaluate in 3d
-                  if (n_blocks[0] == n_rows_n)
-                    {
-                      EvalAniso::template work_normal<0, true, true, 1>(
-                        shape_data[0], values_dofs, values_quad);
-                      EvalAniso::template work_tangential<1, 0, true, true, 1>(
-                        shape_data[1], values_quad, values_quad);
-
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          EvalAniso::template work_normal<0, true, true, 1>(
-                            shape_data[0],
-                            values_dofs + n_blocks[0] * n_blocks[1],
-                            scratch_data);
-                          EvalAniso::
-                            template work_tangential<1, 0, true, true, dim>(
-                              shape_data[1], scratch_data, gradients_quad + 2);
-                        }
-                    }
-                  else if (n_blocks[1] == n_rows_n)
-                    {
-                      EvalAniso::template work_normal<1, true, true, 1>(
-                        shape_data[0], values_dofs, values_quad);
-                      EvalAniso::template work_tangential<0, 1, true, true, 1>(
-                        shape_data[1], values_quad, values_quad);
-
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          EvalAniso::template work_normal<1, true, true, 1>(
-                            shape_data[0],
-                            values_dofs + n_blocks[0] * n_blocks[1],
-                            scratch_data);
-                          EvalAniso::
-                            template work_tangential<0, 1, true, true, dim>(
-                              shape_data[1], scratch_data, gradients_quad + 2);
-                        }
-                    }
-                  else
-                    {
-                      Eval eval(shape_data[1].shape_values_eo.data(), {}, {});
-                      eval.template values<0, true, false>(values_dofs,
-                                                           values_quad);
-                      eval.template values<1, true, false>(values_quad,
-                                                           values_quad);
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          eval.template values<0, true, false>(values_dofs +
-                                                                 n_blocks[0] *
-                                                                   n_blocks[1],
-                                                               scratch_data);
-                          eval.template values<1, true, false, dim>(
-                            scratch_data, gradients_quad + 2);
-                        }
-                    }
-                  if (evaluation_flag & EvaluationFlags::gradients)
-                    {
-                      eval_g.template gradients<0, true, false, dim>(
-                        values_quad, gradients_quad);
-                      eval_g.template gradients<1, true, false, dim>(
-                        values_quad, gradients_quad + 1);
-                    }
+                  for (unsigned int q = 0; q < dofs_per_component; ++q)
+                    vmult<dim>(&inv_coefficient[q * n_coeff_components],
+                               &in[q],
+                               &out[q],
+                               dofs_per_component);
                 }
               else
                 {
-                  // Integrate in 3d
-                  if (evaluation_flag & EvaluationFlags::gradients)
-                    {
-                      if (evaluation_flag & EvaluationFlags::values)
-                        eval_g.template gradients<0, false, true, dim>(
-                          gradients_quad, values_quad);
-                      else
-                        eval_g.template gradients<0, false, false, dim>(
-                          gradients_quad, values_quad);
-                      eval_g.template gradients<1, false, true, dim>(
-                        gradients_quad + 1, values_quad);
-                    }
-                  if (n_blocks[0] == n_rows_n)
-                    {
-                      EvalAniso::template work_tangential<1, 0, false, true, 1>(
-                        shape_data[1], values_quad, values_quad);
-                      EvalAniso::template work_normal<0, false, true, 1>(
-                        shape_data[0], values_quad, values_dofs);
-
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          EvalAniso::
-                            template work_tangential<1, 0, false, true, dim>(
-                              shape_data[1], gradients_quad + 2, scratch_data);
-                          EvalAniso::template work_normal<0, false, true, 1>(
-                            shape_data[0],
-                            scratch_data,
-                            values_dofs + n_blocks[0] * n_blocks[1]);
-                        }
-                    }
-                  else if (n_blocks[1] == n_rows_n)
-                    {
-                      EvalAniso::template work_tangential<0, 1, false, true, 1>(
-                        shape_data[1], values_quad, values_quad);
-                      EvalAniso::template work_normal<1, false, true, 1>(
-                        shape_data[0], values_quad, values_dofs);
-
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          EvalAniso::
-                            template work_tangential<0, 1, false, true, dim>(
-                              shape_data[1], gradients_quad + 2, scratch_data);
-                          EvalAniso::template work_normal<1, false, true, 1>(
-                            shape_data[0],
-                            scratch_data,
-                            values_dofs + n_blocks[0] * n_blocks[1]);
-                        }
-                    }
-                  else
-                    {
-                      Eval eval(shape_data[1].shape_values_eo.data(), {}, {});
-                      eval.template values<1, false, false>(values_quad,
-                                                            values_quad);
-                      eval.template values<0, false, false>(values_quad,
-                                                            values_dofs);
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          eval.template values<1, false, false, dim>(
-                            gradients_quad + 2, scratch_data);
-                          eval.template values<0, false, false>(
-                            scratch_data,
-                            values_dofs + n_blocks[0] * n_blocks[1]);
-                        }
-                    }
-                }
-            }
-          else
-            {
-              using EvalN = EvaluatorTensorProduct<evaluate_evenodd,
-                                                   dim - 1,
-                                                   fe_degree + 1,
-                                                   n_q_points_1d,
-                                                   Number,
-                                                   Number2>;
-              if (!do_integrate)
-                {
-                  // Evaluate in 2d
-                  if (n_blocks[0] == n_rows_n)
-                    {
-                      EvalN eval(shape_data[0].shape_values_eo,
-                                 shape_data[0].shape_gradients_eo,
-                                 {});
-                      eval.template values<0, true, false>(values_dofs,
-                                                           values_quad);
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          eval.template gradients<0, true, false, dim>(
-                            values_dofs, gradients_quad);
-                          eval.template values<0, true, false, dim>(
-                            values_dofs + n_rows_n, gradients_quad + 1);
-                        }
-                    }
-                  else
-                    {
-                      Eval eval(shape_data[1].shape_values_eo,
-                                shape_data[1].shape_gradients_eo,
-                                {});
-                      eval.template values<0, true, false>(values_dofs,
-                                                           values_quad);
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          eval.template gradients<0, true, false, dim>(
-                            values_dofs, gradients_quad);
-                          eval.template values<0, true, false, dim>(
-                            values_dofs + n_rows_t, gradients_quad + 1);
-                        }
-                    }
-                }
-              else
-                {
-                  // Integrate in 2d
-                  if (n_blocks[0] == n_rows_n)
-                    {
-                      EvalN eval(shape_data[0].shape_values_eo,
-                                 shape_data[0].shape_gradients_eo,
-                                 {});
-                      if (evaluation_flag & EvaluationFlags::values)
-                        eval.template values<0, false, false>(values_quad,
-                                                              values_dofs);
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          if (evaluation_flag & EvaluationFlags::values)
-                            eval.template gradients<0, false, true, dim>(
-                              gradients_quad, values_dofs);
-                          else
-                            eval.template gradients<0, false, false, dim>(
-                              gradients_quad, values_dofs);
-                          eval.template values<0, false, false, dim>(
-                            gradients_quad + 1, values_dofs + n_rows_n);
-                        }
-                    }
-                  else
-                    {
-                      Eval eval(shape_data[1].shape_values_eo,
-                                shape_data[1].shape_gradients_eo,
-                                {});
-                      if (evaluation_flag & EvaluationFlags::values)
-                        eval.template values<0, false, false>(values_quad,
-                                                              values_dofs);
-                      if (evaluation_flag & EvaluationFlags::gradients)
-                        {
-                          if (evaluation_flag & EvaluationFlags::values)
-                            eval.template gradients<0, false, true, dim>(
-                              gradients_quad, values_dofs);
-                          else
-                            eval.template gradients<0, false, false, dim>(
-                              gradients_quad, values_dofs);
-                          eval.template values<0, false, false, dim>(
-                            gradients_quad + 1, values_dofs + n_rows_t);
-                        }
-                    }
-                }
-            }
-          values_quad += Utilities::pow(n_q_points_1d, dim - 1);
-          gradients_quad += dim * Utilities::pow(n_q_points_1d, dim - 1);
-        }
-    }
-  };
-
-
-
-  template <int dim, int fe_degree, typename Number>
-  struct FEFaceNormalEvaluationImpl
-  {
-    using Number2 =
-      typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
-
-    template <bool do_evaluate, bool add_into_output>
-    static void
-    interpolate(const unsigned int                             n_components,
-                const EvaluationFlags::EvaluationFlags         flags,
-                const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
-                const Number                                  *input,
-                Number                                        *output,
-                const unsigned int                             face_no)
-    {
-      Assert(static_cast<unsigned int>(fe_degree) ==
-                 shape_info.data.front().fe_degree ||
-               fe_degree == -1,
-             ExcInternalError());
-      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
-        interpolate_raviart_thomas<do_evaluate, add_into_output>(
-          n_components, input, output, flags, face_no, shape_info);
-      else
-        interpolate_generic<do_evaluate, add_into_output>(
-          n_components,
-          input,
-          output,
-          flags,
-          face_no,
-          shape_info.data.front().fe_degree + 1,
-          shape_info.data.front().shape_data_on_face,
-          shape_info.dofs_per_component_on_cell,
-          3 * shape_info.dofs_per_component_on_face);
-    }
-
-    /**
-     * Interpolate the values on the cell quadrature points onto a face.
-     */
-    template <bool do_evaluate, bool add_into_output>
-    static void
-    interpolate_quadrature(
-      const unsigned int                             n_components,
-      const EvaluationFlags::EvaluationFlags         flags,
-      const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
-      const Number                                  *input,
-      Number                                        *output,
-      const unsigned int                             face_no)
-    {
-      Assert(static_cast<unsigned int>(fe_degree + 1) ==
-                 shape_info.data.front().n_q_points_1d ||
-               fe_degree == -1,
-             ExcInternalError());
-
-      interpolate_generic<do_evaluate, add_into_output>(
-        n_components,
-        input,
-        output,
-        flags,
-        face_no,
-        shape_info.data.front().quadrature.size(),
-        shape_info.data.front().quadrature_data_on_face,
-        shape_info.n_q_points,
-        shape_info.n_q_points_face);
-    }
-
-  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 EvaluationFlags::EvaluationFlags flag,
-                        const unsigned int                     face_no,
-                        const unsigned int                     n_points_1d,
-                        const std::array<AlignedVector<Number2>, 2> &shape_data,
-                        const unsigned int dofs_per_component_on_cell,
-                        const unsigned int dofs_per_component_on_face)
-    {
-      if (face_direction == face_no / 2)
-        {
-          constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
-
-          const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
-          const int stride = Utilities::pow(n_rows, face_direction);
-          const std::array<int, 2> n_blocks{
-            {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
-          std::array<int, 2> steps;
-          if constexpr (face_direction == 0)
-            steps = {{n_rows, 0}};
-          else if constexpr (face_direction == 1 && dim == 2)
-            steps = {{1, 0}};
-          else if constexpr (face_direction == 1)
-            // in 3d, the coordinate system is zx, not xz -> switch indices
-            steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
-          else if constexpr (face_direction == 2)
-            steps = {{1, 0}};
-
-          for (unsigned int c = 0; c < n_components; ++c)
-            {
-              if (flag & EvaluationFlags::hessians)
-                interpolate_to_face<fe_degree + 1,
-                                    stride_,
-                                    do_evaluate,
-                                    add_into_output,
-                                    2>(shape_data[face_no % 2].begin(),
-                                       n_blocks,
-                                       steps,
-                                       input,
-                                       output,
-                                       n_rows,
-                                       stride);
-              else if (flag & EvaluationFlags::gradients)
-                interpolate_to_face<fe_degree + 1,
-                                    stride_,
-                                    do_evaluate,
-                                    add_into_output,
-                                    1>(shape_data[face_no % 2].begin(),
-                                       n_blocks,
-                                       steps,
-                                       input,
-                                       output,
-                                       n_rows,
-                                       stride);
-              else
-                interpolate_to_face<fe_degree + 1,
-                                    stride_,
-                                    do_evaluate,
-                                    add_into_output,
-                                    0>(shape_data[face_no % 2].begin(),
-                                       n_blocks,
-                                       steps,
-                                       input,
-                                       output,
-                                       n_rows,
-                                       stride);
-              if (do_evaluate)
-                {
-                  input += dofs_per_component_on_cell;
-                  output += dofs_per_component_on_face;
-                }
-              else
-                {
-                  output += dofs_per_component_on_cell;
-                  input += dofs_per_component_on_face;
-                }
-            }
-        }
-      else if (face_direction < dim)
-        {
-          interpolate_generic<do_evaluate,
-                              add_into_output,
-                              std::min(face_direction + 1, dim - 1)>(
-            n_components,
-            input,
-            output,
-            flag,
-            face_no,
-            n_points_1d,
-            shape_data,
-            dofs_per_component_on_cell,
-            dofs_per_component_on_face);
-        }
-    }
-
-    template <bool do_evaluate,
-              bool add_into_output,
-              int  face_direction = 0,
-              int  max_derivative = 0>
-    static void
-    interpolate_raviart_thomas(
-      const unsigned int                             n_components,
-      const Number                                  *input,
-      Number                                        *output,
-      const EvaluationFlags::EvaluationFlags         flag,
-      const unsigned int                             face_no,
-      const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info)
-    {
-      if (dim == 1)
-        {
-          // This should never happen since the FE_RaviartThomasNodal is not
-          // defined for dim = 1. It prevents compiler warnings of infinite
-          // recursion.
-          Assert(false, ExcInternalError());
-          return;
-        }
-
-      bool increase_max_der = false;
-      if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
-          (flag & EvaluationFlags::gradients && max_derivative < 1))
-        increase_max_der = true;
-
-      if (face_direction == face_no / 2 && !increase_max_der)
-        {
-          constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
-          constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
-          constexpr int stride2 = fe_degree * (fe_degree + 1);
-
-          const int degree =
-            fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
-          const int n_rows_n = degree + 1;
-          const int n_rows_t = degree;
-
-          std::array<int, 3> strides{{1, 1, 1}};
-          if (face_direction > 0)
-            {
-              strides[0] =
-                n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
-              strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
-              strides[2] = Utilities::pow(n_rows_t, face_direction);
-            }
-          const dealii::ndarray<int, 3, 3> dofs_per_direction{
-            {{{n_rows_n, n_rows_t, n_rows_t}},
-             {{n_rows_t, n_rows_n, n_rows_t}},
-             {{n_rows_t, n_rows_t, n_rows_n}}}};
-
-          std::array<int, 2> steps, n_blocks;
-
-          if constexpr (face_direction == 0)
-            steps = {{degree + (face_direction == 0), 0}};
-          else if constexpr (face_direction == 1 && dim == 2)
-            steps = {{1, 0}};
-          else if constexpr (face_direction == 1)
-            // in 3d, the coordinate system is zx, not xz -> switch indices
-            steps = {
-              {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
-          else if constexpr (face_direction == 2)
-            steps = {{1, 0}};
-
-          n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
-          n_blocks[1] =
-            dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
-
-          interpolate_to_face<
-            (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
-            ((face_direction < 2) ? stride1 : stride2),
-            do_evaluate,
-            add_into_output,
-            max_derivative>(shape_info.data[face_direction != 0]
-                              .shape_data_on_face[face_no % 2]
-                              .begin(),
-                            n_blocks,
-                            steps,
-                            input,
-                            output,
-                            degree + (face_direction == 0),
-                            strides[0]);
-
-          if (do_evaluate)
-            {
-              input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
-              output += 3 * n_blocks[0] * n_blocks[1];
-            }
-          else
-            {
-              output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
-              input += 3 * n_blocks[0] * n_blocks[1];
-            }
-
-          // must only change steps only for face direction 0
-          if constexpr (face_direction == 0)
-            steps = {{degree, 0}};
-
-          n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
-          n_blocks[1] =
-            dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
-
-          interpolate_to_face<
-            (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
-            ((face_direction < 2) ? stride0 : stride2),
-            do_evaluate,
-            add_into_output,
-            max_derivative>(shape_info.data[face_direction != 1]
-                              .shape_data_on_face[face_no % 2]
-                              .begin(),
-                            n_blocks,
-                            steps,
-                            input,
-                            output,
-                            degree + (face_direction == 1),
-                            strides[1]);
-
-          if constexpr (dim > 2)
-            {
-              if (do_evaluate)
-                {
-                  input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
-                  output += 3 * n_blocks[0] * n_blocks[1];
-                }
-              else
-                {
-                  output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
-                  input += 3 * n_blocks[0] * n_blocks[1];
-                }
-
-              if constexpr (face_direction == 0)
-                steps = {{degree, 0}};
-              else if constexpr (face_direction == 1)
-                // in 3d, the coordinate system is zx, not xz -> switch indices
-                steps = {
-                  {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
-              else if constexpr (face_direction == 2)
-                steps = {{1, 0}};
-
-              n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
-              n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
-
-              interpolate_to_face<
-                (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
-                stride0,
-                do_evaluate,
-                add_into_output,
-                max_derivative>(shape_info.data[face_direction != 2]
-                                  .shape_data_on_face[face_no % 2]
-                                  .begin(),
-                                n_blocks,
-                                steps,
-                                input,
-                                output,
-                                degree + (face_direction == 2),
-                                strides[2]);
-            }
-        }
-      else if (face_direction == face_no / 2)
-        {
-          // Only increase max_derivative
-          interpolate_raviart_thomas<do_evaluate,
-                                     add_into_output,
-                                     face_direction,
-                                     std::min(max_derivative + 1, 2)>(
-            n_components, input, output, flag, face_no, shape_info);
-        }
-      else if (face_direction < dim)
-        {
-          if (increase_max_der)
-            {
-              interpolate_raviart_thomas<do_evaluate,
-                                         add_into_output,
-                                         std::min(face_direction + 1, dim - 1),
-                                         std::min(max_derivative + 1, 2)>(
-                n_components, input, output, flag, face_no, shape_info);
-            }
-          else
-            {
-              interpolate_raviart_thomas<do_evaluate,
-                                         add_into_output,
-                                         std::min(face_direction + 1, dim - 1),
-                                         max_derivative>(
-                n_components, input, output, flag, face_no, shape_info);
-            }
-        }
-    }
-  };
-
-
-
-  // internal helper function for reading data; base version of different types
-  template <typename VectorizedArrayType, typename Number2>
-  void
-  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
-  {
-    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-      dst[v] = src_ptr[v];
-  }
-
-
-
-  // internal helper function for reading data; specialized version where we
-  // can use a dedicated load function
-  template <typename Number, std::size_t width>
-  void
-  do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
-  {
-    dst.load(src_ptr);
-  }
-
-
-
-  // internal helper function for reading data; base version of different types
-  template <typename VectorizedArrayType, typename Number2>
-  void
-  do_vectorized_gather(const Number2       *src_ptr,
-                       const unsigned int  *indices,
-                       VectorizedArrayType &dst)
-  {
-    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-      dst[v] = src_ptr[indices[v]];
-  }
-
-
-
-  // internal helper function for reading data; specialized version where we
-  // can use a dedicated gather function
-  template <typename Number, std::size_t width>
-  void
-  do_vectorized_gather(const Number                   *src_ptr,
-                       const unsigned int             *indices,
-                       VectorizedArray<Number, width> &dst)
-  {
-    dst.gather(src_ptr, indices);
-  }
-
-
-
-  // internal helper function for reading data; base version of different types
-  template <typename VectorizedArrayType, typename Number2>
-  void
-  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
-  {
-    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-      dst_ptr[v] += src[v];
-  }
-
-
-
-  // internal helper function for reading data; specialized version where we
-  // can use a dedicated load function
-  template <typename Number, std::size_t width>
-  void
-  do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
-  {
-    VectorizedArray<Number, width> tmp;
-    tmp.load(dst_ptr);
-    (tmp + src).store(dst_ptr);
-  }
-
-
-
-  // internal helper function for reading data; base version of different types
-  template <typename VectorizedArrayType, typename Number2>
-  void
-  do_vectorized_scatter_add(const VectorizedArrayType src,
-                            const unsigned int       *indices,
-                            Number2                  *dst_ptr)
-  {
-    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-      dst_ptr[indices[v]] += src[v];
-  }
-
-
-
-  // internal helper function for reading data; specialized version where we
-  // can use a dedicated gather function
-  template <typename Number, std::size_t width>
-  void
-  do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
-                            const unsigned int                  *indices,
-                            Number                              *dst_ptr)
-  {
-#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
-    for (unsigned int v = 0; v < width; ++v)
-      dst_ptr[indices[v]] += src[v];
-#else
-    VectorizedArray<Number, width> tmp;
-    tmp.gather(dst_ptr, indices);
-    (tmp + src).scatter(indices, dst_ptr);
-#endif
-  }
-
-
-
-  template <typename Number>
-  void
-  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)
-  {
-    for (unsigned int c = 0; c < n_components; ++c)
-      {
-        if (flag & EvaluationFlags::values)
-          {
-            if (integrate)
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
-            else
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
-            for (unsigned int q = 0; q < n_q_points; ++q)
-              values_quad[c * n_q_points + q] = tmp_values[q];
-          }
-        if (flag & EvaluationFlags::gradients)
-          for (unsigned int d = 0; d < dim; ++d)
-            {
-              if (integrate)
-                for (unsigned int q = 0; q < n_q_points; ++q)
-                  tmp_values[q] =
-                    gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
-              else
-                for (unsigned int q = 0; q < n_q_points; ++q)
-                  tmp_values[orientation[q]] =
-                    gradients_quad[(c * n_q_points + q) * dim + d];
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
-            }
-        if (flag & EvaluationFlags::hessians)
-          {
-            const unsigned int hdim = (dim * (dim + 1)) / 2;
-            for (unsigned int d = 0; d < hdim; ++d)
-              {
-                if (integrate)
-                  for (unsigned int q = 0; q < n_q_points; ++q)
-                    tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
-                                                  orientation[q]];
-                else
-                  for (unsigned int q = 0; q < n_q_points; ++q)
-                    tmp_values[orientation[q]] =
-                      hessians_quad[(c * hdim + d) * n_q_points + q];
-                for (unsigned int q = 0; q < n_q_points; ++q)
-                  hessians_quad[(c * hdim + d) * n_q_points + q] =
-                    tmp_values[q];
-              }
-          }
-      }
-  }
-
-
-
-  template <typename Number, typename VectorizedArrayType>
-  void
-  adjust_for_face_orientation_per_lane(
-    const unsigned int                     dim,
-    const unsigned int                     n_components,
-    const unsigned int                     v,
-    const EvaluationFlags::EvaluationFlags flag,
-    const unsigned int                    *orientation,
-    const bool                             integrate,
-    const std::size_t                      n_q_points,
-    Number                                *tmp_values,
-    VectorizedArrayType                   *values_quad,
-    VectorizedArrayType                   *gradients_quad = nullptr,
-    VectorizedArrayType                   *hessians_quad  = nullptr)
-  {
-    for (unsigned int c = 0; c < n_components; ++c)
-      {
-        if (flag & EvaluationFlags::values)
-          {
-            if (integrate)
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
-            else
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
-            for (unsigned int q = 0; q < n_q_points; ++q)
-              values_quad[c * n_q_points + q][v] = tmp_values[q];
-          }
-        if (flag & EvaluationFlags::gradients)
-          for (unsigned int d = 0; d < dim; ++d)
-            {
-              Assert(gradients_quad != nullptr, ExcInternalError());
-              if (integrate)
-                for (unsigned int q = 0; q < n_q_points; ++q)
-                  tmp_values[q] =
-                    gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
-                                  [v];
-              else
-                for (unsigned int q = 0; q < n_q_points; ++q)
-                  tmp_values[orientation[q]] =
-                    gradients_quad[(c * n_q_points + q) * dim + d][v];
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                gradients_quad[(c * n_q_points + q) * dim + d][v] =
-                  tmp_values[q];
-            }
-        if (flag & EvaluationFlags::hessians)
-          {
-            Assert(hessians_quad != nullptr, ExcInternalError());
-            const unsigned int hdim = (dim * (dim + 1)) / 2;
-            for (unsigned int d = 0; d < hdim; ++d)
-              {
-                if (integrate)
-                  for (unsigned int q = 0; q < n_q_points; ++q)
-                    tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
-                                                  orientation[q]][v];
-                else
-                  for (unsigned int q = 0; q < n_q_points; ++q)
-                    tmp_values[orientation[q]] =
-                      hessians_quad[(c * hdim + d) * n_q_points + q][v];
-                for (unsigned int q = 0; q < n_q_points; ++q)
-                  hessians_quad[(c * hdim + d) * n_q_points + q][v] =
-                    tmp_values[q];
-              }
-          }
-      }
-  }
-
-
-
-  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 EvaluationFlags::EvaluationFlags evaluation_flag,
-        const Number                          *values_dofs,
-        FEEvaluationData<dim, Number, true>   &fe_eval)
-    {
-      const auto &shape_info = fe_eval.get_shape_info();
-      const auto &shape_data = shape_info.data.front();
-      using Number2 =
-        typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
-
-      if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
-        {
-          Assert((fe_eval.get_dof_access_index() ==
-                    MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-                  fe_eval.is_interior_face() == false) == false,
-                 ExcNotImplemented());
-
-          const unsigned int face_no          = fe_eval.get_face_no();
-          const unsigned int face_orientation = fe_eval.get_face_orientation();
-          const std::size_t  n_dofs     = shape_info.dofs_per_component_on_cell;
-          const std::size_t  n_q_points = shape_info.n_q_points_faces[face_no];
-
-          using Eval =
-            EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
-
-          if (evaluation_flag & EvaluationFlags::values)
-            {
-              const auto *const shape_values =
-                &shape_data.shape_values_face(face_no, face_orientation, 0);
-
-              auto *values_quad_ptr        = fe_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)
-                {
-                  eval.template values<0, true, false>(values_dofs_actual_ptr,
-                                                       values_quad_ptr);
-
-                  values_quad_ptr += n_q_points;
-                  values_dofs_actual_ptr += n_dofs;
-                }
-            }
-
-          if (evaluation_flag & EvaluationFlags::gradients)
-            {
-              auto       *gradients_quad_ptr     = fe_eval.begin_gradients();
-              const auto *values_dofs_actual_ptr = values_dofs;
-
-              std::array<const Number2 *, dim> shape_gradients;
-              for (unsigned int d = 0; d < dim; ++d)
-                shape_gradients[d] = &shape_data.shape_gradients_face(
-                  face_no, face_orientation, d, 0);
-
-              for (unsigned int c = 0; c < n_components; ++c)
-                {
-                  for (unsigned int d = 0; d < dim; ++d)
-                    {
-                      Eval eval(nullptr,
-                                shape_gradients[d],
-                                nullptr,
-                                n_dofs,
-                                n_q_points);
-
-                      eval.template gradients<0, true, false, dim>(
-                        values_dofs_actual_ptr, gradients_quad_ptr + d);
-                    }
-                  gradients_quad_ptr += n_q_points * dim;
-                  values_dofs_actual_ptr += n_dofs;
-                }
-            }
-
-          Assert(!(evaluation_flag & EvaluationFlags::hessians),
-                 ExcNotImplemented());
-
-          return true;
-        }
-
-      const unsigned int dofs_per_face =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                         Utilities::pow(shape_data.fe_degree + 1, dim - 1);
-
-      // Note: we always keep storage of values, 1st and 2nd derivatives in an
-      // array, so reserve space for all three here
-      Number *temp         = fe_eval.get_scratch_data().begin();
-      Number *scratch_data = temp + 3 * n_components * dofs_per_face;
-
-      bool use_vectorization = true;
-
-      if (fe_eval.get_dof_access_index() ==
-            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-          fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
-        for (unsigned int v = 0; v < Number::size(); ++v)
-          if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
-              fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
-            use_vectorization = false;
-
-      if (use_vectorization == false)
-        {
-          for (unsigned int v = 0; v < Number::size(); ++v)
-            {
-              // the loop breaks once an invalid_unsigned_int is hit for
-              // all cases except the exterior faces in the ECL loop (where
-              // some faces might be at the boundaries but others not)
-              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                {
-                  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
-                       ++i)
-                    temp[i][v] = 0;
-                  continue;
-                }
-
-              FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
-                template interpolate<true, false>(n_components,
-                                                  evaluation_flag,
-                                                  shape_info,
-                                                  values_dofs,
-                                                  scratch_data,
-                                                  fe_eval.get_face_no(v));
-
-              for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
-                   ++i)
-                temp[i][v] = scratch_data[i][v];
-            }
-        }
-      else
-        FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
-          template interpolate<true, false>(n_components,
-                                            evaluation_flag,
-                                            shape_info,
-                                            values_dofs,
-                                            temp,
-                                            fe_eval.get_face_no());
-
-      const unsigned int     subface_index = fe_eval.get_subface_index();
-      constexpr unsigned int n_q_points_1d_actual =
-        fe_degree > -1 ? n_q_points_1d : 0;
-
-      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
-        {
-          FEFaceEvaluationImplRaviartThomas<dim,
-                                            fe_degree,
-                                            n_q_points_1d_actual,
-                                            Number>::
-            template evaluate_or_integrate_in_face<false>(
-              evaluation_flag,
-              fe_eval.get_shape_info().data,
-              temp,
-              fe_eval.begin_values(),
-              fe_eval.begin_gradients(),
-              scratch_data,
-              subface_index,
-              fe_eval.get_face_no() / 2);
-        }
-      else if (fe_degree > -1 &&
-               subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-               shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
-        FEFaceEvaluationImpl<true,
-                             dim,
-                             fe_degree,
-                             n_q_points_1d_actual,
-                             Number>::evaluate_in_face(n_components,
-                                                       evaluation_flag,
-                                                       shape_data,
-                                                       temp,
-                                                       fe_eval.begin_values(),
-                                                       fe_eval
-                                                         .begin_gradients(),
-                                                       fe_eval.begin_hessians(),
-                                                       scratch_data,
-                                                       subface_index);
-      else
-        FEFaceEvaluationImpl<false,
-                             dim,
-                             fe_degree,
-                             n_q_points_1d_actual,
-                             Number>::evaluate_in_face(n_components,
-                                                       evaluation_flag,
-                                                       shape_data,
-                                                       temp,
-                                                       fe_eval.begin_values(),
-                                                       fe_eval
-                                                         .begin_gradients(),
-                                                       fe_eval.begin_hessians(),
-                                                       scratch_data,
-                                                       subface_index);
-
-      if (use_vectorization == false)
-        {
-          for (unsigned int v = 0; v < Number::size(); ++v)
-            {
-              // the loop breaks once an invalid_unsigned_int is hit for
-              // all cases except the exterior faces in the ECL loop (where
-              // some faces might be at the boundaries but others not)
-              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                continue;
-
-              if (fe_eval.get_face_orientation(v) != 0)
-                adjust_for_face_orientation_per_lane(
-                  dim,
-                  n_components,
-                  v,
-                  evaluation_flag,
-                  &fe_eval.get_shape_info().face_orientations_quad(
-                    fe_eval.get_face_orientation(v), 0),
-                  false,
-                  shape_info.n_q_points_face,
-                  &temp[0][0],
-                  fe_eval.begin_values(),
-                  fe_eval.begin_gradients(),
-                  fe_eval.begin_hessians());
-            }
-        }
-      else if (fe_eval.get_face_orientation() != 0)
-        adjust_for_face_orientation(
-          dim,
-          n_components,
-          evaluation_flag,
-          &fe_eval.get_shape_info().face_orientations_quad(
-            fe_eval.get_face_orientation(), 0),
-          false,
-          shape_info.n_q_points_face,
-          temp,
-          fe_eval.begin_values(),
-          fe_eval.begin_gradients(),
-          fe_eval.begin_hessians());
-
-      return false;
-    }
-  };
-
-
-
-  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 EvaluationFlags::EvaluationFlags integration_flag,
-        Number                                *values_dofs,
-        FEEvaluationData<dim, Number, true>   &fe_eval)
-    {
-      const auto &shape_info = fe_eval.get_shape_info();
-      const auto &shape_data = shape_info.data.front();
-      using Number2 =
-        typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
-
-      if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
-        {
-          Assert((fe_eval.get_dof_access_index() ==
-                    MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-                  fe_eval.is_interior_face() == false) == false,
-                 ExcNotImplemented());
-
-          const unsigned int face_no          = fe_eval.get_face_no();
-          const unsigned int face_orientation = fe_eval.get_face_orientation();
-          const std::size_t  n_dofs     = shape_info.dofs_per_component_on_cell;
-          const std::size_t  n_q_points = shape_info.n_q_points_faces[face_no];
-
-          using Eval =
-            EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
-
-          if (integration_flag & EvaluationFlags::values)
-            {
-              const auto *const shape_values =
-                &shape_data.shape_values_face(face_no, face_orientation, 0);
-
-              auto *values_quad_ptr        = fe_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)
-                {
-                  eval.template values<0, false, false>(values_quad_ptr,
-                                                        values_dofs_actual_ptr);
-
-                  values_quad_ptr += n_q_points;
-                  values_dofs_actual_ptr += n_dofs;
-                }
-            }
-
-          if (integration_flag & EvaluationFlags::gradients)
-            {
-              auto *gradients_quad_ptr     = fe_eval.begin_gradients();
-              auto *values_dofs_actual_ptr = values_dofs;
-
-              std::array<const Number2 *, dim> shape_gradients;
-              for (unsigned int d = 0; d < dim; ++d)
-                shape_gradients[d] = &shape_data.shape_gradients_face(
-                  face_no, face_orientation, d, 0);
-
-              for (unsigned int c = 0; c < n_components; ++c)
-                {
-                  for (unsigned int d = 0; d < dim; ++d)
-                    {
-                      Eval eval(nullptr,
-                                shape_gradients[d],
-                                nullptr,
-                                n_dofs,
-                                n_q_points);
-
-                      if (!(integration_flag & EvaluationFlags::values) &&
-                          d == 0)
-                        eval.template gradients<0, false, false, dim>(
-                          gradients_quad_ptr + d, values_dofs_actual_ptr);
-                      else
-                        eval.template gradients<0, false, true, dim>(
-                          gradients_quad_ptr + d, values_dofs_actual_ptr);
-                    }
-                  gradients_quad_ptr += n_q_points * dim;
-                  values_dofs_actual_ptr += n_dofs;
-                }
-            }
-
-          Assert(!(integration_flag & EvaluationFlags::hessians),
-                 ExcNotImplemented());
-
-          return true;
-        }
-
-      const unsigned int dofs_per_face =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                         Utilities::pow(shape_data.fe_degree + 1, dim - 1);
-
-      Number *temp         = fe_eval.get_scratch_data().begin();
-      Number *scratch_data = temp + 3 * n_components * dofs_per_face;
-
-      bool use_vectorization = true;
-
-      if (fe_eval.get_dof_access_index() ==
-            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-          fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
-        use_vectorization =
-          fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
-          std::all_of(fe_eval.get_cell_ids().begin() + 1,
-                      fe_eval.get_cell_ids().end(),
-                      [&](const auto &v) {
-                        return v == fe_eval.get_cell_ids()[0] ||
-                               v == numbers::invalid_unsigned_int;
-                      });
-
-      if (use_vectorization == false)
-        {
-          for (unsigned int v = 0; v < Number::size(); ++v)
-            {
-              // the loop breaks once an invalid_unsigned_int is hit for
-              // all cases except the exterior faces in the ECL loop (where
-              // some faces might be at the boundaries but others not)
-              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                continue;
-
-              if (fe_eval.get_face_orientation(v) != 0)
-                adjust_for_face_orientation_per_lane(
-                  dim,
-                  n_components,
-                  v,
-                  integration_flag,
-                  &fe_eval.get_shape_info().face_orientations_quad(
-                    fe_eval.get_face_orientation(v), 0),
-                  true,
-                  shape_info.n_q_points_face,
-                  &temp[0][0],
-                  fe_eval.begin_values(),
-                  fe_eval.begin_gradients(),
-                  fe_eval.begin_hessians());
-            }
-        }
-      else if (fe_eval.get_face_orientation() != 0)
-        adjust_for_face_orientation(
-          dim,
-          n_components,
-          integration_flag,
-          &fe_eval.get_shape_info().face_orientations_quad(
-            fe_eval.get_face_orientation(), 0),
-          true,
-          shape_info.n_q_points_face,
-          temp,
-          fe_eval.begin_values(),
-          fe_eval.begin_gradients(),
-          fe_eval.begin_hessians());
-
-      const unsigned int n_q_points_1d_actual =
-        fe_degree > -1 ? n_q_points_1d : 0;
-      const unsigned int subface_index = fe_eval.get_subface_index();
-
-      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
-        {
-          FEFaceEvaluationImplRaviartThomas<dim,
-                                            fe_degree,
-                                            n_q_points_1d_actual,
-                                            Number>::
-            template evaluate_or_integrate_in_face<true>(
-              integration_flag,
-              fe_eval.get_shape_info().data,
-              temp,
-              fe_eval.begin_values(),
-              fe_eval.begin_gradients(),
-              scratch_data,
-              subface_index,
-              fe_eval.get_face_no() / 2);
-        }
-      else if (fe_degree > -1 &&
-               fe_eval.get_subface_index() >=
-                 GeometryInfo<dim - 1>::max_children_per_cell &&
-               shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
-        FEFaceEvaluationImpl<
-          true,
-          dim,
-          fe_degree,
-          n_q_points_1d_actual,
-          Number>::integrate_in_face(n_components,
-                                     integration_flag,
-                                     shape_data,
-                                     temp,
-                                     fe_eval.begin_values(),
-                                     fe_eval.begin_gradients(),
-                                     fe_eval.begin_hessians(),
-                                     scratch_data,
-                                     subface_index);
-      else
-        FEFaceEvaluationImpl<
-          false,
-          dim,
-          fe_degree,
-          n_q_points_1d_actual,
-          Number>::integrate_in_face(n_components,
-                                     integration_flag,
-                                     shape_data,
-                                     temp,
-                                     fe_eval.begin_values(),
-                                     fe_eval.begin_gradients(),
-                                     fe_eval.begin_hessians(),
-                                     scratch_data,
-                                     subface_index);
-
-      if (use_vectorization == false)
-        {
-          for (unsigned int v = 0; v < Number::size(); ++v)
-            {
-              // the loop breaks once an invalid_unsigned_int is hit for
-              // all cases except the exterior faces in the ECL loop (where
-              // some faces might be at the boundaries but others not)
-              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                continue;
-
-              FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
-                template interpolate<false, false>(n_components,
-                                                   integration_flag,
-                                                   shape_info,
-                                                   values_dofs,
-                                                   scratch_data,
-                                                   fe_eval.get_face_no(v));
-
-              for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
-                   ++i)
-                temp[i][v] = scratch_data[i][v];
-            }
-        }
-      else
-        FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
-          template interpolate<false, false>(n_components,
-                                             integration_flag,
-                                             shape_info,
-                                             temp,
-                                             values_dofs,
-                                             fe_eval.get_face_no());
-      return false;
-    }
-  };
-
-
-
-  template <int n_face_orientations,
-            typename Processor,
-            typename EvaluationData,
-            const bool check_face_orientations = false>
-  void
-  fe_face_evaluation_process_and_io(
-    Processor                             &proc,
-    const unsigned int                     n_components,
-    const EvaluationFlags::EvaluationFlags evaluation_flag,
-    typename Processor::Number2_          *global_vector_ptr,
-    const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
-    const EvaluationData                                             &fe_eval,
-    typename Processor::VectorizedArrayType_                         *temp1)
-  {
-    constexpr int dim         = Processor::dim_;
-    constexpr int fe_degree   = Processor::fe_degree_;
-    using VectorizedArrayType = typename Processor::VectorizedArrayType_;
-    constexpr int n_lanes     = VectorizedArrayType::size();
-
-    using Number   = typename Processor::Number_;
-    using Number2_ = typename Processor::Number2_;
-
-    const auto        &shape_data = fe_eval.get_shape_info().data.front();
-    constexpr bool     integrate  = Processor::do_integrate;
-    const unsigned int face_no    = fe_eval.get_face_no();
-    const auto        &dof_info   = fe_eval.get_dof_info();
-    const unsigned int cell       = fe_eval.get_cell_or_face_batch_id();
-    const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
-      fe_eval.get_dof_access_index();
-    AssertIndexRange(cell,
-                     dof_info.index_storage_variants[dof_access_index].size());
-    constexpr unsigned int dofs_per_face =
-      Utilities::pow(fe_degree + 1, dim - 1);
-    const unsigned int subface_index = fe_eval.get_subface_index();
-
-    const unsigned int n_filled_lanes =
-      dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-    bool all_faces_are_same = n_filled_lanes == n_lanes;
-    if (n_face_orientations == n_lanes)
-      for (unsigned int v = 1; v < n_lanes; ++v)
-        if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
-            fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
-          {
-            all_faces_are_same = false;
-            break;
-          }
-
-    // check for re-orientation ...
-    std::array<const unsigned int *, n_face_orientations> orientation = {};
-
-    if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
-        fe_eval.is_interior_face() == 0)
-      for (unsigned int v = 0; v < n_lanes; ++v)
-        {
-          // the loop breaks once an invalid_unsigned_int is hit for
-          // all cases except the exterior faces in the ECL loop (where
-          // some faces might be at the boundaries but others not)
-          if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-            continue;
-
-          if (shape_data.nodal_at_cell_boundaries &&
-              fe_eval.get_face_orientation(v) != 0)
-            {
-              // ... and in case we detect a re-orientation, go to the other
-              // version of this function that actually allows for this
-              if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
-                  check_face_orientations == false)
-                {
-                  fe_face_evaluation_process_and_io<n_face_orientations,
-                                                    Processor,
-                                                    EvaluationData,
-                                                    true>(proc,
-                                                          n_components,
-                                                          evaluation_flag,
-                                                          global_vector_ptr,
-                                                          sm_ptr,
-                                                          fe_eval,
-                                                          temp1);
-                  return;
-                }
-              orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
-                fe_eval.get_face_orientation(v), 0);
-            }
-        }
-    else if (dim == 3 && fe_eval.get_face_orientation() != 0)
-      {
-        // go to the other version of this function
-        if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
-            check_face_orientations == false)
-          {
-            fe_face_evaluation_process_and_io<n_face_orientations,
-                                              Processor,
-                                              EvaluationData,
-                                              true>(proc,
-                                                    n_components,
-                                                    evaluation_flag,
-                                                    global_vector_ptr,
-                                                    sm_ptr,
-                                                    fe_eval,
-                                                    temp1);
-            return;
-          }
-        for (unsigned int v = 0; v < n_face_orientations; ++v)
-          orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
-            fe_eval.get_face_orientation(), 0);
-      }
-
-    // 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 =
-      shape_data
-        .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
-                                                        (1 + face_no % 2))];
-
-    // face_to_cell_index_hermite
-    std::array<const unsigned int *, n_face_orientations> index_array_hermite =
-      {};
-    if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
-      {
-        if (n_face_orientations == 1)
-          index_array_hermite[0] =
-            &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
-        else
-          {
-            for (unsigned int v = 0; v < n_lanes; ++v)
-              {
-                if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                  continue;
-
-                const auto face_no = fe_eval.get_face_no(v);
-
-                grad_weight[v] =
-                  shape_data.shape_data_on_face[0][fe_degree +
-                                                   (integrate ?
-                                                      (2 - (face_no % 2)) :
-                                                      (1 + (face_no % 2)))][0];
-
-                index_array_hermite[v] =
-                  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
-                                                                       0);
-              }
-          }
-      }
-
-    // face_to_cell_index_nodal
-    std::array<const unsigned int *, n_face_orientations> index_array_nodal =
-      {};
-    if (shape_data.nodal_at_cell_boundaries == true)
-      {
-        if (n_face_orientations == 1)
-          index_array_nodal[0] =
-            &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
-        else
-          {
-            for (unsigned int v = 0; v < n_lanes; ++v)
-              {
-                if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                  continue;
-
-                const auto face_no = fe_eval.get_face_no(v);
-
-                index_array_nodal[v] =
-                  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
-                                                                     0);
-              }
-          }
-      }
-
-
-    const auto reorientate = [&](const unsigned int v, const unsigned int i) {
-      return (!check_face_orientations || orientation[v] == nullptr) ?
-               i :
-               orientation[v][i];
-    };
-
-    const unsigned int cell_index =
-      dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
-        fe_eval.get_cell_ids()[0] :
-        cell * n_lanes;
-    const unsigned int *dof_indices =
-      &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
-
-    for (unsigned int comp = 0; comp < n_components; ++comp)
-      {
-        const std::size_t index_offset =
-          dof_info.component_dof_indices_offset
-            [fe_eval.get_active_fe_index()]
-            [fe_eval.get_first_selected_component()] +
-          comp * Utilities::pow(fe_degree + 1, dim);
-
-        // 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(
-              dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
-              n_lanes);
-            Number2_ *vector_ptr =
-              global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
-
-            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
-              {
-                for (unsigned int i = 0; i < dofs_per_face; ++i)
-                  {
-                    Assert(n_face_orientations == 1, ExcNotImplemented());
-
-                    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(temp1[i_],
-                                                 temp1[i_ + dofs_per_face],
-                                                 vector_ptr + ind1 * n_lanes,
-                                                 vector_ptr + ind2 * n_lanes,
-                                                 grad_weight);
-                  }
-              }
-            else
-              {
-                for (unsigned int i = 0; i < dofs_per_face; ++i)
-                  {
-                    Assert(n_face_orientations == 1, 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 * n_lanes);
-                  }
-              }
-          }
-
-        // 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(
-              dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
-              n_lanes);
-            Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
-            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
-              {
-                for (unsigned int i = 0; i < dofs_per_face; ++i)
-                  {
-                    Assert(n_face_orientations == 1, ExcNotImplemented());
-
-                    const unsigned int i_ = reorientate(0, i);
-                    const unsigned int ind1 =
-                      index_array_hermite[0][2 * i] * n_lanes;
-                    const unsigned int ind2 =
-                      index_array_hermite[0][2 * i + 1] * n_lanes;
-                    proc.hermite_grad_vectorized_indexed(
-                      temp1[i_],
-                      temp1[i_ + dofs_per_face],
-                      vector_ptr + ind1,
-                      vector_ptr + ind2,
-                      grad_weight,
-                      dof_indices,
-                      dof_indices);
-                  }
-              }
-            else
-              {
-                for (unsigned int i = 0; i < dofs_per_face; ++i)
-                  {
-                    Assert(n_face_orientations == 1, ExcNotImplemented());
-
-                    const unsigned int i_  = reorientate(0, i);
-                    const unsigned int ind = index_array_nodal[0][i] * n_lanes;
-                    proc.value_vectorized_indexed(temp1[i_],
-                                                  vector_ptr + ind,
-                                                  dof_indices);
-                  }
-              }
-          }
-
-        // 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)
-          {
-            const unsigned int *strides =
-              &dof_info.dof_indices_interleave_strides[dof_access_index]
-                                                      [cell * n_lanes];
-            unsigned int indices[n_lanes];
-            for (unsigned int v = 0; v < n_lanes; ++v)
-              indices[v] = dof_indices[v] + index_offset * strides[v];
-            const unsigned int n_filled_lanes =
-              dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
-              {
-                if (n_filled_lanes == n_lanes)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      Assert(n_face_orientations == 1, ExcNotImplemented());
-
-                      const unsigned int i_ = reorientate(0, i);
-                      unsigned int       ind1[n_lanes];
-                      DEAL_II_OPENMP_SIMD_PRAGMA
-                      for (unsigned int v = 0; v < n_lanes; ++v)
-                        ind1[v] = indices[v] +
-                                  index_array_hermite[0][2 * i] * strides[v];
-                      unsigned int ind2[n_lanes];
-                      DEAL_II_OPENMP_SIMD_PRAGMA
-                      for (unsigned int v = 0; v < n_lanes; ++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],
-                        global_vector_ptr,
-                        global_vector_ptr,
-                        grad_weight,
-                        ind1,
-                        ind2);
-                    }
-                else
-                  {
-                    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]);
-                        }
-                  }
-              }
-            else
-              {
-                if (n_filled_lanes == n_lanes)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      Assert(n_face_orientations == 1, ExcInternalError());
-                      unsigned int ind[n_lanes];
-                      DEAL_II_OPENMP_SIMD_PRAGMA
-                      for (unsigned int v = 0; v < n_lanes; ++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
-                  {
-                    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]]);
-                  }
-              }
-          }
-
-        // 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)
-          {
-            Number2_ *vector_ptr = global_vector_ptr + index_offset;
-
-            const bool vectorization_possible =
-              all_faces_are_same && (sm_ptr == nullptr);
-
-            std::array<Number2_ *, n_lanes>   vector_ptrs;
-            std::array<unsigned int, n_lanes> reordered_indices;
-
-            if (vectorization_possible == false)
-              {
-                vector_ptrs = {};
-                if (n_face_orientations == 1)
-                  {
-                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                      if (sm_ptr == nullptr)
-                        {
-                          vector_ptrs[v] = vector_ptr + dof_indices[v];
-                        }
-                      else
-                        {
-                          const auto &temp =
-                            dof_info
-                              .dof_indices_contiguous_sm[dof_access_index]
-                                                        [cell * n_lanes + v];
-                          vector_ptrs[v] = const_cast<Number2_ *>(
-                            sm_ptr->operator[](temp.first).data() +
-                            temp.second + index_offset);
-                        }
-                  }
-                else if (n_face_orientations == n_lanes)
-                  {
-                    const auto &cells = fe_eval.get_cell_ids();
-                    for (unsigned int v = 0; v < n_lanes; ++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<Number2_ *>(
-                                sm_ptr->operator[](temp.first).data() +
-                                temp.second + index_offset);
-                            }
-                        }
-                  }
-                else
-                  {
-                    Assert(false, ExcNotImplemented());
-                  }
-              }
-            else if (n_face_orientations == n_lanes)
-              {
-                for (unsigned int v = 0; v < n_lanes; ++v)
-                  reordered_indices[v] =
-                    dof_info.dof_indices_contiguous[dof_access_index]
-                                                   [fe_eval.get_cell_ids()[v]];
-                dof_indices = reordered_indices.data();
-              }
-
-            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
-              {
-                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,
-                        dof_indices,
-                        dof_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 < n_lanes; ++v)
-                          {
-                            temp1[i][v]                 = 0.0;
-                            temp1[i + dofs_per_face][v] = 0.0;
-                          }
-                    }
-                else
-                  {
-                    if (integrate == false && n_filled_lanes < n_lanes)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        temp1[i] = temp1[i + dofs_per_face] = Number();
-
-                    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)
-                  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,
-                                                    dof_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 < n_lanes; ++v)
-                          temp1[i_][v] = 0.0;
-                    }
-                else
-                  {
-                    if (integrate == false && n_filled_lanes < n_lanes)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        temp1[i] = Number();
-
-                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        proc.value(temp1[reorientate(v, i)][v],
-                                   vector_ptrs[v][index_array_nodal[v][i]]);
-                  }
-              }
-          }
-        else
-          {
-            // We should not end up here, this should be caught by
-            // FEFaceEvaluationImplGatherEvaluateSelector::supports()
-            Assert(false, ExcInternalError());
-          }
-        temp1 += 3 * dofs_per_face;
-      }
-  }
-
-
-
-  template <int dim, typename Number2, typename VectorizedArrayType>
-  struct FEFaceEvaluationImplGatherEvaluateSelector
-  {
-    using Number = typename VectorizedArrayType::value_type;
-
-    template <int fe_degree, int n_q_points_1d>
-    static bool
-    run(const unsigned int                                n_components,
-        const EvaluationFlags::EvaluationFlags            evaluation_flag,
-        const Number2                                    *src_ptr,
-        const std::vector<ArrayView<const Number2>>      *sm_ptr,
-        FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
-    {
-      Assert(fe_degree > -1, ExcInternalError());
-      Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric_no_collocation,
-             ExcInternalError());
-
-      const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
-
-      VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
-      VectorizedArrayType *scratch_data =
-        temp + 3 * n_components * dofs_per_face;
-
-      Processor<fe_degree> p;
-
-      if (fe_eval.get_dof_access_index() ==
-            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-          fe_eval.is_interior_face() == false)
-        fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
-          p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
-      else
-        fe_face_evaluation_process_and_io<1>(
-          p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
-
-      const unsigned int subface_index = fe_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,
-                           fe_eval.get_shape_info().data.front(),
-                           temp,
-                           fe_eval.begin_values(),
-                           fe_eval.begin_gradients(),
-                           fe_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,
-                           fe_eval.get_shape_info().data.front(),
-                           temp,
-                           fe_eval.begin_values(),
-                           fe_eval.begin_gradients(),
-                           fe_eval.begin_hessians(),
-                           scratch_data,
-                           subface_index);
-
-      // re-orientation for cases not possible with above algorithm
-      if (subface_index < GeometryInfo<dim>::max_children_per_cell)
-        {
-          if (fe_eval.get_dof_access_index() ==
-                MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-              fe_eval.is_interior_face() == false)
-            {
-              for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-                {
-                  // the loop breaks once an invalid_unsigned_int is hit for
-                  // all cases except the exterior faces in the ECL loop (where
-                  // some faces might be at the boundaries but others not)
-                  if (fe_eval.get_cell_ids()[v] ==
-                      numbers::invalid_unsigned_int)
-                    continue;
-
-                  if (fe_eval.get_face_orientation(v) != 0)
-                    adjust_for_face_orientation_per_lane(
-                      dim,
-                      n_components,
-                      v,
-                      evaluation_flag,
-                      &fe_eval.get_shape_info().face_orientations_quad(
-                        fe_eval.get_face_orientation(v), 0),
-                      false,
-                      Utilities::pow(n_q_points_1d, dim - 1),
-                      &temp[0][0],
-                      fe_eval.begin_values(),
-                      fe_eval.begin_gradients(),
-                      fe_eval.begin_hessians());
-                }
-            }
-          else if (fe_eval.get_face_orientation() != 0)
-            adjust_for_face_orientation(
-              dim,
-              n_components,
-              evaluation_flag,
-              &fe_eval.get_shape_info().face_orientations_quad(
-                fe_eval.get_face_orientation(), 0),
-              false,
-              Utilities::pow(n_q_points_1d, dim - 1),
-              temp,
-              fe_eval.begin_values(),
-              fe_eval.begin_gradients(),
-              fe_eval.begin_hessians());
-        }
-
-      return false;
-    }
-
-    template <typename Number3>
-    static bool
-    supports(const EvaluationFlags::EvaluationFlags             evaluation_flag,
-             const MatrixFreeFunctions::ShapeInfo<Number3>     &shape_info,
-             const Number2                                     *vector_ptr,
-             MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
-    {
-      const unsigned int fe_degree = shape_info.data.front().fe_degree;
-      if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
-          (evaluation_flag & EvaluationFlags::gradients &&
-           (fe_degree < 2 ||
-            shape_info.data.front().element_type !=
-              MatrixFreeFunctions::tensor_symmetric_hermite)) ||
-          (evaluation_flag & EvaluationFlags::hessians) ||
-          vector_ptr == nullptr ||
-          shape_info.data.front().element_type >
-            MatrixFreeFunctions::tensor_symmetric_no_collocation ||
-          storage <
-            MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
-        return false;
-      else
-        return true;
-    }
-
-  private:
-    template <int fe_degree>
-    struct Processor
-    {
-      static const bool do_integrate = false;
-      static const int  dim_         = dim;
-      static const int  fe_degree_   = fe_degree;
-      using VectorizedArrayType_     = VectorizedArrayType;
-      using Number_                  = Number;
-      using Number2_                 = const Number2;
-
-      template <typename T0, typename T1, typename T2>
-      void
-      hermite_grad_vectorized(T0       &temp_1,
-                              T0       &temp_2,
-                              const T1  src_ptr_1,
-                              const T1  src_ptr_2,
-                              const T2 &grad_weight)
-      {
-        do_vectorized_read(src_ptr_1, temp_1);
-        do_vectorized_read(src_ptr_2, temp_2);
-        temp_2 = grad_weight * (temp_1 - temp_2);
-      }
-
-      template <typename T1, typename T2>
-      void
-      value_vectorized(T1 &temp, const T2 src_ptr)
-      {
-        do_vectorized_read(src_ptr, temp);
-      }
-
-      template <typename T0, typename T1, typename T2, typename T3>
-      void
-      hermite_grad_vectorized_indexed(T0       &temp_1,
-                                      T0       &temp_2,
-                                      const T1  src_ptr_1,
-                                      const T1  src_ptr_2,
-                                      const T2 &grad_weight,
-                                      const T3 &indices_1,
-                                      const T3 &indices_2)
-      {
-        do_vectorized_gather(src_ptr_1, indices_1, temp_1);
-        do_vectorized_gather(src_ptr_2, indices_2, temp_2);
-        temp_2 = grad_weight * (temp_1 - temp_2);
-      }
-
-      template <typename T0, typename T1, typename T2>
-      void
-      value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
-      {
-        do_vectorized_gather(src_ptr, indices, temp);
-      }
-
-      template <typename T0, typename T1, typename T2>
-      void
-      hermite_grad(T0       &temp_1,
-                   T0       &temp_2,
-                   const T1 &src_ptr_1,
-                   const T1 &src_ptr_2,
-                   const T2 &grad_weight)
-      {
-        // case 3a)
-        temp_1 = src_ptr_1;
-        temp_2 = grad_weight * (temp_1 - src_ptr_2);
-      }
-
-      template <typename T1, typename T2>
-      void
-      value(T1 &temp, const T2 &src_ptr)
-      {
-        // case 3b)
-        temp = src_ptr;
-      }
-    };
-  };
-
-
-
-  template <int dim, typename Number2, typename VectorizedArrayType>
-  struct FEFaceEvaluationImplIntegrateScatterSelector
-  {
-    using Number = typename VectorizedArrayType::value_type;
-
-    template <int fe_degree, int n_q_points_1d>
-    static bool
-    run(const unsigned int                                n_components,
-        const EvaluationFlags::EvaluationFlags            integration_flag,
-        Number2                                          *dst_ptr,
-        const std::vector<ArrayView<const Number2>>      *sm_ptr,
-        FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
-    {
-      Assert(fe_degree > -1, ExcInternalError());
-      Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric_no_collocation,
-             ExcInternalError());
-
-      const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
-
-      VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
-      VectorizedArrayType *scratch_data =
-        temp + 3 * n_components * dofs_per_face;
-
-      const unsigned int subface_index = fe_eval.get_subface_index();
-
-      // re-orientation for cases not possible with the io function below
-      if (subface_index < GeometryInfo<dim>::max_children_per_cell)
-        {
-          if (fe_eval.get_dof_access_index() ==
-                MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-              fe_eval.is_interior_face() == false)
-            for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-              {
-                // the loop breaks once an invalid_unsigned_int is hit for
-                // all cases except the exterior faces in the ECL loop (where
-                // some faces might be at the boundaries but others not)
-                if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
-                  continue;
-
-                if (fe_eval.get_face_orientation(v) != 0)
-                  adjust_for_face_orientation_per_lane(
-                    dim,
-                    n_components,
-                    v,
-                    integration_flag,
-                    &fe_eval.get_shape_info().face_orientations_quad(
-                      fe_eval.get_face_orientation(v), 0),
-                    true,
-                    Utilities::pow(n_q_points_1d, dim - 1),
-                    &temp[0][0],
-                    fe_eval.begin_values(),
-                    fe_eval.begin_gradients(),
-                    fe_eval.begin_hessians());
-              }
-          else if (fe_eval.get_face_orientation() != 0)
-            adjust_for_face_orientation(
-              dim,
-              n_components,
-              integration_flag,
-              &fe_eval.get_shape_info().face_orientations_quad(
-                fe_eval.get_face_orientation(), 0),
-              true,
-              Utilities::pow(n_q_points_1d, dim - 1),
-              temp,
-              fe_eval.begin_values(),
-              fe_eval.begin_gradients(),
-              fe_eval.begin_hessians());
-        }
-
-      if (fe_degree > -1 && fe_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,
-                            fe_eval.get_shape_info().data.front(),
-                            temp,
-                            fe_eval.begin_values(),
-                            fe_eval.begin_gradients(),
-                            fe_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,
-                            fe_eval.get_shape_info().data.front(),
-                            temp,
-                            fe_eval.begin_values(),
-                            fe_eval.begin_gradients(),
-                            fe_eval.begin_hessians(),
-                            scratch_data,
-                            subface_index);
-
-      Processor<fe_degree> p;
-
-      if (fe_eval.get_dof_access_index() ==
-            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-          fe_eval.is_interior_face() == false)
-        fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
-          p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
-      else
-        fe_face_evaluation_process_and_io<1>(
-          p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
-
-      return false;
-    }
-
-  private:
-    template <int fe_degree>
-    struct Processor
-    {
-      static const bool do_integrate = true;
-      static const int  dim_         = dim;
-      static const int  fe_degree_   = fe_degree;
-      using VectorizedArrayType_     = VectorizedArrayType;
-      using Number_                  = Number;
-      using Number2_                 = Number2;
-
-      template <typename T0, typename T1, typename T2, typename T3, typename T4>
-      void
-      hermite_grad_vectorized(const T0 &temp_1,
-                              const T1 &temp_2,
-                              T2        dst_ptr_1,
-                              T3        dst_ptr_2,
-                              const T4 &grad_weight)
-      {
-        // case 1a)
-        const VectorizedArrayType val  = temp_1 - grad_weight * temp_2;
-        const VectorizedArrayType grad = grad_weight * temp_2;
-        do_vectorized_add(val, dst_ptr_1);
-        do_vectorized_add(grad, dst_ptr_2);
-      }
-
-      template <typename T0, typename T1>
-      void
-      value_vectorized(const T0 &temp, T1 dst_ptr)
-      {
-        // case 1b)
-        do_vectorized_add(temp, dst_ptr);
-      }
-
-      template <typename T0, typename T1, typename T2, typename T3>
-      void
-      hermite_grad_vectorized_indexed(const T0 &temp_1,
-                                      const T0 &temp_2,
-                                      T1        dst_ptr_1,
-                                      T1        dst_ptr_2,
-                                      const T2 &grad_weight,
-                                      const T3 &indices_1,
-                                      const T3 &indices_2)
-      {
-        // case 2a)
-        const VectorizedArrayType val  = temp_1 - grad_weight * temp_2;
-        const VectorizedArrayType grad = grad_weight * temp_2;
-        do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
-        do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
-      }
-
-      template <typename T0, typename T1, typename T2>
-      void
-      value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
-      {
-        // case 2b)
-        do_vectorized_scatter_add(temp, indices, dst_ptr);
-      }
-
-      template <typename T0, typename T1, typename T2>
-      void
-      hermite_grad(const T0 &temp_1,
-                   const T0 &temp_2,
-                   T1       &dst_ptr_1,
-                   T1       &dst_ptr_2,
-                   const T2 &grad_weight)
-      {
-        // case 3a)
-        const Number val  = temp_1 - grad_weight * temp_2;
-        const Number grad = grad_weight * temp_2;
-        dst_ptr_1 += val;
-        dst_ptr_2 += grad;
-      }
-
-      template <typename T0, typename T1>
-      void
-      value(const T0 &temp, T1 &dst_ptr)
-      {
-        // case 3b)
-        dst_ptr += temp;
-      }
-    };
-  };
-
-
-
-  /**
-   * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation,
-   * using an FEEvaluationData argument.
-   */
-  template <int dim, typename Number>
-  struct CellwiseInverseMassMatrixImplBasic
-  {
-    using Number2 =
-      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-
-    template <int fe_degree, int = 0>
-    static bool
-    run(const unsigned int                          n_components,
-        const FEEvaluationData<dim, Number, false> &fe_eval,
-        const Number                               *in_array,
-        Number                                     *out_array)
-    {
-      const unsigned int given_degree =
-        (fe_degree > -1) ? fe_degree :
-                           fe_eval.get_shape_info().data.front().fe_degree;
-
-      const unsigned int dofs_per_component =
-        Utilities::pow(given_degree + 1, dim);
-
-      Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
-      Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric_no_collocation,
-             ExcNotImplemented());
-
-      EvaluatorTensorProduct<evaluate_evenodd,
-                             dim,
-                             fe_degree + 1,
-                             fe_degree + 1,
-                             Number,
-                             Number2>
-        evaluator({},
-                  {},
-                  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
-                  given_degree + 1,
-                  given_degree + 1);
-
-      for (unsigned int d = 0; d < n_components; ++d)
-        {
-          const Number *in  = in_array + d * dofs_per_component;
-          Number       *out = out_array + d * dofs_per_component;
-          // Need to select 'apply' method with hessian slot because values
-          // assume symmetries that do not exist in the inverse shapes
-          evaluator.template hessians<0, true, false>(in, out);
-          if (dim > 1)
-            evaluator.template hessians<1, true, false>(out, out);
-          if (dim > 2)
-            evaluator.template hessians<2, true, false>(out, out);
-        }
-      for (unsigned int q = 0; q < dofs_per_component; ++q)
-        {
-          const Number inverse_JxW_q = Number(1.) / fe_eval.JxW(q);
-          for (unsigned int d = 0; d < n_components; ++d)
-            out_array[q + d * dofs_per_component] *= inverse_JxW_q;
-        }
-      for (unsigned int d = 0; d < n_components; ++d)
-        {
-          Number *out = out_array + d * dofs_per_component;
-          if (dim > 2)
-            evaluator.template hessians<2, false, false>(out, out);
-          if (dim > 1)
-            evaluator.template hessians<1, false, false>(out, out);
-          evaluator.template hessians<0, false, false>(out, out);
-        }
-      return false;
-    }
-  };
-
-
-
-  /**
-   * This struct implements the action of the inverse @ref GlossMassMatrix "mass matrix" operation
-   * with user-provided coefficients at quadrature points (in contrast to
-   * CellwiseInverseMassMatrixImplBasic, which implicitly uses `1/(|J|xW)' as
-   * coefficient).
-   */
-  template <int dim, typename Number>
-  struct CellwiseInverseMassMatrixImplFlexible
-  {
-    using Number2 =
-      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-
-    template <int fe_degree, int = 0>
-    static bool
-    run(const unsigned int                          n_desired_components,
-        const FEEvaluationData<dim, Number, false> &fe_eval,
-        const ArrayView<const Number>              &inverse_coefficients,
-        const bool                                  dyadic_coefficients,
-        const Number                               *in_array,
-        Number                                     *out_array)
-    {
-      const unsigned int given_degree =
-        (fe_degree > -1) ? fe_degree :
-                           fe_eval.get_shape_info().data.front().fe_degree;
-
-      const unsigned int dofs_per_component =
-        Utilities::pow(given_degree + 1, dim);
-
-      Assert(inverse_coefficients.size() > 0 &&
-               inverse_coefficients.size() % dofs_per_component == 0,
-             ExcMessage(
-               "Expected diagonal to be a multiple of scalar dof per cells"));
-
-      if (!dyadic_coefficients)
-        {
-          if (inverse_coefficients.size() != dofs_per_component)
-            AssertDimension(n_desired_components * dofs_per_component,
-                            inverse_coefficients.size());
-        }
-      else
-        {
-          AssertDimension(n_desired_components * n_desired_components *
-                            dofs_per_component,
-                          inverse_coefficients.size());
-        }
-
-      Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
-      Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric_no_collocation,
-             ExcNotImplemented());
-
-      EvaluatorTensorProduct<evaluate_evenodd,
-                             dim,
-                             fe_degree + 1,
-                             fe_degree + 1,
-                             Number,
-                             Number2>
-        evaluator({},
-                  {},
-                  fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
-                  given_degree + 1,
-                  given_degree + 1);
-
-      const Number *in  = in_array;
-      Number       *out = out_array;
-
-      const Number *inv_coefficient = inverse_coefficients.data();
-
-      const unsigned int shift_coefficient =
-        inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
-                                                           0;
-
-      const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components;
-      const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1;
-
-      for (unsigned int d = 0; d < n_comp_outer; ++d)
-        {
-          for (unsigned int di = 0; di < n_comp_inner; ++di)
-            {
-              const Number *in_  = in + di * dofs_per_component;
-              Number       *out_ = out + di * dofs_per_component;
-              evaluator.template hessians<0, true, false>(in_, out_);
-              if (dim > 1)
-                evaluator.template hessians<1, true, false>(out_, out_);
-              if (dim > 2)
-                evaluator.template hessians<2, true, false>(out_, out_);
-            }
-          if (dyadic_coefficients)
-            {
-              const auto n_coeff_components =
-                n_desired_components * n_desired_components;
-              if (n_desired_components == dim)
-                {
-                  for (unsigned int q = 0; q < dofs_per_component; ++q)
-                    vmult<dim>(&inv_coefficient[q * n_coeff_components],
-                               &in[q],
-                               &out[q],
-                               dofs_per_component);
-                }
-              else
-                {
-                  for (unsigned int q = 0; q < dofs_per_component; ++q)
-                    vmult<-1>(&inv_coefficient[q * n_coeff_components],
-                              &in[q],
-                              &out[q],
-                              dofs_per_component,
-                              n_desired_components);
+                  for (unsigned int q = 0; q < dofs_per_component; ++q)
+                    vmult<-1>(&inv_coefficient[q * n_coeff_components],
+                              &in[q],
+                              &out[q],
+                              dofs_per_component,
+                              n_desired_components);
                 }
             }
           else
@@ -5567,21 +2518,6 @@ namespace internal
       return false;
     }
   };
-
-  /**
-   * This struct is used to implement
-   * FEEvaluation::fast_evaluation_supported() and
-   * FEFaceEvaluation::fast_evaluation_supported().
-   */
-  struct FastEvaluationSupported
-  {
-    template <int fe_degree, int n_q_points_1d>
-    static bool
-    run()
-    {
-      return fe_degree != -1;
-    }
-  };
 } // end of namespace internal
 
 
diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h
new file mode 100644 (file)
index 0000000..09a01d4
--- /dev/null
@@ -0,0 +1,2920 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_evaluation_kernels_face_h
+#define dealii_matrix_free_evaluation_kernels_face_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/ndarray.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/dof_info.h>
+#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/fe_evaluation_data.h>
+#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  template <bool symmetric_evaluate,
+            int  dim,
+            int  fe_degree,
+            int  n_q_points_1d,
+            typename Number>
+  struct FEFaceEvaluationImpl
+  {
+    // We enable a transformation to collocation for derivatives if it gives
+    // correct results (first two conditions), if it is the most efficient
+    // choice in terms of operation counts (third condition) and if we were
+    // able to initialize the fields in shape_info.templates.h from the
+    // polynomials (fourth condition).
+    using Number2 =
+      typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
+
+    using Eval = EvaluatorTensorProduct<symmetric_evaluate ? evaluate_evenodd :
+                                                             evaluate_general,
+                                        dim - 1,
+                                        fe_degree + 1,
+                                        n_q_points_1d,
+                                        Number,
+                                        Number2>;
+
+    static Eval
+    create_evaluator_tensor_product(
+      const MatrixFreeFunctions::UnivariateShapeData<Number2> &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 EvaluationFlags::EvaluationFlags                   evaluation_flag,
+      const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+      Number                                                  *values_dofs,
+      Number                                                  *values_quad,
+      Number                                                  *gradients_quad,
+      Number                                                  *hessians_quad,
+      Number                                                  *scratch_data,
+      const unsigned int                                       subface_index)
+    {
+      Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
+      Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
+
+      const std::size_t n_dofs = 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 ((evaluation_flag & EvaluationFlags::values) != 0u &&
+          ((evaluation_flag & EvaluationFlags::gradients) == 0u))
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            switch (dim)
+              {
+                case 3:
+                  eval0.template values<0, true, false>(values_dofs,
+                                                        values_quad);
+                  eval1.template values<1, true, false>(values_quad,
+                                                        values_quad);
+                  break;
+                case 2:
+                  eval0.template values<0, true, false>(values_dofs,
+                                                        values_quad);
+                  break;
+                case 1:
+                  values_quad[0] = values_dofs[0];
+                  break;
+                default:
+                  Assert(false, ExcNotImplemented());
+              }
+            // Note: we always keep storage of values, 1st and 2nd derivatives
+            // in an array
+            values_dofs += 3 * n_dofs;
+            values_quad += n_q_points;
+          }
+      else if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            switch (dim)
+              {
+                case 3:
+                  if (symmetric_evaluate &&
+                      use_collocation_evaluation(fe_degree, n_q_points_1d))
+                    {
+                      eval0.template values<0, true, false>(values_dofs,
+                                                            values_quad);
+                      eval0.template values<1, true, false>(values_quad,
+                                                            values_quad);
+                      EvaluatorTensorProduct<evaluate_evenodd,
+                                             dim - 1,
+                                             n_q_points_1d,
+                                             n_q_points_1d,
+                                             Number,
+                                             Number2>
+                        eval_grad({}, data.shape_gradients_collocation_eo, {});
+                      eval_grad.template gradients<0, true, false, 3>(
+                        values_quad, gradients_quad);
+                      eval_grad.template gradients<1, true, false, 3>(
+                        values_quad, gradients_quad + 1);
+                    }
+                  else
+                    {
+                      // grad x
+                      eval0.template gradients<0, true, false>(values_dofs,
+                                                               scratch_data);
+                      eval1.template values<1, true, false, 3>(scratch_data,
+                                                               gradients_quad);
+
+                      // grad y
+                      eval0.template values<0, true, false>(values_dofs,
+                                                            scratch_data);
+                      eval1.template gradients<1, true, false, 3>(
+                        scratch_data, gradients_quad + 1);
+
+                      if ((evaluation_flag & EvaluationFlags::values) != 0u)
+                        eval1.template values<1, true, false>(scratch_data,
+                                                              values_quad);
+                    }
+                  // grad z
+                  eval0.template values<0, true, false>(values_dofs + n_dofs,
+                                                        scratch_data);
+                  eval1.template values<1, true, false, 3>(scratch_data,
+                                                           gradients_quad + 2);
+
+                  break;
+                case 2:
+                  eval0.template values<0, true, false, 2>(values_dofs + n_dofs,
+                                                           gradients_quad + 1);
+                  eval0.template gradients<0, true, false, 2>(values_dofs,
+                                                              gradients_quad);
+                  if ((evaluation_flag & EvaluationFlags::values) != 0u)
+                    eval0.template values<0, true, false>(values_dofs,
+                                                          values_quad);
+                  break;
+                case 1:
+                  values_quad[0]    = values_dofs[0];
+                  gradients_quad[0] = values_dofs[1];
+                  break;
+                default:
+                  AssertThrow(false, ExcNotImplemented());
+              }
+            values_dofs += 3 * n_dofs;
+            values_quad += n_q_points;
+            gradients_quad += dim * n_q_points;
+          }
+
+      if ((evaluation_flag & EvaluationFlags::hessians) != 0u)
+        {
+          values_dofs = values_dofs_ptr;
+          for (unsigned int c = 0; c < n_components; ++c)
+            {
+              switch (dim)
+                {
+                  case 3:
+                    // grad xx
+                    eval0.template hessians<0, true, false>(values_dofs,
+                                                            scratch_data);
+                    eval1.template values<1, true, false>(scratch_data,
+                                                          hessians_quad);
+
+                    // grad yy
+                    eval0.template values<0, true, false>(values_dofs,
+                                                          scratch_data);
+                    eval1.template hessians<1, true, false>(scratch_data,
+                                                            hessians_quad +
+                                                              n_q_points);
+
+                    // grad zz
+                    eval0.template values<0, true, false>(values_dofs +
+                                                            2 * n_dofs,
+                                                          scratch_data);
+                    eval1.template values<1, true, false>(scratch_data,
+                                                          hessians_quad +
+                                                            2 * n_q_points);
+
+                    // grad xy
+                    eval0.template gradients<0, true, false>(values_dofs,
+                                                             scratch_data);
+                    eval1.template gradients<1, true, false>(scratch_data,
+                                                             hessians_quad +
+                                                               3 * n_q_points);
+
+                    // grad xz
+                    eval0.template gradients<0, true, false>(values_dofs +
+                                                               n_dofs,
+                                                             scratch_data);
+                    eval1.template values<1, true, false>(scratch_data,
+                                                          hessians_quad +
+                                                            4 * n_q_points);
+
+                    // grad yz
+                    eval0.template values<0, true, false>(values_dofs + n_dofs,
+                                                          scratch_data);
+                    eval1.template gradients<1, true, false>(scratch_data,
+                                                             hessians_quad +
+                                                               5 * n_q_points);
+
+                    break;
+                  case 2:
+                    // grad xx
+                    eval0.template hessians<0, true, false>(values_dofs,
+                                                            hessians_quad);
+                    // grad yy
+                    eval0.template values<0, true, false>(
+                      values_dofs + 2 * n_dofs, hessians_quad + n_q_points);
+                    // grad xy
+                    eval0.template gradients<0, true, false>(
+                      values_dofs + n_dofs, hessians_quad + 2 * n_q_points);
+                    break;
+                  case 1:
+                    hessians_quad[0] = values_dofs[2];
+                    break;
+                  default:
+                    AssertThrow(false, ExcNotImplemented());
+                }
+              values_dofs += 3 * n_dofs;
+              hessians_quad += dim * (dim + 1) / 2 * n_q_points;
+            }
+        }
+    }
+
+    static void
+    integrate_in_face(
+      const unsigned int                                       n_components,
+      const EvaluationFlags::EvaluationFlags                   integration_flag,
+      const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+      Number                                                  *values_dofs,
+      Number                                                  *values_quad,
+      Number                                                  *gradients_quad,
+      Number                                                  *hessians_quad,
+      Number                                                  *scratch_data,
+      const unsigned int                                       subface_index)
+    {
+      Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0);
+      Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1);
+
+      const std::size_t n_dofs =
+        fe_degree > -1 ?
+          Utilities::pow(fe_degree + 1, dim - 1) :
+          (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 ((integration_flag & EvaluationFlags::values) != 0u &&
+          (integration_flag & EvaluationFlags::gradients) == 0u)
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            switch (dim)
+              {
+                case 3:
+                  eval1.template values<1, false, false>(values_quad,
+                                                         values_quad);
+                  eval0.template values<0, false, false>(values_quad,
+                                                         values_dofs);
+                  break;
+                case 2:
+                  eval0.template values<0, false, false>(values_quad,
+                                                         values_dofs);
+                  break;
+                case 1:
+                  values_dofs[0] = values_quad[0];
+                  break;
+                default:
+                  Assert(false, ExcNotImplemented());
+              }
+            values_dofs += 3 * n_dofs;
+            values_quad += n_q_points;
+          }
+      else if ((integration_flag & EvaluationFlags::gradients) != 0u)
+        for (unsigned int c = 0; c < n_components; ++c)
+          {
+            switch (dim)
+              {
+                case 3:
+                  // grad z
+                  eval1.template values<1, false, false, 3>(gradients_quad + 2,
+                                                            scratch_data);
+                  eval0.template values<0, false, false>(scratch_data,
+                                                         values_dofs + n_dofs);
+                  if (symmetric_evaluate &&
+                      use_collocation_evaluation(fe_degree, n_q_points_1d))
+                    {
+                      EvaluatorTensorProduct<evaluate_evenodd,
+                                             dim - 1,
+                                             n_q_points_1d,
+                                             n_q_points_1d,
+                                             Number,
+                                             Number2>
+                        eval_grad({}, data.shape_gradients_collocation_eo, {});
+                      if ((integration_flag & EvaluationFlags::values) != 0u)
+                        eval_grad.template gradients<1, false, true, 3>(
+                          gradients_quad + 1, values_quad);
+                      else
+                        eval_grad.template gradients<1, false, false, 3>(
+                          gradients_quad + 1, values_quad);
+                      eval_grad.template gradients<0, false, true, 3>(
+                        gradients_quad, values_quad);
+                      eval0.template values<1, false, false>(values_quad,
+                                                             values_quad);
+                      eval0.template values<0, false, false>(values_quad,
+                                                             values_dofs);
+                    }
+                  else
+                    {
+                      if ((integration_flag & EvaluationFlags::values) != 0u)
+                        {
+                          eval1.template values<1, false, false>(values_quad,
+                                                                 scratch_data);
+                          eval1.template gradients<1, false, true, 3>(
+                            gradients_quad + 1, scratch_data);
+                        }
+                      else
+                        eval1.template gradients<1, false, false, 3>(
+                          gradients_quad + 1, scratch_data);
+
+                      // grad y
+                      eval0.template values<0, false, false>(scratch_data,
+                                                             values_dofs);
+
+                      // grad x
+                      eval1.template values<1, false, false, 3>(gradients_quad,
+                                                                scratch_data);
+                      eval0.template gradients<0, false, true>(scratch_data,
+                                                               values_dofs);
+                    }
+                  break;
+                case 2:
+                  eval0.template values<0, false, false, 2>(gradients_quad + 1,
+                                                            values_dofs +
+                                                              n_dofs);
+                  eval0.template gradients<0, false, false, 2>(gradients_quad,
+                                                               values_dofs);
+                  if ((integration_flag & EvaluationFlags::values) != 0u)
+                    eval0.template values<0, false, true>(values_quad,
+                                                          values_dofs);
+                  break;
+                case 1:
+                  values_dofs[0] = values_quad[0];
+                  values_dofs[1] = gradients_quad[0];
+                  break;
+                default:
+                  AssertThrow(false, ExcNotImplemented());
+              }
+            values_dofs += 3 * n_dofs;
+            values_quad += n_q_points;
+            gradients_quad += dim * n_q_points;
+          }
+
+      if ((integration_flag & EvaluationFlags::hessians) != 0u)
+        {
+          values_dofs = values_dofs_ptr;
+          for (unsigned int c = 0; c < n_components; ++c)
+            {
+              switch (dim)
+                {
+                  case 3:
+                    // grad xx
+                    eval1.template values<1, false, false>(hessians_quad,
+                                                           scratch_data);
+                    if ((integration_flag & (EvaluationFlags::values |
+                                             EvaluationFlags::gradients)) != 0u)
+                      eval0.template hessians<0, false, true>(scratch_data,
+                                                              values_dofs);
+                    else
+                      eval0.template hessians<0, false, false>(scratch_data,
+                                                               values_dofs);
+
+                    // grad yy
+                    eval1.template hessians<1, false, false>(hessians_quad +
+                                                               n_q_points,
+                                                             scratch_data);
+                    eval0.template values<0, false, true>(scratch_data,
+                                                          values_dofs);
+
+                    // grad zz
+                    eval1.template values<1, false, false>(hessians_quad +
+                                                             2 * n_q_points,
+                                                           scratch_data);
+                    eval0.template values<0, false, false>(scratch_data,
+                                                           values_dofs +
+                                                             2 * n_dofs);
+
+                    // grad xy
+                    eval1.template gradients<1, false, false>(hessians_quad +
+                                                                3 * n_q_points,
+                                                              scratch_data);
+                    eval0.template gradients<0, false, true>(scratch_data,
+                                                             values_dofs);
+
+                    // grad xz
+                    eval1.template values<1, false, false>(hessians_quad +
+                                                             4 * n_q_points,
+                                                           scratch_data);
+                    if ((integration_flag & EvaluationFlags::gradients) != 0u)
+                      eval0.template gradients<0, false, true>(scratch_data,
+                                                               values_dofs +
+                                                                 n_dofs);
+                    else
+                      eval0.template gradients<0, false, false>(scratch_data,
+                                                                values_dofs +
+                                                                  n_dofs);
+
+                    // grad yz
+                    eval1.template gradients<1, false, false>(hessians_quad +
+                                                                5 * n_q_points,
+                                                              scratch_data);
+                    eval0.template values<0, false, true>(scratch_data,
+                                                          values_dofs + n_dofs);
+
+                    break;
+                  case 2:
+                    // grad xx
+                    if ((integration_flag & (EvaluationFlags::values |
+                                             EvaluationFlags::gradients)) != 0u)
+                      eval0.template hessians<0, false, true>(hessians_quad,
+                                                              values_dofs);
+                    else
+                      eval0.template hessians<0, false, false>(hessians_quad,
+                                                               values_dofs);
+
+                    // grad yy
+                    eval0.template values<0, false, false>(
+                      hessians_quad + n_q_points, values_dofs + 2 * n_dofs);
+                    // grad xy
+                    if ((integration_flag & EvaluationFlags::gradients) != 0u)
+                      eval0.template gradients<0, false, true>(
+                        hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
+                    else
+                      eval0.template gradients<0, false, false>(
+                        hessians_quad + 2 * n_q_points, values_dofs + n_dofs);
+                    break;
+                  case 1:
+                    values_dofs[2] = hessians_quad[0];
+                    if ((integration_flag & EvaluationFlags::values) == 0u)
+                      values_dofs[0] = 0;
+                    if ((integration_flag & EvaluationFlags::gradients) == 0u)
+                      values_dofs[1] = 0;
+                    break;
+                  default:
+                    AssertThrow(false, ExcNotImplemented());
+                }
+              values_dofs += 3 * n_dofs;
+              hessians_quad += dim * (dim + 1) / 2 * n_q_points;
+            }
+        }
+    }
+  };
+
+
+
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  struct FEFaceEvaluationImplRaviartThomas
+  {
+    using Number2 =
+      typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
+
+    /**
+     * Apply the sum factorization kernels within the face for Raviart-Thomas
+     * elements for either evaluation or integration
+     */
+    template <bool do_integrate>
+    static inline void
+    evaluate_or_integrate_in_face(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      const std::vector<MatrixFreeFunctions::UnivariateShapeData<Number2>>
+                        &shape_data,
+      Number            *values_dofs_in,
+      Number            *values,
+      Number            *gradients,
+      Number            *scratch_data,
+      const unsigned int subface_index,
+      const unsigned int face_direction)
+    {
+      AssertDimension(shape_data.size(), 2);
+
+      const int degree = fe_degree != -1 ? fe_degree : shape_data[0].fe_degree;
+      const int n_rows_n = degree + 1;
+      const int n_rows_t = degree;
+      const dealii::ndarray<int, 3, 3> dofs_per_direction{
+        {{{n_rows_n, n_rows_t, n_rows_t}},
+         {{n_rows_t, n_rows_n, n_rows_t}},
+         {{n_rows_t, n_rows_t, n_rows_n}}}};
+
+      (void)scratch_data;
+      (void)subface_index;
+      // TODO: This is currently not implemented, but the test
+      // matrix_vector_rt_face_03 apparently works without it -> check
+      // if (subface_index < GeometryInfo<dim - 1>::max_children_per_cell)
+      //  Assert(false, ExcNotImplemented());
+
+      using Eval = EvaluatorTensorProduct<evaluate_evenodd,
+                                          dim - 1,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          Number,
+                                          Number2>;
+
+      std::array<int, dim> values_dofs_offsets = {};
+      for (unsigned int comp = 0; comp < dim - 1; ++comp)
+        {
+          if (dim == 2)
+            values_dofs_offsets[comp + 1] =
+              values_dofs_offsets[comp] +
+              3 * dofs_per_direction[comp][(face_direction + 1) % dim];
+          else
+            values_dofs_offsets[comp + 1] =
+              values_dofs_offsets[comp] +
+              3 * dofs_per_direction[comp][(face_direction + 1) % dim] *
+                dofs_per_direction[comp][(face_direction + 2) % dim];
+        }
+
+      // Jacobians on faces are reordered to enable simple access with the
+      // regular evaluators; to get the RT Piola transform right, we need to
+      // pass through the values_dofs array in a permuted right order
+      std::array<unsigned int, dim> components;
+      for (unsigned int comp = 0; comp < dim; ++comp)
+        components[comp] = (face_direction + comp + 1) % dim;
+
+      for (const unsigned int comp : components)
+        {
+          Number *values_dofs = values_dofs_in + values_dofs_offsets[comp];
+
+          std::array<int, 2> n_blocks{
+            {dofs_per_direction[comp][(face_direction + 1) % dim],
+             (dim > 2 ? dofs_per_direction[comp][(face_direction + 2) % dim] :
+                        1)}};
+
+          if constexpr (dim == 3)
+            {
+              EvaluatorTensorProduct<evaluate_evenodd,
+                                     dim - 1,
+                                     n_q_points_1d,
+                                     n_q_points_1d,
+                                     Number,
+                                     Number2>
+                eval_g({},
+                       shape_data[0].shape_gradients_collocation_eo.data(),
+                       {});
+              if (!do_integrate)
+                {
+                  EvaluatorTensorProductAnisotropic<dim - 1,
+                                                    fe_degree,
+                                                    n_q_points_1d,
+                                                    true>
+                    eval;
+                  // Evaluate in 3d
+                  if (n_blocks[0] == n_rows_n)
+                    {
+                      eval.template normal<0>(shape_data[0],
+                                              values_dofs,
+                                              values);
+                      eval.template tangential<1, 0>(shape_data[1],
+                                                     values,
+                                                     values);
+
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template normal<0>(shape_data[0],
+                                                  values_dofs +
+                                                    n_blocks[0] * n_blocks[1],
+                                                  scratch_data);
+                          eval.template tangential<1, 0, dim>(shape_data[1],
+                                                              scratch_data,
+                                                              gradients + 2);
+                        }
+                    }
+                  else if (n_blocks[1] == n_rows_n)
+                    {
+                      eval.template normal<1>(shape_data[0],
+                                              values_dofs,
+                                              values);
+                      eval.template tangential<0, 1>(shape_data[1],
+                                                     values,
+                                                     values);
+
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template normal<1>(shape_data[0],
+                                                  values_dofs +
+                                                    n_blocks[0] * n_blocks[1],
+                                                  scratch_data);
+                          eval.template tangential<0, 1, dim>(shape_data[1],
+                                                              scratch_data,
+                                                              gradients + 2);
+                        }
+                    }
+                  else
+                    {
+                      Eval eval(shape_data[1].shape_values_eo.data(), {}, {});
+                      eval.template values<0, true, false>(values_dofs, values);
+                      eval.template values<1, true, false>(values, values);
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template values<0, true, false>(values_dofs +
+                                                                 n_blocks[0] *
+                                                                   n_blocks[1],
+                                                               scratch_data);
+                          eval.template values<1, true, false, dim>(
+                            scratch_data, gradients + 2);
+                        }
+                    }
+                  if (evaluation_flag & EvaluationFlags::gradients)
+                    {
+                      eval_g.template gradients<0, true, false, dim>(values,
+                                                                     gradients);
+                      eval_g.template gradients<1, true, false, dim>(values,
+                                                                     gradients +
+                                                                       1);
+                    }
+                }
+              else
+                {
+                  EvaluatorTensorProductAnisotropic<dim - 1,
+                                                    fe_degree,
+                                                    n_q_points_1d,
+                                                    false>
+                    eval;
+                  // Integrate in 3d
+                  if (evaluation_flag & EvaluationFlags::gradients)
+                    {
+                      if (evaluation_flag & EvaluationFlags::values)
+                        eval_g.template gradients<0, false, true, dim>(
+                          gradients, values);
+                      else
+                        eval_g.template gradients<0, false, false, dim>(
+                          gradients, values);
+                      eval_g.template gradients<1, false, true, dim>(gradients +
+                                                                       1,
+                                                                     values);
+                    }
+                  if (n_blocks[0] == n_rows_n)
+                    {
+                      eval.template tangential<1, 0>(shape_data[1],
+                                                     values,
+                                                     values);
+                      eval.template normal<0>(shape_data[0],
+                                              values,
+                                              values_dofs);
+
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template tangential<1, 0, dim>(shape_data[1],
+                                                              gradients + 2,
+                                                              scratch_data);
+                          eval.template normal<0>(shape_data[0],
+                                                  scratch_data,
+                                                  values_dofs +
+                                                    n_blocks[0] * n_blocks[1]);
+                        }
+                    }
+                  else if (n_blocks[1] == n_rows_n)
+                    {
+                      eval.template tangential<0, 1>(shape_data[1],
+                                                     values,
+                                                     values);
+                      eval.template normal<1>(shape_data[0],
+                                              values,
+                                              values_dofs);
+
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template tangential<0, 1, dim>(shape_data[1],
+                                                              gradients + 2,
+                                                              scratch_data);
+                          eval.template normal<1>(shape_data[0],
+                                                  scratch_data,
+                                                  values_dofs +
+                                                    n_blocks[0] * n_blocks[1]);
+                        }
+                    }
+                  else
+                    {
+                      Eval eval_iso(shape_data[1].shape_values_eo.data(),
+                                    {},
+                                    {});
+                      eval_iso.template values<1, false, false>(values, values);
+                      eval_iso.template values<0, false, false>(values,
+                                                                values_dofs);
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval_iso.template values<1, false, false, dim>(
+                            gradients + 2, scratch_data);
+                          eval_iso.template values<0, false, false>(
+                            scratch_data,
+                            values_dofs + n_blocks[0] * n_blocks[1]);
+                        }
+                    }
+                }
+            }
+          else
+            {
+              using EvalN = EvaluatorTensorProduct<evaluate_evenodd,
+                                                   dim - 1,
+                                                   fe_degree + 1,
+                                                   n_q_points_1d,
+                                                   Number,
+                                                   Number2>;
+              if (!do_integrate)
+                {
+                  // Evaluate in 2d
+                  if (n_blocks[0] == n_rows_n)
+                    {
+                      EvalN eval(shape_data[0].shape_values_eo,
+                                 shape_data[0].shape_gradients_eo,
+                                 {});
+                      eval.template values<0, true, false>(values_dofs, values);
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template gradients<0, true, false, dim>(
+                            values_dofs, gradients);
+                          eval.template values<0, true, false, dim>(
+                            values_dofs + n_rows_n, gradients + 1);
+                        }
+                    }
+                  else
+                    {
+                      Eval eval(shape_data[1].shape_values_eo,
+                                shape_data[1].shape_gradients_eo,
+                                {});
+                      eval.template values<0, true, false>(values_dofs, values);
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          eval.template gradients<0, true, false, dim>(
+                            values_dofs, gradients);
+                          eval.template values<0, true, false, dim>(
+                            values_dofs + n_rows_t, gradients + 1);
+                        }
+                    }
+                }
+              else
+                {
+                  // Integrate in 2d
+                  if (n_blocks[0] == n_rows_n)
+                    {
+                      EvalN eval(shape_data[0].shape_values_eo,
+                                 shape_data[0].shape_gradients_eo,
+                                 {});
+                      if (evaluation_flag & EvaluationFlags::values)
+                        eval.template values<0, false, false>(values,
+                                                              values_dofs);
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          if (evaluation_flag & EvaluationFlags::values)
+                            eval.template gradients<0, false, true, dim>(
+                              gradients, values_dofs);
+                          else
+                            eval.template gradients<0, false, false, dim>(
+                              gradients, values_dofs);
+                          eval.template values<0, false, false, dim>(
+                            gradients + 1, values_dofs + n_rows_n);
+                        }
+                    }
+                  else
+                    {
+                      Eval eval(shape_data[1].shape_values_eo,
+                                shape_data[1].shape_gradients_eo,
+                                {});
+                      if (evaluation_flag & EvaluationFlags::values)
+                        eval.template values<0, false, false>(values,
+                                                              values_dofs);
+                      if (evaluation_flag & EvaluationFlags::gradients)
+                        {
+                          if (evaluation_flag & EvaluationFlags::values)
+                            eval.template gradients<0, false, true, dim>(
+                              gradients, values_dofs);
+                          else
+                            eval.template gradients<0, false, false, dim>(
+                              gradients, values_dofs);
+                          eval.template values<0, false, false, dim>(
+                            gradients + 1, values_dofs + n_rows_t);
+                        }
+                    }
+                }
+            }
+          values += Utilities::pow(n_q_points_1d, dim - 1);
+          gradients += dim * Utilities::pow(n_q_points_1d, dim - 1);
+        }
+    }
+  };
+
+
+
+  template <int dim, int fe_degree, typename Number>
+  struct FEFaceNormalEvaluationImpl
+  {
+    using Number2 =
+      typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
+
+    template <bool do_evaluate, bool add_into_output>
+    static void
+    interpolate(const unsigned int                             n_components,
+                const EvaluationFlags::EvaluationFlags         flags,
+                const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
+                const Number                                  *input,
+                Number                                        *output,
+                const unsigned int                             face_no)
+    {
+      Assert(static_cast<unsigned int>(fe_degree) ==
+                 shape_info.data.front().fe_degree ||
+               fe_degree == -1,
+             ExcInternalError());
+      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+        interpolate_raviart_thomas<do_evaluate, add_into_output>(
+          n_components, input, output, flags, face_no, shape_info);
+      else
+        interpolate_generic<do_evaluate, add_into_output>(
+          n_components,
+          input,
+          output,
+          flags,
+          face_no,
+          shape_info.data.front().fe_degree + 1,
+          shape_info.data.front().shape_data_on_face,
+          shape_info.dofs_per_component_on_cell,
+          3 * shape_info.dofs_per_component_on_face);
+    }
+
+    /**
+     * Interpolate the values on the cell quadrature points onto a face.
+     */
+    template <bool do_evaluate, bool add_into_output>
+    static void
+    interpolate_quadrature(
+      const unsigned int                             n_components,
+      const EvaluationFlags::EvaluationFlags         flags,
+      const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info,
+      const Number                                  *input,
+      Number                                        *output,
+      const unsigned int                             face_no)
+    {
+      Assert(static_cast<unsigned int>(fe_degree + 1) ==
+                 shape_info.data.front().n_q_points_1d ||
+               fe_degree == -1,
+             ExcInternalError());
+
+      interpolate_generic<do_evaluate, add_into_output>(
+        n_components,
+        input,
+        output,
+        flags,
+        face_no,
+        shape_info.data.front().quadrature.size(),
+        shape_info.data.front().quadrature_data_on_face,
+        shape_info.n_q_points,
+        shape_info.n_q_points_face);
+    }
+
+  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 EvaluationFlags::EvaluationFlags flag,
+                        const unsigned int                     face_no,
+                        const unsigned int                     n_points_1d,
+                        const std::array<AlignedVector<Number2>, 2> &shape_data,
+                        const unsigned int dofs_per_component_on_cell,
+                        const unsigned int dofs_per_component_on_face)
+    {
+      if (face_direction == face_no / 2)
+        {
+          constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
+
+          const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
+          const int stride = Utilities::pow(n_rows, face_direction);
+          const std::array<int, 2> n_blocks{
+            {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
+          std::array<int, 2> steps;
+          if constexpr (face_direction == 0)
+            steps = {{n_rows, 0}};
+          else if constexpr (face_direction == 1 && dim == 2)
+            steps = {{1, 0}};
+          else if constexpr (face_direction == 1)
+            // in 3d, the coordinate system is zx, not xz -> switch indices
+            steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}};
+          else if constexpr (face_direction == 2)
+            steps = {{1, 0}};
+
+          for (unsigned int c = 0; c < n_components; ++c)
+            {
+              if (flag & EvaluationFlags::hessians)
+                interpolate_to_face<fe_degree + 1,
+                                    stride_,
+                                    do_evaluate,
+                                    add_into_output,
+                                    2>(shape_data[face_no % 2].begin(),
+                                       n_blocks,
+                                       steps,
+                                       input,
+                                       output,
+                                       n_rows,
+                                       stride);
+              else if (flag & EvaluationFlags::gradients)
+                interpolate_to_face<fe_degree + 1,
+                                    stride_,
+                                    do_evaluate,
+                                    add_into_output,
+                                    1>(shape_data[face_no % 2].begin(),
+                                       n_blocks,
+                                       steps,
+                                       input,
+                                       output,
+                                       n_rows,
+                                       stride);
+              else
+                interpolate_to_face<fe_degree + 1,
+                                    stride_,
+                                    do_evaluate,
+                                    add_into_output,
+                                    0>(shape_data[face_no % 2].begin(),
+                                       n_blocks,
+                                       steps,
+                                       input,
+                                       output,
+                                       n_rows,
+                                       stride);
+              if (do_evaluate)
+                {
+                  input += dofs_per_component_on_cell;
+                  output += dofs_per_component_on_face;
+                }
+              else
+                {
+                  output += dofs_per_component_on_cell;
+                  input += dofs_per_component_on_face;
+                }
+            }
+        }
+      else if (face_direction < dim)
+        {
+          interpolate_generic<do_evaluate,
+                              add_into_output,
+                              std::min(face_direction + 1, dim - 1)>(
+            n_components,
+            input,
+            output,
+            flag,
+            face_no,
+            n_points_1d,
+            shape_data,
+            dofs_per_component_on_cell,
+            dofs_per_component_on_face);
+        }
+    }
+
+    template <bool do_evaluate,
+              bool add_into_output,
+              int  face_direction = 0,
+              int  max_derivative = 0>
+    static void
+    interpolate_raviart_thomas(
+      const unsigned int                             n_components,
+      const Number                                  *input,
+      Number                                        *output,
+      const EvaluationFlags::EvaluationFlags         flag,
+      const unsigned int                             face_no,
+      const MatrixFreeFunctions::ShapeInfo<Number2> &shape_info)
+    {
+      if (dim == 1)
+        {
+          // This should never happen since the FE_RaviartThomasNodal is not
+          // defined for dim = 1. It prevents compiler warnings of infinite
+          // recursion.
+          Assert(false, ExcInternalError());
+          return;
+        }
+
+      bool increase_max_der = false;
+      if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
+          (flag & EvaluationFlags::gradients && max_derivative < 1))
+        increase_max_der = true;
+
+      if (face_direction == face_no / 2 && !increase_max_der)
+        {
+          constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
+          constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
+          constexpr int stride2 = fe_degree * (fe_degree + 1);
+
+          const int degree =
+            fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
+          const int n_rows_n = degree + 1;
+          const int n_rows_t = degree;
+
+          std::array<int, 3> strides{{1, 1, 1}};
+          if (face_direction > 0)
+            {
+              strides[0] =
+                n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
+              strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
+              strides[2] = Utilities::pow(n_rows_t, face_direction);
+            }
+          const dealii::ndarray<int, 3, 3> dofs_per_direction{
+            {{{n_rows_n, n_rows_t, n_rows_t}},
+             {{n_rows_t, n_rows_n, n_rows_t}},
+             {{n_rows_t, n_rows_t, n_rows_n}}}};
+
+          std::array<int, 2> steps, n_blocks;
+
+          if constexpr (face_direction == 0)
+            steps = {{degree + (face_direction == 0), 0}};
+          else if constexpr (face_direction == 1 && dim == 2)
+            steps = {{1, 0}};
+          else if constexpr (face_direction == 1)
+            // in 3d, the coordinate system is zx, not xz -> switch indices
+            steps = {
+              {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
+          else if constexpr (face_direction == 2)
+            steps = {{1, 0}};
+
+          n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
+          n_blocks[1] =
+            dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
+
+          interpolate_to_face<
+            (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
+            ((face_direction < 2) ? stride1 : stride2),
+            do_evaluate,
+            add_into_output,
+            max_derivative>(shape_info.data[face_direction != 0]
+                              .shape_data_on_face[face_no % 2]
+                              .begin(),
+                            n_blocks,
+                            steps,
+                            input,
+                            output,
+                            degree + (face_direction == 0),
+                            strides[0]);
+
+          if (do_evaluate)
+            {
+              input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+              output += 3 * n_blocks[0] * n_blocks[1];
+            }
+          else
+            {
+              output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+              input += 3 * n_blocks[0] * n_blocks[1];
+            }
+
+          // must only change steps only for face direction 0
+          if constexpr (face_direction == 0)
+            steps = {{degree, 0}};
+
+          n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
+          n_blocks[1] =
+            dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
+
+          interpolate_to_face<
+            (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
+            ((face_direction < 2) ? stride0 : stride2),
+            do_evaluate,
+            add_into_output,
+            max_derivative>(shape_info.data[face_direction != 1]
+                              .shape_data_on_face[face_no % 2]
+                              .begin(),
+                            n_blocks,
+                            steps,
+                            input,
+                            output,
+                            degree + (face_direction == 1),
+                            strides[1]);
+
+          if constexpr (dim > 2)
+            {
+              if (do_evaluate)
+                {
+                  input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+                  output += 3 * n_blocks[0] * n_blocks[1];
+                }
+              else
+                {
+                  output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+                  input += 3 * n_blocks[0] * n_blocks[1];
+                }
+
+              if constexpr (face_direction == 0)
+                steps = {{degree, 0}};
+              else if constexpr (face_direction == 1)
+                // in 3d, the coordinate system is zx, not xz -> switch indices
+                steps = {
+                  {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
+              else if constexpr (face_direction == 2)
+                steps = {{1, 0}};
+
+              n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
+              n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
+
+              interpolate_to_face<
+                (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
+                stride0,
+                do_evaluate,
+                add_into_output,
+                max_derivative>(shape_info.data[face_direction != 2]
+                                  .shape_data_on_face[face_no % 2]
+                                  .begin(),
+                                n_blocks,
+                                steps,
+                                input,
+                                output,
+                                degree + (face_direction == 2),
+                                strides[2]);
+            }
+        }
+      else if (face_direction == face_no / 2)
+        {
+          // Only increase max_derivative
+          interpolate_raviart_thomas<do_evaluate,
+                                     add_into_output,
+                                     face_direction,
+                                     std::min(max_derivative + 1, 2)>(
+            n_components, input, output, flag, face_no, shape_info);
+        }
+      else if (face_direction < dim)
+        {
+          if (increase_max_der)
+            {
+              interpolate_raviart_thomas<do_evaluate,
+                                         add_into_output,
+                                         std::min(face_direction + 1, dim - 1),
+                                         std::min(max_derivative + 1, 2)>(
+                n_components, input, output, flag, face_no, shape_info);
+            }
+          else
+            {
+              interpolate_raviart_thomas<do_evaluate,
+                                         add_into_output,
+                                         std::min(face_direction + 1, dim - 1),
+                                         max_derivative>(
+                n_components, input, output, flag, face_no, shape_info);
+            }
+        }
+    }
+  };
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+      dst[v] = src_ptr[v];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated load function
+  template <typename Number, std::size_t width>
+  void
+  do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
+  {
+    dst.load(src_ptr);
+  }
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_gather(const Number2       *src_ptr,
+                       const unsigned int  *indices,
+                       VectorizedArrayType &dst)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+      dst[v] = src_ptr[indices[v]];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated gather function
+  template <typename Number, std::size_t width>
+  void
+  do_vectorized_gather(const Number                   *src_ptr,
+                       const unsigned int             *indices,
+                       VectorizedArray<Number, width> &dst)
+  {
+    dst.gather(src_ptr, indices);
+  }
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+      dst_ptr[v] += src[v];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated load function
+  template <typename Number, std::size_t width>
+  void
+  do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
+  {
+    VectorizedArray<Number, width> tmp;
+    tmp.load(dst_ptr);
+    (tmp + src).store(dst_ptr);
+  }
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_scatter_add(const VectorizedArrayType src,
+                            const unsigned int       *indices,
+                            Number2                  *dst_ptr)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+      dst_ptr[indices[v]] += src[v];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated gather function
+  template <typename Number, std::size_t width>
+  void
+  do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
+                            const unsigned int                  *indices,
+                            Number                              *dst_ptr)
+  {
+#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
+    for (unsigned int v = 0; v < width; ++v)
+      dst_ptr[indices[v]] += src[v];
+#else
+    VectorizedArray<Number, width> tmp;
+    tmp.gather(dst_ptr, indices);
+    (tmp + src).scatter(indices, dst_ptr);
+#endif
+  }
+
+
+
+  template <typename Number>
+  void
+  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)
+  {
+    for (unsigned int c = 0; c < n_components; ++c)
+      {
+        if (flag & EvaluationFlags::values)
+          {
+            if (integrate)
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
+            else
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              values_quad[c * n_q_points + q] = tmp_values[q];
+          }
+        if (flag & EvaluationFlags::gradients)
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              if (integrate)
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  tmp_values[q] =
+                    gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
+              else
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  tmp_values[orientation[q]] =
+                    gradients_quad[(c * n_q_points + q) * dim + d];
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
+            }
+        if (flag & EvaluationFlags::hessians)
+          {
+            const unsigned int hdim = (dim * (dim + 1)) / 2;
+            for (unsigned int d = 0; d < hdim; ++d)
+              {
+                if (integrate)
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
+                                                  orientation[q]];
+                else
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[orientation[q]] =
+                      hessians_quad[(c * hdim + d) * n_q_points + q];
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  hessians_quad[(c * hdim + d) * n_q_points + q] =
+                    tmp_values[q];
+              }
+          }
+      }
+  }
+
+
+
+  template <typename Number, typename VectorizedArrayType>
+  void
+  adjust_for_face_orientation_per_lane(
+    const unsigned int                     dim,
+    const unsigned int                     n_components,
+    const unsigned int                     v,
+    const EvaluationFlags::EvaluationFlags flag,
+    const unsigned int                    *orientation,
+    const bool                             integrate,
+    const std::size_t                      n_q_points,
+    Number                                *tmp_values,
+    VectorizedArrayType                   *values_quad,
+    VectorizedArrayType                   *gradients_quad = nullptr,
+    VectorizedArrayType                   *hessians_quad  = nullptr)
+  {
+    for (unsigned int c = 0; c < n_components; ++c)
+      {
+        if (flag & EvaluationFlags::values)
+          {
+            if (integrate)
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                tmp_values[q] = values_quad[c * n_q_points + orientation[q]][v];
+            else
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                tmp_values[orientation[q]] = values_quad[c * n_q_points + q][v];
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              values_quad[c * n_q_points + q][v] = tmp_values[q];
+          }
+        if (flag & EvaluationFlags::gradients)
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              Assert(gradients_quad != nullptr, ExcInternalError());
+              if (integrate)
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  tmp_values[q] =
+                    gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
+                                  [v];
+              else
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  tmp_values[orientation[q]] =
+                    gradients_quad[(c * n_q_points + q) * dim + d][v];
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                gradients_quad[(c * n_q_points + q) * dim + d][v] =
+                  tmp_values[q];
+            }
+        if (flag & EvaluationFlags::hessians)
+          {
+            Assert(hessians_quad != nullptr, ExcInternalError());
+            const unsigned int hdim = (dim * (dim + 1)) / 2;
+            for (unsigned int d = 0; d < hdim; ++d)
+              {
+                if (integrate)
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[q] = hessians_quad[(c * hdim + d) * n_q_points +
+                                                  orientation[q]][v];
+                else
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[orientation[q]] =
+                      hessians_quad[(c * hdim + d) * n_q_points + q][v];
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  hessians_quad[(c * hdim + d) * n_q_points + q][v] =
+                    tmp_values[q];
+              }
+          }
+      }
+  }
+
+
+
+  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 EvaluationFlags::EvaluationFlags evaluation_flag,
+        const Number                          *values_dofs,
+        FEEvaluationData<dim, Number, true>   &fe_eval)
+    {
+      const auto &shape_info = fe_eval.get_shape_info();
+      const auto &shape_data = shape_info.data.front();
+      using Number2 =
+        typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
+
+      if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
+        {
+          Assert((fe_eval.get_dof_access_index() ==
+                    MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+                  fe_eval.is_interior_face() == false) == false,
+                 ExcNotImplemented());
+
+          const unsigned int face_no          = fe_eval.get_face_no();
+          const unsigned int face_orientation = fe_eval.get_face_orientation();
+          const std::size_t  n_dofs     = shape_info.dofs_per_component_on_cell;
+          const std::size_t  n_q_points = shape_info.n_q_points_faces[face_no];
+
+          using Eval =
+            EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
+
+          if (evaluation_flag & EvaluationFlags::values)
+            {
+              const auto *const shape_values =
+                &shape_data.shape_values_face(face_no, face_orientation, 0);
+
+              auto *values_quad_ptr        = fe_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)
+                {
+                  eval.template values<0, true, false>(values_dofs_actual_ptr,
+                                                       values_quad_ptr);
+
+                  values_quad_ptr += n_q_points;
+                  values_dofs_actual_ptr += n_dofs;
+                }
+            }
+
+          if (evaluation_flag & EvaluationFlags::gradients)
+            {
+              auto       *gradients_quad_ptr     = fe_eval.begin_gradients();
+              const auto *values_dofs_actual_ptr = values_dofs;
+
+              std::array<const Number2 *, dim> shape_gradients;
+              for (unsigned int d = 0; d < dim; ++d)
+                shape_gradients[d] = &shape_data.shape_gradients_face(
+                  face_no, face_orientation, d, 0);
+
+              for (unsigned int c = 0; c < n_components; ++c)
+                {
+                  for (unsigned int d = 0; d < dim; ++d)
+                    {
+                      Eval eval(nullptr,
+                                shape_gradients[d],
+                                nullptr,
+                                n_dofs,
+                                n_q_points);
+
+                      eval.template gradients<0, true, false, dim>(
+                        values_dofs_actual_ptr, gradients_quad_ptr + d);
+                    }
+                  gradients_quad_ptr += n_q_points * dim;
+                  values_dofs_actual_ptr += n_dofs;
+                }
+            }
+
+          Assert(!(evaluation_flag & EvaluationFlags::hessians),
+                 ExcNotImplemented());
+
+          return true;
+        }
+
+      const unsigned int dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         Utilities::pow(shape_data.fe_degree + 1, dim - 1);
+
+      // Note: we always keep storage of values, 1st and 2nd derivatives in an
+      // array, so reserve space for all three here
+      Number *temp         = fe_eval.get_scratch_data().begin();
+      Number *scratch_data = temp + 3 * n_components * dofs_per_face;
+
+      bool use_vectorization = true;
+
+      if (fe_eval.get_dof_access_index() ==
+            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
+        for (unsigned int v = 0; v < Number::size(); ++v)
+          if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
+              fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
+            use_vectorization = false;
+
+      if (use_vectorization == false)
+        {
+          for (unsigned int v = 0; v < Number::size(); ++v)
+            {
+              // the loop breaks once an invalid_unsigned_int is hit for
+              // all cases except the exterior faces in the ECL loop (where
+              // some faces might be at the boundaries but others not)
+              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                {
+                  for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
+                       ++i)
+                    temp[i][v] = 0;
+                  continue;
+                }
+
+              FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+                template interpolate<true, false>(n_components,
+                                                  evaluation_flag,
+                                                  shape_info,
+                                                  values_dofs,
+                                                  scratch_data,
+                                                  fe_eval.get_face_no(v));
+
+              for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
+                   ++i)
+                temp[i][v] = scratch_data[i][v];
+            }
+        }
+      else
+        FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+          template interpolate<true, false>(n_components,
+                                            evaluation_flag,
+                                            shape_info,
+                                            values_dofs,
+                                            temp,
+                                            fe_eval.get_face_no());
+
+      const unsigned int     subface_index = fe_eval.get_subface_index();
+      constexpr unsigned int n_q_points_1d_actual =
+        fe_degree > -1 ? n_q_points_1d : 0;
+
+      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+        {
+          FEFaceEvaluationImplRaviartThomas<dim,
+                                            fe_degree,
+                                            n_q_points_1d_actual,
+                                            Number>::
+            template evaluate_or_integrate_in_face<false>(
+              evaluation_flag,
+              fe_eval.get_shape_info().data,
+              temp,
+              fe_eval.begin_values(),
+              fe_eval.begin_gradients(),
+              scratch_data,
+              subface_index,
+              fe_eval.get_face_no() / 2);
+        }
+      else if (fe_degree > -1 &&
+               subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+               shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
+        FEFaceEvaluationImpl<true,
+                             dim,
+                             fe_degree,
+                             n_q_points_1d_actual,
+                             Number>::evaluate_in_face(n_components,
+                                                       evaluation_flag,
+                                                       shape_data,
+                                                       temp,
+                                                       fe_eval.begin_values(),
+                                                       fe_eval
+                                                         .begin_gradients(),
+                                                       fe_eval.begin_hessians(),
+                                                       scratch_data,
+                                                       subface_index);
+      else
+        FEFaceEvaluationImpl<false,
+                             dim,
+                             fe_degree,
+                             n_q_points_1d_actual,
+                             Number>::evaluate_in_face(n_components,
+                                                       evaluation_flag,
+                                                       shape_data,
+                                                       temp,
+                                                       fe_eval.begin_values(),
+                                                       fe_eval
+                                                         .begin_gradients(),
+                                                       fe_eval.begin_hessians(),
+                                                       scratch_data,
+                                                       subface_index);
+
+      if (use_vectorization == false)
+        {
+          for (unsigned int v = 0; v < Number::size(); ++v)
+            {
+              // the loop breaks once an invalid_unsigned_int is hit for
+              // all cases except the exterior faces in the ECL loop (where
+              // some faces might be at the boundaries but others not)
+              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                continue;
+
+              if (fe_eval.get_face_orientation(v) != 0)
+                adjust_for_face_orientation_per_lane(
+                  dim,
+                  n_components,
+                  v,
+                  evaluation_flag,
+                  &fe_eval.get_shape_info().face_orientations_quad(
+                    fe_eval.get_face_orientation(v), 0),
+                  false,
+                  shape_info.n_q_points_face,
+                  &temp[0][0],
+                  fe_eval.begin_values(),
+                  fe_eval.begin_gradients(),
+                  fe_eval.begin_hessians());
+            }
+        }
+      else if (fe_eval.get_face_orientation() != 0)
+        adjust_for_face_orientation(
+          dim,
+          n_components,
+          evaluation_flag,
+          &fe_eval.get_shape_info().face_orientations_quad(
+            fe_eval.get_face_orientation(), 0),
+          false,
+          shape_info.n_q_points_face,
+          temp,
+          fe_eval.begin_values(),
+          fe_eval.begin_gradients(),
+          fe_eval.begin_hessians());
+
+      return false;
+    }
+  };
+
+
+
+  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 EvaluationFlags::EvaluationFlags integration_flag,
+        Number                                *values_dofs,
+        FEEvaluationData<dim, Number, true>   &fe_eval)
+    {
+      const auto &shape_info = fe_eval.get_shape_info();
+      const auto &shape_data = shape_info.data.front();
+      using Number2 =
+        typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
+
+      if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
+        {
+          Assert((fe_eval.get_dof_access_index() ==
+                    MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+                  fe_eval.is_interior_face() == false) == false,
+                 ExcNotImplemented());
+
+          const unsigned int face_no          = fe_eval.get_face_no();
+          const unsigned int face_orientation = fe_eval.get_face_orientation();
+          const std::size_t  n_dofs     = shape_info.dofs_per_component_on_cell;
+          const std::size_t  n_q_points = shape_info.n_q_points_faces[face_no];
+
+          using Eval =
+            EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
+
+          if (integration_flag & EvaluationFlags::values)
+            {
+              const auto *const shape_values =
+                &shape_data.shape_values_face(face_no, face_orientation, 0);
+
+              auto *values_quad_ptr        = fe_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)
+                {
+                  eval.template values<0, false, false>(values_quad_ptr,
+                                                        values_dofs_actual_ptr);
+
+                  values_quad_ptr += n_q_points;
+                  values_dofs_actual_ptr += n_dofs;
+                }
+            }
+
+          if (integration_flag & EvaluationFlags::gradients)
+            {
+              auto *gradients_quad_ptr     = fe_eval.begin_gradients();
+              auto *values_dofs_actual_ptr = values_dofs;
+
+              std::array<const Number2 *, dim> shape_gradients;
+              for (unsigned int d = 0; d < dim; ++d)
+                shape_gradients[d] = &shape_data.shape_gradients_face(
+                  face_no, face_orientation, d, 0);
+
+              for (unsigned int c = 0; c < n_components; ++c)
+                {
+                  for (unsigned int d = 0; d < dim; ++d)
+                    {
+                      Eval eval(nullptr,
+                                shape_gradients[d],
+                                nullptr,
+                                n_dofs,
+                                n_q_points);
+
+                      if (!(integration_flag & EvaluationFlags::values) &&
+                          d == 0)
+                        eval.template gradients<0, false, false, dim>(
+                          gradients_quad_ptr + d, values_dofs_actual_ptr);
+                      else
+                        eval.template gradients<0, false, true, dim>(
+                          gradients_quad_ptr + d, values_dofs_actual_ptr);
+                    }
+                  gradients_quad_ptr += n_q_points * dim;
+                  values_dofs_actual_ptr += n_dofs;
+                }
+            }
+
+          Assert(!(integration_flag & EvaluationFlags::hessians),
+                 ExcNotImplemented());
+
+          return true;
+        }
+
+      const unsigned int dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         Utilities::pow(shape_data.fe_degree + 1, dim - 1);
+
+      Number *temp         = fe_eval.get_scratch_data().begin();
+      Number *scratch_data = temp + 3 * n_components * dofs_per_face;
+
+      bool use_vectorization = true;
+
+      if (fe_eval.get_dof_access_index() ==
+            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
+        use_vectorization =
+          fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
+          std::all_of(fe_eval.get_cell_ids().begin() + 1,
+                      fe_eval.get_cell_ids().end(),
+                      [&](const auto &v) {
+                        return v == fe_eval.get_cell_ids()[0] ||
+                               v == numbers::invalid_unsigned_int;
+                      });
+
+      if (use_vectorization == false)
+        {
+          for (unsigned int v = 0; v < Number::size(); ++v)
+            {
+              // the loop breaks once an invalid_unsigned_int is hit for
+              // all cases except the exterior faces in the ECL loop (where
+              // some faces might be at the boundaries but others not)
+              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                continue;
+
+              if (fe_eval.get_face_orientation(v) != 0)
+                adjust_for_face_orientation_per_lane(
+                  dim,
+                  n_components,
+                  v,
+                  integration_flag,
+                  &fe_eval.get_shape_info().face_orientations_quad(
+                    fe_eval.get_face_orientation(v), 0),
+                  true,
+                  shape_info.n_q_points_face,
+                  &temp[0][0],
+                  fe_eval.begin_values(),
+                  fe_eval.begin_gradients(),
+                  fe_eval.begin_hessians());
+            }
+        }
+      else if (fe_eval.get_face_orientation() != 0)
+        adjust_for_face_orientation(
+          dim,
+          n_components,
+          integration_flag,
+          &fe_eval.get_shape_info().face_orientations_quad(
+            fe_eval.get_face_orientation(), 0),
+          true,
+          shape_info.n_q_points_face,
+          temp,
+          fe_eval.begin_values(),
+          fe_eval.begin_gradients(),
+          fe_eval.begin_hessians());
+
+      const unsigned int n_q_points_1d_actual =
+        fe_degree > -1 ? n_q_points_1d : 0;
+      const unsigned int subface_index = fe_eval.get_subface_index();
+
+      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+        {
+          FEFaceEvaluationImplRaviartThomas<dim,
+                                            fe_degree,
+                                            n_q_points_1d_actual,
+                                            Number>::
+            template evaluate_or_integrate_in_face<true>(
+              integration_flag,
+              fe_eval.get_shape_info().data,
+              temp,
+              fe_eval.begin_values(),
+              fe_eval.begin_gradients(),
+              scratch_data,
+              subface_index,
+              fe_eval.get_face_no() / 2);
+        }
+      else if (fe_degree > -1 &&
+               fe_eval.get_subface_index() >=
+                 GeometryInfo<dim - 1>::max_children_per_cell &&
+               shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
+        FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d_actual,
+          Number>::integrate_in_face(n_components,
+                                     integration_flag,
+                                     shape_data,
+                                     temp,
+                                     fe_eval.begin_values(),
+                                     fe_eval.begin_gradients(),
+                                     fe_eval.begin_hessians(),
+                                     scratch_data,
+                                     subface_index);
+      else
+        FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d_actual,
+          Number>::integrate_in_face(n_components,
+                                     integration_flag,
+                                     shape_data,
+                                     temp,
+                                     fe_eval.begin_values(),
+                                     fe_eval.begin_gradients(),
+                                     fe_eval.begin_hessians(),
+                                     scratch_data,
+                                     subface_index);
+
+      if (use_vectorization == false)
+        {
+          for (unsigned int v = 0; v < Number::size(); ++v)
+            {
+              // the loop breaks once an invalid_unsigned_int is hit for
+              // all cases except the exterior faces in the ECL loop (where
+              // some faces might be at the boundaries but others not)
+              if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                continue;
+
+              FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+                template interpolate<false, false>(n_components,
+                                                   integration_flag,
+                                                   shape_info,
+                                                   values_dofs,
+                                                   scratch_data,
+                                                   fe_eval.get_face_no(v));
+
+              for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
+                   ++i)
+                temp[i][v] = scratch_data[i][v];
+            }
+        }
+      else
+        FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+          template interpolate<false, false>(n_components,
+                                             integration_flag,
+                                             shape_info,
+                                             temp,
+                                             values_dofs,
+                                             fe_eval.get_face_no());
+      return false;
+    }
+  };
+
+
+
+  template <int n_face_orientations,
+            typename Processor,
+            typename EvaluationData,
+            const bool check_face_orientations = false>
+  void
+  fe_face_evaluation_process_and_io(
+    Processor                             &proc,
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    typename Processor::Number2_          *global_vector_ptr,
+    const std::vector<ArrayView<const typename Processor::Number2_>> *sm_ptr,
+    const EvaluationData                                             &fe_eval,
+    typename Processor::VectorizedArrayType_                         *temp1)
+  {
+    constexpr int dim         = Processor::dim_;
+    constexpr int fe_degree   = Processor::fe_degree_;
+    using VectorizedArrayType = typename Processor::VectorizedArrayType_;
+    constexpr int n_lanes     = VectorizedArrayType::size();
+
+    using Number   = typename Processor::Number_;
+    using Number2_ = typename Processor::Number2_;
+
+    const auto        &shape_data = fe_eval.get_shape_info().data.front();
+    constexpr bool     integrate  = Processor::do_integrate;
+    const unsigned int face_no    = fe_eval.get_face_no();
+    const auto        &dof_info   = fe_eval.get_dof_info();
+    const unsigned int cell       = fe_eval.get_cell_or_face_batch_id();
+    const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index =
+      fe_eval.get_dof_access_index();
+    AssertIndexRange(cell,
+                     dof_info.index_storage_variants[dof_access_index].size());
+    constexpr unsigned int dofs_per_face =
+      Utilities::pow(fe_degree + 1, dim - 1);
+    const unsigned int subface_index = fe_eval.get_subface_index();
+
+    const unsigned int n_filled_lanes =
+      dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+
+    bool all_faces_are_same = n_filled_lanes == n_lanes;
+    if (n_face_orientations == n_lanes)
+      for (unsigned int v = 1; v < n_lanes; ++v)
+        if (fe_eval.get_face_no(v) != fe_eval.get_face_no(0) ||
+            fe_eval.get_face_orientation(v) != fe_eval.get_face_orientation(0))
+          {
+            all_faces_are_same = false;
+            break;
+          }
+
+    // check for re-orientation ...
+    std::array<const unsigned int *, n_face_orientations> orientation = {};
+
+    if (dim == 3 && n_face_orientations == n_lanes && !all_faces_are_same &&
+        fe_eval.is_interior_face() == 0)
+      for (unsigned int v = 0; v < n_lanes; ++v)
+        {
+          // the loop breaks once an invalid_unsigned_int is hit for
+          // all cases except the exterior faces in the ECL loop (where
+          // some faces might be at the boundaries but others not)
+          if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+            continue;
+
+          if (shape_data.nodal_at_cell_boundaries &&
+              fe_eval.get_face_orientation(v) != 0)
+            {
+              // ... and in case we detect a re-orientation, go to the other
+              // version of this function that actually allows for this
+              if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
+                  check_face_orientations == false)
+                {
+                  fe_face_evaluation_process_and_io<n_face_orientations,
+                                                    Processor,
+                                                    EvaluationData,
+                                                    true>(proc,
+                                                          n_components,
+                                                          evaluation_flag,
+                                                          global_vector_ptr,
+                                                          sm_ptr,
+                                                          fe_eval,
+                                                          temp1);
+                  return;
+                }
+              orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
+                fe_eval.get_face_orientation(v), 0);
+            }
+        }
+    else if (dim == 3 && fe_eval.get_face_orientation() != 0)
+      {
+        // go to the other version of this function
+        if (subface_index == GeometryInfo<dim>::max_children_per_cell &&
+            check_face_orientations == false)
+          {
+            fe_face_evaluation_process_and_io<n_face_orientations,
+                                              Processor,
+                                              EvaluationData,
+                                              true>(proc,
+                                                    n_components,
+                                                    evaluation_flag,
+                                                    global_vector_ptr,
+                                                    sm_ptr,
+                                                    fe_eval,
+                                                    temp1);
+            return;
+          }
+        for (unsigned int v = 0; v < n_face_orientations; ++v)
+          orientation[v] = &fe_eval.get_shape_info().face_orientations_dofs(
+            fe_eval.get_face_orientation(), 0);
+      }
+
+    // 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 =
+      shape_data
+        .shape_data_on_face[0][fe_degree + (integrate ? (2 - face_no % 2) :
+                                                        (1 + face_no % 2))];
+
+    // face_to_cell_index_hermite
+    std::array<const unsigned int *, n_face_orientations> index_array_hermite =
+      {};
+    if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
+      {
+        if (n_face_orientations == 1)
+          index_array_hermite[0] =
+            &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, 0);
+        else
+          {
+            for (unsigned int v = 0; v < n_lanes; ++v)
+              {
+                if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                  continue;
+
+                const auto face_no = fe_eval.get_face_no(v);
+
+                grad_weight[v] =
+                  shape_data.shape_data_on_face[0][fe_degree +
+                                                   (integrate ?
+                                                      (2 - (face_no % 2)) :
+                                                      (1 + (face_no % 2)))][0];
+
+                index_array_hermite[v] =
+                  &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no,
+                                                                       0);
+              }
+          }
+      }
+
+    // face_to_cell_index_nodal
+    std::array<const unsigned int *, n_face_orientations> index_array_nodal =
+      {};
+    if (shape_data.nodal_at_cell_boundaries == true)
+      {
+        if (n_face_orientations == 1)
+          index_array_nodal[0] =
+            &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no, 0);
+        else
+          {
+            for (unsigned int v = 0; v < n_lanes; ++v)
+              {
+                if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                  continue;
+
+                const auto face_no = fe_eval.get_face_no(v);
+
+                index_array_nodal[v] =
+                  &fe_eval.get_shape_info().face_to_cell_index_nodal(face_no,
+                                                                     0);
+              }
+          }
+      }
+
+
+    const auto reorientate = [&](const unsigned int v, const unsigned int i) {
+      return (!check_face_orientations || orientation[v] == nullptr) ?
+               i :
+               orientation[v][i];
+    };
+
+    const unsigned int cell_index =
+      dof_access_index == MatrixFreeFunctions::DoFInfo::dof_access_cell ?
+        fe_eval.get_cell_ids()[0] :
+        cell * n_lanes;
+    const unsigned int *dof_indices =
+      &dof_info.dof_indices_contiguous[dof_access_index][cell_index];
+
+    for (unsigned int comp = 0; comp < n_components; ++comp)
+      {
+        const std::size_t index_offset =
+          dof_info.component_dof_indices_offset
+            [fe_eval.get_active_fe_index()]
+            [fe_eval.get_first_selected_component()] +
+          comp * Utilities::pow(fe_degree + 1, dim);
+
+        // 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(
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+              n_lanes);
+            Number2_ *vector_ptr =
+              global_vector_ptr + dof_indices[0] + index_offset * n_lanes;
+
+            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
+              {
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
+                  {
+                    Assert(n_face_orientations == 1, ExcNotImplemented());
+
+                    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(temp1[i_],
+                                                 temp1[i_ + dofs_per_face],
+                                                 vector_ptr + ind1 * n_lanes,
+                                                 vector_ptr + ind2 * n_lanes,
+                                                 grad_weight);
+                  }
+              }
+            else
+              {
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
+                  {
+                    Assert(n_face_orientations == 1, 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 * n_lanes);
+                  }
+              }
+          }
+
+        // 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(
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+              n_lanes);
+            Number2_ *vector_ptr = global_vector_ptr + index_offset * n_lanes;
+            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
+              {
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
+                  {
+                    Assert(n_face_orientations == 1, ExcNotImplemented());
+
+                    const unsigned int i_ = reorientate(0, i);
+                    const unsigned int ind1 =
+                      index_array_hermite[0][2 * i] * n_lanes;
+                    const unsigned int ind2 =
+                      index_array_hermite[0][2 * i + 1] * n_lanes;
+                    proc.hermite_grad_vectorized_indexed(
+                      temp1[i_],
+                      temp1[i_ + dofs_per_face],
+                      vector_ptr + ind1,
+                      vector_ptr + ind2,
+                      grad_weight,
+                      dof_indices,
+                      dof_indices);
+                  }
+              }
+            else
+              {
+                for (unsigned int i = 0; i < dofs_per_face; ++i)
+                  {
+                    Assert(n_face_orientations == 1, ExcNotImplemented());
+
+                    const unsigned int i_  = reorientate(0, i);
+                    const unsigned int ind = index_array_nodal[0][i] * n_lanes;
+                    proc.value_vectorized_indexed(temp1[i_],
+                                                  vector_ptr + ind,
+                                                  dof_indices);
+                  }
+              }
+          }
+
+        // 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)
+          {
+            const unsigned int *strides =
+              &dof_info.dof_indices_interleave_strides[dof_access_index]
+                                                      [cell * n_lanes];
+            unsigned int indices[n_lanes];
+            for (unsigned int v = 0; v < n_lanes; ++v)
+              indices[v] = dof_indices[v] + index_offset * strides[v];
+            const unsigned int n_filled_lanes =
+              dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+
+            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
+              {
+                if (n_filled_lanes == n_lanes)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      Assert(n_face_orientations == 1, ExcNotImplemented());
+
+                      const unsigned int i_ = reorientate(0, i);
+                      unsigned int       ind1[n_lanes];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0; v < n_lanes; ++v)
+                        ind1[v] = indices[v] +
+                                  index_array_hermite[0][2 * i] * strides[v];
+                      unsigned int ind2[n_lanes];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0; v < n_lanes; ++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],
+                        global_vector_ptr,
+                        global_vector_ptr,
+                        grad_weight,
+                        ind1,
+                        ind2);
+                    }
+                else
+                  {
+                    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]);
+                        }
+                  }
+              }
+            else
+              {
+                if (n_filled_lanes == n_lanes)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      Assert(n_face_orientations == 1, ExcInternalError());
+                      unsigned int ind[n_lanes];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0; v < n_lanes; ++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
+                  {
+                    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]]);
+                  }
+              }
+          }
+
+        // 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)
+          {
+            Number2_ *vector_ptr = global_vector_ptr + index_offset;
+
+            const bool vectorization_possible =
+              all_faces_are_same && (sm_ptr == nullptr);
+
+            std::array<Number2_ *, n_lanes>   vector_ptrs;
+            std::array<unsigned int, n_lanes> reordered_indices;
+
+            if (vectorization_possible == false)
+              {
+                vector_ptrs = {};
+                if (n_face_orientations == 1)
+                  {
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      if (sm_ptr == nullptr)
+                        {
+                          vector_ptrs[v] = vector_ptr + dof_indices[v];
+                        }
+                      else
+                        {
+                          const auto &temp =
+                            dof_info
+                              .dof_indices_contiguous_sm[dof_access_index]
+                                                        [cell * n_lanes + v];
+                          vector_ptrs[v] = const_cast<Number2_ *>(
+                            sm_ptr->operator[](temp.first).data() +
+                            temp.second + index_offset);
+                        }
+                  }
+                else if (n_face_orientations == n_lanes)
+                  {
+                    const auto &cells = fe_eval.get_cell_ids();
+                    for (unsigned int v = 0; v < n_lanes; ++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<Number2_ *>(
+                                sm_ptr->operator[](temp.first).data() +
+                                temp.second + index_offset);
+                            }
+                        }
+                  }
+                else
+                  {
+                    Assert(false, ExcNotImplemented());
+                  }
+              }
+            else if (n_face_orientations == n_lanes)
+              {
+                for (unsigned int v = 0; v < n_lanes; ++v)
+                  reordered_indices[v] =
+                    dof_info.dof_indices_contiguous[dof_access_index]
+                                                   [fe_eval.get_cell_ids()[v]];
+                dof_indices = reordered_indices.data();
+              }
+
+            if (fe_degree > 1 && (evaluation_flag & EvaluationFlags::gradients))
+              {
+                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,
+                        dof_indices,
+                        dof_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 < n_lanes; ++v)
+                          {
+                            temp1[i][v]                 = 0.0;
+                            temp1[i + dofs_per_face][v] = 0.0;
+                          }
+                    }
+                else
+                  {
+                    if (integrate == false && n_filled_lanes < n_lanes)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        temp1[i] = temp1[i + dofs_per_face] = Number();
+
+                    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)
+                  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,
+                                                    dof_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 < n_lanes; ++v)
+                          temp1[i_][v] = 0.0;
+                    }
+                else
+                  {
+                    if (integrate == false && n_filled_lanes < n_lanes)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        temp1[i] = Number();
+
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        proc.value(temp1[reorientate(v, i)][v],
+                                   vector_ptrs[v][index_array_nodal[v][i]]);
+                  }
+              }
+          }
+        else
+          {
+            // We should not end up here, this should be caught by
+            // FEFaceEvaluationImplGatherEvaluateSelector::supports()
+            Assert(false, ExcInternalError());
+          }
+        temp1 += 3 * dofs_per_face;
+      }
+  }
+
+
+
+  template <int dim, typename Number2, typename VectorizedArrayType>
+  struct FEFaceEvaluationImplGatherEvaluateSelector
+  {
+    using Number = typename VectorizedArrayType::value_type;
+
+    template <int fe_degree, int n_q_points_1d>
+    static bool
+    run(const unsigned int                                n_components,
+        const EvaluationFlags::EvaluationFlags            evaluation_flag,
+        const Number2                                    *src_ptr,
+        const std::vector<ArrayView<const Number2>>      *sm_ptr,
+        FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
+    {
+      Assert(fe_degree > -1, ExcInternalError());
+      Assert(fe_eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
+             ExcInternalError());
+
+      const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
+
+      VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
+      VectorizedArrayType *scratch_data =
+        temp + 3 * n_components * dofs_per_face;
+
+      Processor<fe_degree> p;
+
+      if (fe_eval.get_dof_access_index() ==
+            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          fe_eval.is_interior_face() == false)
+        fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
+          p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
+      else
+        fe_face_evaluation_process_and_io<1>(
+          p, n_components, evaluation_flag, src_ptr, sm_ptr, fe_eval, temp);
+
+      const unsigned int subface_index = fe_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,
+                           fe_eval.get_shape_info().data.front(),
+                           temp,
+                           fe_eval.begin_values(),
+                           fe_eval.begin_gradients(),
+                           fe_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,
+                           fe_eval.get_shape_info().data.front(),
+                           temp,
+                           fe_eval.begin_values(),
+                           fe_eval.begin_gradients(),
+                           fe_eval.begin_hessians(),
+                           scratch_data,
+                           subface_index);
+
+      // re-orientation for cases not possible with above algorithm
+      if (subface_index < GeometryInfo<dim>::max_children_per_cell)
+        {
+          if (fe_eval.get_dof_access_index() ==
+                MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+              fe_eval.is_interior_face() == false)
+            {
+              for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+                {
+                  // the loop breaks once an invalid_unsigned_int is hit for
+                  // all cases except the exterior faces in the ECL loop (where
+                  // some faces might be at the boundaries but others not)
+                  if (fe_eval.get_cell_ids()[v] ==
+                      numbers::invalid_unsigned_int)
+                    continue;
+
+                  if (fe_eval.get_face_orientation(v) != 0)
+                    adjust_for_face_orientation_per_lane(
+                      dim,
+                      n_components,
+                      v,
+                      evaluation_flag,
+                      &fe_eval.get_shape_info().face_orientations_quad(
+                        fe_eval.get_face_orientation(v), 0),
+                      false,
+                      Utilities::pow(n_q_points_1d, dim - 1),
+                      &temp[0][0],
+                      fe_eval.begin_values(),
+                      fe_eval.begin_gradients(),
+                      fe_eval.begin_hessians());
+                }
+            }
+          else if (fe_eval.get_face_orientation() != 0)
+            adjust_for_face_orientation(
+              dim,
+              n_components,
+              evaluation_flag,
+              &fe_eval.get_shape_info().face_orientations_quad(
+                fe_eval.get_face_orientation(), 0),
+              false,
+              Utilities::pow(n_q_points_1d, dim - 1),
+              temp,
+              fe_eval.begin_values(),
+              fe_eval.begin_gradients(),
+              fe_eval.begin_hessians());
+        }
+
+      return false;
+    }
+
+    template <typename Number3>
+    static bool
+    supports(const EvaluationFlags::EvaluationFlags             evaluation_flag,
+             const MatrixFreeFunctions::ShapeInfo<Number3>     &shape_info,
+             const Number2                                     *vector_ptr,
+             MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage)
+    {
+      const unsigned int fe_degree = shape_info.data.front().fe_degree;
+      if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries ||
+          (evaluation_flag & EvaluationFlags::gradients &&
+           (fe_degree < 2 ||
+            shape_info.data.front().element_type !=
+              MatrixFreeFunctions::tensor_symmetric_hermite)) ||
+          (evaluation_flag & EvaluationFlags::hessians) ||
+          vector_ptr == nullptr ||
+          shape_info.data.front().element_type >
+            MatrixFreeFunctions::tensor_symmetric_no_collocation ||
+          storage <
+            MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
+        return false;
+      else
+        return true;
+    }
+
+  private:
+    template <int fe_degree>
+    struct Processor
+    {
+      static const bool do_integrate = false;
+      static const int  dim_         = dim;
+      static const int  fe_degree_   = fe_degree;
+      using VectorizedArrayType_     = VectorizedArrayType;
+      using Number_                  = Number;
+      using Number2_                 = const Number2;
+
+      template <typename T0, typename T1, typename T2>
+      void
+      hermite_grad_vectorized(T0       &temp_1,
+                              T0       &temp_2,
+                              const T1  src_ptr_1,
+                              const T1  src_ptr_2,
+                              const T2 &grad_weight)
+      {
+        do_vectorized_read(src_ptr_1, temp_1);
+        do_vectorized_read(src_ptr_2, temp_2);
+        temp_2 = grad_weight * (temp_1 - temp_2);
+      }
+
+      template <typename T1, typename T2>
+      void
+      value_vectorized(T1 &temp, const T2 src_ptr)
+      {
+        do_vectorized_read(src_ptr, temp);
+      }
+
+      template <typename T0, typename T1, typename T2, typename T3>
+      void
+      hermite_grad_vectorized_indexed(T0       &temp_1,
+                                      T0       &temp_2,
+                                      const T1  src_ptr_1,
+                                      const T1  src_ptr_2,
+                                      const T2 &grad_weight,
+                                      const T3 &indices_1,
+                                      const T3 &indices_2)
+      {
+        do_vectorized_gather(src_ptr_1, indices_1, temp_1);
+        do_vectorized_gather(src_ptr_2, indices_2, temp_2);
+        temp_2 = grad_weight * (temp_1 - temp_2);
+      }
+
+      template <typename T0, typename T1, typename T2>
+      void
+      value_vectorized_indexed(T0 &temp, const T1 src_ptr, const T2 &indices)
+      {
+        do_vectorized_gather(src_ptr, indices, temp);
+      }
+
+      template <typename T0, typename T1, typename T2>
+      void
+      hermite_grad(T0       &temp_1,
+                   T0       &temp_2,
+                   const T1 &src_ptr_1,
+                   const T1 &src_ptr_2,
+                   const T2 &grad_weight)
+      {
+        // case 3a)
+        temp_1 = src_ptr_1;
+        temp_2 = grad_weight * (temp_1 - src_ptr_2);
+      }
+
+      template <typename T1, typename T2>
+      void
+      value(T1 &temp, const T2 &src_ptr)
+      {
+        // case 3b)
+        temp = src_ptr;
+      }
+    };
+  };
+
+
+
+  template <int dim, typename Number2, typename VectorizedArrayType>
+  struct FEFaceEvaluationImplIntegrateScatterSelector
+  {
+    using Number = typename VectorizedArrayType::value_type;
+
+    template <int fe_degree, int n_q_points_1d>
+    static bool
+    run(const unsigned int                                n_components,
+        const EvaluationFlags::EvaluationFlags            integration_flag,
+        Number2                                          *dst_ptr,
+        const std::vector<ArrayView<const Number2>>      *sm_ptr,
+        FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
+    {
+      Assert(fe_degree > -1, ExcInternalError());
+      Assert(fe_eval.get_shape_info().element_type <=
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
+             ExcInternalError());
+
+      const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
+
+      VectorizedArrayType *temp = fe_eval.get_scratch_data().begin();
+      VectorizedArrayType *scratch_data =
+        temp + 3 * n_components * dofs_per_face;
+
+      const unsigned int subface_index = fe_eval.get_subface_index();
+
+      // re-orientation for cases not possible with the io function below
+      if (subface_index < GeometryInfo<dim>::max_children_per_cell)
+        {
+          if (fe_eval.get_dof_access_index() ==
+                MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+              fe_eval.is_interior_face() == false)
+            for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+              {
+                // the loop breaks once an invalid_unsigned_int is hit for
+                // all cases except the exterior faces in the ECL loop (where
+                // some faces might be at the boundaries but others not)
+                if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
+                  continue;
+
+                if (fe_eval.get_face_orientation(v) != 0)
+                  adjust_for_face_orientation_per_lane(
+                    dim,
+                    n_components,
+                    v,
+                    integration_flag,
+                    &fe_eval.get_shape_info().face_orientations_quad(
+                      fe_eval.get_face_orientation(v), 0),
+                    true,
+                    Utilities::pow(n_q_points_1d, dim - 1),
+                    &temp[0][0],
+                    fe_eval.begin_values(),
+                    fe_eval.begin_gradients(),
+                    fe_eval.begin_hessians());
+              }
+          else if (fe_eval.get_face_orientation() != 0)
+            adjust_for_face_orientation(
+              dim,
+              n_components,
+              integration_flag,
+              &fe_eval.get_shape_info().face_orientations_quad(
+                fe_eval.get_face_orientation(), 0),
+              true,
+              Utilities::pow(n_q_points_1d, dim - 1),
+              temp,
+              fe_eval.begin_values(),
+              fe_eval.begin_gradients(),
+              fe_eval.begin_hessians());
+        }
+
+      if (fe_degree > -1 && fe_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,
+                            fe_eval.get_shape_info().data.front(),
+                            temp,
+                            fe_eval.begin_values(),
+                            fe_eval.begin_gradients(),
+                            fe_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,
+                            fe_eval.get_shape_info().data.front(),
+                            temp,
+                            fe_eval.begin_values(),
+                            fe_eval.begin_gradients(),
+                            fe_eval.begin_hessians(),
+                            scratch_data,
+                            subface_index);
+
+      Processor<fe_degree> p;
+
+      if (fe_eval.get_dof_access_index() ==
+            MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          fe_eval.is_interior_face() == false)
+        fe_face_evaluation_process_and_io<VectorizedArrayType::size()>(
+          p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
+      else
+        fe_face_evaluation_process_and_io<1>(
+          p, n_components, integration_flag, dst_ptr, sm_ptr, fe_eval, temp);
+
+      return false;
+    }
+
+  private:
+    template <int fe_degree>
+    struct Processor
+    {
+      static const bool do_integrate = true;
+      static const int  dim_         = dim;
+      static const int  fe_degree_   = fe_degree;
+      using VectorizedArrayType_     = VectorizedArrayType;
+      using Number_                  = Number;
+      using Number2_                 = Number2;
+
+      template <typename T0, typename T1, typename T2, typename T3, typename T4>
+      void
+      hermite_grad_vectorized(const T0 &temp_1,
+                              const T1 &temp_2,
+                              T2        dst_ptr_1,
+                              T3        dst_ptr_2,
+                              const T4 &grad_weight)
+      {
+        // case 1a)
+        const VectorizedArrayType val  = temp_1 - grad_weight * temp_2;
+        const VectorizedArrayType grad = grad_weight * temp_2;
+        do_vectorized_add(val, dst_ptr_1);
+        do_vectorized_add(grad, dst_ptr_2);
+      }
+
+      template <typename T0, typename T1>
+      void
+      value_vectorized(const T0 &temp, T1 dst_ptr)
+      {
+        // case 1b)
+        do_vectorized_add(temp, dst_ptr);
+      }
+
+      template <typename T0, typename T1, typename T2, typename T3>
+      void
+      hermite_grad_vectorized_indexed(const T0 &temp_1,
+                                      const T0 &temp_2,
+                                      T1        dst_ptr_1,
+                                      T1        dst_ptr_2,
+                                      const T2 &grad_weight,
+                                      const T3 &indices_1,
+                                      const T3 &indices_2)
+      {
+        // case 2a)
+        const VectorizedArrayType val  = temp_1 - grad_weight * temp_2;
+        const VectorizedArrayType grad = grad_weight * temp_2;
+        do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
+        do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
+      }
+
+      template <typename T0, typename T1, typename T2>
+      void
+      value_vectorized_indexed(const T0 &temp, T1 dst_ptr, const T2 &indices)
+      {
+        // case 2b)
+        do_vectorized_scatter_add(temp, indices, dst_ptr);
+      }
+
+      template <typename T0, typename T1, typename T2>
+      void
+      hermite_grad(const T0 &temp_1,
+                   const T0 &temp_2,
+                   T1       &dst_ptr_1,
+                   T1       &dst_ptr_2,
+                   const T2 &grad_weight)
+      {
+        // case 3a)
+        const Number val  = temp_1 - grad_weight * temp_2;
+        const Number grad = grad_weight * temp_2;
+        dst_ptr_1 += val;
+        dst_ptr_2 += grad;
+      }
+
+      template <typename T0, typename T1>
+      void
+      value(const T0 &temp, T1 &dst_ptr)
+      {
+        // case 3b)
+        dst_ptr += temp;
+      }
+    };
+  };
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index a429df620887e3bc51f4dbd087900c46426473d9..74c7c4899a2df307ad45288c539d64d4fe8c0a35 100644 (file)
@@ -20,8 +20,7 @@
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/matrix_free/evaluation_kernels.h>
-#include <deal.II/matrix_free/evaluation_selector.h>
+#include <deal.II/matrix_free/evaluation_kernels_face.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/evaluation_template_factory_internal.h>
 
index c2378a9812ce78fccce461f48d92b1dcff30664e..ea52c477dd5b1c36fb3c142e1e05cc3ad693817a 100644 (file)
@@ -28,6 +28,21 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
+  /**
+   * This struct is used to implement
+   * FEEvaluation::fast_evaluation_supported() and
+   * FEFaceEvaluation::fast_evaluation_supported().
+   */
+  struct FastEvaluationSupported
+  {
+    template <int fe_degree, int n_q_points_1d>
+    static bool
+    run()
+    {
+      return fe_degree != -1;
+    }
+  };
+
   template <int degree, typename EvaluatorType, typename... Args>
   bool
   instantiation_helper_run(const unsigned int given_degree,
index cf7d1c405730cba6edb27ddc1f2fcf1ad689f40f..9c3bcfdb334aae3403439f47bcb89ba12750515c 100644 (file)
@@ -31,6 +31,7 @@
 
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/evaluation_kernels_face.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/fe_evaluation_data.h>
 #include <deal.II/matrix_free/hanging_nodes_internal.h>
index 696a16c8b2e8ca108b6e54be12659f8fc1b8fb92..9c7a56854236f303c97bc040c0288abd8c46e559 100644 (file)
@@ -29,7 +29,7 @@
 #include <deal.II/fe/mapping.h>
 
 #include <deal.II/matrix_free/evaluation_flags.h>
-#include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/evaluation_kernels_face.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 
index 0f8f660cc68902f5bad6e71cd74247a1b679b540..837510b98638dca3901ad6f740773b0b98e97ac3 100644 (file)
@@ -24,6 +24,8 @@
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/utilities.h>
 
+#include <deal.II/matrix_free/shape_info.h>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -1653,6 +1655,173 @@ namespace internal
 
 
 
+  template <int  dim,
+            int  fe_degree,
+            int  n_q_points_1d,
+            bool contract_over_rows,
+            bool symmetric_evaluate = true>
+  struct EvaluatorTensorProductAnisotropic
+  {
+    template <int direction,
+              int stride       = 1,
+              typename Number  = double,
+              typename Number2 = double>
+    static void
+    normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+           const Number                                            *in,
+           Number                                                  *out,
+           const bool add_into_result  = false,
+           const int  subface_index_1d = 0)
+    {
+      AssertIndexRange(direction, dim);
+      AssertDimension(fe_degree, data.fe_degree);
+      AssertDimension(n_q_points_1d, data.n_q_points_1d);
+      constexpr int  n_rows    = fe_degree + 1;
+      constexpr int  n_columns = n_q_points_1d;
+      constexpr int  mm        = contract_over_rows ? n_rows : n_columns;
+      constexpr int  nn        = contract_over_rows ? n_columns : n_rows;
+      const Number2 *shape_data =
+        symmetric_evaluate ?
+          data.shape_values_eo.data() :
+          data.values_within_subface[subface_index_1d].data();
+      Assert(shape_data != nullptr, ExcNotInitialized());
+      Assert(contract_over_rows == false || !add_into_result,
+             ExcMessage("Cannot add into result if contract_over_rows = true"));
+
+      constexpr int n_blocks1  = Utilities::pow(fe_degree, direction);
+      constexpr int n_blocks2  = Utilities::pow(fe_degree, dim - direction - 1);
+      constexpr int stride_in  = contract_over_rows ? 1 : stride;
+      constexpr int stride_out = contract_over_rows ? stride : 1;
+      constexpr EvaluatorVariant variant =
+        symmetric_evaluate ? evaluate_evenodd : evaluate_general;
+
+      for (int i2 = 0; i2 < n_blocks2; ++i2)
+        {
+          for (int i1 = 0; i1 < n_blocks1; ++i1)
+            {
+              if (contract_over_rows == false && add_into_result)
+                apply_matrix_vector_product<variant,
+                                            EvaluatorQuantity::value,
+                                            n_rows,
+                                            n_columns,
+                                            n_blocks1 * stride_in,
+                                            n_blocks1 * stride_out,
+                                            contract_over_rows,
+                                            true>(shape_data, in, out);
+              else
+                apply_matrix_vector_product<variant,
+                                            EvaluatorQuantity::value,
+                                            n_rows,
+                                            n_columns,
+                                            n_blocks1 * stride_in,
+                                            n_blocks1 * stride_out,
+                                            contract_over_rows,
+                                            false>(shape_data, in, out);
+
+              in += stride_in;
+              out += stride_out;
+            }
+          in += n_blocks1 * (mm - 1) * stride_in;
+          out += n_blocks1 * (nn - 1) * stride_out;
+        }
+    }
+
+    template <int direction,
+              int normal_direction,
+              int stride       = 1,
+              typename Number  = double,
+              typename Number2 = double>
+    static void
+    tangential(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+               const Number                                            *in,
+               Number                                                  *out,
+               const int subface_index_1d = 0)
+    {
+      AssertIndexRange(direction, dim);
+      AssertDimension(fe_degree - 1, data.fe_degree);
+      AssertDimension(n_q_points_1d, data.n_q_points_1d);
+      static_assert(direction != normal_direction,
+                    "Cannot interpolate tangentially in normal direction");
+
+      constexpr int  n_rows    = fe_degree;
+      constexpr int  n_columns = n_q_points_1d;
+      const Number2 *shape_data =
+        symmetric_evaluate ?
+          data.shape_values_eo.data() :
+          data.values_within_subface[subface_index_1d].data();
+      Assert(shape_data != nullptr, ExcNotInitialized());
+
+      constexpr int n_blocks1 =
+        (direction > normal_direction) ?
+          Utilities::pow(n_q_points_1d, direction) :
+          (direction > 0 ?
+             (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
+             1);
+      constexpr int n_blocks2 =
+        (direction > normal_direction) ?
+          Utilities::pow(fe_degree, dim - 1 - direction) :
+          ((direction + 1 < dim) ?
+             (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
+             1);
+
+      constexpr EvaluatorVariant variant =
+        symmetric_evaluate ? evaluate_evenodd : evaluate_general;
+
+      // Since we may perform an in-place interpolation, we must run the step
+      // expanding the size of the basis backward ('contract_over_rows' aka
+      // 'evaluate' case), so shift the pointers and decrement during the loop
+      if (contract_over_rows)
+        {
+          in += (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
+          out +=
+            stride * ((n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1);
+          for (int i2 = 0; i2 < n_blocks2; ++i2)
+            {
+              for (int i1 = 0; i1 < n_blocks1; ++i1)
+                {
+                  apply_matrix_vector_product<variant,
+                                              EvaluatorQuantity::value,
+                                              n_rows,
+                                              n_columns,
+                                              n_blocks1,
+                                              n_blocks1 * stride,
+                                              true,
+                                              false>(shape_data, in, out);
+
+                  --in;
+                  out -= stride;
+                }
+              in -= n_blocks1 * (n_rows - 1);
+              out -= n_blocks1 * (n_columns - 1) * stride;
+            }
+        }
+      else
+        {
+          for (int i2 = 0; i2 < n_blocks2; ++i2)
+            {
+              for (int i1 = 0; i1 < n_blocks1; ++i1)
+                {
+                  apply_matrix_vector_product<variant,
+                                              EvaluatorQuantity::value,
+                                              n_rows,
+                                              n_columns,
+                                              n_blocks1 * stride,
+                                              n_blocks1,
+                                              false,
+                                              false>(shape_data, in, out);
+
+                  in += stride;
+                  ++out;
+                }
+              in += n_blocks1 * (n_columns - 1) * stride;
+              out += n_blocks1 * (n_rows - 1);
+            }
+        }
+    }
+  };
+
+
+
   /**
    * This function applies the tensor product operation to produce face values
    * from cell values. The algorithm involved here can be interpreted as the
@@ -1767,6 +1936,23 @@ namespace internal
 
 
 
+  /**
+   * Helper function to specify whether a transformation to collocation should
+   * be used: It should give correct results (first condition), we need to be
+   * able to initialize the fields in shape_info.templates.h from the
+   * polynomials (second condition), and it should be the most efficient
+   * choice in terms of operation counts (third condition).
+   */
+  constexpr bool
+  use_collocation_evaluation(const unsigned int fe_degree,
+                             const unsigned int n_q_points_1d)
+  {
+    return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
+           (n_q_points_1d <= 3 * fe_degree / 2 + 1);
+  }
+
+
+
   /**
    * This function performs the opposite operation to the interpolate_to_face
    * function, done as the last step in sum factorization to embed face values

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.