]> https://gitweb.dealii.org/ - dealii.git/commitdiff
FEEvaluation: Restructure layout of gradients.
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Wed, 30 Aug 2023 08:21:51 +0000 (10:21 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 11 Sep 2023 19:58:51 +0000 (21:58 +0200)
Generalize and optimize Raviart-Thomas evaluation significantly

include/deal.II/fe/mapping_q_internal.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/tensor_product_kernels.h
tests/matrix_free/matrix_vector_hessians_cells.cc

index edb6e123091d1a0a932253a101331470e97b5917..80250bbdd76c63516a69b50803b6ec39e870978c 100644 (file)
@@ -1165,17 +1165,12 @@ namespace internal
           // We need to reinterpret the data after evaluate has been applied.
           for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
             for (unsigned int point = 0; point < n_q_points; ++point)
-              for (unsigned int j = 0; j < dim; ++j)
-                for (unsigned int in_comp = 0;
-                     in_comp < n_lanes &&
-                     in_comp < spacedim - out_comp * n_lanes;
-                     ++in_comp)
+              for (unsigned int in_comp = 0;
+                   in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
+                   ++in_comp)
+                for (unsigned int j = 0; j < dim; ++j)
                   {
-                    const unsigned int total_number = point * dim + j;
-                    const unsigned int new_comp     = total_number / n_q_points;
-                    const unsigned int new_point    = total_number % n_q_points;
-                    jacobians[new_point][out_comp * n_lanes +
-                                         in_comp][new_comp] =
+                    jacobians[point][out_comp * n_lanes + in_comp][j] =
                       eval.begin_gradients()[(out_comp * n_q_points + point) *
                                                dim +
                                              j][in_comp];
index 95583af217010e56fb12215a765d6768ca660510..1401a041807bde48d97f4846dee00c68e2f686be 100644 (file)
@@ -188,55 +188,6 @@ namespace internal
               const bool                             add_into_values_array);
   };
 
-  /**
-   * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use
-   * specific sum-factorization kernels and with normal/tangential shape_data
-   */
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
-                          dim,
-                          fe_degree,
-                          n_q_points_1d,
-                          Number>
-  {
-    template <bool integrate>
-    static void
-    evaluate_or_integrate(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      Number                                *values_dofs_actual,
-      FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array = false);
-
-  private:
-    template <typename EvalType, typename Number2>
-    static EvalType
-    create_evaluator_tensor_product(
-      const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape_data)
-    {
-      return EvalType(shape_data.shape_values,
-                      shape_data.shape_gradients,
-                      shape_data.shape_hessians);
-    }
-
-    template <int normal_dir>
-    static void
-    evaluate_tensor_product_per_component(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      Number                                *values_dofs_actual,
-      FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array,
-      std::integral_constant<bool, false>);
-
-    template <int normal_dir>
-    static void
-    evaluate_tensor_product_per_component(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      Number                                *values_dofs_actual,
-      FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array,
-      std::integral_constant<bool, true>);
-  };
-
 
 
   template <MatrixFreeFunctions::ElementType type,
@@ -329,7 +280,6 @@ namespace internal
 
     Number *values_quad    = fe_eval.begin_values();
     Number *gradients_quad = fe_eval.begin_gradients();
-    Number *hessians_quad  = fe_eval.begin_hessians();
 
     switch (dim)
       {
@@ -341,15 +291,11 @@ namespace internal
               if (evaluation_flag & EvaluationFlags::gradients)
                 eval0.template gradients<0, true, false>(values_dofs,
                                                          gradients_quad);
-              if (evaluation_flag & EvaluationFlags::hessians)
-                eval0.template hessians<0, true, false>(values_dofs,
-                                                        hessians_quad);
 
               // advance the next component in 1d array
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
               gradients_quad += n_q_points;
-              hessians_quad += n_q_points;
             }
           break;
 
@@ -360,35 +306,15 @@ namespace internal
               if (evaluation_flag & EvaluationFlags::gradients)
                 {
                   eval0.template gradients<0, true, false>(values_dofs, temp1);
-                  eval1.template values<1, true, false>(temp1, gradients_quad);
-                }
-              if (evaluation_flag & EvaluationFlags::hessians)
-                {
-                  // grad xy
-                  if (!(evaluation_flag & EvaluationFlags::gradients))
-                    eval0.template gradients<0, true, false>(values_dofs,
-                                                             temp1);
-                  eval1.template gradients<1, true, false>(temp1,
-                                                           hessians_quad +
-                                                             2 * n_q_points);
-
-                  // grad xx
-                  eval0.template hessians<0, true, false>(values_dofs, temp1);
-                  eval1.template values<1, true, false>(temp1, hessians_quad);
+                  eval1.template values<1, true, false, 2>(temp1,
+                                                           gradients_quad);
                 }
 
               // grad y
               eval0.template values<0, true, false>(values_dofs, temp1);
               if (evaluation_flag & EvaluationFlags::gradients)
-                eval1.template gradients<1, true, false>(temp1,
-                                                         gradients_quad +
-                                                           n_q_points);
-
-              // grad yy
-              if (evaluation_flag & EvaluationFlags::hessians)
-                eval1.template hessians<1, true, false>(temp1,
-                                                        hessians_quad +
-                                                          n_q_points);
+                eval1.template gradients<1, true, false, 2>(temp1,
+                                                            gradients_quad + 1);
 
               // val: can use values applied in x
               if (evaluation_flag & EvaluationFlags::values)
@@ -398,7 +324,6 @@ namespace internal
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
               gradients_quad += 2 * n_q_points;
-              hessians_quad += 3 * n_q_points;
             }
           break;
 
@@ -410,32 +335,8 @@ namespace internal
                   // grad x
                   eval0.template gradients<0, true, false>(values_dofs, temp1);
                   eval1.template values<1, true, false>(temp1, temp2);
-                  eval2.template values<2, true, false>(temp2, gradients_quad);
-                }
-
-              if (evaluation_flag & EvaluationFlags::hessians)
-                {
-                  // grad xz
-                  if (!(evaluation_flag & EvaluationFlags::gradients))
-                    {
-                      eval0.template gradients<0, true, false>(values_dofs,
-                                                               temp1);
-                      eval1.template values<1, true, false>(temp1, temp2);
-                    }
-                  eval2.template gradients<2, true, false>(temp2,
-                                                           hessians_quad +
-                                                             4 * n_q_points);
-
-                  // grad xy
-                  eval1.template gradients<1, true, false>(temp1, temp2);
-                  eval2.template values<2, true, false>(temp2,
-                                                        hessians_quad +
-                                                          3 * n_q_points);
-
-                  // grad xx
-                  eval0.template hessians<0, true, false>(values_dofs, temp1);
-                  eval1.template values<1, true, false>(temp1, temp2);
-                  eval2.template values<2, true, false>(temp2, hessians_quad);
+                  eval2.template values<2, true, false, 3>(temp2,
+                                                           gradients_quad);
                 }
 
               // grad y
@@ -443,41 +344,16 @@ namespace internal
               if (evaluation_flag & EvaluationFlags::gradients)
                 {
                   eval1.template gradients<1, true, false>(temp1, temp2);
-                  eval2.template values<2, true, false>(temp2,
-                                                        gradients_quad +
-                                                          n_q_points);
-                }
-
-              if (evaluation_flag & EvaluationFlags::hessians)
-                {
-                  // grad yz
-                  if (!(evaluation_flag & EvaluationFlags::gradients))
-                    eval1.template gradients<1, true, false>(temp1, temp2);
-                  eval2.template gradients<2, true, false>(temp2,
-                                                           hessians_quad +
-                                                             5 * n_q_points);
-
-                  // grad yy
-                  eval1.template hessians<1, true, false>(temp1, temp2);
-                  eval2.template values<2, true, false>(temp2,
-                                                        hessians_quad +
-                                                          n_q_points);
+                  eval2.template values<2, true, false, 3>(temp2,
+                                                           gradients_quad + 1);
                 }
 
               // grad z: can use the values applied in x direction stored in
               // temp1
               eval1.template values<1, true, false>(temp1, temp2);
               if (evaluation_flag & EvaluationFlags::gradients)
-                eval2.template gradients<2, true, false>(temp2,
-                                                         gradients_quad +
-                                                           2 * n_q_points);
-
-              // grad zz: can use the values applied in x and y direction stored
-              // in temp2
-              if (evaluation_flag & EvaluationFlags::hessians)
-                eval2.template hessians<2, true, false>(temp2,
-                                                        hessians_quad +
-                                                          2 * n_q_points);
+                eval2.template gradients<2, true, false, 3>(temp2,
+                                                            gradients_quad + 2);
 
               // val: can use the values applied in x & y direction stored in
               // temp2
@@ -488,7 +364,6 @@ namespace internal
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
               gradients_quad += 3 * n_q_points;
-              hessians_quad += 6 * n_q_points;
             }
           break;
 
@@ -576,7 +451,6 @@ namespace internal
 
     Number *values_quad    = fe_eval.begin_values();
     Number *gradients_quad = fe_eval.begin_gradients();
-    Number *hessians_quad  = fe_eval.begin_hessians();
 
     switch (dim)
       {
@@ -602,23 +476,11 @@ namespace internal
                     eval0.template gradients<0, false, false>(gradients_quad,
                                                               values_dofs);
                 }
-              if ((integration_flag & EvaluationFlags::hessians) != 0u)
-                {
-                  if ((integration_flag & EvaluationFlags::values) != 0u ||
-                      (integration_flag & EvaluationFlags::gradients) != 0u ||
-                      add_into_values_array == true)
-                    eval0.template hessians<0, false, true>(hessians_quad,
-                                                            values_dofs);
-                  else
-                    eval0.template hessians<0, false, false>(hessians_quad,
-                                                             values_dofs);
-                }
 
               // advance to the next component in 1d array
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
               gradients_quad += n_q_points;
-              hessians_quad += n_q_points;
             }
           break;
 
@@ -636,40 +498,16 @@ namespace internal
                 }
               if (integration_flag & EvaluationFlags::gradients)
                 {
-                  eval1.template gradients<1, false, false>(gradients_quad +
-                                                              n_q_points,
-                                                            temp1);
+                  eval1.template gradients<1, false, false, 2>(gradients_quad +
+                                                                 1,
+                                                               temp1);
                   if (integration_flag & EvaluationFlags::values)
                     eval1.template values<1, false, true>(values_quad, temp1);
                   if (add_into_values_array == false)
                     eval0.template values<0, false, false>(temp1, values_dofs);
                   else
                     eval0.template values<0, false, true>(temp1, values_dofs);
-                  eval1.template values<1, false, false>(gradients_quad, temp1);
-                  eval0.template gradients<0, false, true>(temp1, values_dofs);
-                }
-              if ((integration_flag & EvaluationFlags::hessians) != 0u)
-                {
-                  // grad xx
-                  eval1.template values<1, false, false>(hessians_quad, temp1);
-
-                  if ((integration_flag & EvaluationFlags::values) != 0u ||
-                      (integration_flag & EvaluationFlags::gradients) != 0u ||
-                      add_into_values_array == true)
-                    eval0.template hessians<0, false, true>(temp1, values_dofs);
-                  else
-                    eval0.template hessians<0, false, false>(temp1,
-                                                             values_dofs);
-
-                  // grad yy
-                  eval1.template hessians<1, false, false>(hessians_quad +
-                                                             n_q_points,
-                                                           temp1);
-                  eval0.template values<0, false, true>(temp1, values_dofs);
-
-                  // grad xy
-                  eval1.template gradients<1, false, false>(hessians_quad +
-                                                              2 * n_q_points,
+                  eval1.template values<1, false, false, 2>(gradients_quad,
                                                             temp1);
                   eval0.template gradients<0, false, true>(temp1, values_dofs);
                 }
@@ -678,7 +516,6 @@ namespace internal
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
               gradients_quad += 2 * n_q_points;
-              hessians_quad += 3 * n_q_points;
             }
           break;
 
@@ -697,79 +534,29 @@ namespace internal
                 }
               if (integration_flag & EvaluationFlags::gradients)
                 {
-                  eval2.template gradients<2, false, false>(gradients_quad +
-                                                              2 * n_q_points,
-                                                            temp1);
+                  eval2.template gradients<2, false, false, 3>(gradients_quad +
+                                                                 2,
+                                                               temp1);
                   if (integration_flag & EvaluationFlags::values)
                     eval2.template values<2, false, true>(values_quad, temp1);
                   eval1.template values<1, false, false>(temp1, temp2);
-                  eval2.template values<2, false, false>(gradients_quad +
-                                                           n_q_points,
-                                                         temp1);
+                  eval2.template values<2, false, false, 3>(gradients_quad + 1,
+                                                            temp1);
                   eval1.template gradients<1, false, true>(temp1, temp2);
                   if (add_into_values_array == false)
                     eval0.template values<0, false, false>(temp2, values_dofs);
                   else
                     eval0.template values<0, false, true>(temp2, values_dofs);
-                  eval2.template values<2, false, false>(gradients_quad, temp1);
-                  eval1.template values<1, false, false>(temp1, temp2);
-                  eval0.template gradients<0, false, true>(temp2, values_dofs);
-                }
-              if ((integration_flag & EvaluationFlags::hessians) != 0u)
-                {
-                  // grad xx
-                  eval2.template values<2, false, false>(hessians_quad, temp1);
-                  eval1.template values<1, false, false>(temp1, temp2);
-
-                  if ((integration_flag & EvaluationFlags::values) != 0u ||
-                      (integration_flag & EvaluationFlags::gradients) != 0u ||
-                      add_into_values_array == true)
-                    eval0.template hessians<0, false, true>(temp2, values_dofs);
-                  else
-                    eval0.template hessians<0, false, false>(temp2,
-                                                             values_dofs);
-
-                  // grad yy
-                  eval2.template values<2, false, false>(hessians_quad +
-                                                           n_q_points,
-                                                         temp1);
-                  eval1.template hessians<1, false, false>(temp1, temp2);
-                  eval0.template values<0, false, true>(temp2, values_dofs);
-
-                  // grad zz
-                  eval2.template hessians<2, false, false>(hessians_quad +
-                                                             2 * n_q_points,
-                                                           temp1);
-                  eval1.template values<1, false, false>(temp1, temp2);
-                  eval0.template values<0, false, true>(temp2, values_dofs);
-
-                  // grad xy
-                  eval2.template values<2, false, false>(hessians_quad +
-                                                           3 * n_q_points,
-                                                         temp1);
-                  eval1.template gradients<1, false, false>(temp1, temp2);
-                  eval0.template gradients<0, false, true>(temp2, values_dofs);
-
-                  // grad xz
-                  eval2.template gradients<2, false, false>(hessians_quad +
-                                                              4 * n_q_points,
+                  eval2.template values<2, false, false, 3>(gradients_quad,
                                                             temp1);
                   eval1.template values<1, false, false>(temp1, temp2);
                   eval0.template gradients<0, false, true>(temp2, values_dofs);
-
-                  // grad yz
-                  eval2.template gradients<2, false, false>(hessians_quad +
-                                                              5 * n_q_points,
-                                                            temp1);
-                  eval1.template gradients<1, false, false>(temp1, temp2);
-                  eval0.template values<0, false, true>(temp2, values_dofs);
                 }
 
               // advance to the next component in 1d array
               values_dofs += dofs_per_comp;
               values_quad += n_q_points;
               gradients_quad += 3 * n_q_points;
-              hessians_quad += 6 * n_q_points;
             }
           break;
 
