Use EvaluatorTensorProduct in Simplex code path

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 036b628f597b4ab9e332041a0c9b42d2ee9d31e6..bfa6319cdab0af872f9c96c3d47e2d5fa482bf38 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -696,35 +696,50 @@ namespace internal
     const unsigned int n_dofs     = shape_info.dofs_per_component_on_cell;
     const unsigned int n_q_points = shape_info.n_q_points;
 
+    using Eval =
+      EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
+
     if (evaluation_flag & EvaluationFlags::values)
       {
         const auto shape_values = shape_info.data.front().shape_values.data();
+        auto       values_quad_ptr        = values_quad;
+        auto       values_dofs_actual_ptr = values_dofs_actual;
+
+        Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
         for (unsigned int c = 0; c < n_components; ++c)
-          for (unsigned int q = 0; q < n_q_points; ++q)
-            {
-              values_quad[n_q_points * c + q] = 0.0;
-              for (unsigned int i = 0; i < n_dofs; ++i)
-                values_quad[n_q_points * c + q] +=
-                  shape_values[i * n_q_points + q] *
-                  values_dofs_actual[n_dofs * c + i];
-            }
+          {
+            eval.template values<0, true, false>(values_dofs_actual_ptr,
+                                                 values_quad_ptr);
+
+            values_quad_ptr += n_q_points;
+            values_dofs_actual_ptr += n_dofs;
+          }
       }
 
     if (evaluation_flag & EvaluationFlags::gradients)
       {
         const auto shape_gradients =
           shape_info.data.front().shape_gradients.data();
+        auto gradients_quad_ptr     = gradients_quad;
+        auto values_dofs_actual_ptr = values_dofs_actual;
+
         for (unsigned int c = 0; c < n_components; ++c)
-          for (unsigned int d = 0; d < dim; ++d)
-            for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            for (unsigned int d = 0; d < dim; ++d)
               {
-                gradients_quad[n_q_points * dim * c + n_q_points * d + q] = 0.0;
-                for (unsigned int i = 0; i < n_dofs; ++i)
-                  gradients_quad[n_q_points * dim * c + n_q_points * d + q] +=
-                    shape_gradients[d * n_dofs * n_q_points + i * n_q_points +
-                                    q] *
-                    values_dofs_actual[n_dofs * dim * c + i];
+                Eval eval(nullptr,
+                          shape_gradients + n_q_points * n_dofs * d,
+                          nullptr,
+                          n_dofs,
+                          n_q_points);
+
+                eval.template gradients<0, true, false>(values_dofs_actual_ptr,
+                                                        gradients_quad_ptr);
+
+                gradients_quad_ptr += n_q_points;
               }
+            values_dofs_actual_ptr += n_dofs;
+          }
       }
 
     if (evaluation_flag & EvaluationFlags::hessians)
@@ -757,33 +772,60 @@ namespace internal
     const unsigned int n_dofs     = shape_info.dofs_per_component_on_cell;
     const unsigned int n_q_points = shape_info.n_q_points;
 
-    if (add_into_values_array == false)
-      for (unsigned int i = 0; i < n_components * n_dofs; ++i)
-        values_dofs_actual[i] = 0.0;
+    using Eval =
+      EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number>;
 
     if (integration_flag & EvaluationFlags::values)
       {
         const auto shape_values = shape_info.data.front().shape_values.data();
+        auto       values_quad_ptr        = values_quad;
+        auto       values_dofs_actual_ptr = values_dofs_actual;
+
+        Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
         for (unsigned int c = 0; c < n_components; ++c)
-          for (unsigned int q = 0; q < n_q_points; ++q)
-            for (unsigned int i = 0; i < n_dofs; ++i)
-              values_dofs_actual[n_dofs * c + i] +=
-                shape_values[i * n_q_points + q] *
-                values_quad[n_q_points * c + q];
+          {
+            if (add_into_values_array == false)
+              eval.template values<0, false, false>(values_quad_ptr,
+                                                    values_dofs_actual_ptr);
+            else
+              eval.template values<0, false, true>(values_quad_ptr,
+                                                   values_dofs_actual_ptr);
+
+            values_quad_ptr += n_q_points;
+            values_dofs_actual_ptr += n_dofs;
+          }
       }
 
     if (integration_flag & EvaluationFlags::gradients)
       {
         const auto shape_gradients =
           shape_info.data.front().shape_gradients.data();
+        auto gradients_quad_ptr     = gradients_quad;
+        auto values_dofs_actual_ptr = values_dofs_actual;
+
         for (unsigned int c = 0; c < n_components; ++c)
-          for (unsigned int d = 0; d < dim; ++d)
-            for (unsigned int q = 0; q < n_q_points; ++q)
-              for (unsigned int i = 0; i < n_dofs; ++i)
-                values_dofs_actual[n_dofs * dim * c + i] +=
-                  shape_gradients[d * n_dofs * n_q_points + i * n_q_points +
-                                  q] *
-                  gradients_quad[n_q_points * dim * c + n_q_points * d + q];
+          {
+            for (unsigned int d = 0; d < dim; ++d)
+              {
+                Eval eval(nullptr,
+                          shape_gradients + n_q_points * n_dofs * d,
+                          nullptr,
+                          n_dofs,
+                          n_q_points);
+
+                if ((add_into_values_array == false &&
+                     (integration_flag & EvaluationFlags::values) == false) &&
+                    d == 0)
+                  eval.template gradients<0, false, false>(
+                    gradients_quad_ptr, values_dofs_actual_ptr);
+                else
+                  eval.template gradients<0, false, true>(
+                    gradients_quad_ptr, values_dofs_actual_ptr);
+
+                gradients_quad_ptr += n_q_points;
+              }
+            values_dofs_actual_ptr += n_dofs;
+          }
       }
   }
 

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index 166d7285b7cfb85ece6c2462d4d95991f97df183..21ed078d103038fa80ffa86d386eb3f9d1305c57 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -625,6 +625,21 @@ namespace internal
              ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
     }
 
+    /**
+     * Constructor, taking the data from ShapeInfo
+     */
+    EvaluatorTensorProduct(const Number2 *    shape_values,
+                           const Number2 *    shape_gradients,
+                           const Number2 *    shape_hessians,
+                           const unsigned int n_rows,
+                           const unsigned int n_columns)
+      : shape_values(shape_values)
+      , shape_gradients(shape_gradients)
+      , shape_hessians(shape_hessians)
+      , n_rows(n_rows)
+      , n_columns(n_columns)
+    {}
+
     template <int direction, bool contract_over_rows, bool add>
     void
     values(const Number *in, Number *out) const