Merge pull request #15951 from kronbichler/restructure_gradients

Matrix-free evaluation of Raviart-Thomas: major restructuring
dealii · Sep 12, 2023 · 45ca30f · 45ca30f
2 parents a19c8bc + 8898ee1
commit 45ca30f
Show file tree

Hide file tree

Showing 8 changed files with 2,106 additions and 2,916 deletions.
diff --git a/doc/news/changes/incompatibilities/20230911Kronbichler b/doc/news/changes/incompatibilities/20230911Kronbichler
@@ -0,0 +1,11 @@
+Changed: The class FEEvaluation now uses a different internal data layout for
+the gradients, exposed via FEEvaluation::begin_gradients(). Now, the entries
+of the partial derivatives in the space directions are placed adjacent in
+memory. The entries of different components are still separated by the entries
+on all points. This change has been made to simplify the access in the
+FEEvaluation::get_gradient() and FEEvaluation::submit_gradient() functions,
+which is especially useful for the case with many FE components. For the
+regular use of FEEvaluation apart from the plain pointers mentioned above,
+there is no change in behavior.
+<br>
+(Martin Kronbichler, 2023/09/11)
diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h
@@ -1165,17 +1165,12 @@ namespace internal
           // We need to reinterpret the data after evaluate has been applied.
           for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
             for (unsigned int point = 0; point < n_q_points; ++point)
-              for (unsigned int j = 0; j < dim; ++j)
-                for (unsigned int in_comp = 0;
-                     in_comp < n_lanes &&
-                     in_comp < spacedim - out_comp * n_lanes;
-                     ++in_comp)
+              for (unsigned int in_comp = 0;
+                   in_comp < n_lanes && in_comp < spacedim - out_comp * n_lanes;
+                   ++in_comp)
+                for (unsigned int j = 0; j < dim; ++j)
                   {
-                    const unsigned int total_number = point * dim + j;
-                    const unsigned int new_comp     = total_number / n_q_points;
-                    const unsigned int new_point    = total_number % n_q_points;
-                    jacobians[new_point][out_comp * n_lanes +
-                                         in_comp][new_comp] =
+                    jacobians[point][out_comp * n_lanes + in_comp][j] =
                       eval.begin_gradients()[(out_comp * n_q_points + point) *
                                                dim +
                                              j][in_comp];

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -1230,7 +1230,7 @@ namespace internal
                       for (unsigned int e = 0; e < dim; ++e)
                         jac[d][e] =
                           eval
-                            .begin_gradients()[q + (d * dim + e) * n_q_points];
+                            .begin_gradients()[e + (d * n_q_points + q) * dim];
 
                     // eliminate roundoff errors
                     if (cell_type[cell] == cartesian)
@@ -2180,7 +2180,7 @@ namespace internal
                   for (unsigned int e = 0; e < dim; ++e)
                     for (unsigned int d = 0; d < dim; ++d)
                       jacobi[d][e] =
-                        eval.begin_gradients()[(d * dim + e) * n_q_points + q];
+                        eval.begin_gradients()[(d * n_q_points + q) * dim + e];
                   Tensor<2, dim, VectorizedDouble> inv_transp_jac =
                     transpose(invert(jacobi));
                   Tensor<3, dim, VectorizedDouble> jac_grad;
@@ -2250,7 +2250,7 @@ namespace internal
                       for (unsigned int d = 0; d < dim; ++d)
                         jac[d][ee] =
                           eval_int
-                            .begin_gradients()[(d * dim + e) * n_q_points + q];
+                            .begin_gradients()[(d * n_q_points + q) * dim + e];
                     }
                   Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac);
                   for (unsigned int e = 0; e < dim; ++e)
@@ -2365,8 +2365,8 @@ namespace internal
                           for (unsigned int d = 0; d < dim; ++d)
                             jac[d][ee] =
                               eval_ext
-                                .begin_gradients()[(d * dim + e) * n_q_points +
-                                                   q];
+                                .begin_gradients()[(d * n_q_points + q) * dim +
+                                                   e];
                         }
                       Tensor<2, dim, VectorizedDouble> inv_jac = invert(jac);
                       for (unsigned int e = 0; e < dim; ++e)