@@ -885,11 +672,10 @@ namespace internal
                           n_dofs,
                           n_q_points);
 
-                eval.template gradients<0, true, false>(values_dofs_actual_ptr,
-                                                        gradients_quad_ptr);
-
-                gradients_quad_ptr += 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;
           }
       }
@@ -965,465 +751,19 @@ namespace internal
                 if ((add_into_values_array == false &&
                      !(integration_flag & EvaluationFlags::values)) &&
                     d == 0)
-                  eval.template gradients<0, false, false>(
-                    gradients_quad_ptr, values_dofs_actual_ptr);
+                  eval.template gradients<0, false, false, dim>(
+                    gradients_quad_ptr + d, values_dofs_actual_ptr);
                 else
-                  eval.template gradients<0, false, true>(
-                    gradients_quad_ptr, values_dofs_actual_ptr);
-
-                gradients_quad_ptr += n_q_points;
+                  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;
           }
       }
   }
 
 
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  template <bool integrate>
-  inline void
-  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
-                   dim,
-                   fe_degree,
-                   n_q_points_1d,
-                   Number>::
-    evaluate_or_integrate(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      Number                                *values_dofs_actual,
-      FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array)
-  {
-    if (evaluation_flag == EvaluationFlags::nothing)
-      return;
-
-    AssertDimension(fe_eval.get_shape_info().data.size(), 2);
-    // First component:
-    evaluate_tensor_product_per_component<0>(
-      evaluation_flag,
-      values_dofs_actual,
-      fe_eval,
-      add_into_values_array,
-      std::integral_constant<bool, integrate>());
-    // Second component :
-    evaluate_tensor_product_per_component<1>(
-      evaluation_flag,
-      values_dofs_actual,
-      fe_eval,
-      add_into_values_array,
-      std::integral_constant<bool, integrate>());
-    if (dim == 3)
-      {
-        // Third component
-        evaluate_tensor_product_per_component<2>(
-          evaluation_flag,
-          values_dofs_actual,
-          fe_eval,
-          add_into_values_array,
-          std::integral_constant<bool, integrate>());
-      }
-  }
-
-  // Helper function that applies the 1d evaluation kernels.
-  // std::integral_constant<bool, false> is the interpolation path, and
-  // std::integral_constant<bool, true> below is the integration path.
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  template <int normal_dir>
-  inline void
-  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
-                   dim,
-                   fe_degree,
-                   n_q_points_1d,
-                   Number>::
-    evaluate_tensor_product_per_component(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      Number                                *values_dofs_actual,
-      FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array,
-      std::integral_constant<bool, false>)
-  {
-    (void)add_into_values_array;
-
-    using Number2 =
-      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-    using EvalNormal =
-      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                        dim,
-                                        (fe_degree == -1) ? 1 : fe_degree + 1,
-                                        n_q_points_1d,
-                                        normal_dir,
-                                        Number,
-                                        Number2>;
-
-    using EvalTangent =
-      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                        dim,
-                                        (fe_degree == -1) ? 1 : fe_degree,
-                                        n_q_points_1d,
-                                        normal_dir,
-                                        Number,
-                                        Number2>;
-    using Eval0 = std::conditional_t<normal_dir == 0, EvalNormal, EvalTangent>;
-    using Eval1 = std::conditional_t<normal_dir == 1, EvalNormal, EvalTangent>;
-    using Eval2 = std::conditional_t<normal_dir == 2, EvalNormal, EvalTangent>;
-
-    const auto &shape_info = fe_eval.get_shape_info();
-    Eval0       eval0      = create_evaluator_tensor_product<Eval0>(
-      ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
-    Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
-      ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
-    Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
-      ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
-
-    Number *temp1 = fe_eval.get_scratch_data().begin();
-    Number *temp2;
-
-    temp2 =
-      temp1 +
-      std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
-               Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
-
-    const std::size_t n_q_points    = shape_info.n_q_points;
-    const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
-
-    // Initial shift depending on component (normal_dir)
-    Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
-    Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
-    Number *gradients_quad =
-      fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
-    Number *hessians_quad =
-      (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
-                   fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
-
-    switch (dim)
-      {
-        case 2:
-          if (evaluation_flag & EvaluationFlags::gradients)
-            {
-              eval0.template gradients<0, true, false>(values_dofs, temp1);
-              eval1.template values<1, true, false>(temp1, gradients_quad);
-            }
-          if (evaluation_flag & EvaluationFlags::hessians)
-            {
-              // The evaluation/integration here *should* work, however
-              // the piola transform is not implemented.
-              AssertThrow(false, ExcNotImplemented());
-              // grad xy
-              if (!(evaluation_flag & EvaluationFlags::gradients))
-                eval0.template gradients<0, true, false>(values_dofs, temp1);
-              eval1.template gradients<1, true, false>(temp1,
-                                                       hessians_quad +
-                                                         2 * n_q_points);
-
-              // grad xx
-              eval0.template hessians<0, true, false>(values_dofs, temp1);
-              eval1.template values<1, true, false>(temp1, hessians_quad);
-            }
-
-          // grad y
-          eval0.template values<0, true, false>(values_dofs, temp1);
-          if (evaluation_flag & EvaluationFlags::gradients)
-            eval1.template gradients<1, true, false>(temp1,
-                                                     gradients_quad +
-                                                       n_q_points);
-
-          // grad yy
-          if (evaluation_flag & EvaluationFlags::hessians)
-            eval1.template hessians<1, true, false>(temp1,
-                                                    hessians_quad + n_q_points);
-
-          // val: can use values applied in x
-          if (evaluation_flag & EvaluationFlags::values)
-            eval1.template values<1, true, false>(temp1, values_quad);
-          break;
-        case 3:
-          if (evaluation_flag & EvaluationFlags::gradients)
-            {
-              // grad x
-              eval0.template gradients<0, true, false>(values_dofs, temp1);
-              eval1.template values<1, true, false>(temp1, temp2);
-              eval2.template values<2, true, false>(temp2, gradients_quad);
-            }
-
-          if (evaluation_flag & EvaluationFlags::hessians)
-            {
-              // The evaluation/integration here *should* work, however
-              // the piola transform is not implemented.
-              AssertThrow(false, ExcNotImplemented());
-              // grad xz
-              if (!(evaluation_flag & EvaluationFlags::gradients))
-                {
-                  eval0.template gradients<0, true, false>(values_dofs, temp1);
-                  eval1.template values<1, true, false>(temp1, temp2);
-                }
-              eval2.template gradients<2, true, false>(temp2,
-                                                       hessians_quad +
-                                                         4 * n_q_points);
-
-              // grad xy
-              eval1.template gradients<1, true, false>(temp1, temp2);
-              eval2.template values<2, true, false>(temp2,
-                                                    hessians_quad +
-                                                      3 * n_q_points);
-
-              // grad xx
-              eval0.template hessians<0, true, false>(values_dofs, temp1);
-              eval1.template values<1, true, false>(temp1, temp2);
-              eval2.template values<2, true, false>(temp2, hessians_quad);
-            }
-
-          // grad y
-          eval0.template values<0, true, false>(values_dofs, temp1);
-          if (evaluation_flag & EvaluationFlags::gradients)
-            {
-              eval1.template gradients<1, true, false>(temp1, temp2);
-              eval2.template values<2, true, false>(temp2,
-                                                    gradients_quad +
-                                                      n_q_points);
-            }
-
-          if (evaluation_flag & EvaluationFlags::hessians)
-            {
-              // grad yz
-              if (!(evaluation_flag & EvaluationFlags::gradients))
-                eval1.template gradients<1, true, false>(temp1, temp2);
-              eval2.template gradients<2, true, false>(temp2,
-                                                       hessians_quad +
-                                                         5 * n_q_points);
-
-              // grad yy
-              eval1.template hessians<1, true, false>(temp1, temp2);
-              eval2.template values<2, true, false>(temp2,
-                                                    hessians_quad + n_q_points);
-            }
-
-          // grad z: can use the values applied in x direction stored in
-          // temp1
-          eval1.template values<1, true, false>(temp1, temp2);
-          if (evaluation_flag & EvaluationFlags::gradients)
-            eval2.template gradients<2, true, false>(temp2,
-                                                     gradients_quad +
-                                                       2 * n_q_points);
-
-          // grad zz: can use the values applied in x and y direction stored
-          // in temp2
-          if (evaluation_flag & EvaluationFlags::hessians)
-            eval2.template hessians<2, true, false>(temp2,
-                                                    hessians_quad +
-                                                      2 * n_q_points);
-
-          // val: can use the values applied in x & y direction stored in
-          // temp2
-          if (evaluation_flag & EvaluationFlags::values)
-            eval2.template values<2, true, false>(temp2, values_quad);
-          break;
-        default:
-          AssertThrow(false, ExcNotImplemented());
-      }
-  }
-
-  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  template <int normal_dir>
-  inline void
-  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
-                   dim,
-                   fe_degree,
-                   n_q_points_1d,
-                   Number>::
-    evaluate_tensor_product_per_component(
-      const EvaluationFlags::EvaluationFlags evaluation_flag,
-      Number                                *values_dofs_actual,
-      FEEvaluationData<dim, Number, false>  &fe_eval,
-      const bool                             add_into_values_array,
-      std::integral_constant<bool, true>)
-  {
-    using Number2 =
-      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-    using EvalNormal =
-      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                        dim,
-                                        (fe_degree == -1) ? 1 : fe_degree + 1,
-                                        n_q_points_1d,
-                                        normal_dir,
-                                        Number,
-                                        Number2>;
-
-    using EvalTangent =
-      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                        dim,
-                                        (fe_degree == -1) ? 1 : fe_degree,
-                                        n_q_points_1d,
-                                        normal_dir,
-                                        Number,
-                                        Number2>;
-    using Eval0 = std::conditional_t<normal_dir == 0, EvalNormal, EvalTangent>;
-    using Eval1 = std::conditional_t<normal_dir == 1, EvalNormal, EvalTangent>;
-    using Eval2 = std::conditional_t<normal_dir == 2, EvalNormal, EvalTangent>;
-
-    const auto &shape_info = fe_eval.get_shape_info();
-    Eval0       eval0      = create_evaluator_tensor_product<Eval0>(
-      ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
-    Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
-      ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
-    Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
-      ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
-
-    Number *temp1 = fe_eval.get_scratch_data().begin();
-    Number *temp2;
-
-    temp2 =
-      temp1 +
-      std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
-               Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
-
-    const std::size_t n_q_points    = shape_info.n_q_points;
-    const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
-
-    // Initial shift depending on component (normal_dir)
-    Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
-    Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
-    Number *gradients_quad =
-      fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
-    Number *hessians_quad =
-      (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
-                   fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
-
-    // Integrate path
-    switch (dim)
-      {
-        case 2:
-          if ((evaluation_flag & EvaluationFlags::values) &&
-              !(evaluation_flag & EvaluationFlags::gradients))
-            {
-              eval1.template values<1, false, false>(values_quad, temp1);
-              if (add_into_values_array == false)
-                eval0.template values<0, false, false>(temp1, values_dofs);
-              else
-                eval0.template values<0, false, true>(temp1, values_dofs);
-            }
-          if (evaluation_flag & EvaluationFlags::gradients)
-            {
-              eval1.template gradients<1, false, false>(gradients_quad +
-                                                          n_q_points,
-                                                        temp1);
-              if ((evaluation_flag & EvaluationFlags::values))
-                eval1.template values<1, false, true>(values_quad, temp1);
-              if (add_into_values_array == false)
-                eval0.template values<0, false, false>(temp1, values_dofs);
-              else
-                eval0.template values<0, false, true>(temp1, values_dofs);
-              eval1.template values<1, false, false>(gradients_quad, temp1);
-              eval0.template gradients<0, false, true>(temp1, values_dofs);
-            }
-          if (evaluation_flag & EvaluationFlags::hessians)
-            {
-              // grad xx
-              eval1.template values<1, false, false>(hessians_quad, temp1);
-
-              if ((evaluation_flag & EvaluationFlags::values) ||
-                  (evaluation_flag & EvaluationFlags::gradients) ||
-                  add_into_values_array == true)
-                eval0.template hessians<0, false, true>(temp1, values_dofs);
-              else
-                eval0.template hessians<0, false, false>(temp1, values_dofs);
-
-              // grad yy
-              eval1.template hessians<1, false, false>(hessians_quad +
-                                                         n_q_points,
-                                                       temp1);
-              eval0.template values<0, false, true>(temp1, values_dofs);
-
-              // grad xy
-              eval1.template gradients<1, false, false>(hessians_quad +
-                                                          2 * n_q_points,
-                                                        temp1);
-              eval0.template gradients<0, false, true>(temp1, values_dofs);
-            }
-          break;
-
-        case 3:
-          if ((evaluation_flag & EvaluationFlags::values) &&
-              !(evaluation_flag & EvaluationFlags::gradients))
-            {
-              eval2.template values<2, false, false>(values_quad, temp1);
-              eval1.template values<1, false, false>(temp1, temp2);
-              if (add_into_values_array == false)
-                eval0.template values<0, false, false>(temp2, values_dofs);
-              else
-                eval0.template values<0, false, true>(temp2, values_dofs);
-            }
-          if (evaluation_flag & EvaluationFlags::gradients)
-            {
-              eval2.template gradients<2, false, false>(gradients_quad +
-                                                          2 * n_q_points,
-                                                        temp1);
-              if ((evaluation_flag & EvaluationFlags::values))
-                eval2.template values<2, false, true>(values_quad, temp1);
-              eval1.template values<1, false, false>(temp1, temp2);
-              eval2.template values<2, false, false>(gradients_quad +
-                                                       n_q_points,
-                                                     temp1);
-              eval1.template gradients<1, false, true>(temp1, temp2);
-              if (add_into_values_array == false)
-                eval0.template values<0, false, false>(temp2, values_dofs);
-              else
-                eval0.template values<0, false, true>(temp2, values_dofs);
-              eval2.template values<2, false, false>(gradients_quad, temp1);
-              eval1.template values<1, false, false>(temp1, temp2);
-              eval0.template gradients<0, false, true>(temp2, values_dofs);
-            }
-          if (evaluation_flag & EvaluationFlags::hessians)
-            {
-              // grad xx
-              eval2.template values<2, false, false>(hessians_quad, temp1);
-              eval1.template values<1, false, false>(temp1, temp2);
-
-              if ((evaluation_flag & EvaluationFlags::values) ||
-                  (evaluation_flag & EvaluationFlags::gradients) ||
-                  add_into_values_array == true)
-                eval0.template hessians<0, false, true>(temp2, values_dofs);
-              else
-                eval0.template hessians<0, false, false>(temp2, values_dofs);
-
-              // grad yy
-              eval2.template values<2, false, false>(hessians_quad + n_q_points,
-                                                     temp1);
-              eval1.template hessians<1, false, false>(temp1, temp2);
-              eval0.template values<0, false, true>(temp2, values_dofs);
-
-              // grad zz
-              eval2.template hessians<2, false, false>(hessians_quad +
-                                                         2 * n_q_points,
-                                                       temp1);
-              eval1.template values<1, false, false>(temp1, temp2);
-              eval0.template values<0, false, true>(temp2, values_dofs);
-
-              // grad xy
-              eval2.template values<2, false, false>(hessians_quad +
-                                                       3 * n_q_points,
-                                                     temp1);
-              eval1.template gradients<1, false, false>(temp1, temp2);
-              eval0.template gradients<0, false, true>(temp2, values_dofs);
-
-              // grad xz
-              eval2.template gradients<2, false, false>(hessians_quad +
-                                                          4 * n_q_points,
-                                                        temp1);
-              eval1.template values<1, false, false>(temp1, temp2);
-              eval0.template gradients<0, false, true>(temp2, values_dofs);
-
-              // grad yz
-              eval2.template gradients<2, false, false>(hessians_quad +
-                                                          5 * n_q_points,
-                                                        temp1);
-              eval1.template gradients<1, false, false>(temp1, temp2);
-              eval0.template values<0, false, true>(temp2, values_dofs);
-            }
-
-          break;
-        default:
-          AssertThrow(false, ExcNotImplemented());
-      }
-  }
 
   /**
    * This struct implements the change between two different bases. This is an
@@ -1433,22 +773,12 @@ namespace internal
    *
    * This class allows for dimension-independent application of the operation,
    * implemented by template recursion. It has been tested up to 6d.
-   *
-   * The last two template arguments in this class are unused. They have a
-   * default type set to bool, but any type will compile. These arguments are
-   * present for backward compatibility of this (internal) interface with
-   * previous versions of deal.II, where the template arguments were used to
-   * indicate the number types These are now part of the do_forward() and
-   * do_backward() functions, and they will be removed from a future version
-   * of deal.II.
    */
   template <EvaluatorVariant  variant,
             EvaluatorQuantity quantity,
             int               dim,
             int               basis_size_1,
-            int               basis_size_2,
-            typename = bool,
-            typename = bool>
+            int               basis_size_2>
   struct FEEvaluationImplBasisChange
   {
     static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
@@ -1827,279 +1157,556 @@ namespace internal
 
 
   /**
-   * This struct performs the evaluation of function values, gradients and
-   * Hessians for tensor-product finite elements. This a specialization for
-   * elements where the nodal points coincide with the quadrature points like
-   * FE_Q shape functions on Gauss-Lobatto elements integrated with
-   * Gauss-Lobatto quadrature. The assumption of this class is that the shape
-   * 'values' operation is identity, which allows us to write shorter code.
-   *
-   * In literature, this form of evaluation is often called spectral
-   * evaluation, spectral collocation or simply collocation, meaning the same
-   * location for shape functions and evaluation space (quadrature points).
+   * Internal function that evaluates the gradients of finite element
+   * functions represented by bases in the collocation space, used by
+   * FEEvaluationImplCollocation and FEEvaluationImplTransformToCollocation.
+   * The evaluation strategy uses sum factorization with the even-odd
+   * optimization and fixed loop bounds.
    */
-  template <int dim, int fe_degree, typename Number>
-  struct FEEvaluationImplCollocation
+  template <int n_points_1d, int dim, typename Number, typename Number2>
+  inline void
+  evaluate_gradients_collocation(
+    const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
+    const Number                                            *values,
+    Number                                                  *gradients)
   {
-    using Number2 =
-      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-    using Eval = EvaluatorTensorProduct<evaluate_evenodd,
-                                        dim,
-                                        fe_degree + 1,
-                                        fe_degree + 1,
-                                        Number,
-                                        Number2>;
-
-    static void
-    evaluate(const unsigned int                     n_components,
-             const EvaluationFlags::EvaluationFlags evaluation_flag,
-             const Number                          *values_dofs,
-             FEEvaluationData<dim, Number, false>  &fe_eval);
-
-    static void
-    do_evaluate(const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
-                const EvaluationFlags::EvaluationFlags evaluation_flag,
-                const Number                          *values_dofs,
-                Number                                *gradients_quad,
-                Number                                *hessians_quad);
-
-    static void
-    integrate(const unsigned int                     n_components,
-              const EvaluationFlags::EvaluationFlags integration_flag,
-              Number                                *values_dofs,
-              FEEvaluationData<dim, Number, false>  &fe_eval,
-              const bool                             add_into_values_array);
+    AssertDimension(shape.shape_gradients_collocation_eo.size(),
+                    (n_points_1d + 1) / 2 * n_points_1d);
 
-    static void
-    do_integrate(const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
-                 const EvaluationFlags::EvaluationFlags integration_flag,
-                 Number                                *values_dofs,
-                 Number                                *gradients_quad,
-                 const Number                          *hessians_quad,
-                 const bool                             add_into_values_array);
-  };
+    EvaluatorTensorProduct<evaluate_evenodd,
+                           dim,
+                           n_points_1d,
+                           n_points_1d,
+                           Number,
+                           Number2>
+      eval({}, shape.shape_gradients_collocation_eo, {});
+    EvaluatorTensorProduct<evaluate_evenodd,
+                           2,
+                           n_points_1d,
+                           n_points_1d,
+                           Number,
+                           Number2>
+      eval_2d({}, shape.shape_gradients_collocation_eo, {});
+
+    if (dim == 1)
+      eval.template gradients<0, true, false>(values, gradients);
+    else
+      {
+        if (dim > 2)
+          eval.template gradients<2, true, false, dim>(values, gradients + 2);
+        constexpr unsigned int loop_bound  = (dim > 2 ? n_points_1d : 1);
+        constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
+        const Number          *in = values + (loop_bound - 1) * n_points_2d;
+        Number *out = gradients + (loop_bound - 1) * dim * n_points_2d;
+        for (unsigned int l = 0; l < loop_bound; ++l)
+          {
+            eval_2d.template gradients<0, true, false, dim>(in, out);
+            eval_2d.template gradients<1, true, false, dim>(in, out + 1);
+            in -= n_points_2d;
+            out -= dim * n_points_2d;
+          }
+      }
+  }
 
 
 
-  template <int dim, int fe_degree, typename Number>
+  /**
+   * Internal function that multiplies by the gradients of test functions and
+   * sums over quadrature points for function representations in the
+   * collocation space, used by FEEvaluationImplCollocation and
+   * FEEvaluationImplTransformToCollocation. The evaluation strategy uses sum
+   * factorization with the even-odd optimization and fixed loop bounds.
+   */
+  template <int n_points_1d, int dim, typename Number, typename Number2>
   inline void
-  FEEvaluationImplCollocation<dim, fe_degree, Number>::evaluate(
-    const unsigned int                     n_components,
-    const EvaluationFlags::EvaluationFlags evaluation_flag,
-    const Number                          *values_dofs,
-    FEEvaluationData<dim, Number, false>  &fe_eval)
+  integrate_gradients_collocation(
+    const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
+    Number                                                  *values,
+    const Number                                            *gradients,
+    const bool add_into_values_array)
   {
-    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+    AssertDimension(shape.shape_gradients_collocation_eo.size(),
+                    (n_points_1d + 1) / 2 * n_points_1d);
 
-    for (unsigned int c = 0; c < n_components; ++c)
+    EvaluatorTensorProduct<evaluate_evenodd,
+                           dim,
+                           n_points_1d,
+                           n_points_1d,
+                           Number,
+                           Number2>
+      eval({}, shape.shape_gradients_collocation_eo, {});
+    EvaluatorTensorProduct<evaluate_evenodd,
+                           2,
+                           n_points_1d,
+                           n_points_1d,
+                           Number,
+                           Number2>
+      eval_2d({}, shape.shape_gradients_collocation_eo, {});
+
+    if (dim == 1)
+      {
+        if (add_into_values_array)
+          eval.template gradients<0, false, true>(gradients, values);
+        else
+          eval.template gradients<0, false, false>(gradients, values);
+      }
+    else
       {
-        if ((evaluation_flag & EvaluationFlags::values) != 0u)
-          for (unsigned int i = 0; i < n_points; ++i)
-            fe_eval.begin_values()[n_points * c + i] =
-              values_dofs[n_points * c + i];
-
-        do_evaluate(fe_eval.get_shape_info().data.front(),
-                    evaluation_flag,
-                    values_dofs + c * n_points,
-                    fe_eval.begin_gradients() + c * dim * n_points,
-                    fe_eval.begin_hessians() +
-                      c * dim * (dim + 1) / 2 * n_points);
+        constexpr unsigned int loop_bound  = (dim > 2 ? n_points_1d : 1);
+        constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
+
+        const Number *in  = gradients + (loop_bound - 1) * dim * n_points_2d;
+        Number       *out = values + (loop_bound - 1) * n_points_2d;
+        for (unsigned int l = 0; l < loop_bound; ++l)
+          {
+            if (add_into_values_array)
+              eval_2d.template gradients<0, false, true, dim>(in, out);
+            else
+              eval_2d.template gradients<0, false, false, dim>(in, out);
+            eval_2d.template gradients<1, false, true, dim>(in + 1, out);
+            in -= dim * n_points_2d;
+            out -= n_points_2d;
+          }
       }
+    if (dim > 2)
+      eval.template gradients<2, false, true, dim>(gradients + 2, values);
   }
 
 
 
-  template <int dim, int fe_degree, typename Number>
+  /**
+   * Internal function that evaluates the Hessians of finite element functions
+   * represented by bases in the collocation space, used by
+   * FEEvaluationImplSelector. The evaluation strategy uses sum
+   * factorization with fixed loop bounds.
+   */
+  template <int n_points_1d, int dim, typename Number>
   inline void
-  FEEvaluationImplCollocation<dim, fe_degree, Number>::do_evaluate(
-    const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
-    const EvaluationFlags::EvaluationFlags                   evaluation_flag,
-    const Number                                            *values_dofs,
-    Number                                                  *gradients_quad,
-    Number                                                  *hessians_quad)
+  evaluate_hessians_collocation(const unsigned int n_components,
+                                FEEvaluationData<dim, Number, false> &fe_eval)
   {
-    AssertDimension(shape.shape_gradients_collocation_eo.size(),
-                    (fe_degree + 2) / 2 * (fe_degree + 1));
-    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
-
-    Eval eval({},
-              shape.shape_gradients_collocation_eo,
-              shape.shape_hessians_collocation_eo);
-    if ((evaluation_flag &
-         (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u)
+    using Number2 =
+      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
+
+    // might have non-symmetric quadrature formula, so use the more
+    // conservative 'evaluate_general' scheme rather than 'even_odd' as the
+    // Hessians are not used very often
+    const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+      fe_eval.get_shape_info().data[0];
+    AssertDimension(data.shape_gradients_collocation.size(),
+                    data.n_q_points_1d * data.n_q_points_1d);
+    EvaluatorTensorProduct<evaluate_general,
+                           dim,
+                           n_points_1d,
+                           n_points_1d,
+                           Number,
+                           Number2>
+      eval({},
+           data.shape_gradients_collocation.data(),
+           data.shape_hessians_collocation.data(),
+           data.n_q_points_1d,
+           data.n_q_points_1d);
+
+    const Number     *values   = fe_eval.begin_values();
+    Number           *hessians = fe_eval.begin_hessians();
+    Number           *scratch  = fe_eval.get_scratch_data().begin();
+    const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
+    for (unsigned int comp = 0; comp < n_components; ++comp)
       {
-        eval.template gradients<0, true, false>(values_dofs, gradients_quad);
+        // xx derivative
+        eval.template hessians<0, true, false>(values, hessians);
         if (dim > 1)
-          eval.template gradients<1, true, false>(values_dofs,
-                                                  gradients_quad + n_points);
+          {
+            // xy derivative: we might or might not have the gradients already
+            // computed elsewhere, but we recompute them here since it adds
+            // only moderate extra work (at most 25%)
+            eval.template gradients<0, true, false>(values, scratch);
+            eval.template gradients<1, true, false>(scratch,
+                                                    hessians + dim * n_points);
+            // yy derivative
+            eval.template hessians<1, true, false>(values, hessians + n_points);
+          }
         if (dim > 2)
-          eval.template gradients<2, true, false>(values_dofs,
-                                                  gradients_quad +
-                                                    2 * n_points);
+          {
+            // xz derivative
+            eval.template gradients<2, true, false>(scratch,
+                                                    hessians + 4 * n_points);
+            // yz derivative
+            eval.template gradients<1, true, false>(values, scratch);
+            eval.template gradients<2, true, false>(scratch,
+                                                    hessians + 5 * n_points);
+            // zz derivative
+            eval.template hessians<2, true, false>(values,
+                                                   hessians + 2 * n_points);
+          }
+
+        values += n_points;
+        hessians += (dim * (dim + 1)) / 2 * n_points;
       }
-    if (evaluation_flag & EvaluationFlags::hessians)
+  }
+
+
+
+  /**
+   * Internal function that multiplies by the Hessians of test functions and
+   * sums over quadrature points for function representations in the
+   * collocation space, used by FEEvaluationImplSelector. The evaluation
+   * strategy uses sum factorization with fixed loop bounds.
+   */
+  template <int n_q_points_1d, int dim, typename Number>
+  inline void
+  integrate_hessians_collocation(const unsigned int n_components,
+                                 FEEvaluationData<dim, Number, false> &fe_eval,
+                                 const bool add_into_values_array)
+  {
+    using Number2 =
+      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
+
+    const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+      fe_eval.get_shape_info().data[0];
+    AssertDimension(data.shape_gradients_collocation.size(),
+                    data.n_q_points_1d * data.n_q_points_1d);
+    EvaluatorTensorProduct<evaluate_general,
+                           dim,
+                           n_q_points_1d,
+                           n_q_points_1d,
+                           Number,
+                           Number2>
+                      eval({},
+           data.shape_gradients_collocation.data(),
+           data.shape_hessians_collocation.data(),
+           data.n_q_points_1d,
+           data.n_q_points_1d);
+    Number           *values   = fe_eval.begin_values();
+    const Number     *hessians = fe_eval.begin_hessians();
+    Number           *scratch  = fe_eval.get_scratch_data().begin();
+    const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
+
+    for (unsigned int comp = 0; comp < n_components; ++comp)
       {
-        eval.template hessians<0, true, false>(values_dofs, hessians_quad);
+        // xx derivative
+        if (add_into_values_array == true)
+          eval.template hessians<0, false, true>(hessians, values);
+        else
+          eval.template hessians<0, false, false>(hessians, values);
+
+        // yy derivative
         if (dim > 1)
+          eval.template hessians<1, false, true>(hessians + n_points, values);
+        if (dim > 2)
           {
-            eval.template gradients<1, true, false>(gradients_quad,
-                                                    hessians_quad +
-                                                      dim * n_points);
-            eval.template hessians<1, true, false>(values_dofs,
-                                                   hessians_quad + n_points);
+            // zz derivative
+            eval.template hessians<2, false, true>(hessians + 2 * n_points,
+                                                   values);
+            // yz derivative
+            eval.template gradients<2, false, false>(hessians + 5 * n_points,
+                                                     scratch);
+            eval.template gradients<1, false, true>(scratch, values);
+
+            // xz derivative
+            eval.template gradients<2, false, false>(hessians + 4 * n_points,
+                                                     scratch);
           }
-        if (dim > 2)
+
+        if (dim > 1)
           {
-            eval.template gradients<2, true, false>(gradients_quad,
-                                                    hessians_quad +
-                                                      4 * n_points);
-            eval.template gradients<2, true, false>(gradients_quad + n_points,
-                                                    hessians_quad +
-                                                      5 * n_points);
-            eval.template hessians<2, true, false>(values_dofs,
-                                                   hessians_quad +
-                                                     2 * n_points);
+            // xy derivative, combined with xz in 3d
+            eval.template gradients<1, false, (dim > 2)>(hessians +
+                                                           dim * n_points,
+                                                         scratch);
+            eval.template gradients<0, false, true>(scratch, values);
           }
+
+        values += n_points;
+        hessians += (dim * (dim + 1)) / 2 * n_points;
       }
   }
 
 
 
-  template <int dim, int fe_degree, typename Number>
-  inline void
-  FEEvaluationImplCollocation<dim, fe_degree, Number>::integrate(
-    const unsigned int                     n_components,
-    const EvaluationFlags::EvaluationFlags integration_flag,
-    Number                                *values_dofs,
-    FEEvaluationData<dim, Number, false>  &fe_eval,
-    const bool                             add_into_values_array)
+  /**
+   * Internal function to evaluate the Hessians of finite element functions in
+   * the non-collocation setting as a fall-back. The evaluation strategy uses
+   * sum factorization with run-time loop bounds and is thus slower than the
+   * collocation case, but it is not as widely used and thus uncritical.
+   */
+  template <int dim, typename Number>
+  void
+  evaluate_hessians_slow(const unsigned int                    n_components,
+                         const Number                         *values_dofs,
+                         FEEvaluationData<dim, Number, false> &fe_eval)
+  {
+    const auto &univariate_shape_data = fe_eval.get_shape_info().data;
+    using Impl =
+      FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
+    using Eval = typename Impl::Eval;
+    Eval eval0 =
+      Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
+    Eval eval1 = Impl::create_evaluator_tensor_product(
+      &univariate_shape_data[std::min<int>(1,
+                                           univariate_shape_data.size() - 1)]);
+    Eval eval2 = Impl::create_evaluator_tensor_product(
+      &univariate_shape_data[std::min<int>(2,
+                                           univariate_shape_data.size() - 1)]);
+
+    const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
+    Number            *tmp1     = fe_eval.get_scratch_data().begin();
+    Number            *tmp2 =
+      tmp1 + std::max(Utilities::fixed_power<dim>(
+                        univariate_shape_data.front().fe_degree + 1),
+                      Utilities::fixed_power<dim>(
+                        univariate_shape_data.front().n_q_points_1d));
+    Number *hessians = fe_eval.begin_hessians();
+
+    for (unsigned int comp = 0; comp < n_components;
+         ++comp,
+                      hessians += n_points * dim * (dim + 1) / 2,
+                      values_dofs +=
+                      fe_eval.get_shape_info().dofs_per_component_on_cell)
+      switch (dim)
+        {
+          case 1:
+            eval0.template hessians<0, true, false>(values_dofs, hessians);
+            break;
+          case 2:
+            // xx derivative
+            eval0.template hessians<0, true, false>(values_dofs, tmp1);
+            eval1.template values<1, true, false>(tmp1, hessians);
+            // xy derivative
+            eval0.template gradients<0, true, false>(values_dofs, tmp1);
+            eval1.template gradients<1, true, false>(tmp1,
+                                                     hessians + 2 * n_points);
+            // yy derivative
+            eval0.template values<0, true, false>(values_dofs, tmp1);
+            eval1.template hessians<1, true, false>(tmp1, hessians + n_points);
+            break;
+          case 3:
+            // xx derivative
+            eval0.template hessians<0, true, false>(values_dofs, tmp1);
+            eval1.template values<1, true, false>(tmp1, tmp2);
+            eval2.template values<2, true, false>(tmp2, hessians);
+            // xy derivative
+            eval0.template gradients<0, true, false>(values_dofs, tmp1);
+            eval1.template gradients<1, true, false>(tmp1, tmp2);
+            eval2.template values<2, true, false>(tmp2,
+                                                  hessians + 3 * n_points);
+            // xz derivative
+            eval1.template values<1, true, false>(tmp1, tmp2);
+            eval2.template gradients<2, true, false>(tmp2,
+                                                     hessians + 4 * n_points);
+            // yy derivative
+            eval0.template values<0, true, false>(values_dofs, tmp1);
+            eval1.template hessians<1, true, false>(tmp1, tmp2);
+            eval2.template values<2, true, false>(tmp2, hessians + n_points);
+            // yz derivative
+            eval1.template gradients<1, true, false>(tmp1, tmp2);
+            eval2.template gradients<2, true, false>(tmp2,
+                                                     hessians + 5 * n_points);
+            // zz derivative
+            eval1.template values<1, true, false>(tmp1, tmp2);
+            eval2.template hessians<2, true, false>(tmp2,
+                                                    hessians + 2 * n_points);
+            break;
+
+          default:
+            Assert(false,
+                   ExcNotImplemented(
+                     "Only 1d, 2d and 3d implemented for Hessian"));
+        }
+  }
+
+
+
+  /**
+   * Internal function to multiply by the Hessians of the test functions and
+   * integrate in the non-collocation setting as a fall-back. The evaluation
+   * strategy uses sum factorization with run-time loop bounds and is thus
+   * slower than the collocation case, but it is not as widely used and thus
+   * uncritical.
+   */
+  template <int dim, typename Number>
+  void
+  integrate_hessians_slow(const unsigned int n_components,
+                          const FEEvaluationData<dim, Number, false> &fe_eval,
+                          Number    *values_dofs,
+                          const bool add_into_values_array)
   {
-    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+    const auto &univariate_shape_data = fe_eval.get_shape_info().data;
+    using Impl =
+      FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
+    using Eval = typename Impl::Eval;
+    Eval eval0 =
+      Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
+    Eval eval1 = Impl::create_evaluator_tensor_product(
+      &univariate_shape_data[std::min<int>(1,
+                                           univariate_shape_data.size() - 1)]);
+    Eval eval2 = Impl::create_evaluator_tensor_product(
+      &univariate_shape_data[std::min<int>(2,
+                                           univariate_shape_data.size() - 1)]);
+
+    const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
+    Number            *tmp1     = fe_eval.get_scratch_data().begin();
+    Number            *tmp2 =
+      tmp1 + std::max(Utilities::fixed_power<dim>(
+                        univariate_shape_data.front().fe_degree + 1),
+                      Utilities::fixed_power<dim>(
+                        univariate_shape_data.front().n_q_points_1d));
+    const Number *hessians = fe_eval.begin_hessians();
+
+    for (unsigned int comp = 0; comp < n_components;
+         ++comp,
+                      hessians += n_points * dim * (dim + 1) / 2,
+                      values_dofs +=
+                      fe_eval.get_shape_info().dofs_per_component_on_cell)
+      switch (dim)
+        {
+          case 1:
+            if (add_into_values_array)
+              eval0.template hessians<0, false, true>(hessians, values_dofs);
+            else
+              eval0.template hessians<0, false, false>(hessians, values_dofs);
+            break;
+          case 2:
+            // xx derivative
+            eval1.template values<1, false, false>(hessians, tmp1);
+            if (add_into_values_array)
+              eval0.template hessians<0, false, true>(tmp1, values_dofs);
+            else
+              eval0.template hessians<0, false, false>(tmp1, values_dofs);
+
+            // xy derivative
+            eval1.template gradients<1, false, false>(hessians + 2 * n_points,
+                                                      tmp1);
+            eval0.template gradients<0, false, true>(tmp1, values_dofs);
+            // yy derivative
+            eval1.template hessians<1, false, false>(hessians + n_points, tmp1);
+            eval0.template values<0, false, true>(tmp1, values_dofs);
+            break;
+          case 3:
+            // xx derivative
+            eval2.template values<2, false, false>(hessians, tmp1);
+            eval1.template values<1, false, false>(tmp1, tmp2);
 
-    for (unsigned int c = 0; c < n_components; ++c)
-      {
-        if ((integration_flag & EvaluationFlags::values) != 0u)
-          {
             if (add_into_values_array)
-              for (unsigned int i = 0; i < n_points; ++i)
-                values_dofs[n_points * c + i] +=
-                  fe_eval.begin_values()[n_points * c + i];
+              eval0.template hessians<0, false, true>(tmp2, values_dofs);
             else
-              for (unsigned int i = 0; i < n_points; ++i)
-                values_dofs[n_points * c + i] =
-                  fe_eval.begin_values()[n_points * c + i];
-          }
+              eval0.template hessians<0, false, false>(tmp2, values_dofs);
+
+            // xy derivative
+            eval2.template values<2, false, false>(hessians + 3 * n_points,
+                                                   tmp1);
+            eval1.template gradients<1, false, false>(tmp1, tmp2);
+            // xz derivative
+            eval2.template gradients<2, false, false>(hessians + 4 * n_points,
+                                                      tmp1);
+            eval1.template values<1, false, true>(tmp1, tmp2);
+            eval1.template values<0, false, true>(tmp2, values_dofs);
+
+            // yy derivative
+            eval2.template values<2, false, false>(hessians + n_points, tmp1);
+            eval1.template hessians<1, false, false>(tmp1, tmp2);
+
+            // yz derivative
+            eval2.template gradients<2, false, false>(hessians + 5 * n_points,
+                                                      tmp1);
+            eval1.template gradients<1, false, true>(tmp1, tmp2);
+
+            // zz derivative
+            eval2.template hessians<2, false, false>(hessians + 2 * n_points,
+                                                     tmp1);
+            eval1.template values<1, false, true>(tmp1, tmp2);
+            eval0.template values<0, false, true>(tmp2, values_dofs);
+            break;
 
-        do_integrate(fe_eval.get_shape_info().data.front(),
-                     integration_flag,
-                     values_dofs + c * n_points,
-                     fe_eval.begin_gradients() + c * dim * n_points,
-                     fe_eval.begin_hessians() +
-                       c * dim * (dim + 1) / 2 * n_points,
-                     add_into_values_array ||
-                       ((integration_flag & EvaluationFlags::values) != 0u));
-      }
+          default:
+            Assert(false,
+                   ExcNotImplemented(
+                     "Only 1d, 2d and 3d implemented for Hessian"));
+        }
   }
 
 
 
+  /**
+   * This struct performs the evaluation of function values, gradients and
+   * Hessians for tensor-product finite elements. This a specialization for
+   * elements where the nodal points coincide with the quadrature points like
+   * FE_Q shape functions on Gauss-Lobatto elements integrated with
+   * Gauss-Lobatto quadrature. The assumption of this class is that the shape
+   * 'values' operation is identity, which allows us to write shorter code.
+   *
+   * In literature, this form of evaluation is often called spectral
+   * evaluation, spectral collocation or simply collocation, meaning the same
+   * location for shape functions and evaluation space (quadrature points).
+   */
   template <int dim, int fe_degree, typename Number>
-  inline void
-  FEEvaluationImplCollocation<dim, fe_degree, Number>::do_integrate(
-    const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
-    const EvaluationFlags::EvaluationFlags                   integration_flag,
-    Number                                                  *values_dofs,
-    Number                                                  *gradients_quad,
-    const Number                                            *hessians_quad,
-    const bool add_into_values_array)
+  struct FEEvaluationImplCollocation
   {
-    AssertDimension(shape.shape_gradients_collocation_eo.size(),
-                    (fe_degree + 2) / 2 * (fe_degree + 1));
+    using Number2 =
+      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
+    using Eval = EvaluatorTensorProduct<evaluate_evenodd,
+                                        dim,
+                                        fe_degree + 1,
+                                        fe_degree + 1,
+                                        Number,
+                                        Number2>;
 
-    Eval                  eval({},
-              shape.shape_gradients_collocation_eo,
-              shape.shape_hessians_collocation_eo);
-    constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+    static void
+    evaluate(const unsigned int                     n_components,
+             const EvaluationFlags::EvaluationFlags evaluation_flag,
+             const Number                          *values_dofs,
+             FEEvaluationData<dim, Number, false>  &fe_eval)
+    {
+      constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
 
-    if ((integration_flag & EvaluationFlags::hessians) != 0u)
-      {
-        // diagonal
-        // grad xx
-        if (add_into_values_array == true)
-          eval.template hessians<0, false, true>(hessians_quad, values_dofs);
-        else
-          eval.template hessians<0, false, false>(hessians_quad, values_dofs);
-        // grad yy
-        if (dim > 1)
-          eval.template hessians<1, false, true>(hessians_quad + n_points,
-                                                 values_dofs);
-        // grad zz
-        if (dim > 2)
-          eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
-                                                 values_dofs);
-        // off-diagonal
-        if (dim == 2)
-          {
-            // grad xy, queue into gradient
-            if (integration_flag & EvaluationFlags::gradients)
-              eval.template gradients<1, false, true>(hessians_quad +
-                                                        2 * n_points,
-                                                      gradients_quad);
-            else
-              eval.template gradients<1, false, false>(hessians_quad +
-                                                         2 * n_points,
-                                                       gradients_quad);
-          }
-        if (dim == 3)
-          {
-            // grad xy, queue into gradient
-            if (integration_flag & EvaluationFlags::gradients)
-              eval.template gradients<1, false, true>(hessians_quad +
-                                                        3 * n_points,
-                                                      gradients_quad);
-            else
-              eval.template gradients<1, false, false>(hessians_quad +
-                                                         3 * n_points,
-                                                       gradients_quad);
-
-            // grad xz
-            eval.template gradients<2, false, true>(hessians_quad +
-                                                      4 * n_points,
-                                                    gradients_quad);
-
-            // grad yz
-            if (integration_flag & EvaluationFlags::gradients)
-              eval.template gradients<2, false, true>(
-                hessians_quad + 5 * n_points, gradients_quad + n_points);
-            else
-              eval.template gradients<2, false, false>(
-                hessians_quad + 5 * n_points, gradients_quad + n_points);
-          }
+      for (unsigned int c = 0; c < n_components; ++c)
+        {
+          if ((evaluation_flag & EvaluationFlags::values) != 0u)
+            for (unsigned int i = 0; i < n_points; ++i)
+              fe_eval.begin_values()[n_points * c + i] =
+                values_dofs[n_points * c + i];
+
+          if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+            evaluate_gradients_collocation<fe_degree + 1, dim>(
+              fe_eval.get_shape_info().data.front(),
+              values_dofs + c * n_points,
+              fe_eval.begin_gradients() + c * dim * n_points);
+        }
+    }
 
-        // if we did not integrate gradients, set the last slot to zero
-        // which was not touched before, in order to avoid the if
-        // statement in the gradients loop below
-        if ((integration_flag & EvaluationFlags::gradients) == 0u)
-          for (unsigned int q = 0; q < n_points; ++q)
-            gradients_quad[(dim - 1) * n_points + q] = Number();
-      }
+    static void
+    integrate(const unsigned int                     n_components,
+              const EvaluationFlags::EvaluationFlags integration_flag,
+              Number                                *values_dofs,
+              FEEvaluationData<dim, Number, false>  &fe_eval,
+              const bool                             add_into_values_array)
+    {
+      constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
 
-    if ((integration_flag &
-         (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u)
-      {
-        if (add_into_values_array ||
-            (integration_flag & EvaluationFlags::hessians) != 0u)
-          eval.template gradients<0, false, true>(gradients_quad, values_dofs);
-        else
-          eval.template gradients<0, false, false>(gradients_quad, values_dofs);
-        if (dim > 1)
-          eval.template gradients<1, false, true>(gradients_quad + n_points,
-                                                  values_dofs);
-        if (dim > 2)
-          eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
-                                                  values_dofs);
-      }
-  }
+      for (unsigned int c = 0; c < n_components; ++c)
+        {
+          if ((integration_flag & EvaluationFlags::values) != 0u)
+            {
+              if (add_into_values_array)
+                for (unsigned int i = 0; i < n_points; ++i)
+                  values_dofs[n_points * c + i] +=
+                    fe_eval.begin_values()[n_points * c + i];
+              else
+                for (unsigned int i = 0; i < n_points; ++i)
+                  values_dofs[n_points * c + i] =
+                    fe_eval.begin_values()[n_points * c + i];
+            }
+
+          if ((integration_flag & EvaluationFlags::gradients) != 0u)
+            integrate_gradients_collocation<fe_degree + 1, dim>(
+              fe_eval.get_shape_info().data.front(),
+              values_dofs + c * n_points,
+              fe_eval.begin_gradients() + c * dim * n_points,
+              add_into_values_array ||
+                ((integration_flag & EvaluationFlags::values) != 0u));
+        }
+    }
+  };
 
 
 
@@ -2120,116 +1727,371 @@ namespace internal
     evaluate(const unsigned int                     n_components,
              const EvaluationFlags::EvaluationFlags evaluation_flag,
              const Number                          *values_dofs,
-             FEEvaluationData<dim, Number, false>  &fe_eval);
+             FEEvaluationData<dim, Number, false>  &fe_eval)
+    {
+      const auto &shape_data = fe_eval.get_shape_info().data.front();
+
+      Assert(n_q_points_1d > fe_degree,
+             ExcMessage("You lose information when going to a collocation "
+                        "space of lower degree, so the evaluation results "
+                        "would be wrong. Thus, this class does not permit "
+                        "the chosen operation."));
+      constexpr std::size_t n_dofs     = Utilities::pow(fe_degree + 1, dim);
+      constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+
+      for (unsigned int c = 0; c < n_components; ++c)
+        {
+          FEEvaluationImplBasisChange<
+            evaluate_evenodd,
+            EvaluatorQuantity::value,
+            dim,
+            (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
+            n_q_points_1d>::do_forward(1,
+                                       shape_data.shape_values_eo,
+                                       values_dofs + c * n_dofs,
+                                       fe_eval.begin_values() + c * n_q_points);
+
+          // apply derivatives in the collocation space
+          if (evaluation_flag & EvaluationFlags::gradients)
+            evaluate_gradients_collocation<n_q_points_1d, dim>(
+              shape_data,
+              fe_eval.begin_values() + c * n_q_points,
+              fe_eval.begin_gradients() + c * dim * n_q_points);
+        }
+    }
 
     static void
     integrate(const unsigned int                     n_components,
-              const EvaluationFlags::EvaluationFlags evaluation_flag,
+              const EvaluationFlags::EvaluationFlags integration_flag,
               Number                                *values_dofs,
               FEEvaluationData<dim, Number, false>  &fe_eval,
-              const bool                             add_into_values_array);
+              const bool                             add_into_values_array)
+    {
+      const auto &shape_data = fe_eval.get_shape_info().data.front();
+
+      Assert(n_q_points_1d > fe_degree,
+             ExcMessage("You lose information when going to a collocation "
+                        "space of lower degree, so the evaluation results "
+                        "would be wrong. Thus, this class does not permit "
+                        "the chosen operation."));
+      constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+
+      for (unsigned int c = 0; c < n_components; ++c)
+        {
+          // apply derivatives in collocation space
+          if (integration_flag & EvaluationFlags::gradients)
+            integrate_gradients_collocation<n_q_points_1d, dim>(
+              shape_data,
+              fe_eval.begin_values() + c * n_q_points,
+              fe_eval.begin_gradients() + c * dim * n_q_points,
+              /*add_into_values_array=*/
+              integration_flag & EvaluationFlags::values);
+
+          // transform back to the original space
+          FEEvaluationImplBasisChange<
+            evaluate_evenodd,
+            EvaluatorQuantity::value,
+            dim,
+            (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
+            n_q_points_1d>::do_backward(1,
+                                        shape_data.shape_values_eo,
+                                        add_into_values_array,
+                                        fe_eval.begin_values() + c * n_q_points,
+                                        values_dofs +
+                                          c *
+                                            Utilities::pow(fe_degree + 1, dim));
+        }
+    }
   };
 
 
 
+  /**
+   * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use
+   * specific sum-factorization kernels and with normal/tangential shape_data
+   */
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-  inline void
-  FEEvaluationImplTransformToCollocation<
-    dim,
-    fe_degree,
-    n_q_points_1d,
-    Number>::evaluate(const unsigned int                     n_components,
-                      const EvaluationFlags::EvaluationFlags evaluation_flag,
-                      const Number                          *values_dofs,
-                      FEEvaluationData<dim, Number, false>  &fe_eval)
+  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+                          dim,
+                          fe_degree,
+                          n_q_points_1d,
+                          Number>
   {
-    const auto &shape_data = fe_eval.get_shape_info().data.front();
+    using Number2 =
+      typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
 
-    Assert(n_q_points_1d > fe_degree,
-           ExcMessage("You lose information when going to a collocation space "
-                      "of lower degree, so the evaluation results would be "
-                      "wrong. Thus, this class does not permit the desired "
-                      "operation."));
-    constexpr std::size_t n_dofs     = Utilities::pow(fe_degree + 1, dim);
-    constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+    template <bool integrate>
+    static void
+    evaluate_or_integrate(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number                                *values_dofs_actual,
+      FEEvaluationData<dim, Number, false>  &fe_eval,
+      const bool                             add_into_values_array = false);
 
-    for (unsigned int c = 0; c < n_components; ++c)
-      {
-        FEEvaluationImplBasisChange<
-          evaluate_evenodd,
-          EvaluatorQuantity::value,
-          dim,
-          (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
-          n_q_points_1d>::do_forward(1,
-                                     shape_data.shape_values_eo,
-                                     values_dofs + c * n_dofs,
-                                     fe_eval.begin_values() + c * n_q_points);
-
-        // apply derivatives in the collocation space
-        if (evaluation_flag &
-            (EvaluationFlags::gradients | EvaluationFlags::hessians))
-          FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
-            do_evaluate(shape_data,
-                        evaluation_flag & (EvaluationFlags::gradients |
-                                           EvaluationFlags::hessians),
-                        fe_eval.begin_values() + c * n_q_points,
-                        fe_eval.begin_gradients() + c * dim * n_q_points,
-                        fe_eval.begin_hessians() +
-                          c * dim * (dim + 1) / 2 * n_q_points);
-      }
-  }
+  private:
+    template <int direction, bool contract_over_rows>
+    static void
+    work_normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+                const Number                                            *in,
+                Number                                                  *out,
+                const bool add_into_result = false)
+    {
+      AssertIndexRange(direction, dim);
+      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 = data.shape_values_eo.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);
+
+      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<evaluate_evenodd,
+                                            EvaluatorQuantity::value,
+                                            n_rows,
+                                            n_columns,
+                                            n_blocks1,
+                                            n_blocks1,
+                                            contract_over_rows,
+                                            true>(shape_data, in, out);
+              else
+                apply_matrix_vector_product<evaluate_evenodd,
+                                            EvaluatorQuantity::value,
+                                            n_rows,
+                                            n_columns,
+                                            n_blocks1,
+                                            n_blocks1,
+                                            contract_over_rows,
+                                            false>(shape_data, in, out);
+
+              ++in;
+              ++out;
+            }
+          in += n_blocks1 * (mm - 1);
+          out += n_blocks1 * (nn - 1);
+        }
+    }
+
+    template <int direction, int normal_direction, bool contract_over_rows>
+    static void
+    work_tangential(
+      const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+      Number                                                  *ptr)
+    {
+      AssertIndexRange(direction, dim);
+      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 = data.shape_values_eo.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);
+
+      // Since we perform an in-place interpolation, we must run the step
+      // expanding the size of the basis backward ('contract_over_rows' aka
+      // 'evaluate' case).
+      if (contract_over_rows)
+        {
+          const Number *in =
+            ptr + (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
+          Number *out =
+            ptr + (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<evaluate_evenodd,
+                                              EvaluatorQuantity::value,
+                                              n_rows,
+                                              n_columns,
+                                              n_blocks1,
+                                              n_blocks1,
+                                              true,
+                                              false>(shape_data, in, out);
+
+                  --in;
+                  --out;
+                }
+              in -= n_blocks1 * (n_rows - 1);
+              out -= n_blocks1 * (n_columns - 1);
+            }
+        }
+      else
+        {
+          const Number *in  = ptr;
+          Number       *out = ptr;
+          for (int i2 = 0; i2 < n_blocks2; ++i2)
+            {
+              for (int i1 = 0; i1 < n_blocks1; ++i1)
+                {
+                  apply_matrix_vector_product<evaluate_evenodd,
+                                              EvaluatorQuantity::value,
+                                              n_rows,
+                                              n_columns,
+                                              n_blocks1,
+                                              n_blocks1,
+                                              false,
+                                              false>(shape_data, in, out);
+
+                  ++in;
+                  ++out;
+                }
+              in += n_blocks1 * (n_columns - 1);
+              out += n_blocks1 * (n_rows - 1);
+            }
+        }
+    }
+  };
 
 
 
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <bool integrate>
   inline void
