From: Dominik Still <dominik.still@tum.de>
Date: Fri, 21 Jun 2024 09:57:39 +0000 (+0200)
Subject: Arrange face gradients contiguously
X-Git-Tag: v9.6.0-rc1~162^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6f05a5806e3e47b2bf55a66073a33bb457fb444c;p=dealii.git

Arrange face gradients contiguously
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h
index 9ce0016485..0356dd1c6b 100644
--- a/include/deal.II/matrix_free/evaluation_kernels_face.h
+++ b/include/deal.II/matrix_free/evaluation_kernels_face.h
@@ -1467,50 +1467,46 @@ namespace internal
       const std::size_t  n_dofs     = shape_info.dofs_per_component_on_cell;
       const std::size_t  n_q_points = shape_info.n_q_points_faces[face_no];
 
-      using Eval =
-        EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
-
       if (evaluation_flag & EvaluationFlags::values)
         {
           const auto *const shape_values =
             &shape_data.shape_values_face(face_no, face_orientation, 0);
 
-          auto *values_quad_ptr        = fe_eval.begin_values();
-          auto *values_dofs_actual_ptr = values_dofs;
+          auto *out = fe_eval.begin_values();
+          auto *in  = values_dofs;
 
-          Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
           for (unsigned int c = 0; c < n_components; ++c)
             {
-              eval.template values<0, true, false>(values_dofs_actual_ptr,
-                                                   values_quad_ptr);
-
-              values_quad_ptr += n_q_points;
-              values_dofs_actual_ptr += n_dofs;
+              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;
             }
         }
 
       if (evaluation_flag & EvaluationFlags::gradients)
         {
-          auto       *gradients_quad_ptr     = fe_eval.begin_gradients();
-          const auto *values_dofs_actual_ptr = values_dofs;
+          auto       *out = fe_eval.begin_gradients();
+          const auto *in  = values_dofs;
 
-          std::array<const Number2 *, dim> shape_gradients;
-          for (unsigned int d = 0; d < dim; ++d)
-            shape_gradients[d] =
-              &shape_data.shape_gradients_face(face_no, face_orientation, d, 0);
+          const auto *const shape_gradients =
+            &shape_data.shape_gradients_face(face_no, face_orientation, 0);
 
           for (unsigned int c = 0; c < n_components; ++c)
             {
-              for (unsigned int d = 0; d < dim; ++d)
-                {
-                  Eval eval(
-                    nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points);
-
-                  eval.template gradients<0, true, false, dim>(
-                    values_dofs_actual_ptr, gradients_quad_ptr + d);
-                }
-              gradients_quad_ptr += n_q_points * dim;
-              values_dofs_actual_ptr += n_dofs;
+              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;
             }
         }
 
@@ -1879,58 +1875,63 @@ namespace internal
       const std::size_t  n_dofs     = shape_info.dofs_per_component_on_cell;
       const std::size_t  n_q_points = shape_info.n_q_points_faces[face_no];
 
-      using Eval =
-        EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
 
       if (integration_flag & EvaluationFlags::values)
         {
           const auto *const shape_values =
             &shape_data.shape_values_face(face_no, face_orientation, 0);
 
-          auto *values_quad_ptr        = fe_eval.begin_values();
-          auto *values_dofs_actual_ptr = values_dofs;
+          auto *in  = fe_eval.begin_values();
+          auto *out = values_dofs;
 
-          Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
           for (unsigned int c = 0; c < n_components; ++c)
             {
               if (sum_into_values)
-                eval.template values<0, false, true>(values_quad_ptr,
-                                                     values_dofs_actual_ptr);
+                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);
               else
-                eval.template values<0, false, false>(values_quad_ptr,
-                                                      values_dofs_actual_ptr);
-              values_quad_ptr += n_q_points;
-              values_dofs_actual_ptr += n_dofs;
+                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);
+              in += n_q_points;
+              out += n_dofs;
             }
         }
 
       if (integration_flag & EvaluationFlags::gradients)
         {
-          auto *gradients_quad_ptr     = fe_eval.begin_gradients();
-          auto *values_dofs_actual_ptr = values_dofs;
+          auto *in  = fe_eval.begin_gradients();
+          auto *out = values_dofs;
 
-          std::array<const Number2 *, dim> shape_gradients;
-          for (unsigned int d = 0; d < dim; ++d)
-            shape_gradients[d] =
-              &shape_data.shape_gradients_face(face_no, face_orientation, d, 0);
+          const auto *const shape_gradients =
+            &shape_data.shape_gradients_face(face_no, face_orientation, 0);
 
           for (unsigned int c = 0; c < n_components; ++c)
             {
-              for (unsigned int d = 0; d < dim; ++d)
-                {
-                  Eval eval(
-                    nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points);
-
-                  if (!sum_into_values &&
-                      !(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 (!sum_into_values &&
+                  !(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.h b/include/deal.II/matrix_free/shape_info.h
index b3d5252464..779e1753a6 100644
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -374,7 +374,7 @@ namespace internal
        * quadrature points for all faces, orientations, and directions
        * (no tensor-product structure exploited).
        */
-      Table<4, Number> shape_gradients_face;
+      Table<3, Number> shape_gradients_face;
     };
 
 
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index c4b6b68d18..9858856a05 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -425,10 +425,10 @@ namespace internal
                                           n_max_face_orientations,
                                           n_dofs * n_q_points_face_max});
 
-                shape_gradients_face.reinit({n_faces,
-                                             n_max_face_orientations,
-                                             dim,
-                                             n_dofs * n_q_points_face_max});
+                shape_gradients_face.reinit(
+                  {n_faces,
+                   n_max_face_orientations,
+                   dim * n_dofs * n_q_points_face_max});
 
                 for (unsigned int f = 0; f < n_faces; ++f)
                   {
@@ -464,8 +464,10 @@ namespace internal
                               const auto grad = fe.shape_grad(i, point);
 
                               for (unsigned int d = 0; d < dim; ++d)
-                                shape_gradients_face(
-                                  f, o, d, i * n_q_points_face + q) = grad[d];
+                                shape_gradients_face(f,
+                                                     o,
+                                                     i * dim * n_q_points_face +
+                                                       q * dim + d) = grad[d];
                             }
                       }
                   }