Move template normal_dir from function to struct

Creating new EvaluationTensorProductAnisotropic

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index fe15e2a34f3c7f3ecde946b78b4efc98ce00f8a0..d5d8e57a4da43a8a18a80fe08b799539b43fd28a 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -198,19 +198,6 @@ namespace internal
                           n_q_points_1d,
                           Number>
   {
-    using EvalNormal =
-      EvaluatorTensorProduct<evaluate_raviart_thomas,
-                             dim,
-                             (fe_degree == -1) ? 1 : fe_degree + 1,
-                             n_q_points_1d,
-                             Number>;
-
-    using EvalTangent =
-      EvaluatorTensorProduct<evaluate_raviart_thomas,
-                             dim,
-                             (fe_degree == -1) ? 1 : fe_degree,
-                             n_q_points_1d,
-                             Number>;
     template <bool integrate>
     static void
     evaluate_or_integrate(
@@ -223,22 +210,11 @@ namespace internal
     template <typename EvalType>
     static EvalType
     create_evaluator_tensor_product(
-      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+      const MatrixFreeFunctions::UnivariateShapeData<Number> &shape_data)
     {
-      // Shape_data is of size 2 with normal direction in [0] and tangential
-      // direction in [1]. See shape_info.template.h
-      const MatrixFreeFunctions::UnivariateShapeData<Number> *shape_data;
-      if (std::is_same<EvalType, EvalNormal>::value)
-        {
-          shape_data = &shape_info.data[0];
-        }
-      else
-        {
-          shape_data = &shape_info.data[1];
-        }
-      return EvalType(shape_data->shape_values,
-                      shape_data->shape_gradients,
-                      shape_data->shape_hessians);
+      return EvalType(shape_data.shape_values,
+                      shape_data.shape_gradients,
+                      shape_data.shape_hessians);
     }
 
     template <int normal_dir>
@@ -1061,6 +1037,22 @@ namespace internal
       std::integral_constant<bool, false>)
   {
     (void)add_into_values_array;
+
+    using EvalNormal =
+      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                        dim,
+                                        (fe_degree == -1) ? 1 : fe_degree + 1,
+                                        n_q_points_1d,
+                                        Number,
+                                        normal_dir>;
+
+    using EvalTangent =
+      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                        dim,
+                                        (fe_degree == -1) ? 1 : fe_degree,
+                                        n_q_points_1d,
+                                        Number,
+                                        normal_dir>;
     using Eval0 =
       typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
     using Eval1 =
@@ -1068,24 +1060,25 @@ namespace internal
     using Eval2 =
       typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
 
-    Eval0 eval0 =
-      create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info());
-    Eval1 eval1 =
-      create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info());
-    Eval2 eval2 =
-      create_evaluator_tensor_product<Eval2>(fe_eval.get_shape_info());
+    const MatrixFreeFunctions::ShapeInfo<Number> &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>(
-                               fe_eval.get_shape_info().data[0].fe_degree + 1),
-                             Utilities::fixed_power<dim>(
-                               fe_eval.get_shape_info().data[0].n_q_points_1d));
+    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 = fe_eval.get_shape_info().n_q_points;
-    const std::size_t dofs_per_comp =
-      fe_eval.get_shape_info().dofs_per_component_on_cell;
+    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;
@@ -1101,10 +1094,8 @@ namespace internal
         case 2:
           if (evaluation_flag & EvaluationFlags::gradients)
             {
-              eval0.template gradients<0, true, false, normal_dir>(values_dofs,
-                                                                   temp1);
-              eval1.template values<1, true, false, normal_dir>(temp1,
-                                                                gradients_quad);
+              eval0.template gradients<0, true, false>(values_dofs, temp1);
+              eval1.template values<1, true, false>(temp1, gradients_quad);
             }
           if (evaluation_flag & EvaluationFlags::hessians)
             {
@@ -1113,44 +1104,39 @@ namespace internal
               AssertThrow(false, ExcNotImplemented());
               // grad xy
               if (!(evaluation_flag & EvaluationFlags::gradients))
-                eval0.template gradients<0, true, false, normal_dir>(
-                  values_dofs, temp1);
-              eval1.template gradients<1, true, false, normal_dir>(
-                temp1, hessians_quad + 2 * n_q_points);
+                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, normal_dir>(values_dofs,
-                                                                  temp1);
-              eval1.template values<1, true, false, normal_dir>(temp1,
-                                                                hessians_quad);
+              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, normal_dir>(values_dofs, temp1);
+          eval0.template values<0, true, false>(values_dofs, temp1);
           if (evaluation_flag & EvaluationFlags::gradients)
-            eval1.template gradients<1, true, false, normal_dir>(
-              temp1, gradients_quad + n_q_points);
+            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, normal_dir>(temp1,
-                                                                hessians_quad +
-                                                                  n_q_points);
+            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, normal_dir>(temp1,
-                                                              values_quad);
+            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, normal_dir>(values_dofs,
-                                                                   temp1);
-              eval1.template values<1, true, false, normal_dir>(temp1, temp2);
-              eval2.template values<2, true, false, normal_dir>(temp2,
-                                                                gradients_quad);
+              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)
@@ -1161,73 +1147,69 @@ namespace internal
               // grad xz
               if (!(evaluation_flag & EvaluationFlags::gradients))
                 {
-                  eval0.template gradients<0, true, false, normal_dir>(
-                    values_dofs, temp1);
-                  eval1.template values<1, true, false, normal_dir>(temp1,
-                                                                    temp2);
+                  eval0.template gradients<0, true, false>(values_dofs, temp1);
+                  eval1.template values<1, true, false>(temp1, temp2);
                 }
-              eval2.template gradients<2, true, false, normal_dir>(
-                temp2, hessians_quad + 4 * n_q_points);
+              eval2.template gradients<2, true, false>(temp2,
+                                                       hessians_quad +
+                                                         4 * n_q_points);
 
               // grad xy