-  FEEvaluationImplTransformToCollocation<
-    dim,
-    fe_degree,
-    n_q_points_1d,
-    Number>::integrate(const unsigned int                     n_components,
-                       const EvaluationFlags::EvaluationFlags integration_flag,
-                       Number                                *values_dofs,
-                       FEEvaluationData<dim, Number, false>  &fe_eval,
-                       const bool add_into_values_array)
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+                   dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   Number>::
+    evaluate_or_integrate(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number                                *values_dofs,
+      FEEvaluationData<dim, Number, false>  &fe_eval,
+      const bool                             add_into_values_array)
   {
-    const auto &shape_data = fe_eval.get_shape_info().data.front();
+    Assert(dim == 2 || dim == 3,
+           ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas "
+                      "evaluation/integration"));
 
-    Assert(n_q_points_1d > fe_degree,
-           ExcMessage("You lose information when going to a collocation space "
-                      "of lower degree, so the evaluation results would be "
-                      "wrong. Thus, this class does not permit the desired "
-                      "operation."));
-    constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+    if (evaluation_flag == EvaluationFlags::nothing)
+      return;
 
-    for (unsigned int c = 0; c < n_components; ++c)
+    AssertDimension(fe_eval.get_shape_info().data.size(), 2);
+    AssertDimension(n_q_points_1d,
+                    fe_eval.get_shape_info().data[0].n_q_points_1d);
+    AssertDimension(n_q_points_1d,
+                    fe_eval.get_shape_info().data[1].n_q_points_1d);
+    AssertDimension(fe_degree, fe_eval.get_shape_info().data[0].fe_degree);
+    AssertDimension(fe_degree, fe_eval.get_shape_info().data[1].fe_degree + 1);
+
+    const auto        &shape_data = fe_eval.get_shape_info().data;
+    const unsigned int dofs_per_component =
+      Utilities::pow(fe_degree, dim - 1) * (fe_degree + 1);
+    const unsigned int n_points  = Utilities::pow(n_q_points_1d, dim);
+    Number            *gradients = fe_eval.begin_gradients();
+    Number            *values    = fe_eval.begin_values();
+
+    if (integrate)
       {
-        // apply derivatives in collocation space
-        if (integration_flag &
-            (EvaluationFlags::gradients | EvaluationFlags::hessians))
-          FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
-            do_integrate(shape_data,
-                         integration_flag & (EvaluationFlags::gradients |
-                                             EvaluationFlags::hessians),
-                         fe_eval.begin_values() + c * n_q_points,
-                         fe_eval.begin_gradients() + c * dim * n_q_points,
-                         fe_eval.begin_hessians() +
-                           c * dim * (dim + 1) / 2 * n_q_points,
-                         /*add_into_values_array=*/
-                         integration_flag & EvaluationFlags::values);
-
-        // transform back to the original space
-        FEEvaluationImplBasisChange<
-          evaluate_evenodd,
-          EvaluatorQuantity::value,
-          dim,
-          (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
-          n_q_points_1d>::do_backward(1,
-                                      shape_data.shape_values_eo,
-                                      add_into_values_array,
-                                      fe_eval.begin_values() + c * n_q_points,
-                                      values_dofs +
-                                        c * Utilities::pow(fe_degree + 1, dim));
+        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],
+                                                              values,
+                                                              gradients,
+                                                              do_values);
+        if constexpr (dim > 2)
+          work_tangential<2, 0, false>(shape_data[1], values);
+        work_tangential<1, 0, false>(shape_data[1], values);
+        work_normal<0, false>(shape_data[0],
+                              values,
+                              values_dofs,
+                              add_into_values_array);
+
+        values += n_points;
+        gradients += n_points * dim;
+        values_dofs += dofs_per_component;
+
+        if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+          integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+                                                              values,
+                                                              gradients,
+                                                              do_values);
+        if constexpr (dim > 2)
+          work_tangential<2, 1, false>(shape_data[1], values);
+        work_tangential<0, 1, false>(shape_data[1], values);
+        work_normal<1, false>(shape_data[0],
+                              values,
+                              values_dofs,
+                              add_into_values_array);
+
+        if constexpr (dim > 2)
+          {
+            values += n_points;
+            gradients += n_points * dim;
+            values_dofs += dofs_per_component;
+
+            if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+              integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+                                                                  values,
+                                                                  gradients,
+                                                                  do_values);
+            work_tangential<1, 2, false>(shape_data[1], values);
+            work_tangential<0, 2, false>(shape_data[1], values);
+            work_normal<2, false>(shape_data[0],
+                                  values,
+                                  values_dofs,
+                                  add_into_values_array);
+          }
+      }
+    else
+      {
+        work_normal<0, true>(shape_data[0], values_dofs, values);
+        work_tangential<1, 0, true>(shape_data[1], values);
+        if constexpr (dim > 2)
+          work_tangential<2, 0, true>(shape_data[1], values);
+        if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+          evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+                                                             values,
+                                                             gradients);
+
+        values += n_points;
+        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);
+        if constexpr (dim > 2)
+          work_tangential<2, 1, true>(shape_data[1], values);
+        if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+          evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+                                                             values,
+                                                             gradients);
+
+        if constexpr (dim > 2)
+          {
+            values += n_points;
+            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);
+            work_tangential<1, 2, true>(shape_data[1], values);
+            if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+              evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+                                                                 values,
+                                                                 gradients);
+          }
       }
   }
 
