From 6f05a5806e3e47b2bf55a66073a33bb457fb444c Mon Sep 17 00:00:00 2001 From: Dominik Still Date: Fri, 21 Jun 2024 11:57:39 +0200 Subject: [PATCH] Arrange face gradients contiguously --- .../matrix_free/evaluation_kernels_face.h | 119 +++++++++--------- include/deal.II/matrix_free/shape_info.h | 2 +- .../matrix_free/shape_info.templates.h | 14 ++- 3 files changed, 69 insertions(+), 66 deletions(-) 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; - 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( + 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 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( + 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; 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( + 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( + 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 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( + shape_gradients, in, out, n_dofs, n_q_points * dim, 1, 1); + else + apply_matrix_vector_product( + 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]; } } } -- 2.39.5