From 7f55397a2b8b6254c275f492967a4ce42baca652 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 6 May 2024 13:29:41 +0200 Subject: [PATCH] FEEvaluation for simplices: Arrange gradients contiguously --- .../deal.II/matrix_free/evaluation_kernels.h | 116 +++++++++--------- .../matrix_free/shape_info.templates.h | 3 +- tests/simplex/matrix_free_shape_info_01.cc | 2 +- 3 files changed, 60 insertions(+), 61 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 1a0e39dbdf..c3b08a98ae 100644 --- 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::shape_info_number_type; - using Eval = - EvaluatorTensorProduct; 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( + 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( + 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::shape_info_number_type; - using Eval = - EvaluatorTensorProduct; 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( + 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( + 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( + 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.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 2eafca6dbf..c4b6b68d18 100644 --- 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 602c66c2af..58677743ae 100644 --- 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]; } } -- 2.39.5