@@ -2275,7 +2137,7 @@ namespace internal
         const EvaluationFlags::EvaluationFlags evaluation_flag,
         OtherNumber                           *values_dofs,
         FEEvaluationData<dim, Number, false>  &fe_eval,
-        const bool                             sum_into_values_array = false)
+        const bool                             sum_into_values_array_in = false)
     {
       // `OtherNumber` is either `const Number` (evaluate()) or `Number`
       // (integrate())
@@ -2291,13 +2153,39 @@ namespace internal
                element_type == ElementType::tensor_raviart_thomas,
              ExcNotImplemented());
 
+      EvaluationFlags::EvaluationFlags actual_flag = evaluation_flag;
+      bool sum_into_values_array                   = sum_into_values_array_in;
+      if (evaluation_flag & EvaluationFlags::hessians)
+        {
+          actual_flag |= EvaluationFlags::values;
+          Assert(element_type != MatrixFreeFunctions::tensor_none,
+                 ExcNotImplemented());
+          if constexpr (do_integrate)
+            {
+              if (fe_eval.get_shape_info().data[0].fe_degree <
+                  fe_eval.get_shape_info().data[0].n_q_points_1d)
+                integrate_hessians_collocation<n_q_points_1d>(
+                  n_components,
+                  fe_eval,
+                  evaluation_flag & EvaluationFlags::values);
+              else
+                {
+                  integrate_hessians_slow(n_components,
+                                          fe_eval,
+                                          values_dofs,
+                                          sum_into_values_array);
+                  sum_into_values_array = true;
+                }
+            }
+        }
+
       if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
           element_type == ElementType::tensor_symmetric_collocation)
         {
           evaluate_or_integrate<
             FEEvaluationImplCollocation<dim, fe_degree, Number>>(
             n_components,
-            evaluation_flag,
+            actual_flag,
             values_dofs,
             fe_eval,
             sum_into_values_array);
@@ -2314,7 +2202,7 @@ namespace internal
                                                    n_q_points_1d,
                                                    Number>>(
             n_components,
-            evaluation_flag,
+            actual_flag,
             values_dofs,
             fe_eval,
             sum_into_values_array);
@@ -2328,7 +2216,7 @@ namespace internal
                                                  n_q_points_1d,
                                                  Number>>(
             n_components,
-            evaluation_flag,
+            actual_flag,
             values_dofs,
             fe_eval,
             sum_into_values_array);
@@ -2341,7 +2229,7 @@ namespace internal
                              fe_degree,
                              n_q_points_1d,
                              Number>>(n_components,
-                                      evaluation_flag,
+                                      actual_flag,
                                       values_dofs,
                                       fe_eval,
                                       sum_into_values_array);
@@ -2354,7 +2242,7 @@ namespace internal
                                                  n_q_points_1d,
                                                  Number>>(
             n_components,
-            evaluation_flag,
+            actual_flag,
             values_dofs,
             fe_eval,
             sum_into_values_array);
@@ -2367,23 +2255,32 @@ namespace internal
                                                  n_q_points_1d,
                                                  Number>>(
             n_components,
-            evaluation_flag,
+            actual_flag,
             values_dofs,
             fe_eval,
             sum_into_values_array);
         }
       else if (element_type == ElementType::tensor_raviart_thomas)
         {
-          FEEvaluationImpl<ElementType::tensor_raviart_thomas,
-                           dim,
-                           (fe_degree == -1) ? 1 : fe_degree,
-                           (n_q_points_1d < 1) ? 1 : n_q_points_1d,
-                           Number>::
-            template evaluate_or_integrate<do_integrate>(evaluation_flag,
-                                                         const_cast<Number *>(
-                                                           values_dofs),
-                                                         fe_eval,
-                                                         sum_into_values_array);
+          if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1)
+            {
+              FEEvaluationImpl<ElementType::tensor_raviart_thomas,
+                               dim,
+                               fe_degree,
+                               n_q_points_1d,
+                               Number>::
+                template evaluate_or_integrate<do_integrate>(
+                  actual_flag,
+                  const_cast<Number *>(values_dofs),
+                  fe_eval,
+                  sum_into_values_array);
+            }
+          else
+            {
+              Assert(false,
+                     ExcNotImplemented("Raviart-Thomas currently only possible "
+                                       "in 2d/3d and with templated degree"));
+            }
         }
       else
         {
@@ -2393,12 +2290,23 @@ namespace internal
                                                  n_q_points_1d,
                                                  Number>>(
             n_components,
-            evaluation_flag,
+            actual_flag,
             values_dofs,
             fe_eval,
             sum_into_values_array);
         }
 