-              eval1.template gradients<1, true, false, normal_dir>(temp1,
-                                                                   temp2);
-              eval2.template values<2, true, false, normal_dir>(
-                temp2, hessians_quad + 3 * n_q_points);
+              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, normal_dir>(values_dofs,
-                                                                  temp1);
-              eval1.template values<1, true, false, normal_dir>(temp1, temp2);
-              eval2.template values<2, true, false, normal_dir>(temp2,
-                                                                hessians_quad);
+              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, normal_dir>(values_dofs, temp1);
+          eval0.template values<0, true, false>(values_dofs, temp1);
           if (evaluation_flag & EvaluationFlags::gradients)
             {
-              eval1.template gradients<1, true, false, normal_dir>(temp1,
-                                                                   temp2);
-              eval2.template values<2, true, false, normal_dir>(temp2,
-                                                                gradients_quad +
-                                                                  n_q_points);
+              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, normal_dir>(temp1,
-                                                                     temp2);
-              eval2.template gradients<2, true, false, normal_dir>(
-                temp2, hessians_quad + 5 * n_q_points);
+                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, normal_dir>(temp1, temp2);
-              eval2.template values<2, true, false, normal_dir>(temp2,
-                                                                hessians_quad +
-                                                                  n_q_points);
+              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, normal_dir>(temp1, temp2);
+          eval1.template values<1, true, false>(temp1, temp2);
           if (evaluation_flag & EvaluationFlags::gradients)
-            eval2.template gradients<2, true, false, normal_dir>(
-              temp2, gradients_quad + 2 * n_q_points);
+            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, normal_dir>(
-              temp2, hessians_quad + 2 * n_q_points);
+            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, normal_dir>(temp2,
-                                                              values_quad);
+            eval2.template values<2, true, false>(temp2, values_quad);
           break;
         default:
           AssertThrow(false, ExcNotImplemented());
@@ -1249,6 +1231,21 @@ namespace internal
       const bool                             add_into_values_array,
       std::integral_constant<bool, true>)
   {
+    using EvalNormal =
+      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                        dim,
+                                        (fe_degree == -1) ? 1 : fe_degree + 1,
+                                        n_q_points_1d,
+                                        Number,
+                                        normal_dir>;
+
+    using EvalTangent =
+      EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                        dim,
+                                        (fe_degree == -1) ? 1 : fe_degree,
+                                        n_q_points_1d,
+                                        Number,
+                                        normal_dir>;
     using Eval0 =
       typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
     using Eval1 =
@@ -1256,24 +1253,25 @@ namespace internal
     using Eval2 =
       typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
 
-    Eval0 eval0 =
-      create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info());
-    Eval1 eval1 =
-      create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info());
-    Eval2 eval2 =
-      create_evaluator_tensor_product<Eval2>(fe_eval.get_shape_info());
+    const MatrixFreeFunctions::ShapeInfo<Number> &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>(
-                               fe_eval.get_shape_info().data[0].fe_degree + 1),
-                             Utilities::fixed_power<dim>(
-                               fe_eval.get_shape_info().data[0].n_q_points_1d));
+    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 = fe_eval.get_shape_info().n_q_points;
-    const std::size_t dofs_per_comp =
-      fe_eval.get_shape_info().dofs_per_component_on_cell;
+    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;
@@ -1291,59 +1289,49 @@ namespace internal
           if ((evaluation_flag & EvaluationFlags::values) &&
               !(evaluation_flag & EvaluationFlags::gradients))
             {
-              eval1.template values<1, false, false, normal_dir>(values_quad,
-                                                                 temp1);
+              eval1.template values<1, false, false>(values_quad, temp1);
               if (add_into_values_array == false)
-                eval0.template values<0, false, false, normal_dir>(temp1,
-                                                                   values_dofs);
+                eval0.template values<0, false, false>(temp1, values_dofs);
               else
-                eval0.template values<0, false, true, normal_dir>(temp1,
-                                                                  values_dofs);
+                eval0.template values<0, false, true>(temp1, values_dofs);
             }
           if (evaluation_flag & EvaluationFlags::gradients)
             {
-              eval1.template gradients<1, false, false, normal_dir>(
-                gradients_quad + n_q_points, temp1);
+              eval1.template gradients<1, false, false>(gradients_quad +
+                                                          n_q_points,
+                                                        temp1);
               if ((evaluation_flag & EvaluationFlags::values))
-                eval1.template values<1, false, true, normal_dir>(values_quad,
-                                                                  temp1);
+                eval1.template values<1, false, true>(values_quad, temp1);
               if (add_into_values_array == false)
-                eval0.template values<0, false, false, normal_dir>(temp1,
-                                                                   values_dofs);
+                eval0.template values<0, false, false>(temp1, values_dofs);
               else
-                eval0.template values<0, false, true, normal_dir>(temp1,
-                                                                  values_dofs);
-              eval1.template values<1, false, false, normal_dir>(gradients_quad,
-                                                                 temp1);
-              eval0.template gradients<0, false, true, normal_dir>(temp1,
-                                                                   values_dofs);
+                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, normal_dir>(hessians_quad,
-                                                                 temp1);
+              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, normal_dir>(
-                  temp1, values_dofs);
+                eval0.template hessians<0, false, true>(temp1, values_dofs);
               else
-                eval0.template hessians<0, false, false, normal_dir>(
-                  temp1, values_dofs);
+                eval0.template hessians<0, false, false>(temp1, values_dofs);
 
               // grad yy
-              eval1.template hessians<1, false, false, normal_dir>(
-                hessians_quad + n_q_points, temp1);
-              eval0.template values<0, false, true, normal_dir>(temp1,
-                                                                values_dofs);
+              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, normal_dir>(
-                hessians_quad + 2 * n_q_points, temp1);
-              eval0.template gradients<0, false, true, normal_dir>(temp1,
-                                                                   values_dofs);
+              eval1.template gradients<1, false, false>(hessians_quad +
+                                                          2 * n_q_points,
+                                                        temp1);
+              eval0.template gradients<0, false, true>(temp1, values_dofs);
             }
           break;
 
