FEEvaluation for simplices: Arrange gradients contiguously

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 1a0e39dbdfd9ebe1d323b9d5f5925b6dbe744256..c3b08a98ae86e147a0dca740978696440b04f547 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -630,23 +630,24 @@ namespace internal
 
     using Number2 =
       typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-    using Eval =
-      EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
 
     if (evaluation_flag & EvaluationFlags::values)
       {
         const auto *const shape_values = shape_data.front().shape_values.data();
-        auto             *values_quad_ptr        = fe_eval.begin_values();
-        const auto       *values_dofs_actual_ptr = values_dofs_actual;
+        auto             *out          = fe_eval.begin_values();
+        const auto       *in           = values_dofs_actual;
 
-        Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
         for (unsigned int c = 0; c < n_components; ++c)
           {
-            eval.template values<0, true, false>(values_dofs_actual_ptr,
-                                                 values_quad_ptr);
-
-            values_quad_ptr += n_q_points;
-            values_dofs_actual_ptr += n_dofs;
+            apply_matrix_vector_product<evaluate_general,
+                                        EvaluatorQuantity::value,
+                                        /* transpose_matrix */ true,
+                                        /* add */ false,
+                                        /* consider_strides */ false>(
+              shape_values, in, out, n_dofs, n_q_points, 1, 1);
+
+            out += n_q_points;
+            in += n_dofs;
           }
       }
 
@@ -654,24 +655,19 @@ namespace internal
       {
         const auto *const shape_gradients =
           shape_data.front().shape_gradients.data();
-        auto       *gradients_quad_ptr     = fe_eval.begin_gradients();
-        const auto *values_dofs_actual_ptr = values_dofs_actual;
+        auto       *out = fe_eval.begin_gradients();
+        const auto *in  = values_dofs_actual;
 
         for (unsigned int c = 0; c < n_components; ++c)
           {
-            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);
-
-                eval.template gradients<0, true, false, dim>(
-                  values_dofs_actual_ptr, gradients_quad_ptr + d);
-              }
-            gradients_quad_ptr += n_q_points * dim;
-            values_dofs_actual_ptr += n_dofs;
+            apply_matrix_vector_product<evaluate_general,
+                                        EvaluatorQuantity::value,
+                                        /* transpose_matrix */ true,
+                                        /* add */ false,
+                                        /* consider_strides */ false>(
+              shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+            out += n_q_points * dim;
+            in += n_dofs;
           }
       }
   }
@@ -702,27 +698,32 @@ namespace internal
 
     using Number2 =
       typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
-    using Eval =
-      EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
 
     if (integration_flag & EvaluationFlags::values)
       {
         const auto *const shape_values = shape_data.front().shape_values.data();
-        auto             *values_quad_ptr        = fe_eval.begin_values();
-        auto             *values_dofs_actual_ptr = values_dofs_actual;
+        auto             *in           = fe_eval.begin_values();
+        auto             *out          = values_dofs_actual;
 
-        Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
         for (unsigned int c = 0; c < n_components; ++c)
           {
             if (add_into_values_array == false)
-              eval.template values<0, false, false>(values_quad_ptr,
-                                                    values_dofs_actual_ptr);
+              apply_matrix_vector_product<evaluate_general,
+                                          EvaluatorQuantity::value,
+                                          /* transpose_matrix */ false,
+                                          /* add */ false,
+                                          /* consider_strides */ false>(
+                shape_values, in, out, n_dofs, n_q_points, 1, 1);
             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;
+              apply_matrix_vector_product<evaluate_general,
+                                          EvaluatorQuantity::value,
+                                          /* transpose_matrix */ false,
+                                          /* add */ true,
+                                          /* consider_strides */ false>(
+                shape_values, in, out, n_dofs, n_q_points, 1, 1);
+
+            in += n_q_points;
+            out += n_dofs;
           }
       }
 
@@ -730,30 +731,29 @@ namespace internal
       {
         const auto *const shape_gradients =
           shape_data.front().shape_gradients.data();
-        auto *gradients_quad_ptr     = fe_eval.begin_gradients();
-        auto *values_dofs_actual_ptr = values_dofs_actual;
+        auto *in  = fe_eval.begin_gradients();
+        auto *out = values_dofs_actual;
 
         for (unsigned int c = 0; c < n_components; ++c)
           {
-            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)) &&
-                    d == 0)
-                  eval.template gradients<0, false, false, dim>(
-                    gradients_quad_ptr + d, values_dofs_actual_ptr);
-                else
-                  eval.template gradients<0, false, true, dim>(
-                    gradients_quad_ptr + d, values_dofs_actual_ptr);
-              }
-            gradients_quad_ptr += n_q_points * dim;
-            values_dofs_actual_ptr += n_dofs;
+            if (add_into_values_array == false &&
+                !(integration_flag & EvaluationFlags::values))
+              apply_matrix_vector_product<evaluate_general,
+                                          EvaluatorQuantity::value,
+                                          /* transpose_matrix */ false,
+                                          /* add */ false,
+                                          /* consider_strides */ false>(
+                shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+            else
+              apply_matrix_vector_product<evaluate_general,
+                                          EvaluatorQuantity::value,
+                                          /* transpose_matrix */ false,
+                                          /* add */ true,
+                                          /* consider_strides */ false>(
+                shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1);
+
+            in += n_q_points * dim;
+            out += n_dofs;
           }
       }
   }

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 2eafca6dbf1344d8a414bfca3229eaa3dd86ac46..c4b6b68d18c0de9f47e9a5627f263253adff6d76 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -375,8 +375,7 @@ namespace internal
                 const auto grad = fe.shape_grad(i, quad.point(q));
 
                 for (unsigned int d = 0; d < dim; ++d)
-                  shape_gradients[d * n_dofs * n_q_points + i * n_q_points +
-                                  q] = grad[d];
+                  shape_gradients[i * dim * n_q_points + q * dim + d] = grad[d];
               }
 
           {

diff --git a/tests/simplex/matrix_free_shape_info_01.cc b/tests/simplex/matrix_free_shape_info_01.cc
index 602c66c2afe2c4408a1a13c4e12cbcc2fb8f08fb..58677743ae1e379c68f7ab6f9301e891e97b5ade 100644 (file)
--- a/tests/simplex/matrix_free_shape_info_01.cc
+++ b/tests/simplex/matrix_free_shape_info_01.cc
@@ -83,7 +83,7 @@ public:
           for (unsigned int i = 0; i < n_dofs; ++i)
             gradients[n_q_points * d + q] +=
               shape_info.data[0]
-                .shape_gradients[d * n_dofs * n_q_points + i * n_q_points + q] *
+                .shape_gradients[i * dim * n_q_points + q * dim + d] *
               dof_values[i];
         }
   }