+      if ((evaluation_flag & EvaluationFlags::hessians) && !do_integrate)
+        {
+          Assert(element_type != MatrixFreeFunctions::tensor_none,
+                 ExcNotImplemented());
+          if (fe_eval.get_shape_info().data[0].fe_degree <
+              fe_eval.get_shape_info().data[0].n_q_points_1d)
+            evaluate_hessians_collocation<n_q_points_1d>(n_components, fe_eval);
+          else
+            evaluate_hessians_slow(n_components, values_dofs, fe_eval);
+        }
+
       return false;
     }
 
@@ -2580,25 +2488,24 @@ namespace internal
                                              Number,
                                              Number2>
                         eval_grad({}, data.shape_gradients_collocation_eo, {});
-                      eval_grad.template gradients<0, true, false>(
+                      eval_grad.template gradients<0, true, false, 3>(
                         values_quad, gradients_quad);
-                      eval_grad.template gradients<1, true, false>(
-                        values_quad, gradients_quad + n_q_points);
+                      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>(scratch_data,
-                                                            gradients_quad);
+                      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>(scratch_data,
-                                                               gradients_quad +
-                                                                 n_q_points);
+                      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,
@@ -2607,16 +2514,15 @@ namespace internal
                   // grad z
                   eval0.template values<0, true, false>(values_dofs + n_dofs,
                                                         scratch_data);
-                  eval1.template values<1, true, false>(
-                    scratch_data, gradients_quad + (dim - 1) * n_q_points);
+                  eval1.template values<1, true, false, 3>(scratch_data,
+                                                           gradients_quad + 2);
 
                   break;
                 case 2:
-                  eval0.template values<0, true, false>(values_dofs + n_dofs,
-                                                        gradients_quad +
-                                                          n_q_points);
-                  eval0.template gradients<0, true, false>(values_dofs,
-                                                           gradients_quad);
+                  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);
@@ -2766,12 +2672,9 @@ namespace internal
               {
                 case 3:
                   // grad z
-                  eval1.template values<1, false, false>(gradients_quad +
-                                                           2 * n_q_points,
-                                                         gradients_quad +
-                                                           2 * n_q_points);
-                  eval0.template values<0, false, false>(gradients_quad +
-                                                           2 * n_q_points,
+                  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))
@@ -2784,12 +2687,12 @@ namespace internal
                                              Number2>
                         eval_grad({}, data.shape_gradients_collocation_eo, {});
                       if ((integration_flag & EvaluationFlags::values) != 0u)