@@ -1351,94 +1339,79 @@ namespace internal
           if ((evaluation_flag & EvaluationFlags::values) &&
               !(evaluation_flag & EvaluationFlags::gradients))
             {
-              eval2.template values<2, false, false, normal_dir>(values_quad,
-                                                                 temp1);
-              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              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, normal_dir>(temp2,
-                                                                   values_dofs);
+                eval0.template values<0, false, false>(temp2, values_dofs);
               else
-                eval0.template values<0, false, true, normal_dir>(temp2,
-                                                                  values_dofs);
+                eval0.template values<0, false, true>(temp2, values_dofs);
             }
           if (evaluation_flag & EvaluationFlags::gradients)
             {
-              eval2.template gradients<2, false, false, normal_dir>(
-                gradients_quad + 2 * n_q_points, temp1);
+              eval2.template gradients<2, false, false>(gradients_quad +
+                                                          2 * n_q_points,
+                                                        temp1);
               if ((evaluation_flag & EvaluationFlags::values))
-                eval2.template values<2, false, true, normal_dir>(values_quad,
-                                                                  temp1);
-              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
-              eval2.template values<2, false, false, normal_dir>(
-                gradients_quad + n_q_points, temp1);
-              eval1.template gradients<1, false, true, normal_dir>(temp1,
-                                                                   temp2);
+                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, normal_dir>(temp2,
-                                                                   values_dofs);
+                eval0.template values<0, false, false>(temp2, values_dofs);
               else
-                eval0.template values<0, false, true, normal_dir>(temp2,
-                                                                  values_dofs);
-              eval2.template values<2, false, false, normal_dir>(gradients_quad,
-                                                                 temp1);
-              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
-              eval0.template gradients<0, false, true, normal_dir>(temp2,
-                                                                   values_dofs);
+                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, normal_dir>(hessians_quad,
-                                                                 temp1);
-              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              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, normal_dir>(
-                  temp2, values_dofs);
+                eval0.template hessians<0, false, true>(temp2, values_dofs);
               else
-                eval0.template hessians<0, false, false, normal_dir>(
-                  temp2, values_dofs);
+                eval0.template hessians<0, false, false>(temp2, values_dofs);
 
               // grad yy
-              eval2.template values<2, false, false, normal_dir>(hessians_quad +
-                                                                   n_q_points,
-                                                                 temp1);
-              eval1.template hessians<1, false, false, normal_dir>(temp1,
-                                                                   temp2);
-              eval0.template values<0, false, true, normal_dir>(temp2,
-                                                                values_dofs);
+              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, normal_dir>(
-                hessians_quad + 2 * n_q_points, temp1);
-              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
-              eval0.template values<0, false, true, normal_dir>(temp2,
-                                                                values_dofs);
+              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, normal_dir>(
-                hessians_quad + 3 * n_q_points, temp1);
-              eval1.template gradients<1, false, false, normal_dir>(temp1,
-                                                                    temp2);
-              eval0.template gradients<0, false, true, normal_dir>(temp2,
-                                                                   values_dofs);
+              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, normal_dir>(
-                hessians_quad + 4 * n_q_points, temp1);
-              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
-              eval0.template gradients<0, false, true, normal_dir>(temp2,
-                                                                   values_dofs);
+              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, normal_dir>(
-                hessians_quad + 5 * n_q_points, temp1);
-              eval1.template gradients<1, false, false, normal_dir>(temp1,
-                                                                    temp2);
-              eval0.template values<0, false, true, normal_dir>(temp2,
-                                                                values_dofs);
+              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;
@@ -3035,19 +3008,6 @@ namespace internal
   template <int dim, int fe_degree, int n_q_points_1d, typename Number>
   struct FEFaceEvaluationImplRaviartThomas
   {
-    using EvalNormal =
-      EvaluatorTensorProduct<evaluate_raviart_thomas,
-                             dim - 1,
-                             (fe_degree == -1) ? 1 : fe_degree + 1,
-                             n_q_points_1d,
-                             Number>;
-    using EvalTangent =
-      EvaluatorTensorProduct<evaluate_raviart_thomas,
-                             dim - 1,
-                             (fe_degree == -1) ? 1 : fe_degree,
-                             n_q_points_1d,
-                             Number>;
-
     using EvalGeneral = EvaluatorTensorProduct<evaluate_general,
                                                dim - 1,
                                                fe_degree,
@@ -3056,21 +3016,16 @@ namespace internal
     template <typename EvalType>
     static EvalType
     create_evaluator_tensor_product(
-      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-      const unsigned int                            subface_index,
-      const unsigned int                            direction)
+      const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+      const unsigned int                                      subface_index,
+      const unsigned int                                      direction)
     {
-      const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
-        (std::is_same<EvalType, EvalNormal>::value) ? shape_info.data.front() :
-                                                      shape_info.data.back();
       if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
         return EvalType(data.shape_values,
                         data.shape_gradients,
                         data.shape_hessians);
       else
         {
-          Assert(false, ExcNotImplemented());
-
           const unsigned int index =
             direction == 0 ? subface_index % 2 : subface_index / 2;
           return EvalType(data.values_within_subface[index],
@@ -3091,33 +3046,48 @@ namespace internal
     {
       const unsigned int face_direction = face_no / 2;
 
-      evaluate_in_face_apply<0, EvalNormal, EvalTangent>(
-        values_dofs,
-        fe_eval,
-        scratch_data,
-        evaluation_flag,
-        face_direction,
-        subface_index,
-        std::integral_constant<bool, integrate>());
+      // We first evaluate the anisotropic faces, i.e the faces where
+      // face_direction != component. Note that the call order here is not
+      // important, since the pointers are shifted accordingly within the
+      // function. However, this is the order in which the components will be in
+      // the quadrature points. Furthermore, the isotropic faces have no "normal
+      // direction" but we still pass in normal_dir = 2 since this is used for
+      // the pointers.
+      // -----------------------------------------------------------------------------------
+      // |          |                   Anisotropic faces                 | Isotropic faces|
+      // | Face dir | comp, coords, normal_dir | comp, coords, normal_dir | comp, coords   |
+      // | --------------------------------------------------------------------------------|
+      // |    0     | 1, y, 0                  | -                        | 0, y           |
+      // |    1     | 0, x, 0                  | -                        | 1, x           |
+      // | --------------------------------------------------------------------------------|
+      // |    0     | 1, yz, 0                 | 2, yz, 1                 | 0, yz          |
+      // |    1     | 2, zx, 0                 | 0, zx, 1                 | 1, zx          |
+      // |    2     | 0, xy, 0                 | 1, xy, 1                 | 2, xy          |
+      // -----------------------------------------------------------------------------------
+      evaluate_in_face_apply<0>(values_dofs,
+                                fe_eval,
+                                scratch_data,
+                                evaluation_flag,
+                                face_direction,
+                                subface_index,
+                                std::integral_constant<bool, integrate>());
 
       if (dim == 3)
-        evaluate_in_face_apply<1, EvalTangent, EvalNormal>(
-          values_dofs,
-          fe_eval,
-          scratch_data,
-          evaluation_flag,
-          face_direction,
-          subface_index,
-          std::integral_constant<bool, integrate>());
-
-      evaluate_in_face_apply<2, EvalGeneral, EvalGeneral>(
-        values_dofs,
-        fe_eval,
-        scratch_data,
-        evaluation_flag,
-        face_direction,
-        subface_index,
-        std::integral_constant<bool, integrate>());
+        evaluate_in_face_apply<1>(values_dofs,
+                                  fe_eval,
+                                  scratch_data,
+                                  evaluation_flag,
+                                  face_direction,
+                                  subface_index,
+                                  std::integral_constant<bool, integrate>());
+
+      evaluate_in_face_apply<2>(values_dofs,
+                                fe_eval,
+                                scratch_data,
+                                evaluation_flag,
+                                face_direction,
+                                subface_index,
+                                std::integral_constant<bool, integrate>());
     }
 
     /*
@@ -3128,7 +3098,7 @@ namespace internal
      * integration path. These two functions can be fused together since all
      * offsets and pointers are the exact same.
      */
-    template <int normal_dir, typename Eval0, typename Eval1>
+    template <int normal_dir>
     static inline void
     evaluate_in_face_apply(
       Number *                               values_dofs,
@@ -3139,20 +3109,43 @@ namespace internal
       const unsigned int                     subface_index,
       std::integral_constant<bool, false>)
     {
-      // TODO. Need to make sure hanging nodes work, i.e call
-      // create_eval_tensor_product with correct direction.
-      Eval0 eval0 =
-        create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info(),
-                                               subface_index,
-                                               0);
-      Eval1 eval1 =
-        create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info(),
-                                               subface_index,
-                                               0);
+      using EvalNormal =
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim - 1,
+                                          (fe_degree == -1) ? 1 : fe_degree + 1,
+                                          n_q_points_1d,
+                                          Number,
+                                          normal_dir>;
+      using EvalTangent =
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim - 1,
+                                          (fe_degree == -1) ? 1 : fe_degree,
+                                          n_q_points_1d,
+                                          Number,
+                                          normal_dir>;
+
+      using TempEval0 = typename std::
+        conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
+      using TempEval1 = typename std::
+        conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
+      using Eval0 = typename std::
+        conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
+      using Eval1 = typename std::
+        conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
+
+      const MatrixFreeFunctions::ShapeInfo<Number> &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]),
+        subface_index,
+        0);
+      Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
+        ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
+        subface_index,
+        1);
 
       constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
