From: Martin Kronbichler Date: Wed, 30 Aug 2023 08:21:51 +0000 (+0200) Subject: FEEvaluation: Restructure layout of gradients. X-Git-Tag: relicensing~510^2~8 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5ab7a0f34a082faf6ed99347296d3915ee157b59;p=dealii.git FEEvaluation: Restructure layout of gradients. Generalize and optimize Raviart-Thomas evaluation significantly --- diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index edb6e12309..80250bbdd7 100644 --- 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 index 95583af217..1401a04180 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -188,55 +188,6 @@ namespace internal const bool add_into_values_array); }; - /** - * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use - * specific sum-factorization kernels and with normal/tangential shape_data - */ - template - struct FEEvaluationImpl - { - template - static void - evaluate_or_integrate( - const EvaluationFlags::EvaluationFlags evaluation_flag, - Number *values_dofs_actual, - FEEvaluationData &fe_eval, - const bool add_into_values_array = false); - - private: - template - static EvalType - create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &shape_data) - { - return EvalType(shape_data.shape_values, - shape_data.shape_gradients, - shape_data.shape_hessians); - } - - template - static void - evaluate_tensor_product_per_component( - const EvaluationFlags::EvaluationFlags evaluation_flag, - Number *values_dofs_actual, - FEEvaluationData &fe_eval, - const bool add_into_values_array, - std::integral_constant); - - template - static void - evaluate_tensor_product_per_component( - const EvaluationFlags::EvaluationFlags evaluation_flag, - Number *values_dofs_actual, - FEEvaluationData &fe_eval, - const bool add_into_values_array, - std::integral_constant); - }; - template (values_dofs, gradients_quad); - if (evaluation_flag & EvaluationFlags::hessians) - eval0.template hessians<0, true, false>(values_dofs, - hessians_quad); // advance the next component in 1d array values_dofs += dofs_per_comp; values_quad += n_q_points; gradients_quad += n_q_points; - hessians_quad += n_q_points; } break; @@ -360,35 +306,15 @@ namespace internal if (evaluation_flag & EvaluationFlags::gradients) { eval0.template gradients<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, gradients_quad); - } - if (evaluation_flag & EvaluationFlags::hessians) - { - // grad xy - if (!(evaluation_flag & EvaluationFlags::gradients)) - eval0.template gradients<0, true, false>(values_dofs, - temp1); - eval1.template gradients<1, true, false>(temp1, - hessians_quad + - 2 * n_q_points); - - // grad xx - eval0.template hessians<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, hessians_quad); + eval1.template values<1, true, false, 2>(temp1, + gradients_quad); } // grad y eval0.template values<0, true, false>(values_dofs, temp1); if (evaluation_flag & EvaluationFlags::gradients) - eval1.template gradients<1, true, false>(temp1, - gradients_quad + - n_q_points); - - // grad yy - if (evaluation_flag & EvaluationFlags::hessians) - eval1.template hessians<1, true, false>(temp1, - hessians_quad + - n_q_points); + eval1.template gradients<1, true, false, 2>(temp1, + gradients_quad + 1); // val: can use values applied in x if (evaluation_flag & EvaluationFlags::values) @@ -398,7 +324,6 @@ namespace internal values_dofs += dofs_per_comp; values_quad += n_q_points; gradients_quad += 2 * n_q_points; - hessians_quad += 3 * n_q_points; } break; @@ -410,32 +335,8 @@ namespace internal // grad x eval0.template gradients<0, true, false>(values_dofs, temp1); eval1.template values<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, gradients_quad); - } - - if (evaluation_flag & EvaluationFlags::hessians) - { - // grad xz - if (!(evaluation_flag & EvaluationFlags::gradients)) - { - eval0.template gradients<0, true, false>(values_dofs, - temp1); - eval1.template values<1, true, false>(temp1, temp2); - } - eval2.template gradients<2, true, false>(temp2, - hessians_quad + - 4 * n_q_points); - - // grad xy - eval1.template gradients<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, - hessians_quad + - 3 * n_q_points); - - // grad xx - eval0.template hessians<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, hessians_quad); + eval2.template values<2, true, false, 3>(temp2, + gradients_quad); } // grad y @@ -443,41 +344,16 @@ namespace internal if (evaluation_flag & EvaluationFlags::gradients) { eval1.template gradients<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, - gradients_quad + - n_q_points); - } - - if (evaluation_flag & EvaluationFlags::hessians) - { - // grad yz - if (!(evaluation_flag & EvaluationFlags::gradients)) - eval1.template gradients<1, true, false>(temp1, temp2); - eval2.template gradients<2, true, false>(temp2, - hessians_quad + - 5 * n_q_points); - - // grad yy - eval1.template hessians<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, - hessians_quad + - n_q_points); + eval2.template values<2, true, false, 3>(temp2, + gradients_quad + 1); } // grad z: can use the values applied in x direction stored in // temp1 eval1.template values<1, true, false>(temp1, temp2); if (evaluation_flag & EvaluationFlags::gradients) - eval2.template gradients<2, true, false>(temp2, - gradients_quad + - 2 * n_q_points); - - // grad zz: can use the values applied in x and y direction stored - // in temp2 - if (evaluation_flag & EvaluationFlags::hessians) - eval2.template hessians<2, true, false>(temp2, - hessians_quad + - 2 * n_q_points); + eval2.template gradients<2, true, false, 3>(temp2, + gradients_quad + 2); // val: can use the values applied in x & y direction stored in // temp2 @@ -488,7 +364,6 @@ namespace internal values_dofs += dofs_per_comp; values_quad += n_q_points; gradients_quad += 3 * n_q_points; - hessians_quad += 6 * n_q_points; } break; @@ -576,7 +451,6 @@ namespace internal Number *values_quad = fe_eval.begin_values(); Number *gradients_quad = fe_eval.begin_gradients(); - Number *hessians_quad = fe_eval.begin_hessians(); switch (dim) { @@ -602,23 +476,11 @@ namespace internal eval0.template gradients<0, false, false>(gradients_quad, values_dofs); } - if ((integration_flag & EvaluationFlags::hessians) != 0u) - { - if ((integration_flag & EvaluationFlags::values) != 0u || - (integration_flag & EvaluationFlags::gradients) != 0u || - add_into_values_array == true) - eval0.template hessians<0, false, true>(hessians_quad, - values_dofs); - else - eval0.template hessians<0, false, false>(hessians_quad, - values_dofs); - } // advance to the next component in 1d array values_dofs += dofs_per_comp; values_quad += n_q_points; gradients_quad += n_q_points; - hessians_quad += n_q_points; } break; @@ -636,40 +498,16 @@ namespace internal } if (integration_flag & EvaluationFlags::gradients) { - eval1.template gradients<1, false, false>(gradients_quad + - n_q_points, - temp1); + eval1.template gradients<1, false, false, 2>(gradients_quad + + 1, + temp1); if (integration_flag & EvaluationFlags::values) eval1.template values<1, false, true>(values_quad, temp1); if (add_into_values_array == false) eval0.template values<0, false, false>(temp1, values_dofs); else eval0.template values<0, false, true>(temp1, values_dofs); - eval1.template values<1, false, false>(gradients_quad, temp1); - eval0.template gradients<0, false, true>(temp1, values_dofs); - } - if ((integration_flag & EvaluationFlags::hessians) != 0u) - { - // grad xx - eval1.template values<1, false, false>(hessians_quad, temp1); - - if ((integration_flag & EvaluationFlags::values) != 0u || - (integration_flag & EvaluationFlags::gradients) != 0u || - add_into_values_array == true) - eval0.template hessians<0, false, true>(temp1, values_dofs); - else - eval0.template hessians<0, false, false>(temp1, - values_dofs); - - // grad yy - eval1.template hessians<1, false, false>(hessians_quad + - n_q_points, - temp1); - eval0.template values<0, false, true>(temp1, values_dofs); - - // grad xy - eval1.template gradients<1, false, false>(hessians_quad + - 2 * n_q_points, + eval1.template values<1, false, false, 2>(gradients_quad, temp1); eval0.template gradients<0, false, true>(temp1, values_dofs); } @@ -678,7 +516,6 @@ namespace internal values_dofs += dofs_per_comp; values_quad += n_q_points; gradients_quad += 2 * n_q_points; - hessians_quad += 3 * n_q_points; } break; @@ -697,79 +534,29 @@ namespace internal } if (integration_flag & EvaluationFlags::gradients) { - eval2.template gradients<2, false, false>(gradients_quad + - 2 * n_q_points, - temp1); + eval2.template gradients<2, false, false, 3>(gradients_quad + + 2, + temp1); if (integration_flag & EvaluationFlags::values) eval2.template values<2, false, true>(values_quad, temp1); eval1.template values<1, false, false>(temp1, temp2); - eval2.template values<2, false, false>(gradients_quad + - n_q_points, - temp1); + eval2.template values<2, false, false, 3>(gradients_quad + 1, + temp1); eval1.template gradients<1, false, true>(temp1, temp2); if (add_into_values_array == false) eval0.template values<0, false, false>(temp2, values_dofs); else eval0.template values<0, false, true>(temp2, values_dofs); - eval2.template values<2, false, false>(gradients_quad, temp1); - eval1.template values<1, false, false>(temp1, temp2); - eval0.template gradients<0, false, true>(temp2, values_dofs); - } - if ((integration_flag & EvaluationFlags::hessians) != 0u) - { - // grad xx - eval2.template values<2, false, false>(hessians_quad, temp1); - eval1.template values<1, false, false>(temp1, temp2); - - if ((integration_flag & EvaluationFlags::values) != 0u || - (integration_flag & EvaluationFlags::gradients) != 0u || - add_into_values_array == true) - eval0.template hessians<0, false, true>(temp2, values_dofs); - else - eval0.template hessians<0, false, false>(temp2, - values_dofs); - - // grad yy - eval2.template values<2, false, false>(hessians_quad + - n_q_points, - temp1); - eval1.template hessians<1, false, false>(temp1, temp2); - eval0.template values<0, false, true>(temp2, values_dofs); - - // grad zz - eval2.template hessians<2, false, false>(hessians_quad + - 2 * n_q_points, - temp1); - eval1.template values<1, false, false>(temp1, temp2); - eval0.template values<0, false, true>(temp2, values_dofs); - - // grad xy - eval2.template values<2, false, false>(hessians_quad + - 3 * n_q_points, - temp1); - eval1.template gradients<1, false, false>(temp1, temp2); - eval0.template gradients<0, false, true>(temp2, values_dofs); - - // grad xz - eval2.template gradients<2, false, false>(hessians_quad + - 4 * n_q_points, + eval2.template values<2, false, false, 3>(gradients_quad, temp1); eval1.template values<1, false, false>(temp1, temp2); eval0.template gradients<0, false, true>(temp2, values_dofs); - - // grad yz - eval2.template gradients<2, false, false>(hessians_quad + - 5 * n_q_points, - temp1); - eval1.template gradients<1, false, false>(temp1, temp2); - eval0.template values<0, false, true>(temp2, values_dofs); } // advance to the next component in 1d array values_dofs += dofs_per_comp; values_quad += n_q_points; gradients_quad += 3 * n_q_points; - hessians_quad += 6 * n_q_points; } break; @@ -885,11 +672,10 @@ namespace internal n_dofs, n_q_points); - eval.template gradients<0, true, false>(values_dofs_actual_ptr, - gradients_quad_ptr); - - gradients_quad_ptr += 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; } } @@ -965,465 +751,19 @@ namespace internal if ((add_into_values_array == false && !(integration_flag & EvaluationFlags::values)) && d == 0) - eval.template gradients<0, false, false>( - gradients_quad_ptr, values_dofs_actual_ptr); + eval.template gradients<0, false, false, dim>( + gradients_quad_ptr + d, values_dofs_actual_ptr); else - eval.template gradients<0, false, true>( - gradients_quad_ptr, values_dofs_actual_ptr); - - gradients_quad_ptr += n_q_points; + 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; } } } - template - template - inline void - FEEvaluationImpl:: - evaluate_or_integrate( - const EvaluationFlags::EvaluationFlags evaluation_flag, - Number *values_dofs_actual, - FEEvaluationData &fe_eval, - const bool add_into_values_array) - { - if (evaluation_flag == EvaluationFlags::nothing) - return; - - AssertDimension(fe_eval.get_shape_info().data.size(), 2); - // First component: - evaluate_tensor_product_per_component<0>( - evaluation_flag, - values_dofs_actual, - fe_eval, - add_into_values_array, - std::integral_constant()); - // Second component : - evaluate_tensor_product_per_component<1>( - evaluation_flag, - values_dofs_actual, - fe_eval, - add_into_values_array, - std::integral_constant()); - if (dim == 3) - { - // Third component - evaluate_tensor_product_per_component<2>( - evaluation_flag, - values_dofs_actual, - fe_eval, - add_into_values_array, - std::integral_constant()); - } - } - - // Helper function that applies the 1d evaluation kernels. - // std::integral_constant is the interpolation path, and - // std::integral_constant below is the integration path. - template - template - inline void - FEEvaluationImpl:: - evaluate_tensor_product_per_component( - const EvaluationFlags::EvaluationFlags evaluation_flag, - Number *values_dofs_actual, - FEEvaluationData &fe_eval, - const bool add_into_values_array, - std::integral_constant) - { - (void)add_into_values_array; - - using Number2 = - typename FEEvaluationData::shape_info_number_type; - using EvalNormal = - EvaluatorTensorProductAnisotropic; - - using EvalTangent = - EvaluatorTensorProductAnisotropic; - using Eval0 = std::conditional_t; - using Eval1 = std::conditional_t; - using Eval2 = std::conditional_t; - - const auto &shape_info = fe_eval.get_shape_info(); - Eval0 eval0 = create_evaluator_tensor_product( - ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1])); - Eval1 eval1 = create_evaluator_tensor_product( - ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1])); - Eval2 eval2 = create_evaluator_tensor_product( - ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1])); - - Number *temp1 = fe_eval.get_scratch_data().begin(); - Number *temp2; - - temp2 = - temp1 + - std::max(Utilities::fixed_power(shape_info.data[0].fe_degree + 1), - Utilities::fixed_power(shape_info.data[0].n_q_points_1d)); - - const std::size_t n_q_points = shape_info.n_q_points; - const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell; - - // Initial shift depending on component (normal_dir) - Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir; - Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir; - Number *gradients_quad = - fe_eval.begin_gradients() + dim * n_q_points * normal_dir; - Number *hessians_quad = - (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir : - fe_eval.begin_hessians() + 6 * n_q_points * normal_dir; - - switch (dim) - { - case 2: - if (evaluation_flag & EvaluationFlags::gradients) - { - eval0.template gradients<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, gradients_quad); - } - if (evaluation_flag & EvaluationFlags::hessians) - { - // The evaluation/integration here *should* work, however - // the piola transform is not implemented. - AssertThrow(false, ExcNotImplemented()); - // grad xy - if (!(evaluation_flag & EvaluationFlags::gradients)) - eval0.template gradients<0, true, false>(values_dofs, temp1); - eval1.template gradients<1, true, false>(temp1, - hessians_quad + - 2 * n_q_points); - - // grad xx - eval0.template hessians<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, hessians_quad); - } - - // grad y - eval0.template values<0, true, false>(values_dofs, temp1); - if (evaluation_flag & EvaluationFlags::gradients) - eval1.template gradients<1, true, false>(temp1, - gradients_quad + - n_q_points); - - // grad yy - if (evaluation_flag & EvaluationFlags::hessians) - eval1.template hessians<1, true, false>(temp1, - hessians_quad + n_q_points); - - // val: can use values applied in x - if (evaluation_flag & EvaluationFlags::values) - eval1.template values<1, true, false>(temp1, values_quad); - break; - case 3: - if (evaluation_flag & EvaluationFlags::gradients) - { - // grad x - eval0.template gradients<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, gradients_quad); - } - - if (evaluation_flag & EvaluationFlags::hessians) - { - // The evaluation/integration here *should* work, however - // the piola transform is not implemented. - AssertThrow(false, ExcNotImplemented()); - // grad xz - if (!(evaluation_flag & EvaluationFlags::gradients)) - { - eval0.template gradients<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, temp2); - } - eval2.template gradients<2, true, false>(temp2, - hessians_quad + - 4 * n_q_points); - - // grad xy - eval1.template gradients<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, - hessians_quad + - 3 * n_q_points); - - // grad xx - eval0.template hessians<0, true, false>(values_dofs, temp1); - eval1.template values<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, hessians_quad); - } - - // grad y - eval0.template values<0, true, false>(values_dofs, temp1); - if (evaluation_flag & EvaluationFlags::gradients) - { - eval1.template gradients<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, - gradients_quad + - n_q_points); - } - - if (evaluation_flag & EvaluationFlags::hessians) - { - // grad yz - if (!(evaluation_flag & EvaluationFlags::gradients)) - eval1.template gradients<1, true, false>(temp1, temp2); - eval2.template gradients<2, true, false>(temp2, - hessians_quad + - 5 * n_q_points); - - // grad yy - eval1.template hessians<1, true, false>(temp1, temp2); - eval2.template values<2, true, false>(temp2, - hessians_quad + n_q_points); - } - - // grad z: can use the values applied in x direction stored in - // temp1 - eval1.template values<1, true, false>(temp1, temp2); - if (evaluation_flag & EvaluationFlags::gradients) - eval2.template gradients<2, true, false>(temp2, - gradients_quad + - 2 * n_q_points); - - // grad zz: can use the values applied in x and y direction stored - // in temp2 - if (evaluation_flag & EvaluationFlags::hessians) - eval2.template hessians<2, true, false>(temp2, - hessians_quad + - 2 * n_q_points); - - // val: can use the values applied in x & y direction stored in - // temp2 - if (evaluation_flag & EvaluationFlags::values) - eval2.template values<2, true, false>(temp2, values_quad); - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - } - - template - template - inline void - FEEvaluationImpl:: - evaluate_tensor_product_per_component( - const EvaluationFlags::EvaluationFlags evaluation_flag, - Number *values_dofs_actual, - FEEvaluationData &fe_eval, - const bool add_into_values_array, - std::integral_constant) - { - using Number2 = - typename FEEvaluationData::shape_info_number_type; - using EvalNormal = - EvaluatorTensorProductAnisotropic; - - using EvalTangent = - EvaluatorTensorProductAnisotropic; - using Eval0 = std::conditional_t; - using Eval1 = std::conditional_t; - using Eval2 = std::conditional_t; - - const auto &shape_info = fe_eval.get_shape_info(); - Eval0 eval0 = create_evaluator_tensor_product( - ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1])); - Eval1 eval1 = create_evaluator_tensor_product( - ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1])); - Eval2 eval2 = create_evaluator_tensor_product( - ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1])); - - Number *temp1 = fe_eval.get_scratch_data().begin(); - Number *temp2; - - temp2 = - temp1 + - std::max(Utilities::fixed_power(shape_info.data[0].fe_degree + 1), - Utilities::fixed_power(shape_info.data[0].n_q_points_1d)); - - const std::size_t n_q_points = shape_info.n_q_points; - const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell; - - // Initial shift depending on component (normal_dir) - Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir; - Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir; - Number *gradients_quad = - fe_eval.begin_gradients() + dim * n_q_points * normal_dir; - Number *hessians_quad = - (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir : - fe_eval.begin_hessians() + 6 * n_q_points * normal_dir; - - // Integrate path - switch (dim) - { - case 2: - if ((evaluation_flag & EvaluationFlags::values) && - !(evaluation_flag & EvaluationFlags::gradients)) - { - eval1.template values<1, false, false>(values_quad, temp1); - if (add_into_values_array == false) - eval0.template values<0, false, false>(temp1, values_dofs); - else - eval0.template values<0, false, true>(temp1, values_dofs); - } - if (evaluation_flag & EvaluationFlags::gradients) - { - eval1.template gradients<1, false, false>(gradients_quad + - n_q_points, - temp1); - if ((evaluation_flag & EvaluationFlags::values)) - eval1.template values<1, false, true>(values_quad, temp1); - if (add_into_values_array == false) - eval0.template values<0, false, false>(temp1, values_dofs); - else - eval0.template values<0, false, true>(temp1, values_dofs); - eval1.template values<1, false, false>(gradients_quad, temp1); - eval0.template gradients<0, false, true>(temp1, values_dofs); - } - if (evaluation_flag & EvaluationFlags::hessians) - { - // grad xx - eval1.template values<1, false, false>(hessians_quad, temp1); - - if ((evaluation_flag & EvaluationFlags::values) || - (evaluation_flag & EvaluationFlags::gradients) || - add_into_values_array == true) - eval0.template hessians<0, false, true>(temp1, values_dofs); - else - eval0.template hessians<0, false, false>(temp1, values_dofs); - - // grad yy - eval1.template hessians<1, false, false>(hessians_quad + - n_q_points, - temp1); - eval0.template values<0, false, true>(temp1, values_dofs); - - // grad xy - eval1.template gradients<1, false, false>(hessians_quad + - 2 * n_q_points, - temp1); - eval0.template gradients<0, false, true>(temp1, values_dofs); - } - break; - - case 3: - if ((evaluation_flag & EvaluationFlags::values) && - !(evaluation_flag & EvaluationFlags::gradients)) - { - eval2.template values<2, false, false>(values_quad, temp1); - eval1.template values<1, false, false>(temp1, temp2); - if (add_into_values_array == false) - eval0.template values<0, false, false>(temp2, values_dofs); - else - eval0.template values<0, false, true>(temp2, values_dofs); - } - if (evaluation_flag & EvaluationFlags::gradients) - { - eval2.template gradients<2, false, false>(gradients_quad + - 2 * n_q_points, - temp1); - if ((evaluation_flag & EvaluationFlags::values)) - eval2.template values<2, false, true>(values_quad, temp1); - eval1.template values<1, false, false>(temp1, temp2); - eval2.template values<2, false, false>(gradients_quad + - n_q_points, - temp1); - eval1.template gradients<1, false, true>(temp1, temp2); - if (add_into_values_array == false) - eval0.template values<0, false, false>(temp2, values_dofs); - else - eval0.template values<0, false, true>(temp2, values_dofs); - eval2.template values<2, false, false>(gradients_quad, temp1); - eval1.template values<1, false, false>(temp1, temp2); - eval0.template gradients<0, false, true>(temp2, values_dofs); - } - if (evaluation_flag & EvaluationFlags::hessians) - { - // grad xx - eval2.template values<2, false, false>(hessians_quad, temp1); - eval1.template values<1, false, false>(temp1, temp2); - - if ((evaluation_flag & EvaluationFlags::values) || - (evaluation_flag & EvaluationFlags::gradients) || - add_into_values_array == true) - eval0.template hessians<0, false, true>(temp2, values_dofs); - else - eval0.template hessians<0, false, false>(temp2, values_dofs); - - // grad yy - eval2.template values<2, false, false>(hessians_quad + n_q_points, - temp1); - eval1.template hessians<1, false, false>(temp1, temp2); - eval0.template values<0, false, true>(temp2, values_dofs); - - // grad zz - eval2.template hessians<2, false, false>(hessians_quad + - 2 * n_q_points, - temp1); - eval1.template values<1, false, false>(temp1, temp2); - eval0.template values<0, false, true>(temp2, values_dofs); - - // grad xy - eval2.template values<2, false, false>(hessians_quad + - 3 * n_q_points, - temp1); - eval1.template gradients<1, false, false>(temp1, temp2); - eval0.template gradients<0, false, true>(temp2, values_dofs); - - // grad xz - eval2.template gradients<2, false, false>(hessians_quad + - 4 * n_q_points, - temp1); - eval1.template values<1, false, false>(temp1, temp2); - eval0.template gradients<0, false, true>(temp2, values_dofs); - - // grad yz - eval2.template gradients<2, false, false>(hessians_quad + - 5 * n_q_points, - temp1); - eval1.template gradients<1, false, false>(temp1, temp2); - eval0.template values<0, false, true>(temp2, values_dofs); - } - - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - } /** * This struct implements the change between two different bases. This is an @@ -1433,22 +773,12 @@ namespace internal * * This class allows for dimension-independent application of the operation, * implemented by template recursion. It has been tested up to 6d. - * - * The last two template arguments in this class are unused. They have a - * default type set to bool, but any type will compile. These arguments are - * present for backward compatibility of this (internal) interface with - * previous versions of deal.II, where the template arguments were used to - * indicate the number types These are now part of the do_forward() and - * do_backward() functions, and they will be removed from a future version - * of deal.II. */ template + int basis_size_2> struct FEEvaluationImplBasisChange { static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2, @@ -1827,279 +1157,556 @@ namespace internal /** - * This struct performs the evaluation of function values, gradients and - * Hessians for tensor-product finite elements. This a specialization for - * elements where the nodal points coincide with the quadrature points like - * FE_Q shape functions on Gauss-Lobatto elements integrated with - * Gauss-Lobatto quadrature. The assumption of this class is that the shape - * 'values' operation is identity, which allows us to write shorter code. - * - * In literature, this form of evaluation is often called spectral - * evaluation, spectral collocation or simply collocation, meaning the same - * location for shape functions and evaluation space (quadrature points). + * Internal function that evaluates the gradients of finite element + * functions represented by bases in the collocation space, used by + * FEEvaluationImplCollocation and FEEvaluationImplTransformToCollocation. + * The evaluation strategy uses sum factorization with the even-odd + * optimization and fixed loop bounds. */ - template - struct FEEvaluationImplCollocation + template + inline void + evaluate_gradients_collocation( + const MatrixFreeFunctions::UnivariateShapeData &shape, + const Number *values, + Number *gradients) { - using Number2 = - typename FEEvaluationData::shape_info_number_type; - using Eval = EvaluatorTensorProduct; - - static void - evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - FEEvaluationData &fe_eval); - - static void - do_evaluate(const MatrixFreeFunctions::UnivariateShapeData &shape, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - Number *gradients_quad, - Number *hessians_quad); - - static void - integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, - FEEvaluationData &fe_eval, - const bool add_into_values_array); + AssertDimension(shape.shape_gradients_collocation_eo.size(), + (n_points_1d + 1) / 2 * n_points_1d); - static void - do_integrate(const MatrixFreeFunctions::UnivariateShapeData &shape, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, - Number *gradients_quad, - const Number *hessians_quad, - const bool add_into_values_array); - }; + EvaluatorTensorProduct + eval({}, shape.shape_gradients_collocation_eo, {}); + EvaluatorTensorProduct + eval_2d({}, shape.shape_gradients_collocation_eo, {}); + + if (dim == 1) + eval.template gradients<0, true, false>(values, gradients); + else + { + if (dim > 2) + eval.template gradients<2, true, false, dim>(values, gradients + 2); + constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1); + constexpr unsigned int n_points_2d = n_points_1d * n_points_1d; + const Number *in = values + (loop_bound - 1) * n_points_2d; + Number *out = gradients + (loop_bound - 1) * dim * n_points_2d; + for (unsigned int l = 0; l < loop_bound; ++l) + { + eval_2d.template gradients<0, true, false, dim>(in, out); + eval_2d.template gradients<1, true, false, dim>(in, out + 1); + in -= n_points_2d; + out -= dim * n_points_2d; + } + } + } - template + /** + * Internal function that multiplies by the gradients of test functions and + * sums over quadrature points for function representations in the + * collocation space, used by FEEvaluationImplCollocation and + * FEEvaluationImplTransformToCollocation. The evaluation strategy uses sum + * factorization with the even-odd optimization and fixed loop bounds. + */ + template inline void - FEEvaluationImplCollocation::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - FEEvaluationData &fe_eval) + integrate_gradients_collocation( + const MatrixFreeFunctions::UnivariateShapeData &shape, + Number *values, + const Number *gradients, + const bool add_into_values_array) { - constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); + AssertDimension(shape.shape_gradients_collocation_eo.size(), + (n_points_1d + 1) / 2 * n_points_1d); - for (unsigned int c = 0; c < n_components; ++c) + EvaluatorTensorProduct + eval({}, shape.shape_gradients_collocation_eo, {}); + EvaluatorTensorProduct + eval_2d({}, shape.shape_gradients_collocation_eo, {}); + + if (dim == 1) + { + if (add_into_values_array) + eval.template gradients<0, false, true>(gradients, values); + else + eval.template gradients<0, false, false>(gradients, values); + } + else { - if ((evaluation_flag & EvaluationFlags::values) != 0u) - for (unsigned int i = 0; i < n_points; ++i) - fe_eval.begin_values()[n_points * c + i] = - values_dofs[n_points * c + i]; - - do_evaluate(fe_eval.get_shape_info().data.front(), - evaluation_flag, - values_dofs + c * n_points, - fe_eval.begin_gradients() + c * dim * n_points, - fe_eval.begin_hessians() + - c * dim * (dim + 1) / 2 * n_points); + constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1); + constexpr unsigned int n_points_2d = n_points_1d * n_points_1d; + + const Number *in = gradients + (loop_bound - 1) * dim * n_points_2d; + Number *out = values + (loop_bound - 1) * n_points_2d; + for (unsigned int l = 0; l < loop_bound; ++l) + { + if (add_into_values_array) + eval_2d.template gradients<0, false, true, dim>(in, out); + else + eval_2d.template gradients<0, false, false, dim>(in, out); + eval_2d.template gradients<1, false, true, dim>(in + 1, out); + in -= dim * n_points_2d; + out -= n_points_2d; + } } + if (dim > 2) + eval.template gradients<2, false, true, dim>(gradients + 2, values); } - template + /** + * Internal function that evaluates the Hessians of finite element functions + * represented by bases in the collocation space, used by + * FEEvaluationImplSelector. The evaluation strategy uses sum + * factorization with fixed loop bounds. + */ + template inline void - FEEvaluationImplCollocation::do_evaluate( - const MatrixFreeFunctions::UnivariateShapeData &shape, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - Number *gradients_quad, - Number *hessians_quad) + evaluate_hessians_collocation(const unsigned int n_components, + FEEvaluationData &fe_eval) { - AssertDimension(shape.shape_gradients_collocation_eo.size(), - (fe_degree + 2) / 2 * (fe_degree + 1)); - constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); - - Eval eval({}, - shape.shape_gradients_collocation_eo, - shape.shape_hessians_collocation_eo); - if ((evaluation_flag & - (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u) + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + // might have non-symmetric quadrature formula, so use the more + // conservative 'evaluate_general' scheme rather than 'even_odd' as the + // Hessians are not used very often + const MatrixFreeFunctions::UnivariateShapeData &data = + fe_eval.get_shape_info().data[0]; + AssertDimension(data.shape_gradients_collocation.size(), + data.n_q_points_1d * data.n_q_points_1d); + EvaluatorTensorProduct + eval({}, + data.shape_gradients_collocation.data(), + data.shape_hessians_collocation.data(), + data.n_q_points_1d, + data.n_q_points_1d); + + const Number *values = fe_eval.begin_values(); + Number *hessians = fe_eval.begin_hessians(); + Number *scratch = fe_eval.get_scratch_data().begin(); + const std::size_t n_points = fe_eval.get_shape_info().n_q_points; + for (unsigned int comp = 0; comp < n_components; ++comp) { - eval.template gradients<0, true, false>(values_dofs, gradients_quad); + // xx derivative + eval.template hessians<0, true, false>(values, hessians); if (dim > 1) - eval.template gradients<1, true, false>(values_dofs, - gradients_quad + n_points); + { + // xy derivative: we might or might not have the gradients already + // computed elsewhere, but we recompute them here since it adds + // only moderate extra work (at most 25%) + eval.template gradients<0, true, false>(values, scratch); + eval.template gradients<1, true, false>(scratch, + hessians + dim * n_points); + // yy derivative + eval.template hessians<1, true, false>(values, hessians + n_points); + } if (dim > 2) - eval.template gradients<2, true, false>(values_dofs, - gradients_quad + - 2 * n_points); + { + // xz derivative + eval.template gradients<2, true, false>(scratch, + hessians + 4 * n_points); + // yz derivative + eval.template gradients<1, true, false>(values, scratch); + eval.template gradients<2, true, false>(scratch, + hessians + 5 * n_points); + // zz derivative + eval.template hessians<2, true, false>(values, + hessians + 2 * n_points); + } + + values += n_points; + hessians += (dim * (dim + 1)) / 2 * n_points; } - if (evaluation_flag & EvaluationFlags::hessians) + } + + + + /** + * Internal function that multiplies by the Hessians of test functions and + * sums over quadrature points for function representations in the + * collocation space, used by FEEvaluationImplSelector. The evaluation + * strategy uses sum factorization with fixed loop bounds. + */ + template + inline void + integrate_hessians_collocation(const unsigned int n_components, + FEEvaluationData &fe_eval, + const bool add_into_values_array) + { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + + const MatrixFreeFunctions::UnivariateShapeData &data = + fe_eval.get_shape_info().data[0]; + AssertDimension(data.shape_gradients_collocation.size(), + data.n_q_points_1d * data.n_q_points_1d); + EvaluatorTensorProduct + eval({}, + data.shape_gradients_collocation.data(), + data.shape_hessians_collocation.data(), + data.n_q_points_1d, + data.n_q_points_1d); + Number *values = fe_eval.begin_values(); + const Number *hessians = fe_eval.begin_hessians(); + Number *scratch = fe_eval.get_scratch_data().begin(); + const std::size_t n_points = fe_eval.get_shape_info().n_q_points; + + for (unsigned int comp = 0; comp < n_components; ++comp) { - eval.template hessians<0, true, false>(values_dofs, hessians_quad); + // xx derivative + if (add_into_values_array == true) + eval.template hessians<0, false, true>(hessians, values); + else + eval.template hessians<0, false, false>(hessians, values); + + // yy derivative if (dim > 1) + eval.template hessians<1, false, true>(hessians + n_points, values); + if (dim > 2) { - eval.template gradients<1, true, false>(gradients_quad, - hessians_quad + - dim * n_points); - eval.template hessians<1, true, false>(values_dofs, - hessians_quad + n_points); + // zz derivative + eval.template hessians<2, false, true>(hessians + 2 * n_points, + values); + // yz derivative + eval.template gradients<2, false, false>(hessians + 5 * n_points, + scratch); + eval.template gradients<1, false, true>(scratch, values); + + // xz derivative + eval.template gradients<2, false, false>(hessians + 4 * n_points, + scratch); } - if (dim > 2) + + if (dim > 1) { - eval.template gradients<2, true, false>(gradients_quad, - hessians_quad + - 4 * n_points); - eval.template gradients<2, true, false>(gradients_quad + n_points, - hessians_quad + - 5 * n_points); - eval.template hessians<2, true, false>(values_dofs, - hessians_quad + - 2 * n_points); + // xy derivative, combined with xz in 3d + eval.template gradients<1, false, (dim > 2)>(hessians + + dim * n_points, + scratch); + eval.template gradients<0, false, true>(scratch, values); } + + values += n_points; + hessians += (dim * (dim + 1)) / 2 * n_points; } } - template - inline void - FEEvaluationImplCollocation::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, - FEEvaluationData &fe_eval, - const bool add_into_values_array) + /** + * Internal function to evaluate the Hessians of finite element functions in + * the non-collocation setting as a fall-back. The evaluation strategy uses + * sum factorization with run-time loop bounds and is thus slower than the + * collocation case, but it is not as widely used and thus uncritical. + */ + template + void + evaluate_hessians_slow(const unsigned int n_components, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &univariate_shape_data = fe_eval.get_shape_info().data; + using Impl = + FEEvaluationImpl; + using Eval = typename Impl::Eval; + Eval eval0 = + Impl::create_evaluator_tensor_product(&univariate_shape_data[0]); + Eval eval1 = Impl::create_evaluator_tensor_product( + &univariate_shape_data[std::min(1, + univariate_shape_data.size() - 1)]); + Eval eval2 = Impl::create_evaluator_tensor_product( + &univariate_shape_data[std::min(2, + univariate_shape_data.size() - 1)]); + + const unsigned int n_points = fe_eval.get_shape_info().n_q_points; + Number *tmp1 = fe_eval.get_scratch_data().begin(); + Number *tmp2 = + tmp1 + std::max(Utilities::fixed_power( + univariate_shape_data.front().fe_degree + 1), + Utilities::fixed_power( + univariate_shape_data.front().n_q_points_1d)); + Number *hessians = fe_eval.begin_hessians(); + + for (unsigned int comp = 0; comp < n_components; + ++comp, + hessians += n_points * dim * (dim + 1) / 2, + values_dofs += + fe_eval.get_shape_info().dofs_per_component_on_cell) + switch (dim) + { + case 1: + eval0.template hessians<0, true, false>(values_dofs, hessians); + break; + case 2: + // xx derivative + eval0.template hessians<0, true, false>(values_dofs, tmp1); + eval1.template values<1, true, false>(tmp1, hessians); + // xy derivative + eval0.template gradients<0, true, false>(values_dofs, tmp1); + eval1.template gradients<1, true, false>(tmp1, + hessians + 2 * n_points); + // yy derivative + eval0.template values<0, true, false>(values_dofs, tmp1); + eval1.template hessians<1, true, false>(tmp1, hessians + n_points); + break; + case 3: + // xx derivative + eval0.template hessians<0, true, false>(values_dofs, tmp1); + eval1.template values<1, true, false>(tmp1, tmp2); + eval2.template values<2, true, false>(tmp2, hessians); + // xy derivative + eval0.template gradients<0, true, false>(values_dofs, tmp1); + eval1.template gradients<1, true, false>(tmp1, tmp2); + eval2.template values<2, true, false>(tmp2, + hessians + 3 * n_points); + // xz derivative + eval1.template values<1, true, false>(tmp1, tmp2); + eval2.template gradients<2, true, false>(tmp2, + hessians + 4 * n_points); + // yy derivative + eval0.template values<0, true, false>(values_dofs, tmp1); + eval1.template hessians<1, true, false>(tmp1, tmp2); + eval2.template values<2, true, false>(tmp2, hessians + n_points); + // yz derivative + eval1.template gradients<1, true, false>(tmp1, tmp2); + eval2.template gradients<2, true, false>(tmp2, + hessians + 5 * n_points); + // zz derivative + eval1.template values<1, true, false>(tmp1, tmp2); + eval2.template hessians<2, true, false>(tmp2, + hessians + 2 * n_points); + break; + + default: + Assert(false, + ExcNotImplemented( + "Only 1d, 2d and 3d implemented for Hessian")); + } + } + + + + /** + * Internal function to multiply by the Hessians of the test functions and + * integrate in the non-collocation setting as a fall-back. The evaluation + * strategy uses sum factorization with run-time loop bounds and is thus + * slower than the collocation case, but it is not as widely used and thus + * uncritical. + */ + template + void + integrate_hessians_slow(const unsigned int n_components, + const FEEvaluationData &fe_eval, + Number *values_dofs, + const bool add_into_values_array) { - constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); + const auto &univariate_shape_data = fe_eval.get_shape_info().data; + using Impl = + FEEvaluationImpl; + using Eval = typename Impl::Eval; + Eval eval0 = + Impl::create_evaluator_tensor_product(&univariate_shape_data[0]); + Eval eval1 = Impl::create_evaluator_tensor_product( + &univariate_shape_data[std::min(1, + univariate_shape_data.size() - 1)]); + Eval eval2 = Impl::create_evaluator_tensor_product( + &univariate_shape_data[std::min(2, + univariate_shape_data.size() - 1)]); + + const unsigned int n_points = fe_eval.get_shape_info().n_q_points; + Number *tmp1 = fe_eval.get_scratch_data().begin(); + Number *tmp2 = + tmp1 + std::max(Utilities::fixed_power( + univariate_shape_data.front().fe_degree + 1), + Utilities::fixed_power( + univariate_shape_data.front().n_q_points_1d)); + const Number *hessians = fe_eval.begin_hessians(); + + for (unsigned int comp = 0; comp < n_components; + ++comp, + hessians += n_points * dim * (dim + 1) / 2, + values_dofs += + fe_eval.get_shape_info().dofs_per_component_on_cell) + switch (dim) + { + case 1: + if (add_into_values_array) + eval0.template hessians<0, false, true>(hessians, values_dofs); + else + eval0.template hessians<0, false, false>(hessians, values_dofs); + break; + case 2: + // xx derivative + eval1.template values<1, false, false>(hessians, tmp1); + if (add_into_values_array) + eval0.template hessians<0, false, true>(tmp1, values_dofs); + else + eval0.template hessians<0, false, false>(tmp1, values_dofs); + + // xy derivative + eval1.template gradients<1, false, false>(hessians + 2 * n_points, + tmp1); + eval0.template gradients<0, false, true>(tmp1, values_dofs); + // yy derivative + eval1.template hessians<1, false, false>(hessians + n_points, tmp1); + eval0.template values<0, false, true>(tmp1, values_dofs); + break; + case 3: + // xx derivative + eval2.template values<2, false, false>(hessians, tmp1); + eval1.template values<1, false, false>(tmp1, tmp2); - for (unsigned int c = 0; c < n_components; ++c) - { - if ((integration_flag & EvaluationFlags::values) != 0u) - { if (add_into_values_array) - for (unsigned int i = 0; i < n_points; ++i) - values_dofs[n_points * c + i] += - fe_eval.begin_values()[n_points * c + i]; + eval0.template hessians<0, false, true>(tmp2, values_dofs); else - for (unsigned int i = 0; i < n_points; ++i) - values_dofs[n_points * c + i] = - fe_eval.begin_values()[n_points * c + i]; - } + eval0.template hessians<0, false, false>(tmp2, values_dofs); + + // xy derivative + eval2.template values<2, false, false>(hessians + 3 * n_points, + tmp1); + eval1.template gradients<1, false, false>(tmp1, tmp2); + // xz derivative + eval2.template gradients<2, false, false>(hessians + 4 * n_points, + tmp1); + eval1.template values<1, false, true>(tmp1, tmp2); + eval1.template values<0, false, true>(tmp2, values_dofs); + + // yy derivative + eval2.template values<2, false, false>(hessians + n_points, tmp1); + eval1.template hessians<1, false, false>(tmp1, tmp2); + + // yz derivative + eval2.template gradients<2, false, false>(hessians + 5 * n_points, + tmp1); + eval1.template gradients<1, false, true>(tmp1, tmp2); + + // zz derivative + eval2.template hessians<2, false, false>(hessians + 2 * n_points, + tmp1); + eval1.template values<1, false, true>(tmp1, tmp2); + eval0.template values<0, false, true>(tmp2, values_dofs); + break; - do_integrate(fe_eval.get_shape_info().data.front(), - integration_flag, - values_dofs + c * n_points, - fe_eval.begin_gradients() + c * dim * n_points, - fe_eval.begin_hessians() + - c * dim * (dim + 1) / 2 * n_points, - add_into_values_array || - ((integration_flag & EvaluationFlags::values) != 0u)); - } + default: + Assert(false, + ExcNotImplemented( + "Only 1d, 2d and 3d implemented for Hessian")); + } } + /** + * This struct performs the evaluation of function values, gradients and + * Hessians for tensor-product finite elements. This a specialization for + * elements where the nodal points coincide with the quadrature points like + * FE_Q shape functions on Gauss-Lobatto elements integrated with + * Gauss-Lobatto quadrature. The assumption of this class is that the shape + * 'values' operation is identity, which allows us to write shorter code. + * + * In literature, this form of evaluation is often called spectral + * evaluation, spectral collocation or simply collocation, meaning the same + * location for shape functions and evaluation space (quadrature points). + */ template - inline void - FEEvaluationImplCollocation::do_integrate( - const MatrixFreeFunctions::UnivariateShapeData &shape, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, - Number *gradients_quad, - const Number *hessians_quad, - const bool add_into_values_array) + struct FEEvaluationImplCollocation { - AssertDimension(shape.shape_gradients_collocation_eo.size(), - (fe_degree + 2) / 2 * (fe_degree + 1)); + using Number2 = + typename FEEvaluationData::shape_info_number_type; + using Eval = EvaluatorTensorProduct; - Eval eval({}, - shape.shape_gradients_collocation_eo, - shape.shape_hessians_collocation_eo); - constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); + static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); - if ((integration_flag & EvaluationFlags::hessians) != 0u) - { - // diagonal - // grad xx - if (add_into_values_array == true) - eval.template hessians<0, false, true>(hessians_quad, values_dofs); - else - eval.template hessians<0, false, false>(hessians_quad, values_dofs); - // grad yy - if (dim > 1) - eval.template hessians<1, false, true>(hessians_quad + n_points, - values_dofs); - // grad zz - if (dim > 2) - eval.template hessians<2, false, true>(hessians_quad + 2 * n_points, - values_dofs); - // off-diagonal - if (dim == 2) - { - // grad xy, queue into gradient - if (integration_flag & EvaluationFlags::gradients) - eval.template gradients<1, false, true>(hessians_quad + - 2 * n_points, - gradients_quad); - else - eval.template gradients<1, false, false>(hessians_quad + - 2 * n_points, - gradients_quad); - } - if (dim == 3) - { - // grad xy, queue into gradient - if (integration_flag & EvaluationFlags::gradients) - eval.template gradients<1, false, true>(hessians_quad + - 3 * n_points, - gradients_quad); - else - eval.template gradients<1, false, false>(hessians_quad + - 3 * n_points, - gradients_quad); - - // grad xz - eval.template gradients<2, false, true>(hessians_quad + - 4 * n_points, - gradients_quad); - - // grad yz - if (integration_flag & EvaluationFlags::gradients) - eval.template gradients<2, false, true>( - hessians_quad + 5 * n_points, gradients_quad + n_points); - else - eval.template gradients<2, false, false>( - hessians_quad + 5 * n_points, gradients_quad + n_points); - } + for (unsigned int c = 0; c < n_components; ++c) + { + if ((evaluation_flag & EvaluationFlags::values) != 0u) + for (unsigned int i = 0; i < n_points; ++i) + fe_eval.begin_values()[n_points * c + i] = + values_dofs[n_points * c + i]; + + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + evaluate_gradients_collocation( + fe_eval.get_shape_info().data.front(), + values_dofs + c * n_points, + fe_eval.begin_gradients() + c * dim * n_points); + } + } - // if we did not integrate gradients, set the last slot to zero - // which was not touched before, in order to avoid the if - // statement in the gradients loop below - if ((integration_flag & EvaluationFlags::gradients) == 0u) - for (unsigned int q = 0; q < n_points; ++q) - gradients_quad[(dim - 1) * n_points + q] = Number(); - } + static void + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool add_into_values_array) + { + constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); - if ((integration_flag & - (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u) - { - if (add_into_values_array || - (integration_flag & EvaluationFlags::hessians) != 0u) - eval.template gradients<0, false, true>(gradients_quad, values_dofs); - else - eval.template gradients<0, false, false>(gradients_quad, values_dofs); - if (dim > 1) - eval.template gradients<1, false, true>(gradients_quad + n_points, - values_dofs); - if (dim > 2) - eval.template gradients<2, false, true>(gradients_quad + 2 * n_points, - values_dofs); - } - } + for (unsigned int c = 0; c < n_components; ++c) + { + if ((integration_flag & EvaluationFlags::values) != 0u) + { + if (add_into_values_array) + for (unsigned int i = 0; i < n_points; ++i) + values_dofs[n_points * c + i] += + fe_eval.begin_values()[n_points * c + i]; + else + for (unsigned int i = 0; i < n_points; ++i) + values_dofs[n_points * c + i] = + fe_eval.begin_values()[n_points * c + i]; + } + + if ((integration_flag & EvaluationFlags::gradients) != 0u) + integrate_gradients_collocation( + fe_eval.get_shape_info().data.front(), + values_dofs + c * n_points, + fe_eval.begin_gradients() + c * dim * n_points, + add_into_values_array || + ((integration_flag & EvaluationFlags::values) != 0u)); + } + } + }; @@ -2120,116 +1727,371 @@ namespace internal evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, - FEEvaluationData &fe_eval); + FEEvaluationData &fe_eval) + { + const auto &shape_data = fe_eval.get_shape_info().data.front(); + + Assert(n_q_points_1d > fe_degree, + ExcMessage("You lose information when going to a collocation " + "space of lower degree, so the evaluation results " + "would be wrong. Thus, this class does not permit " + "the chosen operation.")); + constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim); + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim); + + for (unsigned int c = 0; c < n_components; ++c) + { + FEEvaluationImplBasisChange< + evaluate_evenodd, + EvaluatorQuantity::value, + dim, + (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), + n_q_points_1d>::do_forward(1, + shape_data.shape_values_eo, + values_dofs + c * n_dofs, + fe_eval.begin_values() + c * n_q_points); + + // apply derivatives in the collocation space + if (evaluation_flag & EvaluationFlags::gradients) + evaluate_gradients_collocation( + shape_data, + fe_eval.begin_values() + c * n_q_points, + fe_eval.begin_gradients() + c * dim * n_q_points); + } + } static void integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, + const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData &fe_eval, - const bool add_into_values_array); + const bool add_into_values_array) + { + const auto &shape_data = fe_eval.get_shape_info().data.front(); + + Assert(n_q_points_1d > fe_degree, + ExcMessage("You lose information when going to a collocation " + "space of lower degree, so the evaluation results " + "would be wrong. Thus, this class does not permit " + "the chosen operation.")); + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim); + + for (unsigned int c = 0; c < n_components; ++c) + { + // apply derivatives in collocation space + if (integration_flag & EvaluationFlags::gradients) + integrate_gradients_collocation( + shape_data, + fe_eval.begin_values() + c * n_q_points, + fe_eval.begin_gradients() + c * dim * n_q_points, + /*add_into_values_array=*/ + integration_flag & EvaluationFlags::values); + + // transform back to the original space + FEEvaluationImplBasisChange< + evaluate_evenodd, + EvaluatorQuantity::value, + dim, + (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), + n_q_points_1d>::do_backward(1, + shape_data.shape_values_eo, + add_into_values_array, + fe_eval.begin_values() + c * n_q_points, + values_dofs + + c * + Utilities::pow(fe_degree + 1, dim)); + } + } }; + /** + * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use + * specific sum-factorization kernels and with normal/tangential shape_data + */ template - inline void - FEEvaluationImplTransformToCollocation< - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - FEEvaluationData &fe_eval) + struct FEEvaluationImpl { - const auto &shape_data = fe_eval.get_shape_info().data.front(); + using Number2 = + typename FEEvaluationData::shape_info_number_type; - Assert(n_q_points_1d > fe_degree, - ExcMessage("You lose information when going to a collocation space " - "of lower degree, so the evaluation results would be " - "wrong. Thus, this class does not permit the desired " - "operation.")); - constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim); - constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim); + template + static void + evaluate_or_integrate( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number *values_dofs_actual, + FEEvaluationData &fe_eval, + const bool add_into_values_array = false); - for (unsigned int c = 0; c < n_components; ++c) - { - FEEvaluationImplBasisChange< - evaluate_evenodd, - EvaluatorQuantity::value, - dim, - (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), - n_q_points_1d>::do_forward(1, - shape_data.shape_values_eo, - values_dofs + c * n_dofs, - fe_eval.begin_values() + c * n_q_points); - - // apply derivatives in the collocation space - if (evaluation_flag & - (EvaluationFlags::gradients | EvaluationFlags::hessians)) - FEEvaluationImplCollocation:: - do_evaluate(shape_data, - evaluation_flag & (EvaluationFlags::gradients | - EvaluationFlags::hessians), - fe_eval.begin_values() + c * n_q_points, - fe_eval.begin_gradients() + c * dim * n_q_points, - fe_eval.begin_hessians() + - c * dim * (dim + 1) / 2 * n_q_points); - } - } + private: + template + static void + work_normal(const MatrixFreeFunctions::UnivariateShapeData &data, + const Number *in, + Number *out, + const bool add_into_result = false) + { + AssertIndexRange(direction, dim); + constexpr int n_rows = fe_degree + 1; + constexpr int n_columns = n_q_points_1d; + constexpr int mm = contract_over_rows ? n_rows : n_columns; + constexpr int nn = contract_over_rows ? n_columns : n_rows; + const Number2 *shape_data = data.shape_values_eo.data(); + Assert(shape_data != nullptr, ExcNotInitialized()); + Assert(contract_over_rows == false || !add_into_result, + ExcMessage("Cannot add into result if contract_over_rows = true")); + + constexpr int n_blocks1 = Utilities::pow(fe_degree, direction); + constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1); + + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + if (contract_over_rows == false && add_into_result) + apply_matrix_vector_product(shape_data, in, out); + else + apply_matrix_vector_product(shape_data, in, out); + + ++in; + ++out; + } + in += n_blocks1 * (mm - 1); + out += n_blocks1 * (nn - 1); + } + } + + template + static void + work_tangential( + const MatrixFreeFunctions::UnivariateShapeData &data, + Number *ptr) + { + AssertIndexRange(direction, dim); + static_assert(direction != normal_direction, + "Cannot interpolate tangentially in normal direction"); + + constexpr int n_rows = fe_degree; + constexpr int n_columns = n_q_points_1d; + const Number2 *shape_data = data.shape_values_eo.data(); + Assert(shape_data != nullptr, ExcNotInitialized()); + + constexpr int n_blocks1 = + (direction > normal_direction) ? + Utilities::pow(n_q_points_1d, direction) : + (direction > 0 ? + (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) : + 1); + constexpr int n_blocks2 = + (direction > normal_direction) ? + Utilities::pow(fe_degree, dim - 1 - direction) : + ((direction + 1 < dim) ? + (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) : + 1); + + // Since we perform an in-place interpolation, we must run the step + // expanding the size of the basis backward ('contract_over_rows' aka + // 'evaluate' case). + if (contract_over_rows) + { + const Number *in = + ptr + (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1; + Number *out = + ptr + (n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1; + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + apply_matrix_vector_product(shape_data, in, out); + + --in; + --out; + } + in -= n_blocks1 * (n_rows - 1); + out -= n_blocks1 * (n_columns - 1); + } + } + else + { + const Number *in = ptr; + Number *out = ptr; + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + apply_matrix_vector_product(shape_data, in, out); + + ++in; + ++out; + } + in += n_blocks1 * (n_columns - 1); + out += n_blocks1 * (n_rows - 1); + } + } + } + }; template + template inline void - FEEvaluationImplTransformToCollocation< - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, - FEEvaluationData &fe_eval, - const bool add_into_values_array) + FEEvaluationImpl:: + evaluate_or_integrate( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool add_into_values_array) { - const auto &shape_data = fe_eval.get_shape_info().data.front(); + Assert(dim == 2 || dim == 3, + ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas " + "evaluation/integration")); - Assert(n_q_points_1d > fe_degree, - ExcMessage("You lose information when going to a collocation space " - "of lower degree, so the evaluation results would be " - "wrong. Thus, this class does not permit the desired " - "operation.")); - constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim); + if (evaluation_flag == EvaluationFlags::nothing) + return; - for (unsigned int c = 0; c < n_components; ++c) + AssertDimension(fe_eval.get_shape_info().data.size(), 2); + AssertDimension(n_q_points_1d, + fe_eval.get_shape_info().data[0].n_q_points_1d); + AssertDimension(n_q_points_1d, + fe_eval.get_shape_info().data[1].n_q_points_1d); + AssertDimension(fe_degree, fe_eval.get_shape_info().data[0].fe_degree); + AssertDimension(fe_degree, fe_eval.get_shape_info().data[1].fe_degree + 1); + + const auto &shape_data = fe_eval.get_shape_info().data; + const unsigned int dofs_per_component = + Utilities::pow(fe_degree, dim - 1) * (fe_degree + 1); + const unsigned int n_points = Utilities::pow(n_q_points_1d, dim); + Number *gradients = fe_eval.begin_gradients(); + Number *values = fe_eval.begin_values(); + + if (integrate) { - // apply derivatives in collocation space - if (integration_flag & - (EvaluationFlags::gradients | EvaluationFlags::hessians)) - FEEvaluationImplCollocation:: - do_integrate(shape_data, - integration_flag & (EvaluationFlags::gradients | - EvaluationFlags::hessians), - fe_eval.begin_values() + c * n_q_points, - fe_eval.begin_gradients() + c * dim * n_q_points, - fe_eval.begin_hessians() + - c * dim * (dim + 1) / 2 * n_q_points, - /*add_into_values_array=*/ - integration_flag & EvaluationFlags::values); - - // transform back to the original space - FEEvaluationImplBasisChange< - evaluate_evenodd, - EvaluatorQuantity::value, - dim, - (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), - n_q_points_1d>::do_backward(1, - shape_data.shape_values_eo, - add_into_values_array, - fe_eval.begin_values() + c * n_q_points, - values_dofs + - c * Utilities::pow(fe_degree + 1, dim)); + const bool do_values = evaluation_flag & EvaluationFlags::values; + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + integrate_gradients_collocation(shape_data[0], + values, + gradients, + do_values); + if constexpr (dim > 2) + work_tangential<2, 0, false>(shape_data[1], values); + work_tangential<1, 0, false>(shape_data[1], values); + work_normal<0, false>(shape_data[0], + values, + values_dofs, + add_into_values_array); + + values += n_points; + gradients += n_points * dim; + values_dofs += dofs_per_component; + + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + integrate_gradients_collocation(shape_data[0], + values, + gradients, + do_values); + if constexpr (dim > 2) + work_tangential<2, 1, false>(shape_data[1], values); + work_tangential<0, 1, false>(shape_data[1], values); + work_normal<1, false>(shape_data[0], + values, + values_dofs, + add_into_values_array); + + if constexpr (dim > 2) + { + values += n_points; + gradients += n_points * dim; + values_dofs += dofs_per_component; + + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + integrate_gradients_collocation(shape_data[0], + values, + gradients, + do_values); + work_tangential<1, 2, false>(shape_data[1], values); + work_tangential<0, 2, false>(shape_data[1], values); + work_normal<2, false>(shape_data[0], + values, + values_dofs, + add_into_values_array); + } + } + else + { + work_normal<0, true>(shape_data[0], values_dofs, values); + work_tangential<1, 0, true>(shape_data[1], values); + if constexpr (dim > 2) + work_tangential<2, 0, true>(shape_data[1], values); + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + evaluate_gradients_collocation(shape_data[0], + values, + gradients); + + values += n_points; + gradients += n_points * dim; + values_dofs += dofs_per_component; + + work_normal<1, true>(shape_data[0], values_dofs, values); + work_tangential<0, 1, true>(shape_data[1], values); + if constexpr (dim > 2) + work_tangential<2, 1, true>(shape_data[1], values); + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + evaluate_gradients_collocation(shape_data[0], + values, + gradients); + + if constexpr (dim > 2) + { + values += n_points; + gradients += n_points * dim; + values_dofs += dofs_per_component; + + work_normal<2, true>(shape_data[0], values_dofs, values); + work_tangential<0, 2, true>(shape_data[1], values); + work_tangential<1, 2, true>(shape_data[1], values); + if ((evaluation_flag & EvaluationFlags::gradients) != 0u) + evaluate_gradients_collocation(shape_data[0], + values, + gradients); + } } } @@ -2275,7 +2137,7 @@ namespace internal const EvaluationFlags::EvaluationFlags evaluation_flag, OtherNumber *values_dofs, FEEvaluationData &fe_eval, - const bool sum_into_values_array = false) + const bool sum_into_values_array_in = false) { // `OtherNumber` is either `const Number` (evaluate()) or `Number` // (integrate()) @@ -2291,13 +2153,39 @@ namespace internal element_type == ElementType::tensor_raviart_thomas, ExcNotImplemented()); + EvaluationFlags::EvaluationFlags actual_flag = evaluation_flag; + bool sum_into_values_array = sum_into_values_array_in; + if (evaluation_flag & EvaluationFlags::hessians) + { + actual_flag |= EvaluationFlags::values; + Assert(element_type != MatrixFreeFunctions::tensor_none, + ExcNotImplemented()); + if constexpr (do_integrate) + { + if (fe_eval.get_shape_info().data[0].fe_degree < + fe_eval.get_shape_info().data[0].n_q_points_1d) + integrate_hessians_collocation( + n_components, + fe_eval, + evaluation_flag & EvaluationFlags::values); + else + { + integrate_hessians_slow(n_components, + fe_eval, + values_dofs, + sum_into_values_array); + sum_into_values_array = true; + } + } + } + if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && element_type == ElementType::tensor_symmetric_collocation) { evaluate_or_integrate< FEEvaluationImplCollocation>( n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); @@ -2314,7 +2202,7 @@ namespace internal n_q_points_1d, Number>>( n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); @@ -2328,7 +2216,7 @@ namespace internal n_q_points_1d, Number>>( n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); @@ -2341,7 +2229,7 @@ namespace internal fe_degree, n_q_points_1d, Number>>(n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); @@ -2354,7 +2242,7 @@ namespace internal n_q_points_1d, Number>>( n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); @@ -2367,23 +2255,32 @@ namespace internal n_q_points_1d, Number>>( n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); } else if (element_type == ElementType::tensor_raviart_thomas) { - FEEvaluationImpl:: - template evaluate_or_integrate(evaluation_flag, - const_cast( - values_dofs), - fe_eval, - sum_into_values_array); + if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1) + { + FEEvaluationImpl:: + template evaluate_or_integrate( + actual_flag, + const_cast(values_dofs), + fe_eval, + sum_into_values_array); + } + else + { + Assert(false, + ExcNotImplemented("Raviart-Thomas currently only possible " + "in 2d/3d and with templated degree")); + } } else { @@ -2393,12 +2290,23 @@ namespace internal n_q_points_1d, Number>>( n_components, - evaluation_flag, + actual_flag, values_dofs, fe_eval, sum_into_values_array); } + if ((evaluation_flag & EvaluationFlags::hessians) && !do_integrate) + { + Assert(element_type != MatrixFreeFunctions::tensor_none, + ExcNotImplemented()); + if (fe_eval.get_shape_info().data[0].fe_degree < + fe_eval.get_shape_info().data[0].n_q_points_1d) + evaluate_hessians_collocation(n_components, fe_eval); + else + evaluate_hessians_slow(n_components, values_dofs, fe_eval); + } + return false; } @@ -2580,25 +2488,24 @@ namespace internal Number, Number2> eval_grad({}, data.shape_gradients_collocation_eo, {}); - eval_grad.template gradients<0, true, false>( + eval_grad.template gradients<0, true, false, 3>( values_quad, gradients_quad); - eval_grad.template gradients<1, true, false>( - values_quad, gradients_quad + n_q_points); + eval_grad.template gradients<1, true, false, 3>( + values_quad, gradients_quad + 1); } else { // grad x eval0.template gradients<0, true, false>(values_dofs, scratch_data); - eval1.template values<1, true, false>(scratch_data, - gradients_quad); + eval1.template values<1, true, false, 3>(scratch_data, + gradients_quad); // grad y eval0.template values<0, true, false>(values_dofs, scratch_data); - eval1.template gradients<1, true, false>(scratch_data, - gradients_quad + - n_q_points); + eval1.template gradients<1, true, false, 3>( + scratch_data, gradients_quad + 1); if ((evaluation_flag & EvaluationFlags::values) != 0u) eval1.template values<1, true, false>(scratch_data, @@ -2607,16 +2514,15 @@ namespace internal // grad z eval0.template values<0, true, false>(values_dofs + n_dofs, scratch_data); - eval1.template values<1, true, false>( - scratch_data, gradients_quad + (dim - 1) * n_q_points); + eval1.template values<1, true, false, 3>(scratch_data, + gradients_quad + 2); break; case 2: - eval0.template values<0, true, false>(values_dofs + n_dofs, - gradients_quad + - n_q_points); - eval0.template gradients<0, true, false>(values_dofs, - gradients_quad); + eval0.template values<0, true, false, 2>(values_dofs + n_dofs, + gradients_quad + 1); + eval0.template gradients<0, true, false, 2>(values_dofs, + gradients_quad); if ((evaluation_flag & EvaluationFlags::values) != 0u) eval0.template values<0, true, false>(values_dofs, values_quad); @@ -2766,12 +2672,9 @@ namespace internal { case 3: // grad z - eval1.template values<1, false, false>(gradients_quad + - 2 * n_q_points, - gradients_quad + - 2 * n_q_points); - eval0.template values<0, false, false>(gradients_quad + - 2 * n_q_points, + eval1.template values<1, false, false, 3>(gradients_quad + 2, + scratch_data); + eval0.template values<0, false, false>(scratch_data, values_dofs + n_dofs); if (symmetric_evaluate && use_collocation_evaluation(fe_degree, n_q_points_1d)) @@ -2784,12 +2687,12 @@ namespace internal Number2> eval_grad({}, data.shape_gradients_collocation_eo, {}); if ((integration_flag & EvaluationFlags::values) != 0u) - eval_grad.template gradients<1, false, true>( - gradients_quad + n_q_points, values_quad); + eval_grad.template gradients<1, false, true, 3>( + gradients_quad + 1, values_quad); else - eval_grad.template gradients<1, false, false>( - gradients_quad + n_q_points, values_quad); - eval_grad.template gradients<0, false, true>( + eval_grad.template gradients<1, false, false, 3>( + gradients_quad + 1, values_quad); + eval_grad.template gradients<0, false, true, 3>( gradients_quad, values_quad); eval0.template values<1, false, false>(values_quad, values_quad); @@ -2802,30 +2705,30 @@ namespace internal { eval1.template values<1, false, false>(values_quad, scratch_data); - eval1.template gradients<1, false, true>( - gradients_quad + n_q_points, scratch_data); + eval1.template gradients<1, false, true, 3>( + gradients_quad + 1, scratch_data); } else - eval1.template gradients<1, false, false>( - gradients_quad + n_q_points, scratch_data); + eval1.template gradients<1, false, false, 3>( + gradients_quad + 1, scratch_data); // grad y eval0.template values<0, false, false>(scratch_data, values_dofs); // grad x - eval1.template values<1, false, false>(gradients_quad, - scratch_data); + eval1.template values<1, false, false, 3>(gradients_quad, + scratch_data); eval0.template gradients<0, false, true>(scratch_data, values_dofs); } break; case 2: - eval0.template values<0, false, false>(gradients_quad + - n_q_points, - values_dofs + n_dofs); - eval0.template gradients<0, false, false>(gradients_quad, - values_dofs); + eval0.template values<0, false, false, 2>(gradients_quad + 1, + values_dofs + + n_dofs); + eval0.template gradients<0, false, false, 2>(gradients_quad, + values_dofs); if ((integration_flag & EvaluationFlags::values) != 0u) eval0.template values<0, false, true>(values_quad, values_dofs); @@ -2942,6 +2845,8 @@ namespace internal } }; + + template struct FEFaceEvaluationImplRaviartThomas { @@ -3559,44 +3464,71 @@ namespace internal { if (face_direction == face_no / 2) { - EvaluatorTensorProduct - evalf(shape_data[face_no % 2].begin(), - nullptr, - nullptr, - n_points_1d, - 0); - - const unsigned int in_stride = do_evaluate ? - dofs_per_component_on_cell : - dofs_per_component_on_face; - const unsigned int out_stride = do_evaluate ? - dofs_per_component_on_face : - dofs_per_component_on_cell; + constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction); + + const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d; + const int stride = Utilities::pow(n_rows, face_direction); + const std::array n_blocks{ + {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}; + std::array steps; + if constexpr (face_direction == 0) + steps = {{n_rows, 0}}; + else if constexpr (face_direction == 1 && dim == 2) + steps = {{1, 0}}; + else if constexpr (face_direction == 1) + // in 3d, the coordinate system is zx, not xz -> switch indices + steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}}; + else if constexpr (face_direction == 2) + steps = {{1, 0}}; for (unsigned int c = 0; c < n_components; ++c) { if (flag & EvaluationFlags::hessians) - evalf.template apply_face(input, output); + interpolate_to_face(shape_data[face_no % 2].begin(), + n_blocks, + steps, + input, + output, + n_rows, + stride); else if (flag & EvaluationFlags::gradients) - evalf.template apply_face(input, output); + interpolate_to_face(shape_data[face_no % 2].begin(), + n_blocks, + steps, + input, + output, + n_rows, + stride); else - evalf.template apply_face(input, output); - input += in_stride; - output += out_stride; + interpolate_to_face(shape_data[face_no % 2].begin(), + n_blocks, + steps, + input, + output, + n_rows, + stride); + if (do_evaluate) + { + input += dofs_per_component_on_cell; + output += dofs_per_component_on_face; + } + else + { + output += dofs_per_component_on_cell; + input += dofs_per_component_on_face; + } } } else if (face_direction < dim) @@ -3961,13 +3893,13 @@ namespace internal if (integrate) for (unsigned int q = 0; q < n_q_points; ++q) tmp_values[q] = - gradients_quad[(c * dim + d) * n_q_points + orientation[q]]; + gradients_quad[(c * n_q_points + orientation[q]) * dim + d]; else for (unsigned int q = 0; q < n_q_points; ++q) tmp_values[orientation[q]] = - gradients_quad[(c * dim + d) * n_q_points + q]; + gradients_quad[(c * n_q_points + q) * dim + d]; for (unsigned int q = 0; q < n_q_points; ++q) - gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q]; + gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q]; } if (flag & EvaluationFlags::hessians) { @@ -4026,14 +3958,15 @@ namespace internal Assert(gradients_quad != nullptr, ExcInternalError()); if (integrate) for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points + - orientation[q]][v]; + tmp_values[q] = + gradients_quad[(c * n_q_points + orientation[q]) * dim + d] + [v]; else for (unsigned int q = 0; q < n_q_points; ++q) tmp_values[orientation[q]] = - gradients_quad[(c * dim + d) * n_q_points + q][v]; + gradients_quad[(c * n_q_points + q) * dim + d][v]; for (unsigned int q = 0; q < n_q_points; ++q) - gradients_quad[(c * dim + d) * n_q_points + q][v] = + gradients_quad[(c * n_q_points + q) * dim + d][v] = tmp_values[q]; } if (flag & EvaluationFlags::hessians) @@ -4129,11 +4062,10 @@ namespace internal n_dofs, n_q_points); - eval.template gradients<0, true, false>( - values_dofs_actual_ptr, gradients_quad_ptr); - - gradients_quad_ptr += 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; } } @@ -4370,14 +4302,13 @@ namespace internal if (!(integration_flag & EvaluationFlags::values) && d == 0) - eval.template gradients<0, false, false>( - gradients_quad_ptr, values_dofs_actual_ptr); + eval.template gradients<0, false, false, dim>( + gradients_quad_ptr + d, values_dofs_actual_ptr); else - eval.template gradients<0, false, true>( - gradients_quad_ptr, values_dofs_actual_ptr); - - gradients_quad_ptr += n_q_points; + 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; } } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index d1acf28d89..f5288fc65d 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -4526,10 +4526,10 @@ inline DEAL_II_ALWAYS_INLINE // Cartesian cell if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) grad_out[comp][d] = - this->gradients_quad[(comp * dim + d) * nqp + q_point] * + this->gradients_quad[(comp * nqp + q_point) * dim + d] * this->jacobian[0][d][d]; } // cell with general/affine Jacobian @@ -4543,11 +4543,11 @@ inline DEAL_II_ALWAYS_INLINE for (unsigned int d = 0; d < dim; ++d) { grad_out[comp][d] = - jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point]; + jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim]; for (unsigned int e = 1; e < dim; ++e) grad_out[comp][d] += jac[d][e] * - this->gradients_quad[(comp * dim + e) * nqp + q_point]; + this->gradients_quad[(comp * nqp + q_point) * dim + e]; } } return grad_out; @@ -4580,7 +4580,7 @@ FEEvaluationBase:: if (this->cell_type == internal::MatrixFreeFunctions::cartesian) for (unsigned int comp = 0; comp < n_components; ++comp) grad_out[comp] = - this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] * + this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] * (this->normal_x_jacobian[0][dim - 1]); else { @@ -4588,11 +4588,11 @@ FEEvaluationBase:: this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point; for (unsigned int comp = 0; comp < n_components; ++comp) { - grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] * + grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] * this->normal_x_jacobian[index][0]; for (unsigned int d = 1; d < dim; ++d) grad_out[comp] += - this->gradients_quad[(comp * dim + d) * nqp + q_point] * + this->gradients_quad[(comp * nqp + q_point) * dim + d] * this->normal_x_jacobian[index][d]; } } @@ -4770,7 +4770,7 @@ FEEvaluationBase:: for (unsigned int e = 0; e < dim; ++e) hessian_out[comp][d][d] += jac_grad[d][e] * - this->gradients_quad[(comp * dim + e) * nqp + q_point]; + this->gradients_quad[(comp * nqp + q_point) * dim + e]; // add off-diagonal part of J' * grad(u) for (unsigned int d = 0, count = dim; d < dim; ++d) @@ -4778,7 +4778,7 @@ FEEvaluationBase:: for (unsigned int f = 0; f < dim; ++f) hessian_out[comp][d][e] += jac_grad[count][f] * - this->gradients_quad[(comp * dim + f) * nqp + q_point]; + this->gradients_quad[(comp * nqp + q_point) * dim + f]; // take symmetric part for (unsigned int d = 0; d < dim; ++d) @@ -4872,7 +4872,7 @@ FEEvaluationBase:: for (unsigned int e = 0; e < dim; ++e) hessian_out[comp][d] += jac_grad[d][e] * - this->gradients_quad[(comp * dim + e) * nqp + q_point]; + this->gradients_quad[(comp * nqp + q_point) * dim + e]; } } return hessian_out; @@ -5006,7 +5006,7 @@ FEEvaluationBase:: { const VectorizedArrayType factor = jac[d] * JxW; for (unsigned int comp = 0; comp < n_components; ++comp) - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = grad_in[comp][d] * factor; } } @@ -5026,7 +5026,7 @@ FEEvaluationBase:: VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0]; for (unsigned int e = 1; e < dim; ++e) new_val += (jac[e][d] * grad_in[comp][e]); - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = new_val * JxW; } } @@ -5062,9 +5062,9 @@ FEEvaluationBase:: for (unsigned int comp = 0; comp < n_components; ++comp) { for (unsigned int d = 0; d < dim - 1; ++d) - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = VectorizedArrayType(); - this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] = grad_in[comp] * JxW_jac; } } @@ -5081,7 +5081,7 @@ FEEvaluationBase:: for (unsigned int comp = 0; comp < n_components; ++comp) { for (unsigned int d = 0; d < dim; ++d) - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = (grad_in[comp] * JxW) * jac[d]; } } @@ -5233,8 +5233,8 @@ FEEvaluationBase:: for (unsigned int f = e + 1; f < dim; ++f, ++count) sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) * jac_grad[count][d]; - this->gradients_from_hessians_quad[(comp * dim + d) * nqp + - q_point] = sum * JxW; + this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim + + d] = sum * JxW; } } } @@ -5495,7 +5495,7 @@ FEEvaluationAccess:: template inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> FEEvaluationAccess::get_gradient( - const unsigned int q_point) const + const unsigned int q_point_in) const { // could use the base class gradient, but that involves too many expensive // initialization operations on tensors @@ -5504,7 +5504,7 @@ FEEvaluationAccess::get_gradient( Assert(this->gradients_quad_initialized == true, internal::ExcAccessToUninitializedField()); # endif - AssertIndexRange(q_point, this->n_quadrature_points); + AssertIndexRange(q_point_in, this->n_quadrature_points); Assert(this->jacobian != nullptr, internal::ExcMatrixFreeAccessToUninitializedMappingField( @@ -5512,12 +5512,12 @@ FEEvaluationAccess::get_gradient( Tensor<1, dim, VectorizedArrayType> grad_out; - const std::size_t nqp = this->n_quadrature_points; + const std::size_t q_point = q_point_in; if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { for (unsigned int d = 0; d < dim; ++d) grad_out[d] = - this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d]; + this->gradients_quad[dim * q_point + d] * this->jacobian[0][d][d]; } // cell with general/affine Jacobian else @@ -5528,9 +5528,9 @@ FEEvaluationAccess::get_gradient( 0]; for (unsigned int d = 0; d < dim; ++d) { - grad_out[d] = jac[d][0] * this->gradients_quad[q_point]; + grad_out[d] = jac[d][0] * this->gradients_quad[dim * q_point]; for (unsigned int e = 1; e < dim; ++e) - grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point]; + grad_out[d] += jac[d][e] * this->gradients_quad[dim * q_point + e]; } } return grad_out; @@ -5641,12 +5641,12 @@ template inline DEAL_II_ALWAYS_INLINE void FEEvaluationAccess:: submit_gradient(const Tensor<1, dim, VectorizedArrayType> grad_in, - const unsigned int q_point) + const unsigned int q_point_in) { # ifdef DEBUG Assert(this->is_reinitialized, ExcNotInitialized()); # endif - AssertIndexRange(q_point, this->n_quadrature_points); + AssertIndexRange(q_point_in, this->n_quadrature_points); Assert(this->J_value != nullptr, internal::ExcMatrixFreeAccessToUninitializedMappingField( "update_gradients")); @@ -5657,7 +5657,8 @@ FEEvaluationAccess:: this->gradients_quad_submitted = true; # endif - const std::size_t nqp = this->n_quadrature_points; + const std::size_t q_point = q_point_in; + VectorizedArrayType *grad_ptr = this->gradients_quad + dim * q_point; if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArrayType JxW = @@ -5670,7 +5671,7 @@ FEEvaluationAccess:: jac[d] = this->jacobian[0][d][d]; for (unsigned int d = 0; d < dim; ++d) - this->gradients_quad[d * nqp + q_point] = grad_in[d] * jac[d] * JxW; + grad_ptr[d] = grad_in[d] * jac[d] * JxW; } // general/affine cell type else @@ -5688,7 +5689,7 @@ FEEvaluationAccess:: VectorizedArrayType new_val = jac[0][d] * grad_in[0]; for (unsigned int e = 1; e < dim; ++e) new_val += jac[e][d] * grad_in[e]; - this->gradients_quad[d * nqp + q_point] = new_val * JxW; + grad_ptr[d] = new_val * JxW; } } } @@ -5902,7 +5903,7 @@ FEEvaluationAccess:: for (unsigned int d = 0; d < dim; ++d) for (unsigned int comp = 0; comp < n_components; ++comp) grad_out[comp][d] = - this->gradients_quad[(comp * dim + d) * nqp + q_point] * + this->gradients_quad[(comp * nqp + q_point) * dim + d] * inv_t_jac[d][d] * (jac[comp][comp] * inv_det); } else if (this->cell_type <= internal::MatrixFreeFunctions::affine) @@ -5927,7 +5928,7 @@ FEEvaluationAccess:: for (unsigned int f = 0; f < dim; ++f) for (unsigned int e = 0; e < dim; ++e) tmp += jac[comp][f] * inv_t_jac[d][e] * - this->gradients_quad[(f * dim + e) * nqp + q_point]; + this->gradients_quad[(f * nqp + q_point) * dim + e]; grad_out[comp][d] = tmp * inv_det; } @@ -5964,7 +5965,7 @@ FEEvaluationAccess:: for (unsigned int f = 0; f < dim; ++f) for (unsigned int e = 0; e < dim; ++e) tmp += t_jac[f][comp] * inv_t_jac[d][e] * - this->gradients_quad[(f * dim + e) * nqp + q_point]; + this->gradients_quad[(f * nqp + q_point) * dim + e]; grad_out[comp][d] = tmp * inv_det; } @@ -6069,9 +6070,9 @@ FEEvaluationAccess:: this->jacobian[0][2][2]; // div * det(J^-1) - divergence = this->gradients_quad[q_point]; + divergence = this->gradients_quad[q_point * dim]; for (unsigned int d = 1; d < dim; ++d) - divergence += this->gradients_quad[(dim * d + d) * nqp + q_point]; + divergence += this->gradients_quad[(d * nqp + q_point) * dim + d]; divergence *= inv_det; } else @@ -6087,9 +6088,9 @@ FEEvaluationAccess:: Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1); // div * det(J^-1) - divergence = this->gradients_quad[q_point]; + divergence = this->gradients_quad[q_point * dim]; for (unsigned int d = 1; d < dim; ++d) - divergence += this->gradients_quad[(dim * d + d) * nqp + q_point]; + divergence += this->gradients_quad[(d * nqp + q_point) * dim + d]; divergence *= inv_det; } } @@ -6099,9 +6100,10 @@ FEEvaluationAccess:: this->cell_type == internal::MatrixFreeFunctions::cartesian) { // Cartesian cell - divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0]; + divergence = + this->gradients_quad[q_point * dim] * this->jacobian[0][0][0]; for (unsigned int d = 1; d < dim; ++d) - divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] * + divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] * this->jacobian[0][d][d]; } else @@ -6111,13 +6113,13 @@ FEEvaluationAccess:: this->cell_type == internal::MatrixFreeFunctions::general ? this->jacobian[q_point] : this->jacobian[0]; - divergence = jac[0][0] * this->gradients_quad[q_point]; + divergence = jac[0][0] * this->gradients_quad[q_point * dim]; for (unsigned int e = 1; e < dim; ++e) - divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point]; + divergence += jac[0][e] * this->gradients_quad[q_point * dim + e]; for (unsigned int d = 1; d < dim; ++d) for (unsigned int e = 0; e < dim; ++e) divergence += - jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point]; + jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e]; } } return divergence; @@ -6325,7 +6327,7 @@ FEEvaluationAccess:: const VectorizedArrayType weight = this->quadrature_weights[q_point]; for (unsigned int d = 0; d < dim; ++d) for (unsigned int comp = 0; comp < n_components; ++comp) - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight; } else if (this->cell_type <= internal::MatrixFreeFunctions::affine) @@ -6353,7 +6355,7 @@ FEEvaluationAccess:: for (unsigned int e = 0; e < dim; ++e) tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e]; - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = tmp * fac; } } @@ -6385,7 +6387,7 @@ FEEvaluationAccess:: for (unsigned int e = 0; e < dim; ++e) tmp += t_jac[comp][f] * inv_t_jac[e][d] * grad_in[f][e]; - this->gradients_quad[(comp * dim + d) * nqp + q_point] = + this->gradients_quad[(comp * nqp + q_point) * dim + d] = tmp * fac; } @@ -6531,12 +6533,12 @@ FEEvaluationAccess:: for (unsigned int d = 0; d < dim; ++d) { - this->gradients_quad[(dim * d + d) * nqp + q_point] = fac; + this->gradients_quad[(d * nqp + q_point) * dim + d] = fac; for (unsigned int e = d + 1; e < dim; ++e) { - this->gradients_quad[(dim * d + e) * nqp + q_point] = + this->gradients_quad[(d * nqp + q_point) * dim + e] = VectorizedArrayType(); - this->gradients_quad[(dim * e + d) * nqp + q_point] = + this->gradients_quad[(e * nqp + q_point) * dim + d] = VectorizedArrayType(); } } @@ -6551,13 +6553,13 @@ FEEvaluationAccess:: this->J_value[0] * this->quadrature_weights[q_point] * div_in; for (unsigned int d = 0; d < dim; ++d) { - this->gradients_quad[(d * dim + d) * nqp + q_point] = + this->gradients_quad[(d * nqp + q_point) * dim + d] = (fac * this->jacobian[0][d][d]); for (unsigned int e = d + 1; e < dim; ++e) { - this->gradients_quad[(d * dim + e) * nqp + q_point] = + this->gradients_quad[(d * nqp + q_point) * dim + e] = VectorizedArrayType(); - this->gradients_quad[(e * dim + d) * nqp + q_point] = + this->gradients_quad[(e * nqp + q_point) * dim + d] = VectorizedArrayType(); } } @@ -6576,7 +6578,7 @@ FEEvaluationAccess:: for (unsigned int d = 0; d < dim; ++d) { for (unsigned int e = 0; e < dim; ++e) - this->gradients_quad[(d * dim + e) * nqp + q_point] = + this->gradients_quad[(d * nqp + q_point) * dim + e] = jac[d][e] * fac; } } @@ -6620,16 +6622,16 @@ FEEvaluationAccess:: const VectorizedArrayType JxW = this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int d = 0; d < dim; ++d) - this->gradients_quad[(d * dim + d) * nqp + q_point] = + this->gradients_quad[(d * nqp + q_point) * dim + d] = (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]); for (unsigned int e = 0, counter = dim; e < dim; ++e) for (unsigned int d = e + 1; d < dim; ++d, ++counter) { const VectorizedArrayType value = sym_grad.access_raw_entry(counter) * JxW; - this->gradients_quad[(e * dim + d) * nqp + q_point] = + this->gradients_quad[(e * nqp + q_point) * dim + d] = value * this->jacobian[0][d][d]; - this->gradients_quad[(d * dim + e) * nqp + q_point] = + this->gradients_quad[(d * nqp + q_point) * dim + e] = value * this->jacobian[0][e][e]; } } @@ -6661,7 +6663,7 @@ FEEvaluationAccess:: VectorizedArrayType new_val = jac[0][d] * weighted[comp][0]; for (unsigned int e = 1; e < dim; ++e) new_val += jac[e][d] * weighted[comp][e]; - this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val; + this->gradients_quad[(comp * nqp + q_point) * dim + d] = new_val; } } } diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index dfbd342a04..5834f41fb0 100644 --- 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) diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 91a42496c4..11aa7fefc7 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -682,13 +682,13 @@ namespace internal inline DEAL_II_ALWAYS_INLINE #endif std::enable_if_t<(variant == evaluate_evenodd), void> - apply_matrix_vector_product(const Number2 *matrix, - const Number *in, - Number *out, - int n_rows_runtime = 0, - int n_columns_runtime = 0, - int stride_in_runtime = 0, - int stride_out_runtime = 0) + apply_matrix_vector_product(const Number2 *DEAL_II_RESTRICT matrix, + const Number *in, + Number *out, + int n_rows_runtime = 0, + int n_columns_runtime = 0, + int stride_in_runtime = 0, + int stride_out_runtime = 0) { const int n_rows = n_rows_static == 0 ? n_rows_runtime : n_rows_static; const int n_columns = @@ -1165,25 +1165,23 @@ namespace internal (void)dummy2; } - template + template void values(const Number in[], Number out[]) const { - apply(shape_values, in, out); + constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value; + apply( + shape_values, in, out); } - template + template void gradients(const Number in[], Number out[]) const { constexpr EvaluatorQuantity gradient_type = (variant == evaluate_general ? EvaluatorQuantity::value : EvaluatorQuantity::gradient); - apply( + apply( shape_gradients, in, out); } @@ -1263,7 +1261,8 @@ namespace internal bool contract_over_rows, bool add, bool one_line = false, - EvaluatorQuantity = EvaluatorQuantity::value> + EvaluatorQuantity = EvaluatorQuantity::value, + int extra_stride = 1> static void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, @@ -1324,7 +1323,8 @@ namespace internal bool contract_over_rows, bool add, bool one_line, - EvaluatorQuantity quantity> + EvaluatorQuantity quantity, + int extra_stride> inline void EvaluatorTensorProduct:: apply(const Number2 *DEAL_II_RESTRICT shape_data, @@ -1349,6 +1349,8 @@ namespace internal constexpr int n_blocks2 = Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1)); + constexpr int stride_in = !contract_over_rows ? extra_stride : 1; + constexpr int stride_out = contract_over_rows ? extra_stride : 1; for (int i2 = 0; i2 < n_blocks2; ++i2) { for (int i1 = 0; i1 < n_blocks1; ++i1) @@ -1357,22 +1359,146 @@ namespace internal quantity, n_rows, n_columns, - stride, - stride, + stride * stride_in, + stride * stride_out, contract_over_rows, add>(shape_data, in, out); if (one_line == false) { - ++in; - ++out; + in += stride_in; + out += stride_out; } } if (one_line == false) { - in += stride * (mm - 1); - out += stride * (nn - 1); + in += stride * (mm - 1) * stride_in; + out += stride * (nn - 1) * stride_out; + } + } + } + + + + template + inline std::enable_if_t + interpolate_to_face(const Number2 *shape_values, + const std::array &n_blocks, + const std::array &steps, + const Number *input, + Number *DEAL_II_RESTRICT output, + const int n_rows_runtime = 0, + const int stride_runtime = 1) + { + const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime; + const int stride = n_rows_template > 0 ? stride_template : stride_runtime; + + Number *output1 = output + n_blocks[0] * n_blocks[1]; + Number *output2 = output1 + n_blocks[0] * n_blocks[1]; + for (int i2 = 0; i2 < n_blocks[1]; ++i2) + { + for (int i1 = 0; i1 < n_blocks[0]; ++i1) + { + Number res0 = shape_values[0] * input[0]; + Number res1, res2; + if (max_derivative > 0) + res1 = shape_values[n_rows] * input[0]; + if (max_derivative > 1) + res2 = shape_values[2 * n_rows] * input[0]; + for (int ind = 1; ind < n_rows; ++ind) + { + res0 += shape_values[ind] * input[stride * ind]; + if (max_derivative > 0) + res1 += shape_values[ind + n_rows] * input[stride * ind]; + if (max_derivative > 1) + res2 += shape_values[ind + 2 * n_rows] * input[stride * ind]; + } + if (add) + { + output[i1] += res0; + if (max_derivative > 0) + output1[i1] += res1; + if (max_derivative > 1) + output2[i2] += res2; + } + else + { + output[i1] = res0; + if (max_derivative > 0) + output1[i1] = res1; + if (max_derivative > 1) + output2[i1] = res2; + } + input += steps[0]; } + output += n_blocks[0]; + if (max_derivative > 0) + output1 += n_blocks[0]; + if (max_derivative > 1) + output2 += n_blocks[0]; + input += steps[1]; + } + } + + + + template + inline std::enable_if_t + interpolate_to_face(const Number2 *shape_values, + const std::array &n_blocks, + const std::array &steps, + const Number *input, + Number *DEAL_II_RESTRICT output, + const int n_rows_runtime = 0, + const int stride_runtime = 1) + { + const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime; + const int stride = n_rows_template > 0 ? stride_template : stride_runtime; + + const Number *input1 = input + n_blocks[0] * n_blocks[1]; + const Number *input2 = input1 + n_blocks[0] * n_blocks[1]; + for (int i2 = 0; i2 < n_blocks[1]; ++i2) + { + for (int i1 = 0; i1 < n_blocks[0]; ++i1) + { + const Number in = input[i1]; + Number in1, in2; + if (max_derivative > 0) + in1 = input1[i1]; + if (max_derivative > 1) + in2 = input2[i1]; + for (int col = 0; col < n_rows; ++col) + { + Number result = + add ? (output[col * stride] + shape_values[col] * in) : + (shape_values[col] * in); + if (max_derivative > 0) + result += shape_values[col + n_rows] * in1; + if (max_derivative > 1) + result += shape_values[col + 2 * n_rows] * in2; + + output[col * stride] = result; + } + output += steps[0]; + } + input += n_blocks[0]; + if (max_derivative > 0) + input1 += n_blocks[0]; + if (max_derivative > 1) + input2 += n_blocks[0]; + output += steps[1]; } } @@ -1385,7 +1511,7 @@ namespace internal typename Number, typename Number2> template inline void @@ -1400,116 +1526,24 @@ namespace internal ExcMessage( "The given array shape_values must not be the null pointer.")); - constexpr int n_blocks1 = (dim > 1 ? n_rows : 1); - constexpr int n_blocks2 = (dim > 2 ? n_rows : 1); - - AssertIndexRange(face_direction, dim); - constexpr int in_stride = Utilities::pow(n_rows, face_direction); - constexpr int out_stride = Utilities::pow(n_rows, dim - 1); - const Number2 *DEAL_II_RESTRICT shape_values = this->shape_values; - - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - if (contract_onto_face == true) - { - Number res0 = shape_values[0] * in[0]; - Number res1, res2; - if (max_derivative > 0) - res1 = shape_values[n_rows] * in[0]; - if (max_derivative > 1) - res2 = shape_values[2 * n_rows] * in[0]; - for (int ind = 1; ind < n_rows; ++ind) - { - res0 += shape_values[ind] * in[in_stride * ind]; - if (max_derivative > 0) - res1 += shape_values[ind + n_rows] * in[in_stride * ind]; - if (max_derivative > 1) - res2 += - shape_values[ind + 2 * n_rows] * in[in_stride * ind]; - } - if (add) - { - out[0] += res0; - if (max_derivative > 0) - out[out_stride] += res1; - if (max_derivative > 1) - out[2 * out_stride] += res2; - } - else - { - out[0] = res0; - if (max_derivative > 0) - out[out_stride] = res1; - if (max_derivative > 1) - out[2 * out_stride] = res2; - } - } - else - { - for (int col = 0; col < n_rows; ++col) - { - if (add) - out[col * in_stride] += shape_values[col] * in[0]; - else - out[col * in_stride] = shape_values[col] * in[0]; - if (max_derivative > 0) - out[col * in_stride] += - shape_values[col + n_rows] * in[out_stride]; - if (max_derivative > 1) - out[col * in_stride] += - shape_values[col + 2 * n_rows] * in[2 * out_stride]; - } - } - - // increment: in regular case, just go to the next point in - // x-direction. If we are at the end of one chunk in x-dir, need - // to jump over to the next layer in z-direction - switch (face_direction) - { - case 0: - in += contract_onto_face ? n_rows : 1; - out += contract_onto_face ? 1 : n_rows; - break; - case 1: - ++in; - ++out; - // faces 2 and 3 in 3d use local coordinate system zx, which - // is the other way around compared to the tensor - // product. Need to take that into account. - if (dim == 3) - { - if (contract_onto_face) - out += n_rows - 1; - else - in += n_rows - 1; - } - break; - case 2: - ++in; - ++out; - break; - default: - Assert(false, ExcNotImplemented()); - } - } - - // adjust for local coordinate system zx - if (face_direction == 1 && dim == 3) - { - if (contract_onto_face) - { - in += n_rows * (n_rows - 1); - out -= n_rows * n_rows - 1; - } - else - { - out += n_rows * (n_rows - 1); - in -= n_rows * n_rows - 1; - } - } - } + constexpr int stride = Utilities::pow(n_rows, face_direction); + std::array steps; + if constexpr (face_direction == 0) + steps = {{n_rows, 0}}; + else if constexpr (face_direction == 1 && dim == 2) + steps = {{1, 0}}; + else if constexpr (face_direction == 1) + // in 3d, the coordinate system is zx, not xz -> switch indices + steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}}; + else if constexpr (face_direction == 2) + steps = {{1, 0}}; + + interpolate_to_face( + this->shape_values, + {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}, + steps, + in, + out); } @@ -1612,25 +1646,23 @@ namespace internal , n_columns(n_columns) {} - template + template void values(const Number *in, Number *out) const { - apply(shape_values, in, out); + constexpr EvaluatorQuantity value_type = EvaluatorQuantity::value; + apply( + shape_values, in, out); } - template + template void gradients(const Number *in, Number *out) const { constexpr EvaluatorQuantity gradient_type = (variant != evaluate_evenodd ? EvaluatorQuantity::value : EvaluatorQuantity::gradient); - apply( + apply( shape_gradients, in, out); } @@ -1681,8 +1713,9 @@ namespace internal template + bool one_line = false, + EvaluatorQuantity quantity = EvaluatorQuantity::value, + int extra_stride = 1> void apply(const Number2 *DEAL_II_RESTRICT shape_data, const Number *in, @@ -1713,7 +1746,8 @@ namespace internal bool contract_over_rows, bool add, bool one_line, - EvaluatorQuantity quantity> + EvaluatorQuantity quantity, + int extra_stride> inline void EvaluatorTensorProduct::apply( const Number2 *DEAL_II_RESTRICT shape_data, @@ -1741,6 +1775,8 @@ namespace internal Utilities::fixed_power(n_rows); Assert(n_rows <= 128, ExcNotImplemented()); + constexpr int stride_in = !contract_over_rows ? extra_stride : 1; + constexpr int stride_out = contract_over_rows ? extra_stride : 1; for (int i2 = 0; i2 < n_blocks2; ++i2) { for (int i1 = 0; i1 < n_blocks1; ++i1) @@ -1752,19 +1788,24 @@ namespace internal apply_matrix_vector_product( - shape_data, in, out, n_rows, n_columns, stride, stride); + add>(shape_data, + in, + out, + n_rows, + n_columns, + stride * stride_in, + stride * stride_out); if (one_line == false) { - ++in; - ++out; + in += stride_in; + out += stride_out; } } if (one_line == false) { - in += stride * (mm - 1); - out += stride * (nn - 1); + in += stride * (mm - 1) * stride_in; + out += stride * (nn - 1) * stride_out; } } } @@ -1776,7 +1817,7 @@ namespace internal typename Number, typename Number2> template inline void @@ -1788,116 +1829,28 @@ namespace internal ExcMessage( "The given array shape_data must not be the null pointer!")); static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported"); - const int n_blocks1 = dim > 1 ? n_rows : 1; - const int n_blocks2 = dim > 2 ? n_rows : 1; - - AssertIndexRange(face_direction, dim); - const int in_stride = - face_direction > 0 ? Utilities::fixed_power(n_rows) : 1; - const int out_stride = - dim > 1 ? Utilities::fixed_power(n_rows) : 1; - - for (int i2 = 0; i2 < n_blocks2; ++i2) - { - for (int i1 = 0; i1 < n_blocks1; ++i1) - { - if (contract_onto_face == true) - { - Number res0 = shape_values[0] * in[0]; - Number res1, res2; - if (max_derivative > 0) - res1 = shape_values[n_rows] * in[0]; - if (max_derivative > 1) - res2 = shape_values[2 * n_rows] * in[0]; - for (unsigned int ind = 1; ind < n_rows; ++ind) - { - res0 += shape_values[ind] * in[in_stride * ind]; - if (max_derivative > 0) - res1 += shape_values[ind + n_rows] * in[in_stride * ind]; - if (max_derivative > 1) - res2 += - shape_values[ind + 2 * n_rows] * in[in_stride * ind]; - } - if (add) - { - out[0] += res0; - if (max_derivative > 0) - out[out_stride] += res1; - if (max_derivative > 1) - out[2 * out_stride] += res2; - } - else - { - out[0] = res0; - if (max_derivative > 0) - out[out_stride] = res1; - if (max_derivative > 1) - out[2 * out_stride] = res2; - } - } - else - { - for (unsigned int col = 0; col < n_rows; ++col) - { - if (add) - out[col * in_stride] += shape_values[col] * in[0]; - else - out[col * in_stride] = shape_values[col] * in[0]; - if (max_derivative > 0) - out[col * in_stride] += - shape_values[col + n_rows] * in[out_stride]; - if (max_derivative > 1) - out[col * in_stride] += - shape_values[col + 2 * n_rows] * in[2 * out_stride]; - } - } - // increment: in regular case, just go to the next point in - // x-direction. If we are at the end of one chunk in x-dir, need - // to jump over to the next layer in z-direction - switch (face_direction) - { - case 0: - in += contract_onto_face ? n_rows : 1; - out += contract_onto_face ? 1 : n_rows; - break; - case 1: - ++in; - ++out; - // faces 2 and 3 in 3d use local coordinate system zx, which - // is the other way around compared to the tensor - // product. Need to take that into account. - if (dim == 3) - { - if (contract_onto_face) - out += n_rows - 1; - else - in += n_rows - 1; - } - break; - case 2: - ++in; - ++out; - break; - default: - Assert(false, ExcNotImplemented()); - } - } - if (face_direction == 1 && dim == 3) - { - // adjust for local coordinate system zx - if (contract_onto_face) - { - in += n_rows * (n_rows - 1); - out -= n_rows * n_rows - 1; - } - else - { - out += n_rows * (n_rows - 1); - in -= n_rows * n_rows - 1; - } - } - } + const int stride = Utilities::pow(n_rows, face_direction); + const int n_rows = this->n_rows; + std::array steps; + if constexpr (face_direction == 0) + steps = {{n_rows, 0}}; + else if constexpr (face_direction == 1 && dim == 2) + steps = {{1, 0}}; + else if constexpr (face_direction == 1) + // in 3d, the coordinate system is zx, not xz -> switch indices + steps = {{n_rows * n_rows, -n_rows * n_rows * n_rows + 1}}; + else if constexpr (face_direction == 2) + steps = {{1, 0}}; + + interpolate_to_face<0, 0, contract_to_face, add, max_derivative>( + this->shape_values, + {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}}, + steps, + in, + out, + n_rows, + stride); } @@ -2141,35 +2094,14 @@ namespace internal { for (int i1 = 0; i1 < n_blocks1; ++i1) { - Number x[mm]; - for (int i = 0; i < mm; ++i) - x[i] = in[stride * i]; - - for (int col = 0; col < nn; ++col) - { - Number2 val0; - - if (contract_over_rows) - val0 = shape_data[col]; - else - val0 = shape_data[col * n_columns]; - - Number res0 = val0 * x[0]; - for (int i = 1; i < mm; ++i) - { - if (contract_over_rows) - val0 = shape_data[i * n_columns + col]; - else - val0 = shape_data[col * n_columns + i]; - - res0 += val0 * x[i]; - } - if (add) - out[stride * col] += res0; - - else - out[stride * col] = res0; - } + apply_matrix_vector_product(shape_data, in, out); if (one_line == false) { diff --git a/tests/matrix_free/matrix_vector_hessians_cells.cc b/tests/matrix_free/matrix_vector_hessians_cells.cc index 162c917802..5736770329 100644 --- a/tests/matrix_free/matrix_vector_hessians_cells.cc +++ b/tests/matrix_free/matrix_vector_hessians_cells.cc @@ -218,13 +218,13 @@ test_hessians(const dealii::FE_Poly &fe, // compare solutions of matrix vector product { - dst2 -= dst; + dst -= dst2; double error = 0.; - if (dst.l2_norm() > 0) - error = dst2.l2_norm() / dst.l2_norm(); + if (dst2.l2_norm() > 0) + error = dst.l2_norm() / dst2.l2_norm(); else - error = dst2.l2_norm(); + error = dst.l2_norm(); if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) deallog << "FEValues verification: " << error << std::endl << std::endl;