-                        eval_grad.template gradients<1, false, true>(
-                          gradients_quad + n_q_points, values_quad);
+                        eval_grad.template gradients<1, false, true, 3>(
+                          gradients_quad + 1, values_quad);
                       else
-                        eval_grad.template gradients<1, false, false>(
-                          gradients_quad + n_q_points, values_quad);
-                      eval_grad.template gradients<0, false, true>(
+                        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);
@@ -2802,30 +2705,30 @@ namespace internal
                         {
                           eval1.template values<1, false, false>(values_quad,
                                                                  scratch_data);
-                          eval1.template gradients<1, false, true>(
-                            gradients_quad + n_q_points, scratch_data);
+                          eval1.template gradients<1, false, true, 3>(
+                            gradients_quad + 1, scratch_data);
                         }
                       else
-                        eval1.template gradients<1, false, false>(
-                          gradients_quad + n_q_points, scratch_data);
+                        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>(gradients_quad,
-                                                             scratch_data);
+                      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>(gradients_quad +
-                                                           n_q_points,
-                                                         values_dofs + n_dofs);
-                  eval0.template gradients<0, false, false>(gradients_quad,
-                                                            values_dofs);
+                  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);
@@ -2942,6 +2845,8 @@ namespace internal
     }
   };
 
+
+
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
   struct FEFaceEvaluationImplRaviartThomas
   {
@@ -3559,44 +3464,71 @@ namespace internal
     {
       if (face_direction == face_no / 2)
         {
-          EvaluatorTensorProduct<evaluate_general,
-                                 dim,
-                                 fe_degree + 1,
-                                 0,
-                                 Number,
-                                 Number2>
-            evalf(shape_data[face_no % 2].begin(),
-                  nullptr,
-                  nullptr,
-                  n_points_1d,
-                  0);
-
-          const unsigned int in_stride  = do_evaluate ?
-                                            dofs_per_component_on_cell :
-                                            dofs_per_component_on_face;
-          const unsigned int out_stride = do_evaluate ?
-                                            dofs_per_component_on_face :
-                                            dofs_per_component_on_cell;
+          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)
-                evalf.template apply_face<face_direction,
-                                          do_evaluate,
-                                          add_into_output,
-                                          2>(input, output);
+                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)
-                evalf.template apply_face<face_direction,
-                                          do_evaluate,
-                                          add_into_output,
-                                          1>(input, output);
+                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
-                evalf.template apply_face<face_direction,
-                                          do_evaluate,
-                                          add_into_output,
-                                          0>(input, output);
-              input += in_stride;
-              output += out_stride;
+                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)
@@ -3961,13 +3893,13 @@ namespace internal
               if (integrate)
                 for (unsigned int q = 0; q < n_q_points; ++q)
                   tmp_values[q] =
-                    gradients_quad[(c * dim + d) * n_q_points + orientation[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 * dim + d) * n_q_points + q];
+                    gradients_quad[(c * n_q_points + q) * dim + d];
               for (unsigned int q = 0; q < n_q_points; ++q)
-                gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
+                gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
             }
         if (flag & EvaluationFlags::hessians)
           {
@@ -4026,14 +3958,15 @@ namespace internal
               Assert(gradients_quad != nullptr, ExcInternalError());
               if (integrate)
                 for (unsigned int q = 0; q < n_q_points; ++q)
-                  tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
-                                                 orientation[q]][v];
+                  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 * dim + d) * n_q_points + q][v];
+                    gradients_quad[(c * n_q_points + q) * dim + d][v];
               for (unsigned int q = 0; q < n_q_points; ++q)
-                gradients_quad[(c * dim + d) * n_q_points + q][v] =
+                gradients_quad[(c * n_q_points + q) * dim + d][v] =
                   tmp_values[q];
             }
         if (flag & EvaluationFlags::hessians)
@@ -4129,11 +4062,10 @@ namespace internal
                                 n_dofs,
                                 n_q_points);
 
-                      eval.template gradients<0, true, false>(
-                        values_dofs_actual_ptr, gradients_quad_ptr);
-
-                      gradients_quad_ptr += 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;
                 }
             }
@@ -4370,14 +4302,13 @@ namespace internal
 
                       if (!(integration_flag & EvaluationFlags::values) &&
                           d == 0)
-                        eval.template gradients<0, false, false>(
-                          gradients_quad_ptr, values_dofs_actual_ptr);
+                        eval.template gradients<0, false, false, dim>(
+                          gradients_quad_ptr + d, values_dofs_actual_ptr);
                       else
-                        eval.template gradients<0, false, true>(
-                          gradients_quad_ptr, values_dofs_actual_ptr);
-
-                      gradients_quad_ptr += n_q_points;
+                        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;
                 }
             }
index d1acf28d8918804cb1406fd8f75b9ef3822c58c9..f5288fc65d49f9119414eb8631c198334be0e1de 100644 (file)
@@ -4526,10 +4526,10 @@ inline DEAL_II_ALWAYS_INLINE
   // Cartesian cell
   if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
-      for (unsigned int d = 0; d < dim; ++d)
-        for (unsigned int comp = 0; comp < n_components; ++comp)
+      for (unsigned int comp = 0; comp < n_components; ++comp)
+        for (unsigned int d = 0; d < dim; ++d)
           grad_out[comp][d] =
-            this->gradients_quad[(comp * dim + d) * nqp + q_point] *
+            this->gradients_quad[(comp * nqp + q_point) * dim + d] *
             this->jacobian[0][d][d];
     }
   // cell with general/affine Jacobian
@@ -4543,11 +4543,11 @@ inline DEAL_II_ALWAYS_INLINE
         for (unsigned int d = 0; d < dim; ++d)
           {
             grad_out[comp][d] =
-              jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
+              jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
             for (unsigned int e = 1; e < dim; ++e)
               grad_out[comp][d] +=
                 jac[d][e] *
-                this->gradients_quad[(comp * dim + e) * nqp + q_point];
+                this->gradients_quad[(comp * nqp + q_point) * dim + e];
           }
     }
   return grad_out;
@@ -4580,7 +4580,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
     for (unsigned int comp = 0; comp < n_components; ++comp)
       grad_out[comp] =
-        this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
+        this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
         (this->normal_x_jacobian[0][dim - 1]);
   else
     {
@@ -4588,11 +4588,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
       for (unsigned int comp = 0; comp < n_components; ++comp)
         {
-          grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
+          grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
                            this->normal_x_jacobian[index][0];
           for (unsigned int d = 1; d < dim; ++d)
             grad_out[comp] +=
-              this->gradients_quad[(comp * dim + d) * nqp + q_point] *
+              this->gradients_quad[(comp * nqp + q_point) * dim + d] *
               this->normal_x_jacobian[index][d];
         }
     }
@@ -4770,7 +4770,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int e = 0; e < dim; ++e)
               hessian_out[comp][d][d] +=
                 jac_grad[d][e] *
-                this->gradients_quad[(comp * dim + e) * nqp + q_point];
+                this->gradients_quad[(comp * nqp + q_point) * dim + e];
 
           // add off-diagonal part of J' * grad(u)
           for (unsigned int d = 0, count = dim; d < dim; ++d)
@@ -4778,7 +4778,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
               for (unsigned int f = 0; f < dim; ++f)
                 hessian_out[comp][d][e] +=
                   jac_grad[count][f] *
-                  this->gradients_quad[(comp * dim + f) * nqp + q_point];
+                  this->gradients_quad[(comp * nqp + q_point) * dim + f];
 
           // take symmetric part
           for (unsigned int d = 0; d < dim; ++d)
@@ -4872,7 +4872,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int e = 0; e < dim; ++e)
               hessian_out[comp][d] +=
                 jac_grad[d][e] *
-                this->gradients_quad[(comp * dim + e) * nqp + q_point];
+                this->gradients_quad[(comp * nqp + q_point) * dim + e];
         }
     }
   return hessian_out;
@@ -5006,7 +5006,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         {
           const VectorizedArrayType factor = jac[d] * JxW;
           for (unsigned int comp = 0; comp < n_components; ++comp)
-            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+            this->gradients_quad[(comp * nqp + q_point) * dim + d] =
               grad_in[comp][d] * factor;
         }
     }
@@ -5026,7 +5026,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
             for (unsigned int e = 1; e < dim; ++e)
               new_val += (jac[e][d] * grad_in[comp][e]);
-            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+            this->gradients_quad[(comp * nqp + q_point) * dim + d] =
               new_val * JxW;
           }
     }
@@ -5062,9 +5062,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       for (unsigned int comp = 0; comp < n_components; ++comp)
         {
           for (unsigned int d = 0; d < dim - 1; ++d)
-            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+            this->gradients_quad[(comp * nqp + q_point) * dim + d] =
               VectorizedArrayType();
-          this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
+          this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] =
             grad_in[comp] * JxW_jac;
         }
     }
@@ -5081,7 +5081,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       for (unsigned int comp = 0; comp < n_components; ++comp)
         {
           for (unsigned int d = 0; d < dim; ++d)
-            this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+            this->gradients_quad[(comp * nqp + q_point) * dim + d] =
               (grad_in[comp] * JxW) * jac[d];
         }
     }
@@ -5233,8 +5233,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                 for (unsigned int f = e + 1; f < dim; ++f, ++count)
                   sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
                          jac_grad[count][d];
-              this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
-                                                 q_point] = sum * JxW;
+              this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
+                                                 d] = sum * JxW;
             }
         }
     }