-      const std::size_t     n_dofs_tangent =
-        fe_eval.get_shape_info().dofs_per_component_on_face;
+      const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
       const std::size_t n_dofs_normal =
         n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
       const std::size_t dofs_stride =
@@ -3188,14 +3181,11 @@ namespace internal
           switch (dim)
             {
               case 3:
-                eval0.template values<0, true, false, normal_dir>(values_dofs,
-                                                                  values_quad);
-                eval1.template values<1, true, false, normal_dir>(values_quad,
-                                                                  values_quad);
+                eval0.template values<0, true, false>(values_dofs, values_quad);
+                eval1.template values<1, true, false>(values_quad, values_quad);
                 break;
               case 2:
-                eval0.template values<0, true, false, normal_dir>(values_dofs,
-                                                                  values_quad);
+                eval0.template values<0, true, false>(values_dofs, values_quad);
                 break;
               default:
                 Assert(false, ExcNotImplemented());
@@ -3207,37 +3197,39 @@ namespace internal
             {
               case 3:
                 // grad x
-                eval0.template gradients<0, true, false, normal_dir>(
-                  values_dofs, scratch_data);
-                eval1.template values<1, true, false, normal_dir>(
-                  scratch_data, gradients_quad);
+                eval0.template gradients<0, true, false>(values_dofs,
+                                                         scratch_data);
+                eval1.template values<1, true, false>(scratch_data,
+                                                      gradients_quad);
 
                 // grad y
-                eval0.template values<0, true, false, normal_dir>(values_dofs,
-                                                                  scratch_data);
-                eval1.template gradients<1, true, false, normal_dir>(
-                  scratch_data, gradients_quad + n_q_points);
+                eval0.template values<0, true, false>(values_dofs,
+                                                      scratch_data);
+                eval1.template gradients<1, true, false>(scratch_data,
+                                                         gradients_quad +
+                                                           n_q_points);
 
                 if (evaluation_flag & EvaluationFlags::values)
-                  eval1.template values<1, true, false, normal_dir>(
-                    scratch_data, values_quad);
+                  eval1.template values<1, true, false>(scratch_data,
+                                                        values_quad);
 
                 // grad z
-                eval0.template values<0, true, false, normal_dir>(values_dofs +
-                                                                    dofs_stride,
-                                                                  scratch_data);
-                eval1.template values<1, true, false, normal_dir>(
-                  scratch_data, gradients_quad + 2 * n_q_points);
+                eval0.template values<0, true, false>(values_dofs + dofs_stride,
+                                                      scratch_data);
+                eval1.template values<1, true, false>(scratch_data,
+                                                      gradients_quad +
+                                                        2 * n_q_points);
 
                 break;
               case 2:
-                eval0.template values<0, true, false, normal_dir>(
-                  values_dofs + dofs_stride, gradients_quad + n_q_points);
-                eval0.template gradients<0, true, false, normal_dir>(
-                  values_dofs, gradients_quad);
+                eval0.template values<0, true, false>(values_dofs + dofs_stride,
+                                                      gradients_quad +
+                                                        n_q_points);
+                eval0.template gradients<0, true, false>(values_dofs,
+                                                         gradients_quad);
                 if ((evaluation_flag & EvaluationFlags::values))
-                  eval0.template values<0, true, false, normal_dir>(
-                    values_dofs, values_quad);
+                  eval0.template values<0, true, false>(values_dofs,
+                                                        values_quad);
                 break;
               default:
                 AssertThrow(false, ExcNotImplemented());
@@ -3250,52 +3242,58 @@ namespace internal
             {
               case 3:
                 // grad xx
-                eval0.template hessians<0, true, false, normal_dir>(
-                  values_dofs, scratch_data);
-                eval1.template values<1, true, false, normal_dir>(
-                  scratch_data, hessians_quad);
+                eval0.template hessians<0, true, false>(values_dofs,
+                                                        scratch_data);
+                eval1.template values<1, true, false>(scratch_data,
+                                                      hessians_quad);
 
                 // grad yy
-                eval0.template values<0, true, false, normal_dir>(values_dofs,
-                                                                  scratch_data);
-                eval1.template hessians<1, true, false, normal_dir>(
-                  scratch_data, hessians_quad + n_q_points);
+                eval0.template values<0, true, false>(values_dofs,
+                                                      scratch_data);
+                eval1.template hessians<1, true, false>(scratch_data,
+                                                        hessians_quad +
+                                                          n_q_points);
 
                 // grad zz
-                eval0.template values<0, true, false, normal_dir>(
-                  values_dofs + 2 * dofs_stride, scratch_data);
-                eval1.template values<1, true, false, normal_dir>(
-                  scratch_data, hessians_quad + 2 * n_q_points);
+                eval0.template values<0, true, false>(values_dofs +
+                                                        2 * dofs_stride,
+                                                      scratch_data);
+                eval1.template values<1, true, false>(scratch_data,
+                                                      hessians_quad +
+                                                        2 * n_q_points);
 
                 // grad xy
-                eval0.template gradients<0, true, false, normal_dir>(
-                  values_dofs, scratch_data);
-                eval1.template gradients<1, true, false, normal_dir>(
-                  scratch_data, hessians_quad + 3 * n_q_points);
+                eval0.template gradients<0, true, false>(values_dofs,
+                                                         scratch_data);
+                eval1.template gradients<1, true, false>(scratch_data,
+                                                         hessians_quad +
+                                                           3 * n_q_points);
 
                 // grad xz
-                eval0.template gradients<0, true, false, normal_dir>(
-                  values_dofs + dofs_stride, scratch_data);
-                eval1.template values<1, true, false, normal_dir>(
-                  scratch_data, hessians_quad + 4 * n_q_points);
+                eval0.template gradients<0, true, false>(values_dofs +
+                                                           dofs_stride,
+                                                         scratch_data);
+                eval1.template values<1, true, false>(scratch_data,
+                                                      hessians_quad +
+                                                        4 * n_q_points);
 
                 // grad yz
-                eval0.template values<0, true, false, normal_dir>(values_dofs +
-                                                                    dofs_stride,
-                                                                  scratch_data);
-                eval1.template gradients<1, true, false, normal_dir>(
-                  scratch_data, hessians_quad + 5 * n_q_points);
+                eval0.template values<0, true, false>(values_dofs + dofs_stride,
+                                                      scratch_data);
+                eval1.template gradients<1, true, false>(scratch_data,
+                                                         hessians_quad +
+                                                           5 * n_q_points);
 
                 break;
               case 2:
                 // grad xx
-                eval0.template hessians<0, true, false, normal_dir>(
-                  values_dofs, hessians_quad);
+                eval0.template hessians<0, true, false>(values_dofs,
+                                                        hessians_quad);
                 // grad yy
-                eval0.template values<0, true, false, normal_dir>(
+                eval0.template values<0, true, false>(
                   values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
                 // grad xy
-                eval0.template gradients<0, true, false, normal_dir>(
+                eval0.template gradients<0, true, false>(
                   values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
                 break;
               default:
@@ -3304,7 +3302,7 @@ namespace internal
         }
     }
 
-    template <int normal_dir, typename Eval0, typename Eval1>
+    template <int normal_dir>
     static inline void
     evaluate_in_face_apply(
       Number *                               values_dofs,
@@ -3315,20 +3313,43 @@ namespace internal
       const unsigned int                     subface_index,
       std::integral_constant<bool, true>)
     {
-      // TODO. Need to make sure hanging nodes work, i.e call
-      // create_eval_tensor_product with correct direction.
-      Eval0 eval0 =
-        create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info(),
-                                               subface_index,
-                                               0);
-      Eval1 eval1 =
-        create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info(),
-                                               subface_index,
-                                               0);
+      using EvalNormal =
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim - 1,
+                                          (fe_degree == -1) ? 1 : fe_degree + 1,
+                                          n_q_points_1d,
+                                          Number,
+                                          normal_dir>;
+      using EvalTangent =
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim - 1,
+                                          (fe_degree == -1) ? 1 : fe_degree,
+                                          n_q_points_1d,
+                                          Number,
+                                          normal_dir>;
+
+      using TempEval0 = typename std::
+        conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
+      using TempEval1 = typename std::
+        conditional<normal_dir == 0, EvalTangent, EvalNormal>::type;
+      using Eval0 = typename std::
+        conditional<normal_dir == 2, EvalGeneral, TempEval0>::type;
+      using Eval1 = typename std::
+        conditional<normal_dir == 2, EvalGeneral, TempEval1>::type;
+
+      const MatrixFreeFunctions::ShapeInfo<Number> &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]),
+        subface_index,
+        0);
+      Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
+        ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]),
+        subface_index,
+        1);
 
       constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