@@ -5495,7 +5495,7 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
 FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_gradient(
-  const unsigned int q_point) const
+  const unsigned int q_point_in) const
 {
   // could use the base class gradient, but that involves too many expensive
   // initialization operations on tensors
@@ -5504,7 +5504,7 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_gradient(
   Assert(this->gradients_quad_initialized == true,
          internal::ExcAccessToUninitializedField());
 #  endif
-  AssertIndexRange(q_point, this->n_quadrature_points);
+  AssertIndexRange(q_point_in, this->n_quadrature_points);
 
   Assert(this->jacobian != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
@@ -5512,12 +5512,12 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_gradient(
 
   Tensor<1, dim, VectorizedArrayType> grad_out;
 
-  const std::size_t nqp = this->n_quadrature_points;
+  const std::size_t q_point = q_point_in;
   if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       for (unsigned int d = 0; d < dim; ++d)
         grad_out[d] =
-          this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
+          this->gradients_quad[dim * q_point + d] * this->jacobian[0][d][d];
     }
   // cell with general/affine Jacobian
   else
@@ -5528,9 +5528,9 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::get_gradient(
                          0];
       for (unsigned int d = 0; d < dim; ++d)
         {
-          grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
+          grad_out[d] = jac[d][0] * this->gradients_quad[dim * q_point];
           for (unsigned int e = 1; e < dim; ++e)
-            grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
+            grad_out[d] += jac[d][e] * this->gradients_quad[dim * q_point + e];
         }
     }
   return grad_out;
@@ -5641,12 +5641,12 @@ template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline DEAL_II_ALWAYS_INLINE void
 FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
   submit_gradient(const Tensor<1, dim, VectorizedArrayType> grad_in,
-                  const unsigned int                        q_point)
+                  const unsigned int                        q_point_in)
 {
 #  ifdef DEBUG
   Assert(this->is_reinitialized, ExcNotInitialized());
 #  endif
-  AssertIndexRange(q_point, this->n_quadrature_points);
+  AssertIndexRange(q_point_in, this->n_quadrature_points);
   Assert(this->J_value != nullptr,
          internal::ExcMatrixFreeAccessToUninitializedMappingField(
            "update_gradients"));
@@ -5657,7 +5657,8 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
   this->gradients_quad_submitted = true;
 #  endif
 
-  const std::size_t nqp = this->n_quadrature_points;
+  const std::size_t    q_point  = q_point_in;
+  VectorizedArrayType *grad_ptr = this->gradients_quad + dim * q_point;
   if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
       const VectorizedArrayType JxW =
@@ -5670,7 +5671,7 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
         jac[d] = this->jacobian[0][d][d];
 
       for (unsigned int d = 0; d < dim; ++d)
-        this->gradients_quad[d * nqp + q_point] = grad_in[d] * jac[d] * JxW;
+        grad_ptr[d] = grad_in[d] * jac[d] * JxW;
     }
   // general/affine cell type
   else
@@ -5688,7 +5689,7 @@ FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
           VectorizedArrayType new_val = jac[0][d] * grad_in[0];
           for (unsigned int e = 1; e < dim; ++e)
             new_val += jac[e][d] * grad_in[e];
-          this->gradients_quad[d * nqp + q_point] = new_val * JxW;
+          grad_ptr[d] = new_val * JxW;
         }
     }
 }
@@ -5902,7 +5903,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int comp = 0; comp < n_components; ++comp)
               grad_out[comp][d] =
-                this->gradients_quad[(comp * dim + d) * nqp + q_point] *
+                this->gradients_quad[(comp * nqp + q_point) * dim + d] *
                 inv_t_jac[d][d] * (jac[comp][comp] * inv_det);
         }
       else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
@@ -5927,7 +5928,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                 for (unsigned int f = 0; f < dim; ++f)
                   for (unsigned int e = 0; e < dim; ++e)
                     tmp += jac[comp][f] * inv_t_jac[d][e] *
-                           this->gradients_quad[(f * dim + e) * nqp + q_point];
+                           this->gradients_quad[(f * nqp + q_point) * dim + e];
 
                 grad_out[comp][d] = tmp * inv_det;
               }
@@ -5964,7 +5965,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                 for (unsigned int f = 0; f < dim; ++f)
                   for (unsigned int e = 0; e < dim; ++e)
                     tmp += t_jac[f][comp] * inv_t_jac[d][e] *
-                           this->gradients_quad[(f * dim + e) * nqp + q_point];
+                           this->gradients_quad[(f * nqp + q_point) * dim + e];
 
                 grad_out[comp][d] = tmp * inv_det;
               }
@@ -6069,9 +6070,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                            this->jacobian[0][2][2];
 
           // div * det(J^-1)
-          divergence = this->gradients_quad[q_point];
+          divergence = this->gradients_quad[q_point * dim];
           for (unsigned int d = 1; d < dim; ++d)
-            divergence += this->gradients_quad[(dim * d + d) * nqp + q_point];
+            divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
           divergence *= inv_det;
         }
       else
@@ -6087,9 +6088,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
             Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
 
           // div * det(J^-1)
-          divergence = this->gradients_quad[q_point];
+          divergence = this->gradients_quad[q_point * dim];
           for (unsigned int d = 1; d < dim; ++d)
-            divergence += this->gradients_quad[(dim * d + d) * nqp + q_point];
+            divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
           divergence *= inv_det;
         }
     }
@@ -6099,9 +6100,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
           this->cell_type == internal::MatrixFreeFunctions::cartesian)
         {
           // Cartesian cell
-          divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
+          divergence =
+            this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
           for (unsigned int d = 1; d < dim; ++d)
-            divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
+            divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
                           this->jacobian[0][d][d];
         }
       else
@@ -6111,13 +6113,13 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
             this->cell_type == internal::MatrixFreeFunctions::general ?
               this->jacobian[q_point] :
               this->jacobian[0];
-          divergence = jac[0][0] * this->gradients_quad[q_point];
+          divergence = jac[0][0] * this->gradients_quad[q_point * dim];
           for (unsigned int e = 1; e < dim; ++e)
-            divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
+            divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
           for (unsigned int d = 1; d < dim; ++d)
             for (unsigned int e = 0; e < dim; ++e)
               divergence +=
-                jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
+                jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
         }
     }
   return divergence;
@@ -6325,7 +6327,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
           const VectorizedArrayType weight = this->quadrature_weights[q_point];
           for (unsigned int d = 0; d < dim; ++d)
             for (unsigned int comp = 0; comp < n_components; ++comp)
-              this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+              this->gradients_quad[(comp * nqp + q_point) * dim + d] =
                 grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
         }
       else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
@@ -6353,7 +6355,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                   for (unsigned int e = 0; e < dim; ++e)
                     tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e];
 
-                this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+                this->gradients_quad[(comp * nqp + q_point) * dim + d] =
                   tmp * fac;
               }
         }
@@ -6385,7 +6387,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                   for (unsigned int e = 0; e < dim; ++e)
                     tmp += t_jac[comp][f] * inv_t_jac[e][d] * grad_in[f][e];
 
-                this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+                this->gradients_quad[(comp * nqp + q_point) * dim + d] =
                   tmp * fac;
               }
 
@@ -6531,12 +6533,12 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
 
       for (unsigned int d = 0; d < dim; ++d)
         {
-          this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
+          this->gradients_quad[(d * nqp + q_point) * dim + d] = fac;
           for (unsigned int e = d + 1; e < dim; ++e)
             {
-              this->gradients_quad[(dim * d + e) * nqp + q_point] =
+              this->gradients_quad[(d * nqp + q_point) * dim + e] =
                 VectorizedArrayType();
-              this->gradients_quad[(dim * e + d) * nqp + q_point] =
+              this->gradients_quad[(e * nqp + q_point) * dim + d] =
                 VectorizedArrayType();
             }
         }
@@ -6551,13 +6553,13 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
             this->J_value[0] * this->quadrature_weights[q_point] * div_in;
           for (unsigned int d = 0; d < dim; ++d)
             {
-              this->gradients_quad[(d * dim + d) * nqp + q_point] =
+              this->gradients_quad[(d * nqp + q_point) * dim + d] =
                 (fac * this->jacobian[0][d][d]);
               for (unsigned int e = d + 1; e < dim; ++e)
                 {
-                  this->gradients_quad[(d * dim + e) * nqp + q_point] =
+                  this->gradients_quad[(d * nqp + q_point) * dim + e] =
                     VectorizedArrayType();
-                  this->gradients_quad[(e * dim + d) * nqp + q_point] =
+                  this->gradients_quad[(e * nqp + q_point) * dim + d] =
                     VectorizedArrayType();
                 }
             }
@@ -6576,7 +6578,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
           for (unsigned int d = 0; d < dim; ++d)
             {
               for (unsigned int e = 0; e < dim; ++e)
-                this->gradients_quad[(d * dim + e) * nqp + q_point] =
+                this->gradients_quad[(d * nqp + q_point) * dim + e] =
                   jac[d][e] * fac;
             }
         }
@@ -6620,16 +6622,16 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
       const VectorizedArrayType JxW =
         this->J_value[0] * this->quadrature_weights[q_point];
       for (unsigned int d = 0; d < dim; ++d)
-        this->gradients_quad[(d * dim + d) * nqp + q_point] =
+        this->gradients_quad[(d * nqp + q_point) * dim + d] =
           (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
       for (unsigned int e = 0, counter = dim; e < dim; ++e)
         for (unsigned int d = e + 1; d < dim; ++d, ++counter)
           {
             const VectorizedArrayType value =
               sym_grad.access_raw_entry(counter) * JxW;
-            this->gradients_quad[(e * dim + d) * nqp + q_point] =
+            this->gradients_quad[(e * nqp + q_point) * dim + d] =
               value * this->jacobian[0][d][d];
-            this->gradients_quad[(d * dim + e) * nqp + q_point] =
+            this->gradients_quad[(d * nqp + q_point) * dim + e] =
               value * this->jacobian[0][e][e];
           }
     }
@@ -6661,7 +6663,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
             VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
             for (unsigned int e = 1; e < dim; ++e)
               new_val += jac[e][d] * weighted[comp][e];
-            this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
+            this->gradients_quad[(comp * nqp + q_point) * dim + d] = new_val;
           }
     }
 }
index dfbd342a040f52edc5856146a6ea38d6932d6f09..5834f41fb054095d5d0efbb7d228885f2bd6d52b 100644 (file)
@@ -1230,7 +1230,7 @@ namespace internal
                       for (unsigned int e = 0; e < dim; ++e)
                         jac[d][e] =
                           eval
-                            .begin_gradients()[q + (d * dim + e) * n_q_points];
+                            .begin_gradients()[e + (d * n_q_points + q) * dim];
 
                     // eliminate roundoff errors
                     if (cell_type[cell] == cartesian)
@@ -2180,7 +2180,7 @@ namespace internal
                   for (unsigned int e = 0; e < dim; ++e)
                     for (unsigned int d = 0; d < dim; ++d)
                       jacobi[d][e] =
-                        eval.begin_gradients()[(d * dim + e) * n_q_points + q];
+                        eval.begin_gradients()[(d * n_q_points + q) * dim + e];
                   Tensor<2, dim, VectorizedDouble> inv_transp_jac =
                     transpose(invert(jacobi));
                   Tensor<3, dim, VectorizedDouble> jac_grad;
@@ -2250,7 +2250,7 @@ namespace internal
                       for (unsigned int d = 0; d < dim; ++d)
                         jac[d][ee] =
                           eval_int
-                            .begin_gradients()[(d * dim + e) * n_q_points + q];
+                            .begin_gradients()[(d * n_q_points + q) * dim + e];
                     }
                   Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac);
                   for (unsigned int e = 0; e < dim; ++e)
@@ -2365,8 +2365,8 @@ namespace internal
                           for (unsigned int d = 0; d < dim; ++d)
                             jac[d][ee] =
                               eval_ext
-                                .begin_gradients()[(d * dim + e) * n_q_points +
-                                                   q];
+                                .begin_gradients()[(d * n_q_points + q) * dim +
+                                                   e];
                         }
                       Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac);
                       for (unsigned int e = 0; e < dim; ++e)
index 91a42496c4277c1a670fc05d223b8b0b13e30a4b..11aa7fefc738d5eeea65d63a7ae93e229004909c 100644 (file)
@@ -682,13 +682,13 @@ namespace internal
   inline DEAL_II_ALWAYS_INLINE
 #endif
     std::enable_if_t<(variant == evaluate_evenodd), void>
-    apply_matrix_vector_product(const Number2 *matrix,
-                                const Number  *in,
-                                Number        *out,
-                                int            n_rows_runtime     = 0,
-                                int            n_columns_runtime  = 0,
-                                int            stride_in_runtime  = 0,
-                                int            stride_out_runtime = 0)
+    apply_matrix_vector_product(const Number2 *DEAL_II_RESTRICT matrix,
+                                const Number                   *in,
+                                Number                         *out,
+                                int n_rows_runtime     = 0,
+                                int n_columns_runtime  = 0,
+                                int stride_in_runtime  = 0,
+                                int stride_out_runtime = 0)
   {
     const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static;
     const int n_columns =
@@ -1165,25 +1165,23 @@ namespace internal
       (void)dummy2;
     }
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int direction, bool contract_over_rows, bool add, int stride = 1>
     void
     values(const Number in[], Number out[]) const
     {
-      apply<direction,
-            contract_over_rows,
-            add,
-            false,
-            EvaluatorQuantity::value>(shape_values, in, out);
+      constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
+      apply<direction, contract_over_rows, add, false, value_type, stride>(
+        shape_values, in, out);
     }
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int direction, bool contract_over_rows, bool add, int stride = 1>
     void
     gradients(const Number in[], Number out[]) const
     {
       constexpr EvaluatorQuantity gradient_type =
         (variant == evaluate_general ? EvaluatorQuantity::value :
                                        EvaluatorQuantity::gradient);
-      apply<direction, contract_over_rows, add, false, gradient_type>(
+      apply<direction, contract_over_rows, add, false, gradient_type, stride>(
         shape_gradients, in, out);
     }
 
@@ -1263,7 +1261,8 @@ namespace internal
               bool contract_over_rows,
               bool add,
               bool one_line     = false,
-              EvaluatorQuantity = EvaluatorQuantity::value>
+              EvaluatorQuantity = EvaluatorQuantity::value,
+              int extra_stride  = 1>
     static void
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
           const Number                   *in,
@@ -1324,7 +1323,8 @@ namespace internal
             bool              contract_over_rows,
             bool              add,
             bool              one_line,
-            EvaluatorQuantity quantity>
+            EvaluatorQuantity quantity,
+            int               extra_stride>
   inline void
   EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
@@ -1349,6 +1349,8 @@ namespace internal
     constexpr int n_blocks2 =
       Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
 
+    constexpr int stride_in  = !contract_over_rows ? extra_stride : 1;
+    constexpr int stride_out = contract_over_rows ? extra_stride : 1;
     for (int i2 = 0; i2 < n_blocks2; ++i2)
       {
         for (int i1 = 0; i1 < n_blocks1; ++i1)
@@ -1357,22 +1359,146 @@ namespace internal
                                         quantity,
                                         n_rows,
                                         n_columns,
-                                        stride,
-                                        stride,
+                                        stride * stride_in,
+                                        stride * stride_out,
                                         contract_over_rows,
                                         add>(shape_data, in, out);
 
             if (one_line == false)
               {
-                ++in;
-                ++out;
+                in += stride_in;
+                out += stride_out;
               }
           }
         if (one_line == false)
           {
-            in += stride * (mm - 1);
-            out += stride * (nn - 1);
+            in += stride * (mm - 1) * stride_in;
+            out += stride * (nn - 1) * stride_out;
+          }
+      }
+  }
+
+
+
+  template <int  n_rows_template,
+            int  stride_template,
+            bool contract_onto_face,
+            bool add,
+            int  max_derivative,
+            typename Number,
+            typename Number2>
+  inline std::enable_if_t<contract_onto_face, void>
+  interpolate_to_face(const Number2            *shape_values,
+                      const std::array<int, 2> &n_blocks,
+                      const std::array<int, 2> &steps,
+                      const Number             *input,
+                      Number *DEAL_II_RESTRICT  output,
+                      const int                 n_rows_runtime = 0,
+                      const int                 stride_runtime = 1)
+  {
+    const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
+    const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
+
+    Number *output1 = output + n_blocks[0] * n_blocks[1];
+    Number *output2 = output1 + n_blocks[0] * n_blocks[1];
+    for (int i2 = 0; i2 < n_blocks[1]; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks[0]; ++i1)
+          {
+            Number res0 = shape_values[0] * input[0];
+            Number res1, res2;
+            if (max_derivative > 0)
+              res1 = shape_values[n_rows] * input[0];
+            if (max_derivative > 1)
+              res2 = shape_values[2 * n_rows] * input[0];
+            for (int ind = 1; ind < n_rows; ++ind)
+              {
+                res0 += shape_values[ind] * input[stride * ind];
+                if (max_derivative > 0)
+                  res1 += shape_values[ind + n_rows] * input[stride * ind];
+                if (max_derivative > 1)
+                  res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
+              }
+            if (add)
+              {
+                output[i1] += res0;
+                if (max_derivative > 0)
+                  output1[i1] += res1;
+                if (max_derivative > 1)
+                  output2[i2] += res2;
+              }
+            else
+              {
+                output[i1] = res0;
+                if (max_derivative > 0)
+                  output1[i1] = res1;
+                if (max_derivative > 1)
+                  output2[i1] = res2;
+              }
+            input += steps[0];
           }
+        output += n_blocks[0];
+        if (max_derivative > 0)
+          output1 += n_blocks[0];
+        if (max_derivative > 1)
+          output2 += n_blocks[0];
+        input += steps[1];
+      }
+  }
+
+
+
+  template <int  n_rows_template,
+            int  stride_template,
+            bool contract_onto_face,
+            bool add,
+            int  max_derivative,
+            typename Number,
+            typename Number2>
+  inline std::enable_if_t<!contract_onto_face, void>
+  interpolate_to_face(const Number2            *shape_values,
+                      const std::array<int, 2> &n_blocks,
+                      const std::array<int, 2> &steps,
+                      const Number             *input,
+                      Number *DEAL_II_RESTRICT  output,
+                      const int                 n_rows_runtime = 0,
+                      const int                 stride_runtime = 1)
+  {
+    const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
+    const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
+
+    const Number *input1 = input + n_blocks[0] * n_blocks[1];
+    const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
+    for (int i2 = 0; i2 < n_blocks[1]; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks[0]; ++i1)
+          {
+            const Number in = input[i1];
+            Number       in1, in2;
+            if (max_derivative > 0)
+              in1 = input1[i1];
+            if (max_derivative > 1)
+              in2 = input2[i1];
+            for (int col = 0; col < n_rows; ++col)
+              {
+                Number result =
+                  add ? (output[col * stride] + shape_values[col] * in) :
+                        (shape_values[col] * in);
+                if (max_derivative > 0)
+                  result += shape_values[col + n_rows] * in1;
+                if (max_derivative > 1)
+                  result += shape_values[col + 2 * n_rows] * in2;
+
+                output[col * stride] = result;
+              }
+            output += steps[0];
+          }
+        input += n_blocks[0];
+        if (max_derivative > 0)
+          input1 += n_blocks[0];
+        if (max_derivative > 1)
+          input2 += n_blocks[0];
+        output += steps[1];
       }
   }
 