-      const std::size_t     n_dofs_tangent =
-        fe_eval.get_shape_info().dofs_per_component_on_face;
+      const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
       const std::size_t n_dofs_normal =
         n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
       const std::size_t dofs_stride =
@@ -3364,14 +3385,14 @@ namespace internal
           switch (dim)
             {
               case 3:
-                eval1.template values<1, false, false, normal_dir>(values_quad,
-                                                                   values_quad);
-                eval0.template values<0, false, false, normal_dir>(values_quad,
-                                                                   values_dofs);
+                eval1.template values<1, false, false>(values_quad,
+                                                       values_quad);
+                eval0.template values<0, false, false>(values_quad,
+                                                       values_dofs);
                 break;
               case 2:
-                eval0.template values<0, false, false, normal_dir>(values_quad,
-                                                                   values_dofs);
+                eval0.template values<0, false, false>(values_quad,
+                                                       values_dofs);
                 break;
               default:
                 Assert(false, ExcNotImplemented());
@@ -3383,42 +3404,45 @@ namespace internal
             {
               case 3:
                 // grad z
-                eval1.template values<1, false, false, normal_dir>(
-                  gradients_quad + 2 * n_q_points,
-                  gradients_quad + 2 * n_q_points);
-                eval0.template values<0, false, false, normal_dir>(
+                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, values_dofs + dofs_stride);
 
                 if (evaluation_flag & EvaluationFlags::values)
                   {
-                    eval1.template values<1, false, false, normal_dir>(
-                      values_quad, scratch_data);
-                    eval1.template gradients<1, false, true, normal_dir>(
-                      gradients_quad + n_q_points, scratch_data);
+                    eval1.template values<1, false, false>(values_quad,
+                                                           scratch_data);
+                    eval1.template gradients<1, false, true>(gradients_quad +
+                                                               n_q_points,
+                                                             scratch_data);
                   }
                 else
-                  eval1.template gradients<1, false, false, normal_dir>(
-                    gradients_quad + n_q_points, scratch_data);
+                  eval1.template gradients<1, false, false>(gradients_quad +
+                                                              n_q_points,
+                                                            scratch_data);
 
                 // grad y
-                eval0.template values<0, false, false, normal_dir>(scratch_data,
-                                                                   values_dofs);
+                eval0.template values<0, false, false>(scratch_data,
+                                                       values_dofs);
 
                 // grad x
-                eval1.template values<1, false, false, normal_dir>(
-                  gradients_quad, scratch_data);
-                eval0.template gradients<0, false, true, normal_dir>(
-                  scratch_data, values_dofs);
+                eval1.template values<1, false, false>(gradients_quad,
+                                                       scratch_data);
+                eval0.template gradients<0, false, true>(scratch_data,
+                                                         values_dofs);
 
                 break;
               case 2:
-                eval0.template values<0, false, false, normal_dir>(
+                eval0.template values<0, false, false>(
                   gradients_quad + n_q_points, values_dofs + dofs_stride);
-                eval0.template gradients<0, false, false, normal_dir>(
-                  gradients_quad, values_dofs);
+                eval0.template gradients<0, false, false>(gradients_quad,
+                                                          values_dofs);
                 if (evaluation_flag & EvaluationFlags::values)
-                  eval0.template values<0, false, true, normal_dir>(
-                    values_quad, values_dofs);
+                  eval0.template values<0, false, true>(values_quad,
+                                                        values_dofs);
                 break;
               default:
                 AssertThrow(false, ExcNotImplemented());
@@ -3431,70 +3455,79 @@ namespace internal
             {
               case 3:
                 // grad xx
-                eval1.template values<1, false, false, normal_dir>(
-                  hessians_quad, scratch_data);
+                eval1.template values<1, false, false>(hessians_quad,
+                                                       scratch_data);
                 if ((evaluation_flag &
                      (EvaluationFlags::values | EvaluationFlags::gradients)))
-                  eval0.template hessians<0, false, true, normal_dir>(
-                    scratch_data, values_dofs);
+                  eval0.template hessians<0, false, true>(scratch_data,
+                                                          values_dofs);
                 else
-                  eval0.template hessians<0, false, false, normal_dir>(
-                    scratch_data, values_dofs);
+                  eval0.template hessians<0, false, false>(scratch_data,
+                                                           values_dofs);
 
                 // grad yy
-                eval1.template hessians<1, false, false, normal_dir>(
-                  hessians_quad + n_q_points, scratch_data);
-                eval0.template values<0, false, true, normal_dir>(scratch_data,
-                                                                  values_dofs);
+                eval1.template hessians<1, false, false>(hessians_quad +
+                                                           n_q_points,
+                                                         scratch_data);
+                eval0.template values<0, false, true>(scratch_data,
+                                                      values_dofs);
 
                 // grad zz
-                eval1.template values<1, false, false, normal_dir>(
-                  hessians_quad + 2 * n_q_points, scratch_data);
-                eval0.template values<0, false, false, normal_dir>(
-                  scratch_data, values_dofs + 2 * dofs_stride);
+                eval1.template values<1, false, false>(hessians_quad +
+                                                         2 * n_q_points,
+                                                       scratch_data);
+                eval0.template values<0, false, false>(scratch_data,
+                                                       values_dofs +
+                                                         2 * dofs_stride);
 
                 // grad xy
-                eval1.template gradients<1, false, false, normal_dir>(
-                  hessians_quad + 3 * n_q_points, scratch_data);
-                eval0.template gradients<0, false, true, normal_dir>(
-                  scratch_data, values_dofs);
+                eval1.template gradients<1, false, false>(hessians_quad +
+                                                            3 * n_q_points,
+                                                          scratch_data);
+                eval0.template gradients<0, false, true>(scratch_data,
+                                                         values_dofs);
 
                 // grad xz
-                eval1.template values<1, false, false, normal_dir>(
-                  hessians_quad + 4 * n_q_points, scratch_data);
+                eval1.template values<1, false, false>(hessians_quad +
+                                                         4 * n_q_points,
+                                                       scratch_data);
                 if ((evaluation_flag & EvaluationFlags::gradients))
-                  eval0.template gradients<0, false, true, normal_dir>(
-                    scratch_data, values_dofs + dofs_stride);
+                  eval0.template gradients<0, false, true>(scratch_data,
+                                                           values_dofs +
+                                                             dofs_stride);
                 else
-                  eval0.template gradients<0, false, false, normal_dir>(
-                    scratch_data, values_dofs + dofs_stride);
+                  eval0.template gradients<0, false, false>(scratch_data,
+                                                            values_dofs +
+                                                              dofs_stride);
 
                 // grad yz
-                eval1.template gradients<1, false, false, normal_dir>(
-                  hessians_quad + 5 * n_q_points, scratch_data);
-                eval0.template values<0, false, true, normal_dir>(
-                  scratch_data, values_dofs + dofs_stride);
+                eval1.template gradients<1, false, false>(hessians_quad +
+                                                            5 * n_q_points,
+                                                          scratch_data);
+                eval0.template values<0, false, true>(scratch_data,
+                                                      values_dofs +
+                                                        dofs_stride);
 
                 break;
               case 2:
                 // grad xx
                 if (evaluation_flag &
                     (EvaluationFlags::values | EvaluationFlags::gradients))
-                  eval0.template hessians<0, false, true, normal_dir>(
-                    hessians_quad, values_dofs);
+                  eval0.template hessians<0, false, true>(hessians_quad,
+                                                          values_dofs);
                 else
-                  eval0.template hessians<0, false, false, normal_dir>(
-                    hessians_quad, values_dofs);
+                  eval0.template hessians<0, false, false>(hessians_quad,
+                                                           values_dofs);
 
                 // grad yy
-                eval0.template values<0, false, false, normal_dir>(
+                eval0.template values<0, false, false>(
                   hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
                 // grad xy
                 if ((evaluation_flag & EvaluationFlags::gradients))
-                  eval0.template gradients<0, false, true, normal_dir>(
+                  eval0.template gradients<0, false, true>(
                     hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
                 else
-                  eval0.template gradients<0, false, false, normal_dir>(
+                  eval0.template gradients<0, false, false>(
                     hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
                 break;
               default:
@@ -3508,19 +3541,6 @@ namespace internal
   template <int dim, int fe_degree, typename Number, bool lex_faces = false>
   struct FEFaceNormalEvaluationImpl
   {
-    using EvalNormal =
-      EvaluatorTensorProduct<evaluate_raviart_thomas,
-                             dim,
-                             (fe_degree == -1) ? 1 : fe_degree + 1,
-                             0,
-                             Number>;
-    using EvalTangent =
-      EvaluatorTensorProduct<evaluate_raviart_thomas,
-                             dim,
-                             (fe_degree == -1) ? 1 : fe_degree,
-                             0,
-                             Number>;
-
     template <bool do_evaluate, bool add_into_output>
     static void
     interpolate(const unsigned int                            n_components,
@@ -3658,21 +3678,10 @@ namespace internal
     template <typename EvalType>
     static EvalType
     create_evaluator_tensor_product(
-      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
-      const unsigned int                            face_no)
+      const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+      const unsigned int                                      face_no)
     {
-      // Shape_data is of size 2 with normal direction in [0] and tangential
-      // direction in [1]. See shape_info.template.h
-      const MatrixFreeFunctions::UnivariateShapeData<Number> *shape_data;
-      if (std::is_same<EvalType, EvalNormal>::value)
-        {
-          shape_data = &shape_info.data[0];
-        }
-      else
-        {
-          shape_data = &shape_info.data[1];
-        }
-      return EvalType(shape_data->shape_data_on_face[face_no % 2],
+      return EvalType(data.shape_data_on_face[face_no % 2],
                       AlignedVector<Number>(),
                       AlignedVector<Number>());
     }
@@ -3747,19 +3756,67 @@ namespace internal
       const Number *                                input,
       Number *                                      output)
     {
-      using Evalf0 = typename std::
-        conditional<face_direction == 0, EvalNormal, EvalTangent>::type;
-      using Evalf1 = typename std::
-        conditional<face_direction == 1, EvalNormal, EvalTangent>::type;
-      using Evalf2 = typename std::
-        conditional<face_direction == 2, EvalNormal, EvalTangent>::type;
+      // These types are evaluators in either normal or tangential direction
+      // depending on the face direction, with different normal directions for
+      // the different components.
+      using Evalf0 = typename std::conditional<
+        face_direction == 0,
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim,
+                                          (fe_degree == -1) ? 1 : fe_degree + 1,
+                                          0,
+                                          Number,
+                                          0>,
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim,
+                                          (fe_degree == -1) ? 1 : fe_degree,
+                                          0,
+                                          Number,
+                                          0>>::type;
+      using Evalf1 = typename std::conditional<
+        face_direction == 1,
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim,
+                                          (fe_degree == -1) ? 1 : fe_degree + 1,
+                                          0,
+                                          Number,
+                                          1>,
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim,
+                                          (fe_degree == -1) ? 1 : fe_degree,
+                                          0,
+                                          Number,
+                                          1>>::type;
+      using Evalf2 = typename std::conditional<
+        face_direction == 2,
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim,
+                                          (fe_degree == -1) ? 1 : fe_degree + 1,
+                                          0,
+                                          Number,
+                                          2>,
+        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                          dim,
+                                          (fe_degree == -1) ? 1 : fe_degree,
+                                          0,
+                                          Number,
+                                          2>>::type;
 
       Evalf0 evalf0 =
-        create_evaluator_tensor_product<Evalf0>(shape_info, face_no);
+        create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
+                                                  shape_info.data[0] :
+                                                  shape_info.data[1],
+                                                face_no);
       Evalf1 evalf1 =
-        create_evaluator_tensor_product<Evalf1>(shape_info, face_no);
+        create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
+                                                  shape_info.data[0] :
+                                                  shape_info.data[1],
+                                                face_no);
       Evalf2 evalf2 =
-        create_evaluator_tensor_product<Evalf2>(shape_info, face_no);
+        create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
+                                                  shape_info.data[0] :
+                                                  shape_info.data[1],
+                                                face_no);
 
       const unsigned int dofs_per_component_on_cell =
         shape_info.dofs_per_component_on_cell;
@@ -3787,8 +3844,7 @@ namespace internal
                                  do_evaluate,
                                  add_into_output,
                                  max_derivative,
-                                 lex_faces,
-                                 0>(input, output);
+                                 lex_faces>(input, output);
       // stride to next component
       input += (face_direction == 0) ? in_stride_after_normal : in_stride;
       output += (face_direction == 0) ? out_stride_after_normal : out_stride;
@@ -3797,8 +3853,7 @@ namespace internal
                                  do_evaluate,
                                  add_into_output,
                                  max_derivative,
-                                 lex_faces,
-                                 1>(input, output);
+                                 lex_faces>(input, output);
 
       if (dim == 3)
         {
@@ -3811,8 +3866,7 @@ namespace internal
                                      do_evaluate,
                                      add_into_output,
                                      max_derivative,
-                                     lex_faces,
-                                     2>(input, output);
+                                     lex_faces>(input, output);
         }
     }
   };
@@ -4209,12 +4263,13 @@ namespace internal
                                             (n_q_points_1d < 1) ? 1 :
                                                                   n_q_points_1d,
                                             Number>::
-            template evaluate_or_integrate_in_face<false>(evaluation_flag,
-                                                          temp,
-                                                          fe_eval,
-                                                          scratch_data,
-                                                          subface_index,
-                                                          fe_eval.get_face_no());
+            template evaluate_or_integrate_in_face<false>(
+              evaluation_flag,
+              temp,
+              fe_eval,
+              scratch_data,
+              subface_index,
+              fe_eval.get_face_no());
         }
       else if (fe_degree > -1 &&
                subface_index >= GeometryInfo<dim>::max_children_per_cell &&

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index 2e8885f01fa6c8b3a354ab1aa900dec5d889e938..3ff7d54a8156e91bf64c5517cc298bec28500cea 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -122,6 +122,38 @@ namespace internal
   struct EvaluatorTensorProduct
   {};
 
+  /**
+   * Evaluator framework for anisotropic polynomial spaces that valuates the
+   * given shape data in general dimensions using the tensor product form.
+   *
+   * @tparam variant Variant of evaluation used for creating template
+   *                 specializations
+   * @tparam dim Dimension of the function
+   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+   *                to the number of 1d shape functions in the usual tensor
+   *                contraction setting
+   * @tparam n_columns Number of columns in the transformation matrix, which
+   *                   corresponds to the number of 1d shape functions in the
+   *                   usual tensor contraction setting
+   * @tparam Number Abstract number type for input and output arrays
+   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+   *                 same type as the input/output arrays); must implement
+   *                 operator* with Number to be valid
+   * @tparam normal_dir Indicates the direction of the continuous component for the
+   *                    Raviart-Thomas space in terms of the normal onto the
+   * face, e.g 0 if the  is in x-direction, 1 if in y-direction, and 2 if in
+   * z-direction.
+   */
+  template <EvaluatorVariant variant,
+            int              dim,
+            int              n_rows,
+            int              n_columns,
+            typename Number,
+            int normal_dir,
+            typename Number2 = Number>
+  struct EvaluatorTensorProductAnisotropic
+  {};
+
 
 
   /**
@@ -198,30 +230,21 @@ namespace internal
       (void)dummy2;
     }
 
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              int  normal_dir = 0>
+    template <int direction, bool contract_over_rows, bool add>
     void
     values(const Number in[], Number out[]) const
     {
       apply<direction, contract_over_rows, add>(shape_values, in, out);
     }
 
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              int  normal_dir = 0>
+    template <int direction, bool contract_over_rows, bool add>
     void
     gradients(const Number in[], Number out[]) const
     {
       apply<direction, contract_over_rows, add>(shape_gradients, in, out);
     }
 
-    template <int  direction,
-              bool contract_over_rows,
-              bool add,
-              int  normal_dir = 0>
+    template <int direction, bool contract_over_rows, bool add>
     void
     hessians(const Number in[], Number out[]) const
     {
@@ -2364,13 +2387,15 @@ namespace internal
             int n_rows,
             int n_columns,
             typename Number,
+            int normal_dir,
             typename Number2>
-  struct EvaluatorTensorProduct<evaluate_raviart_thomas,
-                                dim,
-                                n_rows,
-                                n_columns,
-                                Number,
-                                Number2>
+  struct EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
+                                           dim,
+                                           n_rows,
+                                           n_columns,
+                                           Number,
+                                           normal_dir,
+                                           Number2>
   {
     static constexpr unsigned int n_rows_of_product =
       numbers::invalid_unsigned_int;
@@ -2381,7 +2406,7 @@ namespace internal
      * Empty constructor. Does nothing. Be careful when using 'values' and
      * related methods because they need to be filled with the other pointer
      */
-    EvaluatorTensorProduct()
+    EvaluatorTensorProductAnisotropic()
       : shape_values(nullptr)
       , shape_gradients(nullptr)
       , shape_hessians(nullptr)
@@ -2390,11 +2415,12 @@ namespace internal
     /**
      * Constructor, taking the data from ShapeInfo
      */
-    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
-                           const AlignedVector<Number2> &shape_gradients,
-                           const AlignedVector<Number2> &shape_hessians,
-                           const unsigned int            dummy1 = 0,
-                           const unsigned int            dummy2 = 0)
+    EvaluatorTensorProductAnisotropic(
+      const AlignedVector<Number2> &shape_values,
+      const AlignedVector<Number2> &shape_gradients,
+      const AlignedVector<Number2> &shape_hessians,
+      const unsigned int            dummy1 = 0,
+      const unsigned int            dummy2 = 0)
       : shape_values(shape_values.begin())
       , shape_gradients(shape_gradients.begin())
       , shape_hessians(shape_hessians.begin())
@@ -2417,31 +2443,25 @@ namespace internal
       (void)dummy2;
     }
 
-    template <int direction, bool contract_over_rows, bool add, int normal_dir>
+    template <int direction, bool contract_over_rows, bool add>
     void
     values(const Number in[], Number out[]) const
     {
-      apply<direction, contract_over_rows, add, normal_dir>(shape_values,
-                                                            in,
-                                                            out);
+      apply<direction, contract_over_rows, add>(shape_values, in, out);
     }
 
-    template <int direction, bool contract_over_rows, bool add, int normal_dir>
+    template <int direction, bool contract_over_rows, bool add>
     void
     gradients(const Number in[], Number out[]) const
     {
-      apply<direction, contract_over_rows, add, normal_dir>(shape_gradients,
-                                                            in,
-                                                            out);
+      apply<direction, contract_over_rows, add>(shape_gradients, in, out);
     }
 
-    template <int direction, bool contract_over_rows, bool add, int normal_dir>
+    template <int direction, bool contract_over_rows, bool add>
     void
     hessians(const Number in[], Number out[]) const
     {
-      apply<direction, contract_over_rows, add, normal_dir>(shape_hessians,
-                                                            in,
-                                                            out);
+      apply<direction, contract_over_rows, add>(shape_hessians, in, out);
     }
 
     /**
@@ -2475,7 +2495,6 @@ namespace internal
     template <int  direction,
               bool contract_over_rows,
               bool add,
-              int  normal_dir,
               bool one_line = false>
     static void
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
@@ -2486,8 +2505,7 @@ namespace internal
               bool contract_onto_face,
               bool add,
               int  max_derivative,
-              bool lex_faces = false,
-              int  normal_direction>
+              bool lex_faces = false>
     void
     apply_face(const Number *DEAL_II_RESTRICT in,
                Number *DEAL_II_RESTRICT       out) const;
@@ -2501,22 +2519,20 @@ namespace internal
             int n_rows,
             int n_columns,
             typename Number,
+            int normal_dir,
             typename Number2>
-  template <int  direction,
-            bool contract_over_rows,
-            bool add,
-            int  normal_dir,
-            bool one_line>
+  template <int direction, bool contract_over_rows, bool add, bool one_line>
   inline void
-  EvaluatorTensorProduct<evaluate_raviart_thomas,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::apply(const Number2 *DEAL_II_RESTRICT
-                                                       shape_data,
-                                         const Number *in,
-                                         Number *      out)
+  EvaluatorTensorProductAnisotropic<
+    evaluate_raviart_thomas,
+    dim,
+    n_rows,
+    n_columns,
+    Number,
+    normal_dir,
+    Number2>::apply(const Number2 *DEAL_II_RESTRICT shape_data,
+                    const Number *                  in,
+                    Number *                        out)
   {
     static_assert(one_line == false || direction == dim - 1,
                   "Single-line evaluation only works for direction=dim-1.");
@@ -2596,22 +2612,23 @@ namespace internal
             int n_rows,
             int n_columns,
             typename Number,
+            int normal_dir,
             typename Number2>
   template <int  face_direction,
             bool contract_onto_face,
             bool add,
             int  max_derivative,
-            bool lex_faces,
-            int  normal_direction>
+            bool lex_faces>
   inline void
-  EvaluatorTensorProduct<evaluate_raviart_thomas,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
-                                              Number *DEAL_II_RESTRICT
-                                                out) const
+  EvaluatorTensorProductAnisotropic<
+    evaluate_raviart_thomas,
+    dim,
+    n_rows,
+    n_columns,
+    Number,
+    normal_dir,
+    Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
+                         Number *DEAL_II_RESTRICT       out) const
   {
     Assert(dim > 1 && (lex_faces || dim < 4),
            ExcMessage("Only dim=2,3 supported"));
@@ -2625,34 +2642,33 @@ namespace internal
 
     // Determine the number of blocks depending on the face and normaldirection,
     // as well as dimension.
-    constexpr int n_blocks1 = (face_direction == normal_direction) ?
-                                (n_rows - 1) :
-                              ((face_direction == 0 && normal_direction == 2) ||
-                               (face_direction == 1 && normal_direction == 2) ||
-                               (face_direction == 2 && normal_direction == 1)) ?
-                                n_rows :
-                                (n_rows + 1);
-    constexpr int n_blocks2 =
-      (dim == 2) ? 1 :
-                   ((face_direction == normal_direction) ?
-                      (n_rows - 1) :
-                      (((face_direction == 0 && normal_direction == 1) ||
-                        (face_direction == 1 && normal_direction == 0) ||
-                        (face_direction == 2 && normal_direction == 0)) ?
-                         n_rows :
-                         (n_rows + 1)));
+    constexpr int n_blocks1 = (face_direction == normal_dir) ? (n_rows - 1) :
+                              ((face_direction == 0 && normal_dir == 2) ||
+                               (face_direction == 1 && normal_dir == 2) ||
+                               (face_direction == 2 && normal_dir == 1)) ?
+                                                               n_rows :
+                                                               (n_rows + 1);
+    constexpr int n_blocks2 = (dim == 2) ?
+                                1 :
+                                ((face_direction == normal_dir) ?
+                                   (n_rows - 1) :
+                                   (((face_direction == 0 && normal_dir == 1) ||
+                                     (face_direction == 1 && normal_dir == 0) ||
+                                     (face_direction == 2 && normal_dir == 0)) ?
+                                      n_rows :
+                                      (n_rows + 1)));
 
     AssertIndexRange(face_direction, dim);
 
     constexpr int stride =
-      (face_direction == normal_direction) ?
+      (face_direction == normal_dir) ?
         Utilities::pow(n_rows - 1, face_direction) :
         ((face_direction == 0) ?
            1 :
            ((face_direction == 2) ?
               n_rows * (n_rows + 1) :
-              ((face_direction == 1 && normal_direction == 0) ? (n_rows + 1) :
-                                                                n_rows)));
+              ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) :
+                                                          n_rows)));
     constexpr int out_stride = n_blocks1 * n_blocks2;
 
     const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
@@ -2739,21 +2755,21 @@ namespace internal
                   // product. Need to take that into account.
                   if (dim == 3)
                     {
-                      if (normal_direction == 0)
+                      if (normal_dir == 0)
                         {
                           if (contract_onto_face)
                             out += n_rows - 1;
                           else
                             in += n_rows - 1;
                         }
-                      if (normal_direction == 1)
+                      if (normal_dir == 1)
                         {
                           if (contract_onto_face)
                             out += n_rows - 2;
                           else
                             in += n_rows - 2;
                         }
-                      if (normal_direction == 2)
+                      if (normal_dir == 2)
                         {
                           if (contract_onto_face)
                             out += n_rows;
@@ -2777,17 +2793,17 @@ namespace internal
             // adjust for local coordinate system zx
             if (contract_onto_face)
               {
-                if (normal_direction == 0)
+                if (normal_dir == 0)
                   {
                     in += (n_rows + 1) * (n_rows - 1);
                     out -= n_rows * (n_rows + 1) - 1;
                   }
-                if (normal_direction == 1)
+                if (normal_dir == 1)
                   {
                     in += (n_rows - 1) * (n_rows - 1);
                     out -= (n_rows - 1) * (n_rows - 1) - 1;
                   }
-                if (normal_direction == 2)
+                if (normal_dir == 2)
                   {
                     in += (n_rows - 1) * (n_rows);
                     out -= (n_rows) * (n_rows + 1) - 1;
@@ -2795,17 +2811,17 @@ namespace internal
               }
             else
               {
-                if (normal_direction == 0)
+                if (normal_dir == 0)
                   {
                     out += (n_rows + 1) * (n_rows - 1);
                     in -= n_rows * (n_rows + 1) - 1;
                   }
-                if (normal_direction == 1)
+                if (normal_dir == 1)
                   {
                     out += (n_rows - 1) * (n_rows - 1);
                     in -= (n_rows - 1) * (n_rows - 1) - 1;
                   }
-                if (normal_direction == 2)
+                if (normal_dir == 2)
                   {
                     out += (n_rows - 1) * (n_rows);
                     in -= (n_rows) * (n_rows + 1) - 1;