@@ -1385,7 +1511,7 @@ namespace internal
             typename Number,
             typename Number2>
   template <int  face_direction,
-            bool contract_onto_face,
+            bool contract_to_face,
             bool add,
             int  max_derivative>
   inline void
@@ -1400,116 +1526,24 @@ namespace internal
            ExcMessage(
              "The given array shape_values must not be the null pointer."));
 
-    constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
-    constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
-
-    AssertIndexRange(face_direction, dim);
-    constexpr int in_stride  = Utilities::pow(n_rows, face_direction);
-    constexpr int out_stride = Utilities::pow(n_rows, dim - 1);
-    const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values;
-
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
-          {
-            if (contract_onto_face == true)
-              {
-                Number res0 = shape_values[0] * in[0];
-                Number res1, res2;
-                if (max_derivative > 0)
-                  res1 = shape_values[n_rows] * in[0];
-                if (max_derivative > 1)
-                  res2 = shape_values[2 * n_rows] * in[0];
-                for (int ind = 1; ind < n_rows; ++ind)
-                  {
-                    res0 += shape_values[ind] * in[in_stride * ind];
-                    if (max_derivative > 0)
-                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
-                    if (max_derivative > 1)
-                      res2 +=
-                        shape_values[ind + 2 * n_rows] * in[in_stride * ind];
-                  }
-                if (add)
-                  {
-                    out[0] += res0;
-                    if (max_derivative > 0)
-                      out[out_stride] += res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] += res2;
-                  }
-                else
-                  {
-                    out[0] = res0;
-                    if (max_derivative > 0)
-                      out[out_stride] = res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] = res2;
-                  }
-              }
-            else
-              {
-                for (int col = 0; col < n_rows; ++col)
-                  {
-                    if (add)
-                      out[col * in_stride] += shape_values[col] * in[0];
-                    else
-                      out[col * in_stride] = shape_values[col] * in[0];
-                    if (max_derivative > 0)
-                      out[col * in_stride] +=
-                        shape_values[col + n_rows] * in[out_stride];
-                    if (max_derivative > 1)
-                      out[col * in_stride] +=
-                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
-                  }
-              }
-
-            // increment: in regular case, just go to the next point in
-            // x-direction. If we are at the end of one chunk in x-dir, need
-            // to jump over to the next layer in z-direction
-            switch (face_direction)
-              {
-                case 0:
-                  in += contract_onto_face ? n_rows : 1;
-                  out += contract_onto_face ? 1 : n_rows;
-                  break;
-                case 1:
-                  ++in;
-                  ++out;
-                  // faces 2 and 3 in 3d use local coordinate system zx, which
-                  // is the other way around compared to the tensor
-                  // product. Need to take that into account.
-                  if (dim == 3)
-                    {
-                      if (contract_onto_face)
-                        out += n_rows - 1;
-                      else
-                        in += n_rows - 1;
-                    }
-                  break;
-                case 2:
-                  ++in;
-                  ++out;
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
-          }
-
-        // adjust for local coordinate system zx
-        if (face_direction == 1 && dim == 3)
-          {
-            if (contract_onto_face)
-              {
-                in += n_rows * (n_rows - 1);
-                out -= n_rows * n_rows - 1;
-              }
-            else
-              {
-                out += n_rows * (n_rows - 1);
-                in -= n_rows * n_rows - 1;
-              }
-          }
-      }
+    constexpr int      stride = Utilities::pow(n_rows, face_direction);
+    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}};
+
+    interpolate_to_face<n_rows, stride, contract_to_face, add, max_derivative>(
+      this->shape_values,
+      {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}},
+      steps,
+      in,
+      out);
   }
 
 
@@ -1612,25 +1646,23 @@ namespace internal
       , n_columns(n_columns)
     {}
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int direction, bool contract_over_rows, bool add, int stride = 1>
     void
     values(const Number *in, Number *out) const
     {
-      apply<direction,
-            contract_over_rows,
-            add,
-            false,
-            EvaluatorQuantity::value>(shape_values, in, out);
+      constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value;
+      apply<direction, contract_over_rows, add, false, value_type, stride>(
+        shape_values, in, out);
     }
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int direction, bool contract_over_rows, bool add, int stride = 1>
     void
     gradients(const Number *in, Number *out) const
     {
       constexpr EvaluatorQuantity gradient_type =
         (variant != evaluate_evenodd ? EvaluatorQuantity::value :
                                        EvaluatorQuantity::gradient);
-      apply<direction, contract_over_rows, add, false, gradient_type>(
+      apply<direction, contract_over_rows, add, false, gradient_type, stride>(
         shape_gradients, in, out);
     }
 
@@ -1681,8 +1713,9 @@ namespace internal
     template <int               direction,
               bool              contract_over_rows,
               bool              add,
-              bool              one_line = false,
-              EvaluatorQuantity quantity = EvaluatorQuantity::value>
+              bool              one_line     = false,
+              EvaluatorQuantity quantity     = EvaluatorQuantity::value,
+              int               extra_stride = 1>
     void
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
           const Number                   *in,
@@ -1713,7 +1746,8 @@ namespace internal
             bool              contract_over_rows,
             bool              add,
             bool              one_line,
-            EvaluatorQuantity quantity>
+            EvaluatorQuantity quantity,
+            int               extra_stride>
   inline void
   EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply(
     const Number2 *DEAL_II_RESTRICT shape_data,
@@ -1741,6 +1775,8 @@ namespace internal
                             Utilities::fixed_power<dim - direction - 1>(n_rows);
     Assert(n_rows <= 128, ExcNotImplemented());
 
+    constexpr int stride_in  = !contract_over_rows ? extra_stride : 1;
+    constexpr int stride_out = contract_over_rows ? extra_stride : 1;
     for (int i2 = 0; i2 < n_blocks2; ++i2)
       {
         for (int i1 = 0; i1 < n_blocks1; ++i1)
@@ -1752,19 +1788,24 @@ namespace internal
             apply_matrix_vector_product<restricted_variant,
                                         quantity,
                                         contract_over_rows,
-                                        add>(
-              shape_data, in, out, n_rows, n_columns, stride, stride);
+                                        add>(shape_data,
+                                             in,
+                                             out,
+                                             n_rows,
+                                             n_columns,
+                                             stride * stride_in,
+                                             stride * stride_out);
 
             if (one_line == false)
               {
-                ++in;
-                ++out;
+                in += stride_in;
+                out += stride_out;
               }
           }
         if (one_line == false)
           {
-            in += stride * (mm - 1);
-            out += stride * (nn - 1);
+            in += stride * (mm - 1) * stride_in;
+            out += stride * (nn - 1) * stride_out;
           }
       }
   }
@@ -1776,7 +1817,7 @@ namespace internal
             typename Number,
             typename Number2>
   template <int  face_direction,
-            bool contract_onto_face,
+            bool contract_to_face,
             bool add,
             int  max_derivative>
   inline void
@@ -1788,116 +1829,28 @@ namespace internal
            ExcMessage(
              "The given array shape_data must not be the null pointer!"));
     static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
-    const int n_blocks1 = dim > 1 ? n_rows : 1;
-    const int n_blocks2 = dim > 2 ? n_rows : 1;
-
-    AssertIndexRange(face_direction, dim);
-    const int in_stride =
-      face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
-    const int out_stride =
-      dim > 1 ? Utilities::fixed_power<dim - 1>(n_rows) : 1;
-
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
-          {
-            if (contract_onto_face == true)
-              {
-                Number res0 = shape_values[0] * in[0];
-                Number res1, res2;
-                if (max_derivative > 0)
-                  res1 = shape_values[n_rows] * in[0];
-                if (max_derivative > 1)
-                  res2 = shape_values[2 * n_rows] * in[0];
-                for (unsigned int ind = 1; ind < n_rows; ++ind)
-                  {
-                    res0 += shape_values[ind] * in[in_stride * ind];
-                    if (max_derivative > 0)
-                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
-                    if (max_derivative > 1)
-                      res2 +=
-                        shape_values[ind + 2 * n_rows] * in[in_stride * ind];
-                  }
-                if (add)
-                  {
-                    out[0] += res0;
-                    if (max_derivative > 0)
-                      out[out_stride] += res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] += res2;
-                  }
-                else
-                  {
-                    out[0] = res0;
-                    if (max_derivative > 0)
-                      out[out_stride] = res1;
-                    if (max_derivative > 1)
-                      out[2 * out_stride] = res2;
-                  }
-              }
-            else
-              {
-                for (unsigned int col = 0; col < n_rows; ++col)
-                  {
-                    if (add)
-                      out[col * in_stride] += shape_values[col] * in[0];
-                    else
-                      out[col * in_stride] = shape_values[col] * in[0];
-                    if (max_derivative > 0)
-                      out[col * in_stride] +=
-                        shape_values[col + n_rows] * in[out_stride];
-                    if (max_derivative > 1)
-                      out[col * in_stride] +=
-                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
-                  }
-              }
 
-            // increment: in regular case, just go to the next point in
-            // x-direction. If we are at the end of one chunk in x-dir, need
-            // to jump over to the next layer in z-direction
-            switch (face_direction)
-              {
-                case 0:
-                  in += contract_onto_face ? n_rows : 1;
-                  out += contract_onto_face ? 1 : n_rows;
-                  break;
-                case 1:
-                  ++in;
-                  ++out;
-                  // faces 2 and 3 in 3d use local coordinate system zx, which
-                  // is the other way around compared to the tensor
-                  // product. Need to take that into account.
-                  if (dim == 3)
-                    {
-                      if (contract_onto_face)
-                        out += n_rows - 1;
-                      else
-                        in += n_rows - 1;
-                    }
-                  break;
-                case 2:
-                  ++in;
-                  ++out;
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
-          }
-        if (face_direction == 1 && dim == 3)
-          {
-            // adjust for local coordinate system zx
-            if (contract_onto_face)
-              {
-                in += n_rows * (n_rows - 1);
-                out -= n_rows * n_rows - 1;
-              }
-            else
-              {
-                out += n_rows * (n_rows - 1);
-                in -= n_rows * n_rows - 1;
-              }
-          }
-      }
+    const int          stride = Utilities::pow(n_rows, face_direction);
+    const int          n_rows = this->n_rows;
+    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}};
+
+    interpolate_to_face<0, 0, contract_to_face, add, max_derivative>(
+      this->shape_values,
+      {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}},
+      steps,
+      in,
+      out,
+      n_rows,
+      stride);
   }
 
 
@@ -2141,35 +2094,14 @@ namespace internal
       {
         for (int i1 = 0; i1 < n_blocks1; ++i1)
           {
-            Number x[mm];
-            for (int i = 0; i < mm; ++i)
-              x[i] = in[stride * i];
-
-            for (int col = 0; col < nn; ++col)
-              {
-                Number2 val0;
-
-                if (contract_over_rows)
-                  val0 = shape_data[col];
-                else
-                  val0 = shape_data[col * n_columns];
-
-                Number res0 = val0 * x[0];
-                for (int i = 1; i < mm; ++i)
-                  {
-                    if (contract_over_rows)
-                      val0 = shape_data[i * n_columns + col];
-                    else
-                      val0 = shape_data[col * n_columns + i];
-
-                    res0 += val0 * x[i];
-                  }
-                if (add)
-                  out[stride * col] += res0;
-
-                else
-                  out[stride * col] = res0;
-              }
+            apply_matrix_vector_product<evaluate_general,
+                                        EvaluatorQuantity::value,
+                                        n_rows,
+                                        n_columns,
+                                        stride,
+                                        stride,
+                                        contract_over_rows,
+                                        add>(shape_data, in, out);
 
             if (one_line == false)
               {
index 162c917802e5d8f1de20ac8b1fdcfcab7d3ff7d5..57367703291b72a003b6e537432a98bf950d51ba 100644 (file)
@@ -218,13 +218,13 @@ test_hessians(const dealii::FE_Poly<dim>                    &fe,
 
   // compare solutions of matrix vector product
   {
-    dst2 -= dst;
+    dst -= dst2;
 
     double error = 0.;
-    if (dst.l2_norm() > 0)
-      error = dst2.l2_norm() / dst.l2_norm();
+    if (dst2.l2_norm() > 0)
+      error = dst.l2_norm() / dst2.l2_norm();
     else
-      error = dst2.l2_norm();
+      error = dst.l2_norm();
 
     if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
       deallog << "FEValues verification: " << error << std::endl << std::endl;

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.