From: Maximilian Bergbauer Date: Fri, 3 Nov 2023 13:39:54 +0000 (+0100) Subject: FCL path X-Git-Tag: relicensing~108^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a7418f049dff9421f70e086838dabcda645d6424;p=dealii.git FCL path --- diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h index b1cf07adce..bc1fcd49d2 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_face.h +++ b/include/deal.II/matrix_free/evaluation_kernels_face.h @@ -1447,107 +1447,102 @@ namespace internal template struct FEFaceEvaluationImplEvaluateSelector { - template static bool - run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number *values_dofs, - FEEvaluationData &fe_eval) + evaluate_tensor_none(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) { const auto &shape_info = fe_eval.get_shape_info(); const auto &shape_data = shape_info.data.front(); using Number2 = typename FEEvaluationData::shape_info_number_type; - if (shape_info.element_type == MatrixFreeFunctions::tensor_none) - { - Assert((fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) == false, - ExcNotImplemented()); + Assert((fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) == false, + ExcNotImplemented()); - const unsigned int face_no = fe_eval.get_face_no(); - const unsigned int face_orientation = fe_eval.get_face_orientation(); - const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; - const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; + const unsigned int face_no = fe_eval.get_face_no(); + const unsigned int face_orientation = fe_eval.get_face_orientation(); + const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; + const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; - using Eval = - EvaluatorTensorProduct; + using Eval = + EvaluatorTensorProduct; - if (evaluation_flag & EvaluationFlags::values) - { - const auto *const shape_values = - &shape_data.shape_values_face(face_no, face_orientation, 0); + if (evaluation_flag & EvaluationFlags::values) + { + const auto *const shape_values = + &shape_data.shape_values_face(face_no, face_orientation, 0); - auto *values_quad_ptr = fe_eval.begin_values(); - auto *values_dofs_actual_ptr = values_dofs; + auto *values_quad_ptr = fe_eval.begin_values(); + auto *values_dofs_actual_ptr = values_dofs; - Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); - for (unsigned int c = 0; c < n_components; ++c) - { - eval.template values<0, true, false>(values_dofs_actual_ptr, - values_quad_ptr); + Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); + for (unsigned int c = 0; c < n_components; ++c) + { + eval.template values<0, true, false>(values_dofs_actual_ptr, + values_quad_ptr); - values_quad_ptr += n_q_points; - values_dofs_actual_ptr += n_dofs; - } + values_quad_ptr += n_q_points; + values_dofs_actual_ptr += n_dofs; } + } - if (evaluation_flag & EvaluationFlags::gradients) - { - auto *gradients_quad_ptr = fe_eval.begin_gradients(); - const auto *values_dofs_actual_ptr = values_dofs; + if (evaluation_flag & EvaluationFlags::gradients) + { + auto *gradients_quad_ptr = fe_eval.begin_gradients(); + const auto *values_dofs_actual_ptr = values_dofs; - std::array shape_gradients; - for (unsigned int d = 0; d < dim; ++d) - shape_gradients[d] = &shape_data.shape_gradients_face( - face_no, face_orientation, d, 0); + std::array shape_gradients; + for (unsigned int d = 0; d < dim; ++d) + shape_gradients[d] = + &shape_data.shape_gradients_face(face_no, face_orientation, d, 0); - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; ++c) + { + for (unsigned int d = 0; d < dim; ++d) { - for (unsigned int d = 0; d < dim; ++d) - { - Eval eval(nullptr, - shape_gradients[d], - nullptr, - n_dofs, - n_q_points); - - eval.template gradients<0, true, false, dim>( - values_dofs_actual_ptr, gradients_quad_ptr + d); - } - gradients_quad_ptr += n_q_points * dim; - values_dofs_actual_ptr += n_dofs; + Eval eval( + nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points); + + eval.template gradients<0, true, false, dim>( + values_dofs_actual_ptr, gradients_quad_ptr + d); } + gradients_quad_ptr += n_q_points * dim; + values_dofs_actual_ptr += n_dofs; } - - Assert(!(evaluation_flag & EvaluationFlags::hessians), - ExcNotImplemented()); - - return true; } - const unsigned int dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(shape_data.fe_degree + 1, dim - 1); + Assert(!(evaluation_flag & EvaluationFlags::hessians), + ExcNotImplemented()); - // Note: we always keep storage of values, 1st and 2nd derivatives in an - // array, so reserve space for all three here - Number *temp = fe_eval.get_scratch_data().begin(); - Number *scratch_data = temp + 3 * n_components * dofs_per_face; - - bool use_vectorization = true; + return true; + } - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) // exterior faces in the ECL loop - for (unsigned int v = 0; v < Number::size(); ++v) - if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int && - fe_eval.get_face_no(v) != fe_eval.get_face_no(0)) - use_vectorization = false; + template + static void + project_to_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval, + const bool use_vectorization, + Number *temp, + Number *scratch_data) + { + const auto &shape_info = fe_eval.get_shape_info(); if (use_vectorization == false) { + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + const unsigned int dofs_per_face = n_components * dofs_per_comp_face; + for (unsigned int v = 0; v < Number::size(); ++v) { // the loop breaks once an invalid_unsigned_int is hit for @@ -1555,8 +1550,7 @@ namespace internal // some faces might be at the boundaries but others not) if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int) { - for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; - ++i) + for (unsigned int i = 0; i < 3 * dofs_per_face; ++i) temp[i][v] = 0; continue; } @@ -1569,8 +1563,7 @@ namespace internal scratch_data, fe_eval.get_face_no(v)); - for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; - ++i) + for (unsigned int i = 0; i < 3 * dofs_per_face; ++i) temp[i][v] = scratch_data[i][v]; } } @@ -1582,6 +1575,19 @@ namespace internal values_dofs, temp, fe_eval.get_face_no()); + } + + + template + static void + evaluate_in_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + FEEvaluationData &fe_eval, + Number *temp, + Number *scratch_data) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); const unsigned int subface_index = fe_eval.get_subface_index(); constexpr unsigned int n_q_points_1d_actual = @@ -1635,6 +1641,17 @@ namespace internal fe_eval.begin_hessians(), scratch_data, subface_index); + } + + static void + adjust_quadrature_for_face_orientation( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + FEEvaluationData &fe_eval, + const bool use_vectorization, + Number *temp) + { + const auto &shape_info = fe_eval.get_shape_info(); if (use_vectorization == false) { @@ -1652,7 +1669,7 @@ namespace internal n_components, v, evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( + &shape_info.face_orientations_quad( fe_eval.get_face_orientation(v), 0), false, shape_info.n_q_points_face, @@ -1667,129 +1684,263 @@ namespace internal dim, n_components, evaluation_flag, - &fe_eval.get_shape_info().face_orientations_quad( - fe_eval.get_face_orientation(), 0), + &shape_info.face_orientations_quad(fe_eval.get_face_orientation(), 0), false, shape_info.n_q_points_face, temp, fe_eval.begin_values(), fe_eval.begin_gradients(), fe_eval.begin_hessians()); + } + + + + template + static bool + evaluate_tensor(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + + // Note: we always keep storage of values, 1st and 2nd derivatives in an + // array, so reserve space for all three here + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + + bool use_vectorization = true; + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) // exterior faces in the ECL loop + for (unsigned int v = 0; v < Number::size(); ++v) + if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int && + fe_eval.get_face_no(v) != fe_eval.get_face_no(0)) + use_vectorization = false; + + project_to_face(n_components, + evaluation_flag, + values_dofs, + fe_eval, + use_vectorization, + temp, + scratch_data); + + evaluate_in_face( + n_components, evaluation_flag, fe_eval, temp, scratch_data); + + if (dim == 3) + adjust_quadrature_for_face_orientation( + n_components, evaluation_flag, fe_eval, use_vectorization, temp); return false; } + + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + + if (shape_info.element_type == MatrixFreeFunctions::tensor_none) + return evaluate_tensor_none(n_components, + evaluation_flag, + values_dofs, + fe_eval); + else + return evaluate_tensor(n_components, + evaluation_flag, + values_dofs, + fe_eval); + } }; template - struct FEFaceEvaluationImplIntegrateSelector + struct FEFaceEvaluationImplProjectToFaceSelector + { + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + + // Note: we always keep storage of values, 1st and 2nd derivatives in an + // array, so reserve space for all three here + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + + bool use_vectorization = true; + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) // exterior faces in the ECL loop + for (unsigned int v = 0; v < Number::size(); ++v) + if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int && + fe_eval.get_face_no(v) != fe_eval.get_face_no(0)) + use_vectorization = false; + + FEFaceEvaluationImplEvaluateSelector:: + template project_to_face(n_components, + evaluation_flag, + values_dofs, + fe_eval, + use_vectorization, + temp, + scratch_data); + + return false; + } + }; + + + + template + struct FEFaceEvaluationImplEvaluateInFaceSelector { template static bool run(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number *values_dofs, + const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData &fe_eval) { const auto &shape_info = fe_eval.get_shape_info(); const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + + // Note: we always keep storage of values, 1st and 2nd derivatives in an + // array, so reserve space for all three here + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + + FEFaceEvaluationImplEvaluateSelector:: + template evaluate_in_face( + n_components, evaluation_flag, fe_eval, temp, scratch_data); + + return false; + } + }; + + + + template + struct FEFaceEvaluationImplIntegrateSelector + { + static bool + integrate_tensor_none( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool sum_into_values) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); using Number2 = typename FEEvaluationData::shape_info_number_type; - if (shape_info.element_type == MatrixFreeFunctions::tensor_none) - { - Assert((fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) == false, - ExcNotImplemented()); - - const unsigned int face_no = fe_eval.get_face_no(); - const unsigned int face_orientation = fe_eval.get_face_orientation(); - const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; - const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; + Assert((fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) == false, + ExcNotImplemented()); - using Eval = - EvaluatorTensorProduct; + const unsigned int face_no = fe_eval.get_face_no(); + const unsigned int face_orientation = fe_eval.get_face_orientation(); + const std::size_t n_dofs = shape_info.dofs_per_component_on_cell; + const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; - if (integration_flag & EvaluationFlags::values) - { - const auto *const shape_values = - &shape_data.shape_values_face(face_no, face_orientation, 0); + using Eval = + EvaluatorTensorProduct; - auto *values_quad_ptr = fe_eval.begin_values(); - auto *values_dofs_actual_ptr = values_dofs; + if (integration_flag & EvaluationFlags::values) + { + const auto *const shape_values = + &shape_data.shape_values_face(face_no, face_orientation, 0); - Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); - for (unsigned int c = 0; c < n_components; ++c) - { - eval.template values<0, false, false>(values_quad_ptr, - values_dofs_actual_ptr); + auto *values_quad_ptr = fe_eval.begin_values(); + auto *values_dofs_actual_ptr = values_dofs; - values_quad_ptr += n_q_points; - values_dofs_actual_ptr += n_dofs; - } + Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points); + for (unsigned int c = 0; c < n_components; ++c) + { + if (sum_into_values) + eval.template values<0, false, true>(values_quad_ptr, + values_dofs_actual_ptr); + else + eval.template values<0, false, false>(values_quad_ptr, + values_dofs_actual_ptr); + values_quad_ptr += n_q_points; + values_dofs_actual_ptr += n_dofs; } + } - if (integration_flag & EvaluationFlags::gradients) - { - auto *gradients_quad_ptr = fe_eval.begin_gradients(); - auto *values_dofs_actual_ptr = values_dofs; + if (integration_flag & EvaluationFlags::gradients) + { + auto *gradients_quad_ptr = fe_eval.begin_gradients(); + auto *values_dofs_actual_ptr = values_dofs; - std::array shape_gradients; - for (unsigned int d = 0; d < dim; ++d) - shape_gradients[d] = &shape_data.shape_gradients_face( - face_no, face_orientation, d, 0); + std::array shape_gradients; + for (unsigned int d = 0; d < dim; ++d) + shape_gradients[d] = + &shape_data.shape_gradients_face(face_no, face_orientation, d, 0); - for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int c = 0; c < n_components; ++c) + { + for (unsigned int d = 0; d < dim; ++d) { - for (unsigned int d = 0; d < dim; ++d) - { - Eval eval(nullptr, - shape_gradients[d], - nullptr, - n_dofs, - n_q_points); - - if (!(integration_flag & EvaluationFlags::values) && - d == 0) - eval.template gradients<0, false, false, dim>( - gradients_quad_ptr + d, values_dofs_actual_ptr); - else - eval.template gradients<0, false, true, dim>( - gradients_quad_ptr + d, values_dofs_actual_ptr); - } - gradients_quad_ptr += n_q_points * dim; - values_dofs_actual_ptr += n_dofs; + Eval eval( + nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points); + + if (!sum_into_values && + !(integration_flag & EvaluationFlags::values) && d == 0) + eval.template gradients<0, false, false, dim>( + gradients_quad_ptr + d, values_dofs_actual_ptr); + else + eval.template gradients<0, false, true, dim>( + gradients_quad_ptr + d, values_dofs_actual_ptr); } + gradients_quad_ptr += n_q_points * dim; + values_dofs_actual_ptr += n_dofs; } - - Assert(!(integration_flag & EvaluationFlags::hessians), - ExcNotImplemented()); - - return true; } - const unsigned int dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(shape_data.fe_degree + 1, dim - 1); - - Number *temp = fe_eval.get_scratch_data().begin(); - Number *scratch_data = temp + 3 * n_components * dofs_per_face; + Assert(!(integration_flag & EvaluationFlags::hessians), + ExcNotImplemented()); - bool use_vectorization = true; + return true; + } - if (fe_eval.get_dof_access_index() == - MatrixFreeFunctions::DoFInfo::dof_access_cell && - fe_eval.is_interior_face() == false) // exterior faces in the ECL loop - use_vectorization = - fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int && - std::all_of(fe_eval.get_cell_ids().begin() + 1, - fe_eval.get_cell_ids().end(), - [&](const auto &v) { - return v == fe_eval.get_cell_ids()[0] || - v == numbers::invalid_unsigned_int; - }); + static void + adjust_quadrature_for_face_orientation( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + FEEvaluationData &fe_eval, + const bool use_vectorization, + Number *temp) + { + const auto &shape_info = fe_eval.get_shape_info(); if (use_vectorization == false) { @@ -1830,6 +1981,18 @@ namespace internal fe_eval.begin_values(), fe_eval.begin_gradients(), fe_eval.begin_hessians()); + } + + template + static void + integrate_in_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + FEEvaluationData &fe_eval, + Number *temp, + Number *scratch_data) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; @@ -1884,6 +2047,27 @@ namespace internal fe_eval.begin_hessians(), scratch_data, subface_index); + } + + template + static void + collect_from_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool use_vectorization, + const Number *temp, + Number *scratch_data, + const bool sum_into_values) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + const unsigned int dofs_per_face = n_components * dofs_per_comp_face; if (use_vectorization == false) { @@ -1899,23 +2083,195 @@ namespace internal template interpolate(n_components, integration_flag, shape_info, - values_dofs, + temp, scratch_data, fe_eval.get_face_no(v)); - for (unsigned int i = 0; i < 3 * n_components * dofs_per_face; - ++i) - temp[i][v] = scratch_data[i][v]; + if (sum_into_values) + for (unsigned int i = 0; i < 3 * dofs_per_face; ++i) + values_dofs[i][v] += scratch_data[i][v]; + else + for (unsigned int i = 0; i < 3 * dofs_per_face; ++i) + values_dofs[i][v] = scratch_data[i][v]; } } else - FEFaceNormalEvaluationImpl:: - template interpolate(n_components, - integration_flag, - shape_info, - temp, - values_dofs, - fe_eval.get_face_no()); + { + if (sum_into_values) + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flag, + shape_info, + temp, + values_dofs, + fe_eval.get_face_no()); + else + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flag, + shape_info, + temp, + values_dofs, + fe_eval.get_face_no()); + } + } + + template + static bool + integrate_tensor(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool sum_into_values) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + + bool use_vectorization = true; + + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) // exterior faces in the ECL loop + use_vectorization = + fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int && + std::all_of(fe_eval.get_cell_ids().begin() + 1, + fe_eval.get_cell_ids().end(), + [&](const auto &v) { + return v == fe_eval.get_cell_ids()[0] || + v == numbers::invalid_unsigned_int; + }); + + if (dim == 3) + adjust_quadrature_for_face_orientation( + n_components, integration_flag, fe_eval, use_vectorization, temp); + + integrate_in_face( + n_components, integration_flag, fe_eval, temp, scratch_data); + + collect_from_face(n_components, + integration_flag, + values_dofs, + fe_eval, + use_vectorization, + temp, + scratch_data, + sum_into_values); + + return false; + } + + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool sum_into_values) + { + const auto &shape_info = fe_eval.get_shape_info(); + + if (shape_info.element_type == MatrixFreeFunctions::tensor_none) + return integrate_tensor_none(n_components, + integration_flag, + values_dofs, + fe_eval, + sum_into_values); + else + return integrate_tensor(n_components, + integration_flag, + values_dofs, + fe_eval, + sum_into_values); + } + }; + + + + template + struct FEFaceEvaluationImplCollectFromFaceSelector + { + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool sum_into_values) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + + bool use_vectorization = true; + + if (fe_eval.get_dof_access_index() == + MatrixFreeFunctions::DoFInfo::dof_access_cell && + fe_eval.is_interior_face() == false) // exterior faces in the ECL loop + use_vectorization = + fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int && + std::all_of(fe_eval.get_cell_ids().begin() + 1, + fe_eval.get_cell_ids().end(), + [&](const auto &v) { + return v == fe_eval.get_cell_ids()[0] || + v == numbers::invalid_unsigned_int; + }); + + FEFaceEvaluationImplIntegrateSelector:: + template collect_from_face(n_components, + integration_flag, + values_dofs, + fe_eval, + use_vectorization, + temp, + scratch_data, + sum_into_values); + + return false; + } + }; + + + + template + struct FEFaceEvaluationImplIntegrateInFaceSelector + { + template + static bool + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + + FEEvaluationData &fe_eval) + { + const auto &shape_info = fe_eval.get_shape_info(); + const auto &shape_data = shape_info.data.front(); + + const unsigned int dofs_per_comp_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::fixed_power(shape_data.fe_degree + 1); + + Number *temp = fe_eval.get_scratch_data().begin(); + Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face; + + FEFaceEvaluationImplIntegrateSelector:: + template integrate_in_face( + n_components, integration_flag, fe_eval, temp, scratch_data); + return false; } }; diff --git a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h index 74c7c4899a..79519e39a8 100644 --- a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h @@ -47,6 +47,44 @@ namespace internal } + template + void + FEFaceEvaluationFactory::project_to_face( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval) + { + instantiation_helper_degree_run< + 1, + FEFaceEvaluationImplProjectToFaceSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + n_components, + evaluation_flag, + values_dofs, + fe_eval); + } + + + + template + void + FEFaceEvaluationFactory::evaluate_in_face( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + FEEvaluationData &fe_eval) + { + instantiation_helper_run< + 1, + FEFaceEvaluationImplEvaluateInFaceSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + fe_eval.get_shape_info().data[0].n_q_points_1d, + n_components, + evaluation_flag, + fe_eval); + } + + template void @@ -54,7 +92,8 @@ namespace internal const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, - FEEvaluationData &fe_eval) + FEEvaluationData &fe_eval, + const bool sum_into_values) { instantiation_helper_run< 1, @@ -64,6 +103,48 @@ namespace internal n_components, integration_flag, values_dofs, + fe_eval, + sum_into_values); + } + + + + template + void + FEFaceEvaluationFactory::collect_from_face( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool sum_into_values) + { + instantiation_helper_degree_run< + 1, + FEFaceEvaluationImplCollectFromFaceSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + n_components, + integration_flag, + values_dofs, + fe_eval, + sum_into_values); + } + + + + template + void + FEFaceEvaluationFactory::integrate_in_face( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + FEEvaluationData &fe_eval) + { + instantiation_helper_run< + 1, + FEFaceEvaluationImplIntegrateInFaceSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + fe_eval.get_shape_info().data[0].n_q_points_1d, + n_components, + integration_flag, fe_eval); } diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index 9e04bc977a..ac2d4f09d9 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -66,11 +66,35 @@ namespace internal const Number *values_dofs, FEEvaluationData &fe_eval); + static void + project_to_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number *values_dofs, + FEEvaluationData &fe_eval); + + static void + evaluate_in_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + FEEvaluationData &fe_eval); + static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, - FEEvaluationData &fe_eval); + FEEvaluationData &fe_eval, + const bool sum_into_values); + + static void + collect_from_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number *values_dofs, + FEEvaluationData &fe_eval, + const bool sum_into_values); + + static void + integrate_in_face(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + FEEvaluationData &fe_eval); static bool fast_evaluation_supported(const unsigned int given_degree, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index cf135ebd54..1ed360ce38 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2604,6 +2604,28 @@ public: evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag); + /** + * Projects the values, the gradients, and the Hessians into the face DoFs of + * the current face using the internally stored cell DoFs. + */ + void + project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag); + + /** + * Projects the values, the gradients, and the Hessians into the face DoFs of + * the current face using the cell DoFs provided via `values_array`. + */ + void + project_to_face(const VectorizedArrayType *values_array, + const EvaluationFlags::EvaluationFlags evaluation_flag); + + /** + * Evaluates the values, the gradients, and the Hessians in-face, + * interpolating into the face quadrature points. + */ + void + evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag); + /** * Reads from the input vector and evaluates the function values, the * gradients, and the Laplacians of the FE function at the quadrature points @@ -2630,7 +2652,8 @@ public: * distribute_local_to_global() or set_dof_values() methods). */ void - integrate(const EvaluationFlags::EvaluationFlags integration_flag); + integrate(const EvaluationFlags::EvaluationFlags integration_flag, + const bool sum_into_values = false); /** * This function takes the values and/or gradients that are stored on @@ -2642,7 +2665,34 @@ public: */ void integrate(const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType *values_array); + VectorizedArrayType *values_array, + const bool sum_into_values = false); + + /** + * This function tests the values, gradients and Hessians submitted on the + * face quadrature points by multiplying with the in-face basis function + * values, gradients and Hessians and accumulating to the respective face + * DoFs. + */ + void + integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag); + + /** + * Collects the contributions from the face DoFs of values, normal gradients + * and normal Hessians to the internal cell DoFs. + */ + void + collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, + const bool sum_into_values = false); + + /** + * Collects the contributions from the face DoFs of values, normal gradients + * and normal Hessians to the cell DoFs specified via `values_array`. + */ + void + collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, + VectorizedArrayType *values_array, + const bool sum_into_values = false); /** * This function takes the values and/or gradients that are stored on @@ -2705,7 +2755,7 @@ public: * static_dofs_per_component, but the number depends on the actual element * selected and is thus not static. */ - const unsigned int dofs_per_component; + const unsigned int dofs_per_component_on_cell; /** * The number of degrees of freedom on the cell accumulated over all @@ -2715,6 +2765,22 @@ public: */ const unsigned int dofs_per_cell; + /** + * The number of degrees of freedom of a single component on the cell for + * the underlying evaluation object. Usually close to + * static_dofs_per_component, but the number depends on the actual element + * selected and is thus not static. + */ + const unsigned int dofs_per_component_on_face; + + /** + * The number of degrees of freedom on the cell accumulated over all + * components in the current evaluation object. Usually close to + * static_dofs_per_cell = static_dofs_per_component*n_components, but the + * number depends on the actual element selected and is thus not static. + */ + const unsigned int dofs_per_face; + /** * The number of quadrature points in use. If the number of quadrature * points in 1d is given as a template, this number is simply the @@ -8244,8 +8310,10 @@ inline FEFaceEvaluationdata->dofs_per_component_on_cell) + , dofs_per_component_on_cell(this->data->dofs_per_component_on_cell) , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_) + , dofs_per_component_on_face(this->data->dofs_per_component_on_face) + , dofs_per_face(this->data->dofs_per_component_on_face * n_components_) , n_q_points(this->n_quadrature_points) {} @@ -8646,9 +8714,144 @@ FEFaceEvaluation:: - integrate(const EvaluationFlags::EvaluationFlags integration_flag) + project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag) { - integrate(integration_flag, this->values_dofs); +# ifdef DEBUG + Assert(this->dof_values_initialized, ExcNotInitialized()); +# endif + + project_to_face(this->values_dofs, evaluation_flag); +} + + + +template +inline void +FEFaceEvaluation:: + project_to_face(const VectorizedArrayType *values_array, + const EvaluationFlags::EvaluationFlags evaluation_flag) +{ + Assert((evaluation_flag & + ~(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians)) == 0, + ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, " + "and EvaluationFlags::hessians are supported.")); + + const bool hessians_on_general_cells = + evaluation_flag & EvaluationFlags::hessians && + (this->cell_type > internal::MatrixFreeFunctions::affine); + EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag; + if (hessians_on_general_cells) + evaluation_flag_actual |= EvaluationFlags::gradients; + + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas && + evaluation_flag & EvaluationFlags::gradients && + (this->cell_type > internal::MatrixFreeFunctions::affine)) + evaluation_flag_actual |= EvaluationFlags::values; + + if constexpr (fe_degree > -1) + internal::FEFaceEvaluationImplProjectToFaceSelector< + dim, + VectorizedArrayType>::template run(n_components, + evaluation_flag_actual, + values_array, + *this); + else + internal::FEFaceEvaluationFactory:: + project_to_face(n_components, + evaluation_flag_actual, + values_array, + *this); + + // face dofs initialized +} + + + +template +inline void +FEFaceEvaluation:: + evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag) +{ + Assert((evaluation_flag & + ~(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians)) == 0, + ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, " + "and EvaluationFlags::hessians are supported.")); + + const bool hessians_on_general_cells = + evaluation_flag & EvaluationFlags::hessians && + (this->cell_type > internal::MatrixFreeFunctions::affine); + EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag; + if (hessians_on_general_cells) + evaluation_flag_actual |= EvaluationFlags::gradients; + + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas && + evaluation_flag & EvaluationFlags::gradients && + (this->cell_type > internal::MatrixFreeFunctions::affine)) + evaluation_flag_actual |= EvaluationFlags::values; + + if constexpr (fe_degree > -1) + internal::FEFaceEvaluationImplEvaluateInFaceSelector< + dim, + VectorizedArrayType>::template run(n_components, + evaluation_flag_actual, + *this); + else + internal::FEFaceEvaluationFactory:: + evaluate_in_face(n_components, evaluation_flag_actual, *this); + +# ifdef DEBUG + if (evaluation_flag_actual & EvaluationFlags::values) + this->values_quad_initialized = true; + if (evaluation_flag_actual & EvaluationFlags::gradients) + this->gradients_quad_initialized = true; + if ((evaluation_flag_actual & EvaluationFlags::hessians) != 0u) + this->hessians_quad_initialized = true; +# endif +} + + + +template +inline void +FEFaceEvaluation:: + integrate(const EvaluationFlags::EvaluationFlags integration_flag, + const bool sum_into_values) +{ + integrate(integration_flag, this->values_dofs, sum_into_values); # ifdef DEBUG this->dof_values_initialized = true; @@ -8671,7 +8874,8 @@ FEFaceEvaluation:: integrate(const EvaluationFlags::EvaluationFlags integration_flag, - VectorizedArrayType *values_array) + VectorizedArrayType *values_array, + const bool sum_into_values) { Assert((integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients | @@ -8722,10 +8926,168 @@ FEFaceEvaluation(n_components, integration_flag_actual, values_array, - *this); + *this, + sum_into_values); else internal::FEFaceEvaluationFactory::integrate( - n_components, integration_flag_actual, values_array, *this); + n_components, + integration_flag_actual, + values_array, + *this, + sum_into_values); +} + + + +template +inline void +FEFaceEvaluation:: + integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag) +{ + Assert((integration_flag & + ~(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians)) == 0, + ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, " + "and EvaluationFlags::hessians are supported.")); + + EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag; + if (integration_flag & EvaluationFlags::hessians && + (this->cell_type > internal::MatrixFreeFunctions::affine)) + { + unsigned int size = n_components * dim * n_q_points; + if ((integration_flag & EvaluationFlags::gradients) != 0u) + { + for (unsigned int i = 0; i < size; ++i) + this->gradients_quad[i] += this->gradients_from_hessians_quad[i]; + } + else + { + for (unsigned int i = 0; i < size; ++i) + this->gradients_quad[i] = this->gradients_from_hessians_quad[i]; + integration_flag_actual |= EvaluationFlags::gradients; + } + } + + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas && + integration_flag & EvaluationFlags::gradients && + this->cell_type > internal::MatrixFreeFunctions::affine && + this->divergence_is_requested == false) + { + unsigned int size = n_components * n_q_points; + if ((integration_flag & EvaluationFlags::values) != 0u) + { + for (unsigned int i = 0; i < size; ++i) + this->values_quad[i] += this->values_from_gradients_quad[i]; + } + else + { + for (unsigned int i = 0; i < size; ++i) + this->values_quad[i] = this->values_from_gradients_quad[i]; + integration_flag_actual |= EvaluationFlags::values; + } + } + + if constexpr (fe_degree > -1) + internal::FEFaceEvaluationImplIntegrateInFaceSelector< + dim, + VectorizedArrayType>::template run(n_components, + integration_flag_actual, + *this); + else + internal::FEFaceEvaluationFactory:: + integrate_in_face(n_components, integration_flag_actual, *this); + + // face dofs initialized +} + + + +template +inline void +FEFaceEvaluation:: + collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, + const bool sum_into_values) +{ + collect_from_face(integration_flag, this->values_dofs, sum_into_values); + +# ifdef DEBUG + this->dof_values_initialized = true; +# endif +} + + + +template +inline void +FEFaceEvaluation:: + collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, + VectorizedArrayType *values_array, + const bool sum_into_values) +{ + Assert((integration_flag & + ~(EvaluationFlags::values | EvaluationFlags::gradients | + EvaluationFlags::hessians)) == 0, + ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, " + "and EvaluationFlags::hessians are supported.")); + + EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag; + if (integration_flag & EvaluationFlags::hessians && + (this->cell_type > internal::MatrixFreeFunctions::affine)) + integration_flag_actual |= EvaluationFlags::gradients; + + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas && + integration_flag & EvaluationFlags::gradients && + this->cell_type > internal::MatrixFreeFunctions::affine && + this->divergence_is_requested == false) + integration_flag_actual |= EvaluationFlags::values; + + if constexpr (fe_degree > -1) + internal::FEFaceEvaluationImplCollectFromFaceSelector< + dim, + VectorizedArrayType>::template run(n_components, + integration_flag_actual, + values_array, + *this, + sum_into_values); + else + internal::FEFaceEvaluationFactory:: + collect_from_face(n_components, + integration_flag_actual, + values_array, + *this, + sum_into_values); } diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 57f5971f6a..57161ae20b 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -45,7 +45,7 @@ namespace internal DeclException1( ExcFEPointEvaluationAccessToUninitializedMappingField, std::string, - << "You are requesting information from an FEPointEvaluation " + << "You are requesting information from an FEPointEvaluationBase " << "object for which this kind of information has not been computed. " << "What information these objects compute is determined by the update_* " << "flags you pass to MappingInfo() in the Constructor. Here, " @@ -682,40 +682,16 @@ namespace internal /** - * This class provides an interface to the evaluation of interpolated solution - * values and gradients on cells on arbitrary reference point positions. These - * points can change from cell to cell, both with respect to their quantity as - * well to the location. The two typical use cases are evaluations on - * non-matching grids and particle simulations. - * - * The use of this class is similar to FEValues or FEEvaluation: The class is - * first initialized to a cell by calling `FEPointEvaluation::reinit(cell, - * unit_points)`, with the main difference to the other concepts that the - * underlying points in reference coordinates need to be passed along. Then, - * upon call to evaluate() or integrate(), the user can compute information at - * the give points. Eventually, the access functions get_value() or - * get_gradient() allow to query this information at a specific point index. - * - * The functionality is similar to creating an FEValues object with a - * Quadrature object on the `unit_points` on every cell separately and then - * calling FEValues::get_function_values or FEValues::get_function_gradients, - * and for some elements and mappings this is what actually happens - * internally. For specific combinations of Mapping and FiniteElement - * realizations, however, there is a much more efficient implementation that - * avoids the memory allocation and other expensive start-up cost of - * FEValues. Currently, the functionality is specialized for mappings derived - * from MappingQ and MappingCartesian and for finite elements with tensor - * product structure that work with the - * @ref matrixfree - * module. In those cases, the cost implied - * by this class is similar (or sometimes even somewhat lower) than using - * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`. + * Base class of FEPointEvaluation and FEFacePointEvaluation. This class needs + * usually not be called in user code and does not have any public + * constructor. The usage is through the class + * FEPointEvaluation/FEFacePointEvaluation instead. */ template -class FEPointEvaluation +class FEPointEvaluationBase { public: static constexpr unsigned int dimension = dim; @@ -736,6 +712,7 @@ public: using interface_vectorized_gradient_type = typename ETT::interface_vectorized_gradient_type; +protected: /** * Constructor. * @@ -754,14 +731,14 @@ public: * objects, this parameter allows to select a range of `n_components` * components starting from this parameter. */ - FEPointEvaluation(const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component = 0); + FEPointEvaluationBase(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component = 0); /** * Constructor to make the present class able to re-use the geometry - * data also used by other `FEPointEvaluation` objects. + * data also used by other `FEPointEvaluationBase` objects. * * @param mapping_info The MappingInfo class describes the geometry-related * data for evaluating finite-element solutions. This object enables to @@ -774,219 +751,37 @@ public: * @param first_selected_component For multi-component FiniteElement * objects, this parameter allows to select a range of `n_components` * components starting from this parameter. + * + * @param is_interior Defines if interior or exterior. Only makes sense for + * faces. */ - FEPointEvaluation( + FEPointEvaluationBase( NonMatching::MappingInfo &mapping_info, const FiniteElement &fe, - const unsigned int first_selected_component = 0); + const unsigned int first_selected_component = 0, + const bool is_interior = true); /** * Copy constructor. */ - FEPointEvaluation(FEPointEvaluation &other) noexcept; + FEPointEvaluationBase(FEPointEvaluationBase &other) noexcept; /** * Move constructor. */ - FEPointEvaluation(FEPointEvaluation &&other) noexcept; + FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept; /** * Destructor. */ - ~FEPointEvaluation(); - - /** - * Set up the mapping information for the given cell, e.g., by computing the - * Jacobian of the mapping for the given points if gradients of the functions - * are requested. - * - * @param[in] cell An iterator to the current cell - * - * @param[in] unit_points List of points in the reference locations of the - * current cell where the FiniteElement object should be - * evaluated/integrated in the evaluate() and integrate() functions. - */ - void - reinit(const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points); - - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the single cell in the MappingInfo object. - */ - void - reinit(); - - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the cell in the MappingInfo object. - */ - void - reinit(const unsigned int cell_index); - - /** - * Reinitialize the evaluator to point to the correct precomputed mapping of - * the face in the MappingInfo object. - */ - void - reinit(const unsigned int cell_index, const unsigned int face_number); - - /** - * This function interpolates the finite element solution, represented by - * `solution_values`, on the cell and `unit_points` passed to reinit(). - * - * @param[in] solution_values This array is supposed to contain the unknown - * values on the element read out by - * `FEEvaluation::read_dof_values(global_vector)`. - * - * @param[in] evaluation_flags Flags specifying which quantities should be - * evaluated at the points. - */ - template - void - evaluate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * This function interpolates the finite element solution, represented by - * `solution_values`, on the cell and `unit_points` passed to reinit(). - * - * @param[in] solution_values This array is supposed to contain the unknown - * values on the element as returned by `cell->get_dof_values(global_vector, - * solution_values)`. - * - * @param[in] evaluation_flags Flags specifying which quantities should be - * evaluated at the points. - */ - void - evaluate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `FEEvaluation::set_dof_values(global_vector)` or - * `FEEvaluation::distribute_local_to_global(global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. For the default value - * `sum_into_values=false` every value of @p solution_values is zeroed out. - * - */ - template - void - integrate(const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points multiplied be - * the Jacobian determinant times the quadrature weight (JxW). - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used to during - * `cell->set_dof_values(solution_values, global_vector)` or - * `cell->distribute_local_to_global(solution_values, global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. For the default value - * `sum_into_values=false` every value of @p solution_values is zeroed out. - * - */ - void - integrate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points. This is - * similar to the integration of a bilinear form in terms of the test - * function, with the difference that this formula does not include a `JxW` - * factor (in contrast to the integrate function of this class). This allows - * the class to naturally embed point information (e.g. particles) into a - * finite element formulation. - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `FEEvaluation::set_dof_values(global_vector)` or - * `FEEvaluation::distribute_local_to_global(global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. For the default value - * `sum_into_values=false` every value of @p solution_values is zeroed out. - * - */ - template - void - test_and_sum( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); - - /** - * This function multiplies the quantities passed in by previous - * submit_value() or submit_gradient() calls by the value or gradient of the - * test functions, and performs summation over all given points. This is - * similar to the integration of a bilinear form in terms of the test - * function, with the difference that this formula does not include a `JxW` - * factor (in contrast to the integrate function of this class). This allows - * the class to naturally embed point information (e.g. particles) into a - * finite element formulation. - * - * @param[out] solution_values This array will contain the result of the - * integral, which can be used during - * `cell->set_dof_values(solution_values, global_vector)` or - * `cell->distribute_local_to_global(solution_values, global_vector)`. Note - * that for multi-component systems where only some of the components are - * selected by the present class, the entries in `solution_values` not touched - * by this class will be set to zero. - * - * @param[in] integration_flags Flags specifying which quantities should be - * integrated at the points. - * - * @param[in] sum_into_values Flag specifying if the integrated values - * should be summed into the solution values. For the default value - * `sum_into_values=false` every value of @p solution_values is zeroed out. - * - */ - void - test_and_sum(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values = false); + ~FEPointEvaluationBase(); +public: /** * Return the value at quadrature point number @p point_index after a call to * FEPointEvaluation::evaluate() with EvaluationFlags::values set, or * the value that has been stored there with a call to - * FEPointEvaluation::submit_value(). If the object is vector-valued, a + * FEPointEvaluationBase::submit_value(). If the object is vector-valued, a * vector-valued return argument is given. */ const value_type & @@ -995,10 +790,10 @@ public: /** * Write a value to the field containing the values on points * with component point_index. Access to the same field as through - * get_value(). If applied before the function FEPointEvaluation::integrate() - * with EvaluationFlags::values set is called, this specifies the value - * which is tested by all basis function on the current cell and - * integrated over. + * get_value(). If applied before the function + * FEPointEvaluation::integrate() with EvaluationFlags::values set is + * called, this specifies the value which is tested by all basis function on + * the current cell and integrated over. */ void submit_value(const value_type &value, const unsigned int point_index); @@ -1007,10 +802,10 @@ public: * Return the gradient in real coordinates at the point with index * `point_index` after a call to FEPointEvaluation::evaluate() with * EvaluationFlags::gradients set, or the gradient that has been stored there - * with a call to FEPointEvaluation::submit_gradient(). The gradient in real - * coordinates is obtained by taking the unit gradient (also accessible via - * get_unit_gradient()) and applying the inverse Jacobian of the mapping. If - * the object is vector-valued, a vector-valued return argument is given. + * with a call to FEPointEvaluationBase::submit_gradient(). The gradient in + * real coordinates is obtained by taking the unit gradient (also accessible + * via get_unit_gradient()) and applying the inverse Jacobian of the mapping. + * If the object is vector-valued, a vector-valued return argument is given. */ const gradient_type & get_gradient(const unsigned int point_index) const; @@ -1019,9 +814,9 @@ public: * Write a contribution that is tested by the gradient to the field * containing the values on points with the given `point_index`. Access to * the same field as through get_gradient(). If applied before the function - * FEPointEvaluation::integrate(EvaluationFlags::gradients) is called, this - * specifies what is tested by all basis function gradients on the current - * cell and integrated over. + * FEPointEvaluation::integrate(EvaluationFlags::gradients) is called, + * this specifies what is tested by all basis function gradients on the + * current cell and integrated over. */ void submit_gradient(const gradient_type &, const unsigned int point_index); @@ -1051,14 +846,6 @@ public: Number JxW(const unsigned int point_index) const; - /** - * Return the normal vector. This class or the MappingInfo object passed to - * this function needs to be constructed with UpdateFlags containing - * `update_normal_vectors`. - */ - Tensor<1, spacedim, Number> - normal_vector(const unsigned int point_index) const; - /** * Return the position in real coordinates of the given point index among * the points passed to reinit(). @@ -1081,7 +868,7 @@ public: inline std_cxx20::ranges::iota_view quadrature_point_indices() const; -private: +protected: static constexpr std::size_t n_lanes_user_interface = internal::VectorizedArrayTrait::width(); static constexpr std::size_t n_lanes_internal = @@ -1110,125 +897,22 @@ private: do_reinit(); /** - * Resizes necessary data fields, reads in and renumbers solution values. - * Interpolates onto face if face path is selected. + * Number of quadrature batches of the current cell/face. */ - template - void - prepare_evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + const unsigned int n_q_batches; /** - * Evaluates the actual interpolation on the cell or face for a quadrature - * batch. + * Number of quadrature points/batches of the current cell/face. */ - template - void - compute_evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags, - const unsigned int n_shapes, - const unsigned int qb, - vectorized_value_type &value, - interface_vectorized_gradient_type &gradient); + const unsigned int n_q_points; /** - * Fast path of the evaluate function. + * Number of quadrature points of the current cell/face. */ - template - void - evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + const unsigned int n_q_points_scalar; /** - * Slow path of the evaluate function using FEValues. - */ - template - void - evaluate_slow( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); - - /** - * Integrates the product of the data passed in by submit_value() and - * submit_gradient() with the values or gradients of test functions on the - * cell or face for a given quadrature batch. - */ - template - void - compute_integrate_fast( - const EvaluationFlags::EvaluationFlags &integration_flags, - const unsigned int n_shapes, - const unsigned int qb, - const vectorized_value_type value, - const interface_vectorized_gradient_type gradient, - vectorized_value_type *solution_values_vectorized_linear); - - /** - * Addition across the lanes of VectorizedArray as accumulated by the - * compute_integrate_fast_function(), writing the sum into the result vector. - * Applies face contributions to cell contributions for face path. - */ - template - void - finish_integrate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - vectorized_value_type *solution_values_vectorized_linear, - const bool sum_into_values); - - /** - * Fast path of the integrate function. - */ - template - void - integrate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); - - /** - * Slow path of the integrate function using FEValues. - */ - template - void - integrate_slow( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); - - /** - * Implementation of the integrate/test_and_sum function. - */ - template - void - do_integrate( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values); - - /** - * Number of quadrature batches of the current cell/face. - */ - const unsigned int n_q_batches; - - /** - * Number of quadrature points/batches of the current cell/face. - */ - const unsigned int n_q_points; - - /** - * Number of quadrature points of the current cell/face. - */ - const unsigned int n_q_points_scalar; - - /** - * Pointer to the Mapping object passed to the constructor. + * Pointer to the Mapping object passed to the constructor. */ SmartPointer> mapping; @@ -1272,7 +956,7 @@ private: AlignedVector solution_renumbered_vectorized; /** - * Temporary array for the use_face_path path (scalar). + * Temporary array for the face path (scalar). */ AlignedVector scratch_data_scalar; @@ -1294,7 +978,7 @@ private: /** * Pointer to first unit point batch of current face from MappingInfo, - * set internally during do_reinit(). Needed for use_face_path path. + * set internally during do_reinit(). Needed for face path. */ const Point *unit_point_faces_ptr; @@ -1348,13 +1032,7 @@ private: unsigned int dofs_per_component_face; /** - * Bool indicating if use_face_path path should be chosen. Set during - * do_reinit(). - */ - bool use_face_path; - - /** - * Scalar ShapeInfo object needed for use_face_path path. + * Scalar ShapeInfo object needed for face path. */ internal::MatrixFreeFunctions::ShapeInfo shape_info; @@ -1429,24 +1107,25 @@ private: * reinit()) at the vectorized unit points on faces. */ AlignedVector> shapes_faces; + + const bool is_interior; }; // ----------------------- template and inline function ---------------------- template -FEPointEvaluation::FEPointEvaluation( - const Mapping &mapping, - const FiniteElement &fe, - const UpdateFlags update_flags, - const unsigned int first_selected_component) +FEPointEvaluationBase:: + FEPointEvaluationBase(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component) : n_q_batches(numbers::invalid_unsigned_int) , n_q_points(numbers::invalid_unsigned_int) , n_q_points_scalar(numbers::invalid_unsigned_int) , mapping(&mapping) , fe(&fe) , JxW_ptr(nullptr) - , use_face_path(false) , update_flags(update_flags) , mapping_info_on_the_fly( std::make_unique>( @@ -1456,6 +1135,7 @@ FEPointEvaluation::FEPointEvaluation( , current_cell_index(numbers::invalid_unsigned_int) , current_face_number(numbers::invalid_unsigned_int) , is_reinitialized(false) + , is_interior(true) { setup(first_selected_component); } @@ -1463,22 +1143,24 @@ FEPointEvaluation::FEPointEvaluation( template -FEPointEvaluation::FEPointEvaluation( - NonMatching::MappingInfo &mapping_info, - const FiniteElement &fe, - const unsigned int first_selected_component) +FEPointEvaluationBase:: + FEPointEvaluationBase( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const unsigned int first_selected_component, + const bool is_interior) : n_q_batches(numbers::invalid_unsigned_int) , n_q_points(numbers::invalid_unsigned_int) , n_q_points_scalar(numbers::invalid_unsigned_int) , mapping(&mapping_info.get_mapping()) , fe(&fe) , JxW_ptr(nullptr) - , use_face_path(false) , update_flags(mapping_info.get_update_flags()) , mapping_info(&mapping_info) , current_cell_index(numbers::invalid_unsigned_int) , current_face_number(numbers::invalid_unsigned_int) , is_reinitialized(false) + , is_interior(is_interior) { setup(first_selected_component); connection_is_reinitialized = mapping_info.connect_is_reinitialized( @@ -1488,8 +1170,9 @@ FEPointEvaluation::FEPointEvaluation( template -FEPointEvaluation::FEPointEvaluation( - FEPointEvaluation &other) noexcept +FEPointEvaluationBase:: + FEPointEvaluationBase( + FEPointEvaluationBase &other) noexcept : n_q_batches(other.n_q_batches) , n_q_points(other.n_q_points) , n_q_points_scalar(other.n_q_points_scalar) @@ -1504,7 +1187,6 @@ FEPointEvaluation::FEPointEvaluation( , gradients(other.gradients) , dofs_per_component(other.dofs_per_component) , dofs_per_component_face(other.dofs_per_component_face) - , use_face_path(false) , component_in_base_element(other.component_in_base_element) , nonzero_shape_function_component(other.nonzero_shape_function_component) , update_flags(other.update_flags) @@ -1522,6 +1204,7 @@ FEPointEvaluation::FEPointEvaluation( , is_reinitialized(false) , shapes(other.shapes) , shapes_faces(other.shapes_faces) + , is_interior(other.is_interior) { connection_is_reinitialized = mapping_info->connect_is_reinitialized( [this]() { this->is_reinitialized = false; }); @@ -1530,8 +1213,10 @@ FEPointEvaluation::FEPointEvaluation( template -FEPointEvaluation::FEPointEvaluation( - FEPointEvaluation &&other) noexcept +FEPointEvaluationBase:: + FEPointEvaluationBase( + FEPointEvaluationBase + &&other) noexcept : n_q_batches(other.n_q_batches) , n_q_points(other.n_q_points) , n_q_points_scalar(other.n_q_points_scalar) @@ -1546,7 +1231,6 @@ FEPointEvaluation::FEPointEvaluation( , gradients(other.gradients) , dofs_per_component(other.dofs_per_component) , dofs_per_component_face(other.dofs_per_component_face) - , use_face_path(false) , component_in_base_element(other.component_in_base_element) , nonzero_shape_function_component(other.nonzero_shape_function_component) , update_flags(other.update_flags) @@ -1559,6 +1243,7 @@ FEPointEvaluation::FEPointEvaluation( , is_reinitialized(false) , shapes(other.shapes) , shapes_faces(other.shapes_faces) + , is_interior(other.is_interior) { connection_is_reinitialized = mapping_info->connect_is_reinitialized( [this]() { this->is_reinitialized = false; }); @@ -1567,7 +1252,8 @@ FEPointEvaluation::FEPointEvaluation( template -FEPointEvaluation::~FEPointEvaluation() +FEPointEvaluationBase:: + ~FEPointEvaluationBase() { connection_is_reinitialized.disconnect(); } @@ -1576,7 +1262,7 @@ FEPointEvaluation::~FEPointEvaluation() template void -FEPointEvaluation::setup( +FEPointEvaluationBase::setup( const unsigned int first_selected_component) { AssertIndexRange(first_selected_component + n_components, @@ -1626,8 +1312,8 @@ FEPointEvaluation::setup( poly[1].value(1.) == 1.) && (fe->n_components() == n_components); - const unsigned int size_face = 2 * dofs_per_component_face; - const unsigned int size_cell = dofs_per_component; + const unsigned int size_face = 3 * dofs_per_component_face * n_components; + const unsigned int size_cell = dofs_per_component * n_components; scratch_data_scalar.resize(size_face + size_cell); solution_renumbered.resize(dofs_per_component); @@ -1660,107 +1346,10 @@ FEPointEvaluation::setup( -template -inline void -FEPointEvaluation::reinit( - const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points) -{ - // reinit is only allowed for mapping computation on the fly - AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented()); - - mapping_info->reinit(cell, unit_points); - - if (!fast_path) - { - fe_values = std::make_shared>( - *mapping, - *fe, - Quadrature( - std::vector>(unit_points.begin(), unit_points.end())), - update_flags); - fe_values->reinit(cell); - } - - if (use_linear_path) - do_reinit(); - else - do_reinit(); -} - - - -template -inline void -FEPointEvaluation::reinit() -{ - current_cell_index = numbers::invalid_unsigned_int; - current_face_number = numbers::invalid_unsigned_int; - - if (use_linear_path) - do_reinit(); - else - do_reinit(); -} - - - -template -inline void -FEPointEvaluation::reinit( - const unsigned int cell_index) -{ - current_cell_index = cell_index; - current_face_number = numbers::invalid_unsigned_int; - - if (use_linear_path) - do_reinit(); - else - do_reinit(); - - if (!fast_path) - { - std::vector> unit_points(n_q_points_scalar); - - for (unsigned int v = 0; v < n_q_points_scalar; ++v) - for (unsigned int d = 0; d < dim; ++d) - unit_points[v][d] = - unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal]; - - fe_values = std::make_shared>( - *mapping, - *fe, - Quadrature( - std::vector>(unit_points.begin(), unit_points.end())), - update_flags); - - fe_values->reinit(mapping_info->get_cell_iterator(current_cell_index)); - } -} - - - -template -inline void -FEPointEvaluation::reinit( - const unsigned int cell_index, - const unsigned int face_number) -{ - current_cell_index = cell_index; - current_face_number = face_number; - - if (use_linear_path) - do_reinit(); - else - do_reinit(); -} - - - template template inline void -FEPointEvaluation::do_reinit() +FEPointEvaluationBase::do_reinit() { const unsigned int geometry_index = mapping_info->template compute_geometry_index_offset( @@ -1794,16 +1383,11 @@ FEPointEvaluation::do_reinit() return; } - // use face path if mapping_info in face state and number of quadrature points - // is large enough - use_face_path = - is_face && (mapping_info->is_face_state() && n_q_points_scalar >= 6); - // set unit point pointer const unsigned int unit_point_offset = mapping_info->compute_unit_point_index_offset(geometry_index); - if (use_face_path) + if (is_face) unit_point_faces_ptr = mapping_info->get_unit_point_faces(unit_point_offset); else @@ -1820,19 +1404,21 @@ FEPointEvaluation::do_reinit() if (update_flags_mapping & UpdateFlags::update_quadrature_points) real_point_ptr = mapping_info->get_real_point(data_offset); if (update_flags_mapping & UpdateFlags::update_jacobians) - jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset); + jacobian_ptr = + mapping_info->get_jacobian(compressed_data_offset, is_interior); if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) inverse_jacobian_ptr = - mapping_info->get_inverse_jacobian(compressed_data_offset); + mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior); if (update_flags_mapping & UpdateFlags::update_normal_vectors) normal_ptr = mapping_info->get_normal_vector(data_offset); if (update_flags_mapping & UpdateFlags::update_JxW_values) JxW_ptr = mapping_info->get_JxW(data_offset); #else real_point_ptr = mapping_info->get_real_point(data_offset); - jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset); + jacobian_ptr = + mapping_info->get_jacobian(compressed_data_offset, is_interior); inverse_jacobian_ptr = - mapping_info->get_inverse_jacobian(compressed_data_offset); + mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior); normal_ptr = mapping_info->get_normal_vector(data_offset); JxW_ptr = mapping_info->get_JxW(data_offset); #endif @@ -1842,7 +1428,7 @@ FEPointEvaluation::do_reinit() const std::size_t n_shapes = poly.size(); for (unsigned int qb = 0; qb < n_q_batches; ++qb) - if (use_face_path) + if (is_face) { if (dim > 1) { @@ -1871,220 +1457,1244 @@ FEPointEvaluation::do_reinit() template -template -inline void -FEPointEvaluation::prepare_evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) +inline const typename FEPointEvaluationBase::value_type & +FEPointEvaluationBase::get_value( + const unsigned int point_index) const { - const unsigned int dofs_per_comp = - is_linear ? Utilities::pow(2, dim) : dofs_per_component; + AssertIndexRange(point_index, values.size()); + return values[point_index]; +} - for (unsigned int comp = 0; comp < n_components; ++comp) - { - const std::size_t offset = - (component_in_base_element + comp) * dofs_per_comp; - if (is_face_path) - { - const ScalarNumber *input; - if (is_linear || renumber.empty()) - { - for (unsigned int i = 0; i < dofs_per_comp; ++i) - scratch_data_scalar[i] = solution_values[i + offset]; - input = scratch_data_scalar.data(); - } - else - { - const unsigned int *renumber_ptr = renumber.data() + offset; - for (unsigned int i = 0; i < dofs_per_comp; ++i) - scratch_data_scalar[i] = solution_values[renumber_ptr[i]]; - input = scratch_data_scalar.data(); - } - ScalarNumber *output = scratch_data_scalar.begin() + dofs_per_comp; +template +inline const typename FEPointEvaluationBase::gradient_type & +FEPointEvaluationBase::get_gradient( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, gradients.size()); + return gradients[point_index]; +} - internal::FEFaceNormalEvaluationImpl:: - template interpolate(1, - evaluation_flags, - shape_info, - input, - output, - current_face_number); - const unsigned int dofs_per_comp_face = - is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face; - for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i) - ETT::read_value(output[i], comp, solution_renumbered[i]); - } - else - { - if ((is_linear && n_components == 1) || renumber.empty()) - { - for (unsigned int i = 0; i < dofs_per_comp; ++i) - ETT::read_value(solution_values[i + offset], - comp, - solution_renumbered[i]); - } - else - { - const unsigned int *renumber_ptr = renumber.data() + offset; - for (unsigned int i = 0; i < dofs_per_comp; ++i) - ETT::read_value(solution_values[renumber_ptr[i]], - comp, - solution_renumbered[i]); - } - } - } + +template +inline void +FEPointEvaluationBase::submit_value( + const value_type &value, + const unsigned int point_index) +{ + AssertIndexRange(point_index, n_q_points); + values[point_index] = value; } template -template inline void -FEPointEvaluation::compute_evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags, - const unsigned int n_shapes, - const unsigned int qb, - vectorized_value_type &value, - interface_vectorized_gradient_type &gradient) +FEPointEvaluationBase::submit_gradient( + const gradient_type &gradient, + const unsigned int point_index) { - if (is_face_path) - { - if (evaluation_flags & EvaluationFlags::gradients) - { - const std::array interpolated_value = - is_linear ? - internal::evaluate_tensor_product_value_and_gradient_linear< - dim - 1, - scalar_value_type, - VectorizedArrayType, - 2>(solution_renumbered.data(), unit_point_faces_ptr[qb]) : - internal::evaluate_tensor_product_value_and_gradient_shapes< - dim - 1, - scalar_value_type, - VectorizedArrayType, - 2, - false>(shapes_faces.data() + qb * n_shapes, - n_shapes, - solution_renumbered.data()); + AssertIndexRange(point_index, n_q_points); + gradients[point_index] = gradient; +} - value = interpolated_value[dim - 1]; - // reorder derivative from tangential/normal derivatives into tensor - // in physical coordinates - if (current_face_number / 2 == 0) - { - gradient[0] = interpolated_value[dim]; - if (dim > 1) - gradient[1] = interpolated_value[0]; - if (dim > 2) - gradient[2] = interpolated_value[1]; - } - else if (current_face_number / 2 == 1) - { - if (dim > 1) - gradient[1] = interpolated_value[dim]; - if (dim == 3) - { - gradient[0] = interpolated_value[1]; - gradient[2] = interpolated_value[0]; - } - else if (dim == 2) - gradient[0] = interpolated_value[0]; - else - Assert(false, ExcInternalError()); - } - else if (current_face_number / 2 == 2) - { - if (dim > 2) - { - gradient[0] = interpolated_value[0]; - gradient[1] = interpolated_value[1]; - gradient[2] = interpolated_value[dim]; - } - else - Assert(false, ExcInternalError()); - } - else - Assert(false, ExcInternalError()); - } - else - { - value = is_linear ? internal::evaluate_tensor_product_value_linear< - dim - 1, - scalar_value_type, - VectorizedArrayType>(solution_renumbered.data(), - unit_point_faces_ptr[qb]) : - internal::evaluate_tensor_product_value_shapes< - dim - 1, - scalar_value_type, - VectorizedArrayType, - false>(shapes_faces.data() + qb * n_shapes, - n_shapes, - solution_renumbered.data()); + + +template +inline DerivativeForm<1, dim, spacedim, Number> +FEPointEvaluationBase::jacobian( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(jacobian_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_jacobians")); + return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions:: + GeometryType::affine ? + 0 : + point_index]; +} + + + +template +inline DerivativeForm<1, spacedim, dim, Number> +FEPointEvaluationBase::inverse_jacobian( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(inverse_jacobian_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_inverse_jacobians")); + return inverse_jacobian_ptr + [cell_type <= + ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ? + 0 : + point_index]; +} + + + +template +inline Number +FEPointEvaluationBase::JxW( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(JxW_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_JxW_values")); + return JxW_ptr[point_index]; +} + + + +template +inline Point +FEPointEvaluationBase::real_point( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(real_point_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_quadrature_points")); + return real_point_ptr[point_index]; +} + + + +template +inline Point +FEPointEvaluationBase::unit_point( + const unsigned int point_index) const +{ + AssertIndexRange(point_index, n_q_points); + Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!")); + Point unit_point; + for (unsigned int d = 0; d < dim; ++d) + unit_point[d] = internal::VectorizedArrayTrait::get_from_vectorized( + unit_point_ptr[point_index / stride][d], point_index % stride); + return unit_point; +} + + + +template +inline std_cxx20::ranges::iota_view +FEPointEvaluationBase:: + quadrature_point_indices() const +{ + return {0U, n_q_points}; +} + + + +/** + * This class provides an interface to the evaluation of interpolated solution + * values and gradients on cells on arbitrary reference point positions. These + * points can change from cell to cell, both with respect to their quantity as + * well to the location. The two typical use cases are evaluations on + * non-matching grids and particle simulations. + * + * The use of this class is similar to FEValues or FEEvaluation: The class is + * first initialized to a cell by calling `FEPointEvaluation::reinit(cell, + * unit_points)`, with the main difference to the other concepts that the + * underlying points in reference coordinates need to be passed along. Then, + * upon call to evaluate() or integrate(), the user can compute information at + * the give points. Eventually, the access functions get_value() or + * get_gradient() allow to query this information at a specific point index. + * + * The functionality is similar to creating an FEValues object with a + * Quadrature object on the `unit_points` on every cell separately and then + * calling FEValues::get_function_values or FEValues::get_function_gradients, + * and for some elements and mappings this is what actually happens + * internally. For specific combinations of Mapping and FiniteElement + * realizations, however, there is a much more efficient implementation that + * avoids the memory allocation and other expensive start-up cost of + * FEValues. Currently, the functionality is specialized for mappings derived + * from MappingQ and MappingCartesian and for finite elements with tensor + * product structure that work with the + * @ref matrixfree + * module. In those cases, the cost implied + * by this class is similar (or sometimes even somewhat lower) than using + * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`. + */ +template +class FEPointEvaluation + : public FEPointEvaluationBase +{ +public: + static constexpr unsigned int dimension = dim; + static constexpr unsigned int n_components = n_components_; + + using number_type = Number; + + using ScalarNumber = + typename internal::VectorizedArrayTrait::value_type; + using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; + using ETT = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits; + using value_type = typename ETT::value_type; + using scalar_value_type = typename ETT::scalar_value_type; + using vectorized_value_type = typename ETT::vectorized_value_type; + using gradient_type = typename ETT::gradient_type; + using interface_vectorized_gradient_type = + typename ETT::interface_vectorized_gradient_type; + + /** + * Constructor. + * + * @param mapping The Mapping class describing the actual geometry of a cell + * passed to the evaluate() function. + * + * @param fe The FiniteElement object that is used for the evaluation, which + * is typically the same on all cells to be evaluated. + * + * @param update_flags Specify the quantities to be computed by the mapping + * during the call of reinit(). During evaluate() or integrate(), this data + * is queried to produce the desired result (e.g., the gradient of a finite + * element solution). + * + * @param first_selected_component For multi-component FiniteElement + * objects, this parameter allows to select a range of `n_components` + * components starting from this parameter. + */ + FEPointEvaluation(const Mapping &mapping, + const FiniteElement &fe, + const UpdateFlags update_flags, + const unsigned int first_selected_component = 0) + : FEPointEvaluationBase( + mapping, + fe, + update_flags, + first_selected_component) + {} + + /** + * Constructor to make the present class able to re-use the geometry + * data also used by other `FEPointEvaluation` objects. + * + * @param mapping_info The MappingInfo class describes the geometry-related + * data for evaluating finite-element solutions. This object enables to + * construct such an object on the outside, possibly re-using it between + * several objects or between several calls to the same cell and unit points. + * + * @param fe The FiniteElement object that is used for the evaluation, which + * is typically the same on all cells to be evaluated. + * + * @param first_selected_component For multi-component FiniteElement + * objects, this parameter allows to select a range of `n_components` + * components starting from this parameter. + */ + FEPointEvaluation( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const unsigned int first_selected_component = 0) + : FEPointEvaluationBase( + mapping_info, + fe, + first_selected_component) + {} + + /** + * Set up the mapping information for the given cell, e.g., by computing the + * Jacobian of the mapping for the given points if gradients of the functions + * are requested. + * + * @param[in] cell An iterator to the current cell + * + * @param[in] unit_points List of points in the reference locations of the + * current cell where the FiniteElement object should be + * evaluated/integrated in the evaluate() and integrate() functions. + */ + void + reinit(const typename Triangulation::cell_iterator &cell, + const ArrayView> &unit_points); + + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the single cell in the MappingInfo object. + */ + void + reinit(); + + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the cell in the MappingInfo object. + */ + void + reinit(const unsigned int cell_index); + + + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element read out by + * `FEEvaluation::read_dof_values(global_vector)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + template + void + evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element as returned by `cell->get_dof_values(global_vector, + * solution_values)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + void + evaluate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + template + void + integrate(const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used to during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + void + integrate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points. This is + * similar to the integration of a bilinear form in terms of the test + * function, with the difference that this formula does not include a `JxW` + * factor (in contrast to the integrate function of this class). This allows + * the class to naturally embed point information (e.g. particles) into a + * finite element formulation. + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + template + void + test_and_sum( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points. This is + * similar to the integration of a bilinear form in terms of the test + * function, with the difference that this formula does not include a `JxW` + * factor (in contrast to the integrate function of this class). This allows + * the class to naturally embed point information (e.g. particles) into a + * finite element formulation. + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + void + test_and_sum(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * Return the normal vector. This class or the MappingInfo object passed to + * this function needs to be constructed with UpdateFlags containing + * `update_normal_vectors`. + */ + Tensor<1, spacedim, Number> + normal_vector(const unsigned int point_index) const; + +private: + static constexpr std::size_t n_lanes_user_interface = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t n_lanes_internal = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t stride = + internal::VectorizedArrayTrait::stride(); + + /** + * Resizes necessary data fields, reads in and renumbers solution values. + * Interpolates onto face if face path is selected. + */ + template + void + prepare_evaluate_fast( + const StridedArrayView &solution_values); + + /** + * Evaluates the actual interpolation on the cell or face for a quadrature + * batch. + */ + template + void + compute_evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags, + const unsigned int n_shapes, + const unsigned int qb, + vectorized_value_type &value, + interface_vectorized_gradient_type &gradient); + + /** + * Fast path of the evaluate function. + */ + template + void + evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * Slow path of the evaluate function using FEValues. + */ + template + void + evaluate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * Integrates the product of the data passed in by submit_value() and + * submit_gradient() with the values or gradients of test functions on the + * cell or face for a given quadrature batch. + */ + template + void + compute_integrate_fast( + const EvaluationFlags::EvaluationFlags &integration_flags, + const unsigned int n_shapes, + const unsigned int qb, + const vectorized_value_type value, + const interface_vectorized_gradient_type gradient, + vectorized_value_type *solution_values_vectorized_linear); + + /** + * Addition across the lanes of VectorizedArray as accumulated by the + * compute_integrate_fast_function(), writing the sum into the result vector. + * Applies face contributions to cell contributions for face path. + */ + template + void + finish_integrate_fast( + const StridedArrayView &solution_values, + vectorized_value_type *solution_values_vectorized_linear, + const bool sum_into_values); + + /** + * Fast path of the integrate function. + */ + template + void + integrate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); + + /** + * Slow path of the integrate function using FEValues. + */ + template + void + integrate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); + + /** + * Implementation of the integrate/test_and_sum function. + */ + template + void + do_integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); +}; + + + +template +inline void +FEPointEvaluation::reinit() +{ + this->current_cell_index = numbers::invalid_unsigned_int; + this->current_face_number = numbers::invalid_unsigned_int; + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); +} + + + +template +inline void +FEPointEvaluation::reinit( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &unit_points) +{ + // reinit is only allowed for mapping computation on the fly + AssertThrow(this->mapping_info_on_the_fly.get() != nullptr, + ExcNotImplemented()); + + this->mapping_info->reinit(cell, unit_points); + + if (!this->fast_path) + { + this->fe_values = std::make_shared>( + *this->mapping, + *this->fe, + Quadrature( + std::vector>(unit_points.begin(), unit_points.end())), + this->update_flags); + this->fe_values->reinit(cell); + } + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); +} + + + +template +inline void +FEPointEvaluation::reinit( + const unsigned int cell_index) +{ + this->current_cell_index = cell_index; + this->current_face_number = numbers::invalid_unsigned_int; + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); + + if (!this->fast_path) + { + std::vector> unit_points(this->n_q_points_scalar); + + for (unsigned int v = 0; v < this->n_q_points_scalar; ++v) + for (unsigned int d = 0; d < dim; ++d) + unit_points[v][d] = + this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal]; + + this->fe_values = std::make_shared>( + *this->mapping, + *this->fe, + Quadrature( + std::vector>(unit_points.begin(), unit_points.end())), + this->update_flags); + + this->fe_values->reinit( + this->mapping_info->get_cell_iterator(this->current_cell_index)); + } +} + + + +template +template +void +FEPointEvaluation::evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) +{ + if (!this->is_reinitialized) + reinit(); + + if (this->n_q_points == 0) + return; + + Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented()); + + if (!((evaluation_flags & EvaluationFlags::values) || + (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags + return; + + AssertDimension(solution_values.size(), this->fe->dofs_per_cell); + if (this->fast_path) + { + if (this->use_linear_path) + evaluate_fast(solution_values, evaluation_flags); + else + evaluate_fast(solution_values, evaluation_flags); + } + else + evaluate_slow(solution_values, evaluation_flags); +} + + + +template +void +FEPointEvaluation::evaluate( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) +{ + evaluate(StridedArrayView(solution_values.data(), + solution_values.size()), + evaluation_flags); +} + + + +template +template +void +FEPointEvaluation::integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + do_integrate(solution_values, integration_flags, sum_into_values); +} + + + +template +void +FEPointEvaluation::integrate( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + integrate(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags, + sum_into_values); +} + + + +template +template +void +FEPointEvaluation::test_and_sum( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + do_integrate(solution_values, integration_flags, sum_into_values); +} + + + +template +void +FEPointEvaluation::test_and_sum( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + test_and_sum(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags, + sum_into_values); +} + + + +template +template +inline void +FEPointEvaluation::prepare_evaluate_fast( + const StridedArrayView &solution_values) +{ + const unsigned int dofs_per_comp = + is_linear ? Utilities::pow(2, dim) : this->dofs_per_component; + + for (unsigned int comp = 0; comp < n_components; ++comp) + { + const std::size_t offset = + (this->component_in_base_element + comp) * dofs_per_comp; + + if ((is_linear && n_components == 1) || this->renumber.empty()) + { + for (unsigned int i = 0; i < dofs_per_comp; ++i) + ETT::read_value(solution_values[i + offset], + comp, + this->solution_renumbered[i]); + } + else + { + const unsigned int *renumber_ptr = this->renumber.data() + offset; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + ETT::read_value(solution_values[renumber_ptr[i]], + comp, + this->solution_renumbered[i]); + } + } +} + + + +template +template +inline void +FEPointEvaluation::compute_evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags, + const unsigned int n_shapes, + const unsigned int qb, + vectorized_value_type &value, + interface_vectorized_gradient_type &gradient) +{ + if (evaluation_flags & EvaluationFlags::gradients) + { + std::array result; + if constexpr (is_linear) + { + if constexpr (n_components == 1) + result = + internal::evaluate_tensor_product_value_and_gradient_linear< + dim, + scalar_value_type, + VectorizedArrayType, + 1, + stride_view>(solution_values.data(), this->unit_point_ptr[qb]); + else + result = + internal::evaluate_tensor_product_value_and_gradient_linear( + this->solution_renumbered.data(), this->unit_point_ptr[qb]); + } + else + result = internal::evaluate_tensor_product_value_and_gradient_shapes< + dim, + scalar_value_type, + VectorizedArrayType, + 1, + false>(this->shapes.data() + qb * n_shapes, + n_shapes, + this->solution_renumbered.data()); + gradient[0] = result[0]; + if (dim > 1) + gradient[1] = result[1]; + if (dim > 2) + gradient[2] = result[2]; + value = result[dim]; + } + else + { + if constexpr (is_linear) + { + if constexpr (n_components == 1) + value = internal::evaluate_tensor_product_value_linear< + dim, + scalar_value_type, + VectorizedArrayType, + stride_view>(solution_values.data(), this->unit_point_ptr[qb]); + else + value = internal::evaluate_tensor_product_value_linear( + this->solution_renumbered.data(), this->unit_point_ptr[qb]); + } + else + value = + internal::evaluate_tensor_product_value_shapes( + this->shapes.data() + qb * n_shapes, + n_shapes, + this->solution_renumbered.data()); + } +} + + + +template +template +inline void +FEPointEvaluation::evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) +{ + if (!(is_linear && n_components == 1)) + prepare_evaluate_fast(solution_values); + + // loop over quadrature batches qb + const unsigned int n_shapes = is_linear ? 2 : this->poly.size(); + + for (unsigned int qb = 0; qb < this->n_q_batches; ++qb) + { + vectorized_value_type value; + interface_vectorized_gradient_type gradient; + + compute_evaluate_fast( + solution_values, evaluation_flags, n_shapes, qb, value, gradient); + + if (evaluation_flags & EvaluationFlags::values) + { + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + ETT::set_value(value, v, this->values[offset]); + } + if (evaluation_flags & EvaluationFlags::gradients) + { + Assert(this->update_flags & update_gradients || + this->update_flags & update_inverse_jacobians, + ExcNotInitialized()); + + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + { + gradient_type unit_gradient; + ETT::set_gradient(gradient, v, unit_gradient); + this->gradients[offset] = + this->cell_type <= + internal::MatrixFreeFunctions::GeometryType::cartesian ? + apply_diagonal_transformation(this->inverse_jacobian_ptr[0], + unit_gradient) : + apply_transformation( + this + ->inverse_jacobian_ptr[this->cell_type <= + internal::MatrixFreeFunctions:: + GeometryType::affine ? + 0 : + offset] + .transpose(), + unit_gradient); + } + } + } +} + + + +template +template +inline void +FEPointEvaluation::evaluate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) +{ + // slow path with FEValues + Assert(this->fe_values.get() != nullptr, + ExcMessage( + "Not initialized. Please call FEPointEvaluation::reinit()!")); + + const std::size_t n_points = this->fe_values->get_quadrature().size(); + + if (evaluation_flags & EvaluationFlags::values) + { + this->values.resize(this->n_q_points); + std::fill(this->values.begin(), this->values.end(), value_type()); + for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i) + { + const ScalarNumber value = solution_values[i]; + for (unsigned int d = 0; d < n_components; ++d) + if (this->nonzero_shape_function_component[i][d] && + (this->fe->is_primitive(i) || this->fe->is_primitive())) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + ETT::access(this->values[qb], + v, + d, + this->fe_values->shape_value(i, q + v) * value); + else if (this->nonzero_shape_function_component[i][d]) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + ETT::access(this->values[qb], + v, + d, + this->fe_values->shape_value_component(i, + q + v, + d) * + value); + } + } + + if (evaluation_flags & EvaluationFlags::gradients) + { + this->gradients.resize(this->n_q_points); + std::fill(this->gradients.begin(), + this->gradients.end(), + gradient_type()); + for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i) + { + const ScalarNumber value = solution_values[i]; + for (unsigned int d = 0; d < n_components; ++d) + if (this->nonzero_shape_function_component[i][d] && + (this->fe->is_primitive(i) || this->fe->is_primitive())) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + ETT::access(this->gradients[qb], + v, + d, + this->fe_values->shape_grad(i, q + v) * value); + else if (this->nonzero_shape_function_component[i][d]) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + ETT::access( + this->gradients[qb], + v, + d, + this->fe_values->shape_grad_component(i, q + v, d) * value); } } +} + + + +template +template +inline void +FEPointEvaluation::compute_integrate_fast( + const EvaluationFlags::EvaluationFlags &integration_flags, + const unsigned int n_shapes, + const unsigned int qb, + const vectorized_value_type value, + const interface_vectorized_gradient_type gradient, + vectorized_value_type *solution_values_vectorized_linear) +{ + if (integration_flags & EvaluationFlags::gradients) + internal::integrate_tensor_product_value_and_gradient< + is_linear, + dim, + VectorizedArrayType, + vectorized_value_type>(this->shapes.data() + qb * n_shapes, + n_shapes, + &value, + gradient, + is_linear ? + solution_values_vectorized_linear : + this->solution_renumbered_vectorized.data(), + this->unit_point_ptr[qb], + qb != 0); else + internal::integrate_tensor_product_value( + this->shapes.data() + qb * n_shapes, + n_shapes, + value, + is_linear ? solution_values_vectorized_linear : + this->solution_renumbered_vectorized.data(), + this->unit_point_ptr[qb], + qb != 0); +} + + + +template +template +inline void +FEPointEvaluation::finish_integrate_fast( + const StridedArrayView &solution_values, + vectorized_value_type *solution_values_vectorized_linear, + const bool sum_into_values) +{ + if (!sum_into_values && this->fe->n_components() > n_components) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; + + const unsigned int dofs_per_comp = + is_linear ? Utilities::pow(2, dim) : this->dofs_per_component; + + for (unsigned int comp = 0; comp < n_components; ++comp) { - if (evaluation_flags & EvaluationFlags::gradients) + const std::size_t offset = + (this->component_in_base_element + comp) * dofs_per_comp; + + if (is_linear || this->renumber.empty()) { - std::array result; - if constexpr (is_linear) - { - if constexpr (n_components == 1) - result = - internal::evaluate_tensor_product_value_and_gradient_linear< - dim, - scalar_value_type, - VectorizedArrayType, - 1, - stride_view>(solution_values.data(), unit_point_ptr[qb]); - else - result = - internal::evaluate_tensor_product_value_and_gradient_linear( - solution_renumbered.data(), unit_point_ptr[qb]); - } - else - result = - internal::evaluate_tensor_product_value_and_gradient_shapes< - dim, - scalar_value_type, - VectorizedArrayType, - 1, - false>(shapes.data() + qb * n_shapes, - n_shapes, - solution_renumbered.data()); - gradient[0] = result[0]; - if (dim > 1) - gradient[1] = result[1]; - if (dim > 2) - gradient[2] = result[2]; - value = result[dim]; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + if (sum_into_values) + solution_values[i + offset] += + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear + i) : + this->solution_renumbered_vectorized[i]); + else + solution_values[i + offset] = + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear + i) : + this->solution_renumbered_vectorized[i]); } else { - if constexpr (is_linear) - { - if constexpr (n_components == 1) - value = internal::evaluate_tensor_product_value_linear< - dim, - scalar_value_type, - VectorizedArrayType, - stride_view>(solution_values.data(), unit_point_ptr[qb]); - else - value = internal::evaluate_tensor_product_value_linear( - solution_renumbered.data(), unit_point_ptr[qb]); - } - else - value = internal::evaluate_tensor_product_value_shapes< - dim, - scalar_value_type, - VectorizedArrayType, - false>(shapes.data() + qb * n_shapes, - n_shapes, - solution_renumbered.data()); + const unsigned int *renumber_ptr = this->renumber.data() + offset; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + if (sum_into_values) + solution_values[renumber_ptr[i]] += + ETT::sum_value(comp, this->solution_renumbered_vectorized[i]); + else + solution_values[renumber_ptr[i]] = + ETT::sum_value(comp, this->solution_renumbered_vectorized[i]); + } + } +} + + + +template +template +inline void +FEPointEvaluation::integrate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + // zero out lanes of incomplete last quadrature point batch + if constexpr (stride == 1) + if (const unsigned int n_filled_lanes = + this->n_q_points_scalar & (n_lanes_internal - 1); + n_filled_lanes > 0) + { + if (integration_flags & EvaluationFlags::values) + for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) + ETT::set_zero_value(this->values.back(), v); + if (integration_flags & EvaluationFlags::gradients) + for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) + ETT::set_zero_gradient(this->gradients.back(), v); + } + + std::array + solution_values_vectorized_linear = {}; + + // loop over quadrature batches qb + const unsigned int n_shapes = is_linear ? 2 : this->poly.size(); + + const bool cartesian_cell = + this->cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian; + const bool affine_cell = + this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine; + for (unsigned int qb = 0; qb < this->n_q_batches; ++qb) + { + vectorized_value_type value = {}; + Tensor<1, dim, vectorized_value_type> gradient; + + if (integration_flags & EvaluationFlags::values) + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + ETT::get_value(value, + v, + do_JxW ? this->values[offset] * this->JxW_ptr[offset] : + this->values[offset]); + + if (integration_flags & EvaluationFlags::gradients) + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + { + const auto grad_w = + do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] : + this->gradients[offset]; + ETT::get_gradient( + gradient, + v, + cartesian_cell ? + apply_diagonal_transformation(this->inverse_jacobian_ptr[0], + grad_w) : + apply_transformation( + this->inverse_jacobian_ptr[affine_cell ? 0 : offset], + grad_w)); + } + + compute_integrate_fast( + integration_flags, + n_shapes, + qb, + value, + gradient, + solution_values_vectorized_linear.data()); + } + + // add between the lanes and write into the result + finish_integrate_fast(solution_values, + solution_values_vectorized_linear.data(), + sum_into_values); +} + + + +template +template +inline void +FEPointEvaluation::integrate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + // slow path with FEValues + Assert(this->fe_values.get() != nullptr, + ExcMessage( + "Not initialized. Please call FEPointEvaluation::reinit()!")); + if (!sum_into_values) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; + + const std::size_t n_points = this->fe_values->get_quadrature().size(); + + if (integration_flags & EvaluationFlags::values) + { + AssertIndexRange(this->n_q_points, this->values.size() + 1); + for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i) + { + for (unsigned int d = 0; d < n_components; ++d) + if (this->nonzero_shape_function_component[i][d] && + (this->fe->is_primitive(i) || this->fe->is_primitive())) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + solution_values[i] += + this->fe_values->shape_value(i, q + v) * + ETT::access(this->values[qb], v, d) * + (do_JxW ? this->fe_values->JxW(q + v) : 1.); + else if (this->nonzero_shape_function_component[i][d]) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + solution_values[i] += + this->fe_values->shape_value_component(i, q + v, d) * + ETT::access(this->values[qb], v, d) * + (do_JxW ? this->fe_values->JxW(q + v) : 1.); + } + } + + if (integration_flags & EvaluationFlags::gradients) + { + AssertIndexRange(this->n_q_points, this->gradients.size() + 1); + for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i) + { + for (unsigned int d = 0; d < n_components; ++d) + if (this->nonzero_shape_function_component[i][d] && + (this->fe->is_primitive(i) || this->fe->is_primitive())) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + solution_values[i] += + this->fe_values->shape_grad(i, q + v) * + ETT::access(this->gradients[qb], v, d) * + (do_JxW ? this->fe_values->JxW(q + v) : 1.); + else if (this->nonzero_shape_function_component[i][d]) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) + for (unsigned int v = 0; + v < n_lanes_user_interface && q + v < n_points; + ++v) + solution_values[i] += + this->fe_values->shape_grad_component(i, q + v, d) * + ETT::access(this->gradients[qb], v, d) * + (do_JxW ? this->fe_values->JxW(q + v) : 1.); } } } @@ -2092,145 +2702,396 @@ FEPointEvaluation::compute_evaluate_fast( template -template -inline void -FEPointEvaluation::evaluate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) +template +void +FEPointEvaluation::do_integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + if (!this->is_reinitialized) + reinit(); + + Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented()); + + if (this->n_q_points == 0 || // no evaluation points provided + !((integration_flags & EvaluationFlags::values) || + (integration_flags & + EvaluationFlags::gradients))) // no integration flags + { + if (!sum_into_values) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; + return; + } + + Assert( + !do_JxW || this->JxW_ptr != nullptr, + ExcMessage( + "JxW pointer is not set! If you do not want to integrate() use test_and_sum()")); + + AssertDimension(solution_values.size(), this->fe->dofs_per_cell); + if (this->fast_path) + { + if (this->use_linear_path) + integrate_fast(solution_values, + integration_flags, + sum_into_values); + else + integrate_fast(solution_values, + integration_flags, + sum_into_values); + } + else + integrate_slow(solution_values, integration_flags, sum_into_values); +} + + + +template +inline Tensor<1, spacedim, Number> +FEPointEvaluation::normal_vector( + const unsigned int point_index) const { - if (!(is_linear && n_components == 1) || is_face_path) - prepare_evaluate_fast(solution_values, - evaluation_flags); + AssertIndexRange(point_index, this->n_q_points); + Assert(this->normal_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_normal_vectors")); + if (this->is_interior) + return this->normal_ptr[point_index]; + else + return -this->normal_ptr[point_index]; +} - // loop over quadrature batches qb - const unsigned int n_shapes = is_linear ? 2 : poly.size(); - for (unsigned int qb = 0; qb < n_q_batches; ++qb) - { - vectorized_value_type value; - interface_vectorized_gradient_type gradient; - compute_evaluate_fast( - solution_values, evaluation_flags, n_shapes, qb, value, gradient); +/** + * This class provides an interface to the evaluation of interpolated solution + * values and gradients on faces on arbitrary reference point positions. These + * points can change from face to face, both with respect to their quantity as + * well to the location. A typical use case is evaluations on non-matching + * grids. + * + * The use of this class is similar to FEEvaluation: In the constructor, a + * reference to a NonMatching::MappingInfo object is passed, where the + * quadrature points in reference position is stored together with the mapping + * information. The class is then reinitialized to a cell by calling + * `FEFacePointEvaluation::reinit(face_index)` or + * `FEFacePointEvaluation::reinit(cell_index, face_number)`. Then, upon call to + * evaluate() or integrate(), the user can compute information at the given + * points. Eventually, the access functions get_value() or get_gradient() allow + * to query this information at a specific point index. + */ +template +class FEFacePointEvaluation + : public FEPointEvaluationBase +{ +public: + static constexpr unsigned int dimension = dim; + static constexpr unsigned int n_components = n_components_; - if (evaluation_flags & EvaluationFlags::values) - { - for (unsigned int v = 0, offset = qb * stride; - v < stride && (stride == 1 || offset < n_q_points_scalar); - ++v, ++offset) - ETT::set_value(value, v, values[offset]); - } - if (evaluation_flags & EvaluationFlags::gradients) - { - Assert(update_flags & update_gradients || - update_flags & update_inverse_jacobians, - ExcNotInitialized()); + using number_type = Number; - for (unsigned int v = 0, offset = qb * stride; - v < stride && (stride == 1 || offset < n_q_points_scalar); - ++v, ++offset) - { - gradient_type unit_gradient; - ETT::set_gradient(gradient, v, unit_gradient); - gradients[offset] = - cell_type <= - internal::MatrixFreeFunctions::GeometryType::cartesian ? - apply_diagonal_transformation(inverse_jacobian_ptr[0], - unit_gradient) : - apply_transformation( - inverse_jacobian_ptr - [cell_type <= - internal::MatrixFreeFunctions::GeometryType::affine ? - 0 : - offset] - .transpose(), - unit_gradient); - } - } - } -} + using ScalarNumber = + typename internal::VectorizedArrayTrait::value_type; + using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; + using ETT = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits; + using value_type = typename ETT::value_type; + using scalar_value_type = typename ETT::scalar_value_type; + using vectorized_value_type = typename ETT::vectorized_value_type; + using gradient_type = typename ETT::gradient_type; + using interface_vectorized_gradient_type = + typename ETT::interface_vectorized_gradient_type; + + /** + * Constructor. Allows to select if interior or exterior face is selected. + */ + FEFacePointEvaluation( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const bool is_interior = true, + const unsigned int first_selected_component = 0); + + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the face in the MappingInfo object. Used in element-centric loops (ECL). + */ + void + reinit(const unsigned int cell_index, const unsigned int face_number); + + /** + * Reinitialize the evaluator to point to the correct precomputed mapping of + * the face in the MappingInfo object. Used in face-centric loops (FCL). + */ + void + reinit(const unsigned int face_index); + + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element read out by + * `FEEvaluation::read_dof_values(global_vector)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + template + void + evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element as returned by `cell->get_dof_values(global_vector, + * solution_values)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + void + evaluate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + template + void + integrate(const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used to during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + void + integrate(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + template + void + test_and_sum( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used to during + * `cell->set_dof_values(solution_values, global_vector)` or + * `cell->distribute_local_to_global(solution_values, global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + * @param[in] sum_into_values Flag specifying if the integrated values + * should be summed into the solution values. Defaults to false. + * + */ + void + test_and_sum(const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * Evaluate values and gradients in face for the selected face (lane) of the + * batch. Default stride into the face dofs is width of + * VectorizedArray which is the default + * vectorization over faces for FEFaceEvaluation. + */ + template + void + evaluate_in_face(const ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * Integrate values and gradients in face for the selected face (lane) of the + * batch. Default stride into the face dofs is width of + * VectorizedArray which is the default + * vectorization over faces for FEFaceEvaluation. + */ + template + void + integrate_in_face(ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values = false); + + /** + * Return the normal vector. This class or the MappingInfo object passed to + * this function needs to be constructed with UpdateFlags containing + * `update_normal_vectors`. + */ + Tensor<1, spacedim, Number> + normal_vector(const unsigned int point_index) const; + +private: + static constexpr std::size_t n_lanes_user_interface = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t n_lanes_internal = + internal::VectorizedArrayTrait::width(); + static constexpr std::size_t stride = + internal::VectorizedArrayTrait::stride(); + + template + void + do_evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + template + void + do_integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); + + /** + * Actually does the evaluation templated on the chosen code path (linear or + * higher order). + */ + template + void + do_evaluate_in_face(const ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + + /** + * Actually does the integration templated on the chosen code path (linear or + * higher order). + */ + template + void + do_integrate_in_face( + ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values); +}; + + + +template +FEFacePointEvaluation:: + FEFacePointEvaluation( + NonMatching::MappingInfo &mapping_info, + const FiniteElement &fe, + const bool is_interior, + const unsigned int first_selected_component) + : FEPointEvaluationBase( + mapping_info, + fe, + first_selected_component, + is_interior) +{} template -template inline void -FEPointEvaluation::evaluate_slow( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) +FEFacePointEvaluation::reinit( + const unsigned int cell_index, + const unsigned int face_number) { - // slow path with FEValues - Assert(fe_values.get() != nullptr, - ExcMessage( - "Not initialized. Please call FEPointEvaluation::reinit()!")); + this->current_cell_index = cell_index; + this->current_face_number = face_number; - const std::size_t n_points = fe_values->get_quadrature().size(); + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); +} - if (evaluation_flags & EvaluationFlags::values) - { - values.resize(n_q_points); - std::fill(values.begin(), values.end(), value_type()); - for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) - { - const ScalarNumber value = solution_values[i]; - for (unsigned int d = 0; d < n_components; ++d) - if (nonzero_shape_function_component[i][d] && - (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - ETT::access(values[qb], - v, - d, - fe_values->shape_value(i, q + v) * value); - else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - ETT::access(values[qb], - v, - d, - fe_values->shape_value_component(i, q + v, d) * - value); - } - } - if (evaluation_flags & EvaluationFlags::gradients) - { - gradients.resize(n_q_points); - std::fill(gradients.begin(), gradients.end(), gradient_type()); - for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) - { - const ScalarNumber value = solution_values[i]; - for (unsigned int d = 0; d < n_components; ++d) - if (nonzero_shape_function_component[i][d] && - (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - ETT::access(gradients[qb], - v, - d, - fe_values->shape_grad(i, q + v) * value); - else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - ETT::access(gradients[qb], - v, - d, - fe_values->shape_grad_component(i, q + v, d) * - value); - } - } + +template +inline void +FEFacePointEvaluation::reinit( + const unsigned int face_index) +{ + this->current_cell_index = face_index; + this->current_face_number = + this->mapping_info->get_face_number(face_index, this->is_interior); + + if (this->use_linear_path) + this->template do_reinit(); + else + this->template do_reinit(); } @@ -2238,14 +3099,13 @@ FEPointEvaluation::evaluate_slow( template template void -FEPointEvaluation::evaluate( +FEFacePointEvaluation::evaluate( const StridedArrayView &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags) { - if (!is_reinitialized) - reinit(); + Assert(this->is_reinitialized, ExcMessage("Is not reinitalized!")); - if (n_q_points == 0) + if (this->n_q_points == 0) return; Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented()); @@ -2254,33 +3114,19 @@ FEPointEvaluation::evaluate( (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags return; - AssertDimension(solution_values.size(), fe->dofs_per_cell); - if (fast_path) - { - if (use_face_path) - { - if (use_linear_path) - evaluate_fast(solution_values, evaluation_flags); - else - evaluate_fast(solution_values, evaluation_flags); - } - else - { - if (use_linear_path) - evaluate_fast(solution_values, evaluation_flags); - else - evaluate_fast(solution_values, evaluation_flags); - } - } + AssertDimension(solution_values.size(), this->fe->dofs_per_cell); + + if (this->use_linear_path) + do_evaluate(solution_values, evaluation_flags); else - evaluate_slow(solution_values, evaluation_flags); + do_evaluate(solution_values, evaluation_flags); } template void -FEPointEvaluation::evaluate( +FEFacePointEvaluation::evaluate( const ArrayView &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags) { @@ -2292,400 +3138,124 @@ FEPointEvaluation::evaluate( template -template -inline void -FEPointEvaluation::compute_integrate_fast( - const EvaluationFlags::EvaluationFlags &integration_flags, - const unsigned int n_shapes, - const unsigned int qb, - const vectorized_value_type value, - const interface_vectorized_gradient_type gradient, - vectorized_value_type *solution_values_vectorized_linear) -{ - if (is_face_path) - { - if (integration_flags & EvaluationFlags::gradients) - { - std::array value_face = {}; - Tensor<1, dim - 1, vectorized_value_type> gradient_in_face; - - value_face[0] = value; - // fill derivative in physical coordinates into tangential/normal - // derivatives - if (current_face_number / 2 == 0) - { - value_face[1] = gradient[0]; - if (dim > 1) - gradient_in_face[0] = gradient[1]; - if (dim > 2) - gradient_in_face[1] = gradient[2]; - } - else if (current_face_number / 2 == 1) - { - if (dim > 1) - value_face[1] = gradient[1]; - if (dim == 3) - { - gradient_in_face[0] = gradient[2]; - gradient_in_face[1] = gradient[0]; - } - else if (dim == 2) - gradient_in_face[0] = gradient[0]; - else - Assert(false, ExcInternalError()); - } - else if (current_face_number / 2 == 2) - { - if (dim > 2) - { - value_face[1] = gradient[2]; - gradient_in_face[0] = gradient[0]; - gradient_in_face[1] = gradient[1]; - } - else - Assert(false, ExcInternalError()); - } - else - Assert(false, ExcInternalError()); - - internal::integrate_tensor_product_value_and_gradient< - is_linear, - dim - 1, - VectorizedArrayType, - vectorized_value_type, - 2>(shapes_faces.data() + qb * n_shapes, - n_shapes, - value_face.data(), - gradient_in_face, - is_linear ? solution_values_vectorized_linear : - solution_renumbered_vectorized.data(), - unit_point_faces_ptr[qb], - qb != 0); - } - else - internal::integrate_tensor_product_value( - shapes_faces.data() + qb * n_shapes, - n_shapes, - value, - is_linear ? solution_values_vectorized_linear : - solution_renumbered_vectorized.data(), - unit_point_faces_ptr[qb], - qb != 0); - } - else - { - if (integration_flags & EvaluationFlags::gradients) - internal::integrate_tensor_product_value_and_gradient< - is_linear, - dim, - VectorizedArrayType, - vectorized_value_type>(shapes.data() + qb * n_shapes, - n_shapes, - &value, - gradient, - is_linear ? - solution_values_vectorized_linear : - solution_renumbered_vectorized.data(), - unit_point_ptr[qb], - qb != 0); - else - internal::integrate_tensor_product_value( - shapes.data() + qb * n_shapes, - n_shapes, - value, - is_linear ? solution_values_vectorized_linear : - solution_renumbered_vectorized.data(), - unit_point_ptr[qb], - qb != 0); - } -} - - - -template -template -inline void -FEPointEvaluation::finish_integrate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - vectorized_value_type *solution_values_vectorized_linear, - const bool sum_into_values) +template +void +FEFacePointEvaluation::do_evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) { - if (!sum_into_values && fe->n_components() > n_components) - for (unsigned int i = 0; i < solution_values.size(); ++i) - solution_values[i] = 0; - const unsigned int dofs_per_comp = - is_linear ? Utilities::pow(2, dim) : dofs_per_component; - - for (unsigned int comp = 0; comp < n_components; ++comp) - { - const std::size_t offset = - (component_in_base_element + comp) * dofs_per_comp; - - if (is_face_path) - { - const unsigned int dofs_per_comp_face = - is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face; - - const unsigned int size_input = 2 * dofs_per_comp_face; - ScalarNumber *input = scratch_data_scalar.begin(); - ScalarNumber *output = input + size_input; - - for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i) - input[i] = - ETT::sum_value(comp, - is_linear ? - *(solution_values_vectorized_linear + i) : - solution_renumbered_vectorized[i]); - - internal::FEFaceNormalEvaluationImpl:: - template interpolate(1, - integration_flags, - shape_info, - input, - output, - current_face_number); - - if (is_linear || renumber.empty()) - { - for (unsigned int i = 0; i < dofs_per_comp; ++i) - if (sum_into_values) - solution_values[i + offset] += output[i]; - else - solution_values[i + offset] = output[i]; - } - else - { - const unsigned int *renumber_ptr = renumber.data() + offset; - for (unsigned int i = 0; i < dofs_per_comp; ++i) - if (sum_into_values) - solution_values[renumber_ptr[i]] += output[i]; - else - solution_values[renumber_ptr[i]] = output[i]; - } - } - else - { - if (is_linear || renumber.empty()) - { - for (unsigned int i = 0; i < dofs_per_comp; ++i) - if (sum_into_values) - solution_values[i + offset] += - ETT::sum_value(comp, - is_linear ? - *(solution_values_vectorized_linear + i) : - solution_renumbered_vectorized[i]); - else - solution_values[i + offset] = - ETT::sum_value(comp, - is_linear ? - *(solution_values_vectorized_linear + i) : - solution_renumbered_vectorized[i]); - } - else - { - const unsigned int *renumber_ptr = renumber.data() + offset; - for (unsigned int i = 0; i < dofs_per_comp; ++i) - if (sum_into_values) - solution_values[renumber_ptr[i]] += - ETT::sum_value(comp, solution_renumbered_vectorized[i]); - else - solution_values[renumber_ptr[i]] = - ETT::sum_value(comp, solution_renumbered_vectorized[i]); - } - } - } -} - - - -template -template -inline void -FEPointEvaluation::integrate_fast( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values) -{ - // zero out lanes of incomplete last quadrature point batch - if constexpr (stride == 1) - if (const unsigned int n_filled_lanes = - n_q_points_scalar & (n_lanes_internal - 1); - n_filled_lanes > 0) - { - if (integration_flags & EvaluationFlags::values) - for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) - ETT::set_zero_value(values.back(), v); - if (integration_flags & EvaluationFlags::gradients) - for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) - ETT::set_zero_gradient(gradients.back(), v); - } - - std::array - solution_values_vectorized_linear = {}; - - // loop over quadrature batches qb - const unsigned int n_shapes = is_linear ? 2 : poly.size(); + is_linear ? Utilities::pow(2, dim) : this->dofs_per_component; - const bool cartesian_cell = - cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian; - const bool affine_cell = - cell_type <= internal::MatrixFreeFunctions::GeometryType::affine; - for (unsigned int qb = 0; qb < n_q_batches; ++qb) + const ScalarNumber *input; + if (stride_view == 1 && this->component_in_base_element == 0 && + (is_linear || this->renumber.empty())) + input = solution_values.data(); + else { - vectorized_value_type value = {}; - Tensor<1, dim, vectorized_value_type> gradient; + for (unsigned int comp = 0; comp < n_components; ++comp) + { + const std::size_t offset = + (this->component_in_base_element + comp) * dofs_per_comp; - if (integration_flags & EvaluationFlags::values) - for (unsigned int v = 0, offset = qb * stride; - v < stride && (stride == 1 || offset < n_q_points_scalar); - ++v, ++offset) - ETT::get_value(value, - v, - do_JxW ? values[offset] * JxW_ptr[offset] : - values[offset]); + if (is_linear || this->renumber.empty()) + { + for (unsigned int i = 0; i < dofs_per_comp; ++i) + this->scratch_data_scalar[i + comp * dofs_per_comp] = + solution_values[i + offset]; + } + else + { + const unsigned int *renumber_ptr = this->renumber.data() + offset; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + this->scratch_data_scalar[i + comp * dofs_per_comp] = + solution_values[renumber_ptr[i]]; + } + } + input = this->scratch_data_scalar.data(); + } - if (integration_flags & EvaluationFlags::gradients) - for (unsigned int v = 0, offset = qb * stride; - v < stride && (stride == 1 || offset < n_q_points_scalar); - ++v, ++offset) - { - const auto grad_w = - do_JxW ? gradients[offset] * JxW_ptr[offset] : gradients[offset]; - ETT::get_gradient( - gradient, - v, - cartesian_cell ? - apply_diagonal_transformation(inverse_jacobian_ptr[0], grad_w) : - apply_transformation( - inverse_jacobian_ptr[affine_cell ? 0 : offset], grad_w)); - } + ScalarNumber *output = + this->scratch_data_scalar.begin() + dofs_per_comp * n_components; - compute_integrate_fast( - integration_flags, - n_shapes, - qb, - value, - gradient, - solution_values_vectorized_linear.data()); - } + internal::FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + evaluation_flags, + this->shape_info, + input, + output, + this->current_face_number); - // add between the lanes and write into the result - finish_integrate_fast( - solution_values, - integration_flags, - solution_values_vectorized_linear.data(), - sum_into_values); + do_evaluate_in_face(output, evaluation_flags); } template -template -inline void -FEPointEvaluation::integrate_slow( +template +void +FEFacePointEvaluation::integrate( const StridedArrayView &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values) { - // slow path with FEValues - Assert(fe_values.get() != nullptr, - ExcMessage( - "Not initialized. Please call FEPointEvaluation::reinit()!")); - if (!sum_into_values) - for (unsigned int i = 0; i < solution_values.size(); ++i) - solution_values[i] = 0; + Assert(this->is_reinitialized, ExcMessage("Is not reinitalized!")); - const std::size_t n_points = fe_values->get_quadrature().size(); + Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented()); - if (integration_flags & EvaluationFlags::values) + if (this->n_q_points == 0 || // no evaluation points provided + !((integration_flags & EvaluationFlags::values) || + (integration_flags & + EvaluationFlags::gradients))) // no integration flags { - AssertIndexRange(n_q_points, values.size() + 1); - for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) - { - for (unsigned int d = 0; d < n_components; ++d) - if (nonzero_shape_function_component[i][d] && - (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - solution_values[i] += fe_values->shape_value(i, q + v) * - ETT::access(values[qb], v, d) * - (do_JxW ? fe_values->JxW(q + v) : 1.); - else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - solution_values[i] += - fe_values->shape_value_component(i, q + v, d) * - ETT::access(values[qb], v, d) * - (do_JxW ? fe_values->JxW(q + v) : 1.); - } + if (!sum_into_values) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; + return; } - if (integration_flags & EvaluationFlags::gradients) - { - AssertIndexRange(n_q_points, gradients.size() + 1); - for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) - { - for (unsigned int d = 0; d < n_components; ++d) - if (nonzero_shape_function_component[i][d] && - (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - solution_values[i] += fe_values->shape_grad(i, q + v) * - ETT::access(gradients[qb], v, d) * - (do_JxW ? fe_values->JxW(q + v) : 1.); - else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; - ++qb, q += n_lanes_user_interface) - for (unsigned int v = 0; - v < n_lanes_user_interface && q + v < n_points; - ++v) - solution_values[i] += - fe_values->shape_grad_component(i, q + v, d) * - ETT::access(gradients[qb], v, d) * - (do_JxW ? fe_values->JxW(q + v) : 1.); - } - } + AssertDimension(solution_values.size(), this->fe->dofs_per_cell); + + if (this->use_linear_path) + do_integrate(solution_values, + integration_flags, + sum_into_values); + else + do_integrate(solution_values, + integration_flags, + sum_into_values); } template -template void -FEPointEvaluation::do_integrate( +FEFacePointEvaluation::integrate( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + integrate(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags, + sum_into_values); +} + + + +template +template +void +FEFacePointEvaluation::test_and_sum( const StridedArrayView &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values) { - if (!is_reinitialized) - reinit(); + Assert(this->is_reinitialized, ExcMessage("Is not reinitalized!")); Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented()); - if (n_q_points == 0 || // no evaluation points provided + if (this->n_q_points == 0 || // no evaluation points provided !((integration_flags & EvaluationFlags::values) || (integration_flags & EvaluationFlags::gradients))) // no integration flags @@ -2696,251 +3266,505 @@ FEPointEvaluation::do_integrate( return; } - Assert( - !do_JxW || JxW_ptr != nullptr, - ExcMessage( - "JxW pointer is not set! If you do not want to integrate() use test_and_sum()")); + AssertDimension(solution_values.size(), this->fe->dofs_per_cell); - AssertDimension(solution_values.size(), fe->dofs_per_cell); - if (fast_path) - { - if (use_face_path) - { - if (use_linear_path) - integrate_fast(solution_values, - integration_flags, - sum_into_values); - else - integrate_fast(solution_values, - integration_flags, - sum_into_values); - } - else - { - if (use_linear_path) - integrate_fast(solution_values, - integration_flags, - sum_into_values); - else - integrate_fast(solution_values, - integration_flags, - sum_into_values); - } - } + if (this->use_linear_path) + do_integrate(solution_values, + integration_flags, + sum_into_values); else - integrate_slow(solution_values, integration_flags, sum_into_values); + do_integrate(solution_values, + integration_flags, + sum_into_values); } template -template void -FEPointEvaluation::integrate( +FEFacePointEvaluation::test_and_sum( + const ArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) +{ + test_and_sum(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags, + sum_into_values); +} + + + +template +template +void +FEFacePointEvaluation::do_integrate( const StridedArrayView &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values) { - do_integrate(solution_values, integration_flags, sum_into_values); + if (!sum_into_values && this->fe->n_components() > n_components) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; + + do_integrate_in_face(this->scratch_data_scalar.begin(), + integration_flags, + false); + + ScalarNumber *input = this->scratch_data_scalar.begin(); + + if (stride_view == 1 && this->component_in_base_element == 0 && + (is_linear || this->renumber.empty())) + { + if (sum_into_values) + internal:: + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flags, + this->shape_info, + input, + solution_values.data(), + this->current_face_number); + else + internal:: + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flags, + this->shape_info, + input, + solution_values.data(), + this->current_face_number); + } + else + { + const unsigned int dofs_per_comp_face = + is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face; + + const unsigned int size_input = 3 * dofs_per_comp_face * n_components; + ScalarNumber *output = input + size_input; + + internal:: + FEFaceNormalEvaluationImpl:: + template interpolate(n_components, + integration_flags, + this->shape_info, + input, + output, + this->current_face_number); + + const unsigned int dofs_per_comp = + is_linear ? Utilities::pow(2, dim) : this->dofs_per_component; + + for (unsigned int comp = 0; comp < n_components; ++comp) + { + const std::size_t offset = + (this->component_in_base_element + comp) * dofs_per_comp; + + if (is_linear || this->renumber.empty()) + { + for (unsigned int i = 0; i < dofs_per_comp; ++i) + if (sum_into_values) + solution_values[i + offset] += + output[i + comp * dofs_per_comp]; + else + solution_values[i + offset] = + output[i + comp * dofs_per_comp]; + } + else + { + const unsigned int *renumber_ptr = this->renumber.data() + offset; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + if (sum_into_values) + solution_values[renumber_ptr[i]] += + output[i + comp * dofs_per_comp]; + else + solution_values[renumber_ptr[i]] = + output[i + comp * dofs_per_comp]; + } + } + } } template +template void -FEPointEvaluation::integrate( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values) +FEFacePointEvaluation::evaluate_in_face( + const ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) { - integrate(StridedArrayView(solution_values.data(), - solution_values.size()), - integration_flags, - sum_into_values); + if (this->use_linear_path) + do_evaluate_in_face(face_dof_values, + evaluation_flags); + else + do_evaluate_in_face(face_dof_values, + evaluation_flags); } -template -template -void -FEPointEvaluation::test_and_sum( - const StridedArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags, - const bool sum_into_values) -{ - do_integrate(solution_values, integration_flags, sum_into_values); +template +template +inline void +FEFacePointEvaluation:: + do_evaluate_in_face(const ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) +{ + const scalar_value_type *face_dof_values_ptr; + if constexpr (n_components == 1) + face_dof_values_ptr = face_dof_values; + else + { + const unsigned int dofs_per_comp_face = + is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face; + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i) + ETT::read_value(face_dof_values[i + comp * 3 * dofs_per_comp_face], + comp, + this->solution_renumbered[i]); + + face_dof_values_ptr = this->solution_renumbered.data(); + } + + // loop over quadrature batches qb + const unsigned int n_shapes = is_linear ? 2 : this->poly.size(); + + for (unsigned int qb = 0; qb < this->n_q_batches; ++qb) + { + vectorized_value_type value; + interface_vectorized_gradient_type gradient; + + if (evaluation_flags & EvaluationFlags::gradients) + { + const std::array interpolated_value = + is_linear ? + internal::evaluate_tensor_product_value_and_gradient_linear< + dim - 1, + scalar_value_type, + VectorizedArrayType, + 2, + stride_face_dof>(face_dof_values_ptr, + this->unit_point_faces_ptr[qb]) : + internal::evaluate_tensor_product_value_and_gradient_shapes< + dim - 1, + scalar_value_type, + VectorizedArrayType, + 2, + false, + stride_face_dof>(this->shapes_faces.data() + qb * n_shapes, + n_shapes, + face_dof_values_ptr); + + value = interpolated_value[dim - 1]; + // reorder derivative from tangential/normal derivatives into tensor + // in physical coordinates + if (this->current_face_number / 2 == 0) + { + gradient[0] = interpolated_value[dim]; + if (dim > 1) + gradient[1] = interpolated_value[0]; + if (dim > 2) + gradient[2] = interpolated_value[1]; + } + else if (this->current_face_number / 2 == 1) + { + if (dim > 1) + gradient[1] = interpolated_value[dim]; + if (dim == 3) + { + gradient[0] = interpolated_value[1]; + gradient[2] = interpolated_value[0]; + } + else if (dim == 2) + gradient[0] = interpolated_value[0]; + else + Assert(false, ExcInternalError()); + } + else if (this->current_face_number / 2 == 2) + { + if (dim > 2) + { + gradient[0] = interpolated_value[0]; + gradient[1] = interpolated_value[1]; + gradient[2] = interpolated_value[dim]; + } + else + Assert(false, ExcInternalError()); + } + else + Assert(false, ExcInternalError()); + } + else + { + value = + is_linear ? + internal::evaluate_tensor_product_value_linear< + dim - 1, + scalar_value_type, + VectorizedArrayType, + stride_face_dof>(face_dof_values_ptr, + this->unit_point_faces_ptr[qb]) : + internal::evaluate_tensor_product_value_shapes< + dim - 1, + scalar_value_type, + VectorizedArrayType, + false, + stride_face_dof>(this->shapes_faces.data() + qb * n_shapes, + n_shapes, + face_dof_values_ptr); + } + + if (evaluation_flags & EvaluationFlags::values) + { + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + ETT::set_value(value, v, this->values[offset]); + } + if (evaluation_flags & EvaluationFlags::gradients) + { + Assert(this->update_flags & update_gradients || + this->update_flags & update_inverse_jacobians, + ExcNotInitialized()); + + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + { + gradient_type unit_gradient; + ETT::set_gradient(gradient, v, unit_gradient); + this->gradients[offset] = + this->cell_type <= + internal::MatrixFreeFunctions::GeometryType::cartesian ? + apply_diagonal_transformation(this->inverse_jacobian_ptr[0], + unit_gradient) : + apply_transformation( + this + ->inverse_jacobian_ptr[this->cell_type <= + internal::MatrixFreeFunctions:: + GeometryType::affine ? + 0 : + offset] + .transpose(), + unit_gradient); + } + } + } } template +template void -FEPointEvaluation::test_and_sum( - const ArrayView &solution_values, +FEFacePointEvaluation::integrate_in_face( + ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values) { - test_and_sum(StridedArrayView(solution_values.data(), - solution_values.size()), - integration_flags, - sum_into_values); -} - - - -template -inline const typename FEPointEvaluation:: - value_type & - FEPointEvaluation::get_value( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, values.size()); - return values[point_index]; -} - - - -template -inline const typename FEPointEvaluation:: - gradient_type & - FEPointEvaluation::get_gradient( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, gradients.size()); - return gradients[point_index]; + if (this->use_linear_path) + do_integrate_in_face(face_dof_values, + integration_flags, + sum_into_values); + else + do_integrate_in_face(face_dof_values, + integration_flags, + sum_into_values); } template +template inline void -FEPointEvaluation::submit_value( - const value_type &value, - const unsigned int point_index) +FEFacePointEvaluation:: + do_integrate_in_face( + ScalarNumber *face_dof_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + const bool sum_into_values) { - AssertIndexRange(point_index, n_q_points); - values[point_index] = value; -} - - + // zero out lanes of incomplete last quadrature point batch + if constexpr (stride == 1) + if (const unsigned int n_filled_lanes = + this->n_q_points_scalar & (n_lanes_internal - 1); + n_filled_lanes > 0) + { + if (integration_flags & EvaluationFlags::values) + for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) + ETT::set_zero_value(this->values.back(), v); + if (integration_flags & EvaluationFlags::gradients) + for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) + ETT::set_zero_gradient(this->gradients.back(), v); + } -template -inline void -FEPointEvaluation::submit_gradient( - const gradient_type &gradient, - const unsigned int point_index) -{ - AssertIndexRange(point_index, n_q_points); - gradients[point_index] = gradient; -} + std::array + solution_values_vectorized_linear = {}; + // loop over quadrature batches qb + const unsigned int n_shapes = is_linear ? 2 : this->poly.size(); + const bool cartesian_cell = + this->cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian; + const bool affine_cell = + this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine; + for (unsigned int qb = 0; qb < this->n_q_batches; ++qb) + { + vectorized_value_type value = {}; + Tensor<1, dim, vectorized_value_type> gradient; -template -inline DerivativeForm<1, dim, spacedim, Number> -FEPointEvaluation::jacobian( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, n_q_points); - Assert(jacobian_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_jacobians")); - return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions:: - GeometryType::affine ? - 0 : - point_index]; -} + if (integration_flags & EvaluationFlags::values) + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + ETT::get_value(value, + v, + do_JxW ? this->values[offset] * this->JxW_ptr[offset] : + this->values[offset]); + if (integration_flags & EvaluationFlags::gradients) + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < this->n_q_points_scalar); + ++v, ++offset) + { + const auto grad_w = + do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] : + this->gradients[offset]; + ETT::get_gradient( + gradient, + v, + cartesian_cell ? + apply_diagonal_transformation(this->inverse_jacobian_ptr[0], + grad_w) : + apply_transformation( + this->inverse_jacobian_ptr[affine_cell ? 0 : offset], + grad_w)); + } + if (integration_flags & EvaluationFlags::gradients) + { + std::array value_face = {}; + Tensor<1, dim - 1, vectorized_value_type> gradient_in_face; -template -inline DerivativeForm<1, spacedim, dim, Number> -FEPointEvaluation::inverse_jacobian( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, n_q_points); - Assert(inverse_jacobian_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_inverse_jacobians")); - return inverse_jacobian_ptr - [cell_type <= - ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ? - 0 : - point_index]; -} + value_face[0] = value; + // fill derivative in physical coordinates into tangential/normal + // derivatives + if (this->current_face_number / 2 == 0) + { + value_face[1] = gradient[0]; + if (dim > 1) + gradient_in_face[0] = gradient[1]; + if (dim > 2) + gradient_in_face[1] = gradient[2]; + } + else if (this->current_face_number / 2 == 1) + { + if (dim > 1) + value_face[1] = gradient[1]; + if (dim == 3) + { + gradient_in_face[0] = gradient[2]; + gradient_in_face[1] = gradient[0]; + } + else if (dim == 2) + gradient_in_face[0] = gradient[0]; + else + Assert(false, ExcInternalError()); + } + else if (this->current_face_number / 2 == 2) + { + if (dim > 2) + { + value_face[1] = gradient[2]; + gradient_in_face[0] = gradient[0]; + gradient_in_face[1] = gradient[1]; + } + else + Assert(false, ExcInternalError()); + } + else + Assert(false, ExcInternalError()); + internal::integrate_tensor_product_value_and_gradient< + is_linear, + dim - 1, + VectorizedArrayType, + vectorized_value_type, + 2>(this->shapes_faces.data() + qb * n_shapes, + n_shapes, + value_face.data(), + gradient_in_face, + is_linear ? solution_values_vectorized_linear.data() : + this->solution_renumbered_vectorized.data(), + this->unit_point_faces_ptr[qb], + qb != 0); + } + else + internal::integrate_tensor_product_value( + this->shapes_faces.data() + qb * n_shapes, + n_shapes, + value, + is_linear ? solution_values_vectorized_linear.data() : + this->solution_renumbered_vectorized.data(), + this->unit_point_faces_ptr[qb], + qb != 0); + } + const unsigned int dofs_per_comp_face = + is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face; -template -inline Number -FEPointEvaluation::JxW( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, n_q_points); - Assert(JxW_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_JxW_values")); - return JxW_ptr[point_index]; + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i) + if (sum_into_values) + face_dof_values[(i + comp * 3 * dofs_per_comp_face) * + stride_face_dof] += + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear.data() + i) : + this->solution_renumbered_vectorized[i]); + else + face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] = + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear.data() + i) : + this->solution_renumbered_vectorized[i]); } template inline Tensor<1, spacedim, Number> -FEPointEvaluation::normal_vector( +FEFacePointEvaluation::normal_vector( const unsigned int point_index) const { - AssertIndexRange(point_index, n_q_points); - Assert(normal_ptr != nullptr, + AssertIndexRange(point_index, this->n_q_points); + Assert(this->normal_ptr != nullptr, internal::FEPointEvaluation:: ExcFEPointEvaluationAccessToUninitializedMappingField( "update_normal_vectors")); - return normal_ptr[point_index]; -} - - - -template -inline Point -FEPointEvaluation::real_point( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, n_q_points); - Assert(real_point_ptr != nullptr, - internal::FEPointEvaluation:: - ExcFEPointEvaluationAccessToUninitializedMappingField( - "update_quadrature_points")); - return real_point_ptr[point_index]; -} - - - -template -inline Point -FEPointEvaluation::unit_point( - const unsigned int point_index) const -{ - AssertIndexRange(point_index, n_q_points); - Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!")); - Point unit_point; - for (unsigned int d = 0; d < dim; ++d) - unit_point[d] = internal::VectorizedArrayTrait::get_from_vectorized( - unit_point_ptr[point_index / stride][d], point_index % stride); - return unit_point; -} - - - -template -inline std_cxx20::ranges::iota_view -FEPointEvaluation:: - quadrature_point_indices() const -{ - return {0U, n_q_points}; + if (this->cell_type <= dealii::internal::MatrixFreeFunctions::affine) + { + Tensor<1, spacedim, Number> normal; + for (unsigned int d = 0; d < dim; ++d) + normal[d] = + internal::VectorizedArrayTrait::get(this->normal_ptr[0][d], + 0); + if (this->is_interior) + return normal; + else + return -normal; + } + else + { + if (this->is_interior) + return this->normal_ptr[point_index]; + else + return -(this->normal_ptr[point_index]); + } } DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 8832ec3809..93fb94d5da 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2084,7 +2084,8 @@ namespace internal typename Number2, typename Number, int n_values = 1, - bool do_renumber = true> + bool do_renumber = true, + int stride = 1> inline #ifndef DEBUG DEAL_II_ALWAYS_INLINE @@ -2108,8 +2109,9 @@ namespace internal // is used to interpolate normal derivatives onto faces. const Number *values_2 = n_values > 1 ? - values + (length > 0 ? Utilities::pow(length, dim) : - Utilities::fixed_power(n_shapes_runtime)) : + values + stride * (length > 0 ? + Utilities::pow(length, dim) : + Utilities::fixed_power(n_shapes_runtime)) : nullptr; using Number3 = typename ProductTypeNoPoint::type; std::array result = {}; @@ -2124,21 +2126,24 @@ namespace internal for (int i0 = 0; i0 < n_shapes; ++i0, ++i) { // gradient - inner_result[0] += shapes[i0][1][0] * values[renumber[i]]; + inner_result[0] += + shapes[i0][1][0] * values[renumber[i] * stride]; // values - inner_result[1] += shapes[i0][0][0] * values[renumber[i]]; + inner_result[1] += + shapes[i0][0][0] * values[renumber[i] * stride]; if (n_values > 1) - inner_result[2] += shapes[i0][0][0] * values_2[renumber[i]]; + inner_result[2] += + shapes[i0][0][0] * values_2[renumber[i] * stride]; } else for (int i0 = 0; i0 < n_shapes; ++i0, ++i) { // gradient - inner_result[0] += shapes[i0][1][0] * values[i]; + inner_result[0] += shapes[i0][1][0] * values[i * stride]; // values - inner_result[1] += shapes[i0][0][0] * values[i]; + inner_result[1] += shapes[i0][0][0] * values[i * stride]; if (n_values > 1) - inner_result[2] += shapes[i0][0][0] * values_2[i]; + inner_result[2] += shapes[i0][0][0] * values_2[i * stride]; } if (dim > 1) @@ -2175,7 +2180,8 @@ namespace internal typename Number, typename Number2, int n_values = 1, - bool do_renumber = true> + bool do_renumber = true, + int stride = 1> inline std::array::type, dim + n_values> evaluate_tensor_product_value_and_gradient_shapes( @@ -2198,7 +2204,7 @@ namespace internal // point, simply set the result vector accordingly. result[0] = values[0]; if (n_values > 1) - result[1] = values[1]; + result[1] = values[1 * stride]; return result; } @@ -2211,28 +2217,58 @@ namespace internal // cases if (n_shapes == 2) inner_result = - do_interpolate_xy( - values, renumber, shapes, n_shapes, i); + do_interpolate_xy(values, renumber, shapes, n_shapes, i); else if (n_shapes == 3) inner_result = - do_interpolate_xy( - values, renumber, shapes, n_shapes, i); + do_interpolate_xy(values, renumber, shapes, n_shapes, i); else if (n_shapes == 4) inner_result = - do_interpolate_xy( - values, renumber, shapes, n_shapes, i); + do_interpolate_xy(values, renumber, shapes, n_shapes, i); else if (n_shapes == 5) inner_result = - do_interpolate_xy( - values, renumber, shapes, n_shapes, i); + do_interpolate_xy(values, renumber, shapes, n_shapes, i); else if (n_shapes == 6) inner_result = - do_interpolate_xy( - values, renumber, shapes, n_shapes, i); + do_interpolate_xy(values, renumber, shapes, n_shapes, i); else inner_result = - do_interpolate_xy( - values, renumber, shapes, n_shapes, i); + do_interpolate_xy(values, renumber, shapes, n_shapes, i); if (dim == 3) { // derivative + interpolation in z direction @@ -2292,9 +2328,6 @@ namespace internal using Number3 = typename ProductTypeNoPoint::type; - static_assert( - n_values == 1 || stride == 1, - "Either n_values or stride has to be one for correct data access!"); // If n_values > 1, we want to interpolate from a second array, // placed in the same array immediately after the main data. This // is used to interpolate normal derivatives onto faces. @@ -2305,7 +2338,7 @@ namespace internal // we only need the value on faces of a 1d element result[0] = values[0]; if (n_values > 1) - result[1] = values[1]; + result[1] = values[1 * stride]; } else if (dim == 1) { @@ -2314,7 +2347,8 @@ namespace internal // values result[1] = Number3(values[0]) + p[0] * result[0]; if (n_values > 1) - result[2] = Number3(values[2]) + p[0] * (values[3] - values[2]); + result[2] = Number3(values[2 * stride]) + + p[0] * (values[3 * stride] - values[2 * stride]); } else if (dim == 2) { @@ -2333,9 +2367,11 @@ namespace internal if (n_values > 1) { const Number3 tmp0_2 = - Number3(values[4]) + p[0] * (values[5] - values[4]); + Number3(values[4 * stride]) + + p[0] * (values[5 * stride] - values[4 * stride]); const Number3 tmp1_2 = - Number3(values[6]) + p[0] * (values[7] - values[6]); + Number3(values[6 * stride]) + + p[0] * (values[7 * stride] - values[6 * stride]); result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2); } } @@ -2445,7 +2481,8 @@ namespace internal int length, typename Number2, typename Number, - bool do_renumber = true> + bool do_renumber = true, + int stride = 1> inline #ifndef DEBUG DEAL_II_ALWAYS_INLINE @@ -2469,10 +2506,10 @@ namespace internal // renumbering or not if (do_renumber && !renumber.empty()) for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - value += shapes[i0][0][0] * values[renumber[i]]; + value += shapes[i0][0][0] * values[renumber[i] * stride]; else for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - value += shapes[i0][0][0] * values[i]; + value += shapes[i0][0][0] * values[i * stride]; if (dim > 1) result += value * shapes[i1][0][1]; @@ -2484,7 +2521,11 @@ namespace internal - template + template inline typename ProductTypeNoPoint::type evaluate_tensor_product_value_shapes( const dealii::ndarray *shapes, @@ -2511,29 +2552,53 @@ namespace internal // Generate separate code with known loop bounds for the most common // cases if (n_shapes == 2) - inner_result = - do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); else if (n_shapes == 3) - inner_result = - do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); else if (n_shapes == 4) - inner_result = - do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); else if (n_shapes == 5) - inner_result = - do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); else if (n_shapes == 6) - inner_result = - do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); else - inner_result = - do_interpolate_xy_value( - values, renumber, shapes, n_shapes, i); + inner_result = do_interpolate_xy_value( + values, renumber, shapes, n_shapes, i); if (dim == 3) { // Interpolation + derivative in z direction diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index c7d091a01e..8663c0f4f6 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -378,6 +378,16 @@ namespace NonMatching const std::vector>> &quadrature_vector, const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int); + /** + * Compute the mapping information incoming vector of faces and + * corresponding vector of quadratures. + */ + template + void + reinit_faces(const std::vector> + &face_iterator_range_interior, + const std::vector> &quadrature_vector); + /** * Return if this MappingInfo object is reinitialized for faces (by * reinit_faces()) or not. @@ -385,6 +395,12 @@ namespace NonMatching bool is_face_state() const; + /** + * Returns the face number of the interior/exterior face. + */ + unsigned int + get_face_number(const unsigned int offset, const bool is_interior) const; + /** * Getter function for unit points. The offset can be obtained with * compute_unit_point_index_offset(). @@ -404,14 +420,16 @@ namespace NonMatching * compute_data_index_offset(). */ const DerivativeForm<1, dim, spacedim, Number> * - get_jacobian(const unsigned int offset) const; + get_jacobian(const unsigned int offset, + const bool is_interior = true) const; /** * Getter function for inverse Jacobians. The offset can be obtained with * compute_data_index_offset(). */ const DerivativeForm<1, spacedim, dim, Number> * - get_inverse_jacobian(const unsigned int offset) const; + get_inverse_jacobian(const unsigned int offset, + const bool is_interior = true) const; /** * Getter function for normal vectors. The offset can be obtained with @@ -545,7 +563,8 @@ namespace NonMatching * Resize the mapping data fields. */ void - resize_data_fields(const unsigned int n_data_point_batches); + resize_data_fields(const unsigned int n_data_point_batches, + const bool is_face_centric = false); /** * Store the unit points. @@ -575,7 +594,8 @@ namespace NonMatching const MappingData &mapping_data, const std::vector &weights, const unsigned int compressed_unit_point_index_offset, - const bool affine_cell); + const bool affine_cell, + const bool is_interior = true); /** * Compute the compressed cell index. @@ -605,7 +625,8 @@ namespace NonMatching invalid, single_cell, cell_vector, - faces_on_cells_in_vector + faces_on_cells_in_vector, + face_vector }; /** @@ -698,7 +719,8 @@ namespace NonMatching * * Indexed by @p compressed_data_index_offsets. */ - AlignedVector> jacobians; + std::array>, 2> + jacobians; /** * The storage of covariant transformation on quadrature points, i.e., @@ -707,7 +729,8 @@ namespace NonMatching * * Indexed by @p compressed_data_index_offsets. */ - AlignedVector> inverse_jacobians; + std::array>, 2> + inverse_jacobians; /** * The mapped real points. @@ -756,8 +779,20 @@ namespace NonMatching * vector is only filled if AdditionalData::store_cells is enabled. */ std::vector> cell_level_and_indices; + + std::vector> face_number; }; + template + inline unsigned int + MappingInfo::get_face_number( + const unsigned int offset, + const bool is_interior) const + { + const auto &face_pair = face_number[offset]; + return is_interior ? face_pair.first : face_pair.second; + } + // ----------------------- template functions ---------------------- @@ -1116,13 +1151,13 @@ namespace NonMatching if (update_flags_mapping & UpdateFlags::update_jacobians) { - jacobians.resize(size_compressed_data); - jacobians.shrink_to_fit(); + jacobians[0].resize(size_compressed_data); + jacobians[0].shrink_to_fit(); } if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) { - inverse_jacobians.resize(size_compressed_data); - inverse_jacobians.shrink_to_fit(); + inverse_jacobians[0].resize(size_compressed_data); + inverse_jacobians[0].shrink_to_fit(); } state = State::cell_vector; @@ -1449,13 +1484,13 @@ namespace NonMatching if (update_flags_mapping & UpdateFlags::update_jacobians) { - jacobians.resize(size_compressed_data); - jacobians.shrink_to_fit(); + jacobians[0].resize(size_compressed_data); + jacobians[0].shrink_to_fit(); } if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) { - inverse_jacobians.resize(size_compressed_data); - inverse_jacobians.shrink_to_fit(); + inverse_jacobians[0].resize(size_compressed_data); + inverse_jacobians[0].shrink_to_fit(); } state = State::faces_on_cells_in_vector; @@ -1464,6 +1499,289 @@ namespace NonMatching + template + template + void + MappingInfo::reinit_faces( + const std::vector> + &face_iterator_range_interior, + const std::vector> &quadrature_vector) + { + clear(); + + do_cell_index_compression = false; + + Assert(additional_data.store_cells == false, ExcNotImplemented()); + + + const unsigned int n_faces = quadrature_vector.size(); + AssertDimension(n_faces, + std::distance(face_iterator_range_interior.begin(), + face_iterator_range_interior.end())); + + n_q_points_unvectorized.reserve(n_faces); + + cell_type.reserve(n_faces); + face_number.reserve(n_faces); + + // fill unit points index offset vector + unit_points_index.reserve(n_faces + 1); + unit_points_index.push_back(0); + data_index_offsets.reserve(n_faces + 1); + data_index_offsets.push_back(0); + for (const auto &quadrature : quadrature_vector) + { + const unsigned int n_points = quadrature.size(); + n_q_points_unvectorized.push_back(n_points); + + const unsigned int n_q_points = + compute_n_q_points(n_points); + unit_points_index.push_back(unit_points_index.back() + n_q_points); + + const unsigned int n_q_points_data = + compute_n_q_points(n_points); + data_index_offsets.push_back(data_index_offsets.back() + + n_q_points_data); + } + + const unsigned int n_unit_points = unit_points_index.back(); + const unsigned int n_data_points = data_index_offsets.back(); + + // resize data vectors + resize_unit_points(n_unit_points); + resize_unit_points_faces(n_unit_points); + resize_data_fields(n_data_points, true); + + std::array mapping_data; + std::array mapping_data_previous_cell; + std::array mapping_data_first; + bool first_set = false; + unsigned int size_compressed_data = 0; + unsigned int face_index = 0; + QProjector q_projector; + for (const auto &cell_and_f : face_iterator_range_interior) + { + const auto &quadrature_on_face = quadrature_vector[face_index]; + const bool empty = quadrature_on_face.empty(); + + // get interior cell and face number + const auto &cell_m = cell_and_f.first; + const auto f_m = cell_and_f.second; + + // get exterior cell and face number + const auto &cell_p = + cell_m->at_boundary(f_m) ? cell_m : cell_m->neighbor(f_m); + const auto f_p = + cell_m->at_boundary(f_m) ? f_m : cell_m->neighbor_of_neighbor(f_m); + + face_number.emplace_back(f_m, f_p); + + Assert( + cell_m->combined_face_orientation(f_m) == 1 && + cell_p->combined_face_orientation(f_p) == 1, + ExcMessage( + "Non standard face orientation is currently not implemented.")); + + const auto quadrature_on_cell_m = + q_projector.project_to_face(cell_m->reference_cell(), + quadrature_on_face, + f_m); + + // store unit points + const unsigned int n_q_points = compute_n_q_points( + n_q_points_unvectorized[face_index]); + store_unit_points(unit_points_index[face_index], + n_q_points, + n_q_points_unvectorized[face_index], + quadrature_on_cell_m.get_points()); + + store_unit_points_faces(unit_points_index[face_index], + n_q_points, + n_q_points_unvectorized[face_index], + quadrature_on_face.get_points()); + + // compute mapping for interior face + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_face_quadrature(mapping, + update_flags_mapping, + cell_m, + f_m, + quadrature_on_face, + internal_mapping_data, + mapping_data[0]); + + // compute mapping for exterior face + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_face_quadrature(mapping, + update_flags_mapping, + cell_p, + f_p, + quadrature_on_face, + internal_mapping_data, + mapping_data[1]); + + // check for cartesian/affine cell + if (!empty && + update_flags_mapping & UpdateFlags::update_inverse_jacobians) + { + // select more general type of interior and exterior cell + cell_type.push_back(std::max( + internal::compute_geometry_type( + cell_m->diameter(), mapping_data[0].inverse_jacobians), + internal::compute_geometry_type( + cell_m->diameter(), mapping_data[1].inverse_jacobians))); + + // cache mapping data of first cell pair with non-empty quadrature + // on the face + if (!first_set) + { + mapping_data_first = mapping_data; + first_set = true; + } + } + else + cell_type.push_back( + dealii::internal::MatrixFreeFunctions::GeometryType::general); + + if (face_index > 0) + { + // check if current and previous cell pairs are affine + const bool affine_cells = + cell_type[face_index] <= + dealii::internal::MatrixFreeFunctions::affine && + cell_type[face_index - 1] <= + dealii::internal::MatrixFreeFunctions::affine; + + // create a comparator to compare inverse Jacobian of current + // and previous cell pair + FloatingPointComparator comparator( + 1e4 / cell_m->diameter() * + std::numeric_limits::epsilon() * 1024.); + + // we can only compare if current and previous cell have at + // least one quadrature point and both cells are at least affine + const auto comparison_result_m = + (!affine_cells || mapping_data[0].inverse_jacobians.empty() || + mapping_data_previous_cell[0].inverse_jacobians.empty()) ? + FloatingPointComparator::ComparisonResult::less : + comparator.compare( + mapping_data[0].inverse_jacobians[0], + mapping_data_previous_cell[0].inverse_jacobians[0]); + + const auto comparison_result_p = + (!affine_cells || mapping_data[1].inverse_jacobians.empty() || + mapping_data_previous_cell[1].inverse_jacobians.empty()) ? + FloatingPointComparator::ComparisonResult::less : + comparator.compare( + mapping_data[1].inverse_jacobians[0], + mapping_data_previous_cell[1].inverse_jacobians[0]); + + // we can compress the Jacobians and inverse Jacobians if + // inverse Jacobians are equal and cells are affine + if (affine_cells && + comparison_result_m == + FloatingPointComparator::ComparisonResult::equal && + comparison_result_p == + FloatingPointComparator::ComparisonResult::equal) + { + compressed_data_index_offsets.push_back( + compressed_data_index_offsets.back()); + } + else if (first_set && + (cell_type[face_index] <= + dealii::internal::MatrixFreeFunctions::affine) && + (comparator.compare( + mapping_data[0].inverse_jacobians[0], + mapping_data_first[0].inverse_jacobians[0]) == + FloatingPointComparator< + double>::ComparisonResult::equal) && + (comparator.compare( + mapping_data[1].inverse_jacobians[0], + mapping_data_first[1].inverse_jacobians[0]) == + FloatingPointComparator::ComparisonResult::equal)) + { + compressed_data_index_offsets.push_back(0); + } + else + { + const unsigned int n_compressed_data_last_cell = + cell_type[face_index - 1] <= + dealii::internal::MatrixFreeFunctions::affine ? + 1 : + compute_n_q_points( + n_q_points_unvectorized[face_index - 1]); + + compressed_data_index_offsets.push_back( + compressed_data_index_offsets.back() + + n_compressed_data_last_cell); + } + } + else + compressed_data_index_offsets.push_back(0); + + // cache mapping_data from previous cell pair + mapping_data_previous_cell = mapping_data; + + const unsigned int n_q_points_data = + compute_n_q_points(n_q_points_unvectorized[face_index]); + + // store mapping data of interior face + store_mapping_data(data_index_offsets[face_index], + n_q_points_data, + n_q_points_unvectorized[face_index], + mapping_data[0], + quadrature_on_face.get_weights(), + data_index_offsets[face_index], + cell_type[face_index] <= + dealii::internal::MatrixFreeFunctions::affine, + true); + + // store only necessary mapping data for exterior face (Jacobians and + // inverse Jacobians) + store_mapping_data(data_index_offsets[face_index], + n_q_points_data, + n_q_points_unvectorized[face_index], + mapping_data[1], + quadrature_on_face.get_weights(), + data_index_offsets[face_index], + cell_type[face_index] <= + dealii::internal::MatrixFreeFunctions::affine, + false); + + // update size of compressed data depending on cell type and handle + // empty quadratures + if (cell_type[face_index] <= + dealii::internal::MatrixFreeFunctions::affine) + size_compressed_data = compressed_data_index_offsets.back() + 1; + else + size_compressed_data = + std::max(size_compressed_data, + compressed_data_index_offsets.back() + n_q_points_data); + + ++face_index; + } + + if (update_flags_mapping & UpdateFlags::update_jacobians) + { + jacobians[0].resize(size_compressed_data); + jacobians[0].shrink_to_fit(); + jacobians[1].resize(size_compressed_data); + jacobians[1].shrink_to_fit(); + } + if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) + { + inverse_jacobians[0].resize(size_compressed_data); + inverse_jacobians[0].shrink_to_fit(); + inverse_jacobians[1].resize(size_compressed_data); + inverse_jacobians[1].shrink_to_fit(); + } + + state = State::face_vector; + is_reinitialized(); + } + + + template bool MappingInfo::is_face_state() const @@ -1550,10 +1868,14 @@ namespace NonMatching Assert(cell_index != numbers::invalid_unsigned_int, ExcMessage( "cell_index has to be set if face_number is specified!")); - Assert(state == State::faces_on_cells_in_vector, + Assert(state == State::faces_on_cells_in_vector || + state == State::face_vector, ExcMessage("This mapping info is not reinitialized for faces" " on cells in a vector!")); - return cell_index_offset[compressed_cell_index] + face_number; + if (state == State::faces_on_cells_in_vector) + return cell_index_offset[compressed_cell_index] + face_number; + else if (state == State::face_vector) + return cell_index; } } @@ -1636,7 +1958,8 @@ namespace NonMatching const MappingInfo::MappingData &mapping_data, const std::vector &weights, const unsigned int compressed_unit_point_index_offset, - const bool affine_cell) + const bool affine_cell, + const bool is_interior) { const unsigned int n_lanes = dealii::internal::VectorizedArrayTrait::width(); @@ -1656,40 +1979,47 @@ namespace NonMatching for (unsigned int d = 0; d < dim; ++d) for (unsigned int s = 0; s < spacedim; ++s) dealii::internal::VectorizedArrayTrait::get( - jacobians[compressed_offset][d][s], v) = - mapping_data.jacobians[q * n_lanes + v][d][s]; + jacobians[is_interior ? 0 : 1][compressed_offset][d][s], + v) = mapping_data.jacobians[q * n_lanes + v][d][s]; if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) for (unsigned int d = 0; d < dim; ++d) for (unsigned int s = 0; s < spacedim; ++s) dealii::internal::VectorizedArrayTrait::get( - inverse_jacobians[compressed_offset][s][d], v) = + inverse_jacobians[is_interior ? 0 : 1] + [compressed_offset][s][d], + v) = mapping_data.inverse_jacobians[q * n_lanes + v][s][d]; } - if (update_flags_mapping & UpdateFlags::update_JxW_values) + + if (is_interior) { - if (additional_data.use_global_weights) + if (update_flags_mapping & UpdateFlags::update_JxW_values) { - dealii::internal::VectorizedArrayTrait::get( - JxW_values[offset], v) = weights[q * n_lanes + v]; + if (additional_data.use_global_weights) + { + dealii::internal::VectorizedArrayTrait::get( + JxW_values[offset], v) = weights[q * n_lanes + v]; + } + else + { + dealii::internal::VectorizedArrayTrait::get( + JxW_values[offset], v) = + mapping_data.JxW_values[q * n_lanes + v]; + } } - else - { + if (update_flags_mapping & UpdateFlags::update_normal_vectors) + for (unsigned int s = 0; s < spacedim; ++s) dealii::internal::VectorizedArrayTrait::get( - JxW_values[offset], v) = - mapping_data.JxW_values[q * n_lanes + v]; - } + normal_vectors[offset][s], v) = + mapping_data.normal_vectors[q * n_lanes + v][s]; + if (update_flags_mapping & + UpdateFlags::update_quadrature_points) + for (unsigned int s = 0; s < spacedim; ++s) + dealii::internal::VectorizedArrayTrait::get( + real_points[offset][s], v) = + mapping_data.quadrature_points[q * n_lanes + v][s]; } - if (update_flags_mapping & UpdateFlags::update_normal_vectors) - for (unsigned int s = 0; s < spacedim; ++s) - dealii::internal::VectorizedArrayTrait::get( - normal_vectors[offset][s], v) = - mapping_data.normal_vectors[q * n_lanes + v][s]; - if (update_flags_mapping & UpdateFlags::update_quadrature_points) - for (unsigned int s = 0; s < spacedim; ++s) - dealii::internal::VectorizedArrayTrait::get( - real_points[offset][s], v) = - mapping_data.quadrature_points[q * n_lanes + v][s]; } } } @@ -1719,12 +2049,21 @@ namespace NonMatching template void MappingInfo::resize_data_fields( - const unsigned int n_data_point_batches) + const unsigned int n_data_point_batches, + const bool is_face_centric) { if (update_flags_mapping & UpdateFlags::update_jacobians) - jacobians.resize(n_data_point_batches); + { + jacobians[0].resize(n_data_point_batches); + if (is_face_centric) + jacobians[1].resize(n_data_point_batches); + } if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) - inverse_jacobians.resize(n_data_point_batches); + { + inverse_jacobians[0].resize(n_data_point_batches); + if (is_face_centric) + inverse_jacobians[1].resize(n_data_point_batches); + } if (update_flags_mapping & UpdateFlags::update_JxW_values) JxW_values.resize(n_data_point_batches); if (update_flags_mapping & UpdateFlags::update_normal_vectors) @@ -1800,10 +2139,10 @@ namespace NonMatching template inline const DerivativeForm<1, dim, spacedim, Number> * - MappingInfo::get_jacobian( - const unsigned int offset) const + MappingInfo::get_jacobian(const unsigned int offset, + const bool is_interior) const { - return jacobians.data() + offset; + return jacobians[is_interior ? 0 : 1].data() + offset; } @@ -1811,12 +2150,14 @@ namespace NonMatching template inline const DerivativeForm<1, spacedim, dim, Number> * MappingInfo::get_inverse_jacobian( - const unsigned int offset) const + const unsigned int offset, + const bool is_interior) const { - return inverse_jacobians.data() + offset; + return inverse_jacobians[is_interior ? 0 : 1].data() + offset; } + template inline const Tensor<1, spacedim, Number> * MappingInfo::get_normal_vector( diff --git a/tests/matrix_free/point_evaluation_23.cc b/tests/matrix_free/point_evaluation_23.cc index 3aa622b03f..1bcca19668 100644 --- a/tests/matrix_free/point_evaluation_23.cc +++ b/tests/matrix_free/point_evaluation_23.cc @@ -105,7 +105,7 @@ test(const unsigned int degree) mapping_info.reinit_faces(tria.active_cell_iterators(), quad_vec_faces); - FEPointEvaluation<1, dim> evaluator(mapping_info, fe); + FEFacePointEvaluation<1, dim> evaluator(mapping_info, fe); VectorTools::interpolate(mapping, dof_handler, MyFunction(), vector); diff --git a/tests/matrix_free/point_evaluation_24.cc b/tests/matrix_free/point_evaluation_24.cc index c21330e1b4..e432fa42f4 100644 --- a/tests/matrix_free/point_evaluation_24.cc +++ b/tests/matrix_free/point_evaluation_24.cc @@ -105,7 +105,7 @@ test(const unsigned int degree) mapping_info.reinit_faces(tria.active_cell_iterators(), quad_vec_faces); - FEPointEvaluation evaluator(mapping_info, fe); + FEFacePointEvaluation evaluator(mapping_info, fe); VectorTools::interpolate(mapping, dof_handler, MyFunction(), vector); diff --git a/tests/matrix_free/point_evaluation_27.cc b/tests/matrix_free/point_evaluation_27.cc index 1d37406611..836c7fad1e 100644 --- a/tests/matrix_free/point_evaluation_27.cc +++ b/tests/matrix_free/point_evaluation_27.cc @@ -80,8 +80,8 @@ test(const FiniteElement &fe) std::vector>(4, QGauss(7))); dealii::NonMatching::MappingInfo mapping_info(mapping, update_values); mapping_info.reinit_faces(tria.active_cell_iterators(), quad_vec); - FEPointEvaluation - fe_point_eval(mapping_info, fe, first_selected_component); + FEFacePointEvaluation + fe_point_eval(mapping_info, fe, true, first_selected_component); std::vector buffer(fe.dofs_per_cell); diff --git a/tests/non_matching/mapping_info.cc b/tests/non_matching/mapping_info.cc index 912fd13602..55375f4027 100644 --- a/tests/non_matching/mapping_info.cc +++ b/tests/non_matching/mapping_info.cc @@ -183,10 +183,10 @@ test(const bool filtered_compression) FEPointEvaluation<1, dim, dim, double> fe_point_cell(mapping_info_cell, fe_q); FEPointEvaluation<1, dim, dim, double> fe_point_surface(mapping_info_surface, fe_q); - FEPointEvaluation<1, dim, dim, double> fe_point_faces_m(mapping_info_faces, - fe_q); - FEPointEvaluation<1, dim, dim, double> fe_point_faces_p(mapping_info_faces, - fe_q); + FEFacePointEvaluation<1, dim, dim, double> fe_point_faces_m( + mapping_info_faces, fe_q); + FEFacePointEvaluation<1, dim, dim, double> fe_point_faces_p( + mapping_info_faces, fe_q); std::vector solution_values_in(fe_q.dofs_per_cell); std::vector solution_values_neighbor_in(fe_q.dofs_per_cell); diff --git a/tests/non_matching/mapping_info_03.cc b/tests/non_matching/mapping_info_03.cc index ecce078503..0531a6b276 100644 --- a/tests/non_matching/mapping_info_03.cc +++ b/tests/non_matching/mapping_info_03.cc @@ -78,7 +78,7 @@ test() // 2) reinit mapping info mapping_info.reinit_faces(tria.active_cell_iterators(), quad_vec); - FEPointEvaluation fe_point_eval( + FEFacePointEvaluation fe_point_eval( mapping_info, fe); // 3) print JxW diff --git a/tests/non_matching/mapping_info_04.cc b/tests/non_matching/mapping_info_04.cc new file mode 100644 index 0000000000..bee712570e --- /dev/null +++ b/tests/non_matching/mapping_info_04.cc @@ -0,0 +1,431 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the NonMatching::MappingInfo class together with FEPointEvaluation and + * compare to FEEvaluation + */ + +#include + +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +do_flux_term(Integrator &evaluator_m, + Integrator &evaluator_p, + const Number2 &tau, + const unsigned int q) +{ + const auto gradient_m = evaluator_m.get_gradient(q); + const auto gradient_p = evaluator_p.get_gradient(q); + + const auto value_m = evaluator_m.get_value(q); + const auto value_p = evaluator_p.get_value(q); + + const auto normal = evaluator_m.normal_vector(q); + + const auto jump_value = (value_m - value_p) * normal; + + const auto central_flux_gradient = 0.5 * (gradient_m + gradient_p); + + const auto value_terms = normal * (central_flux_gradient - tau * jump_value); + + evaluator_m.submit_value(-value_terms, q); + evaluator_p.submit_value(value_terms, q); + + const auto gradient_terms = -0.5 * jump_value; + + evaluator_m.submit_gradient(gradient_terms, q); + evaluator_p.submit_gradient(gradient_terms, q); +} + +template +void +test_dg_fcl(const unsigned int degree, const bool curved_mesh) +{ + constexpr unsigned int n_lanes = VectorizedArray::size(); + + const unsigned int n_q_points = degree + 1; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + if (curved_mesh && dim > 1) + GridGenerator::hyper_shell(tria, Point(), 0.5, 1, 6); + else + GridGenerator::subdivided_hyper_cube(tria, 2, 0, 1); + + tria.refine_global(1); + + FE_DGQ fe(degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + MappingQGeneric mapping(degree); + + AffineConstraints constraints; + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, Functions::ZeroFunction(), constraints); + constraints.close(); + + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags_inner_faces = + update_values | update_gradients; + additional_data.mapping_update_flags_boundary_faces = + update_values | update_gradients; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, QGauss<1>(n_q_points), additional_data); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "Working with " << fe.get_name() << " and " + << dof_handler.n_dofs() << " dofs" << std::endl; + + LinearAlgebra::distributed::Vector src, dst, dst2; + matrix_free.initialize_dof_vector(src); + for (auto &v : src) + v = static_cast(rand()) / RAND_MAX; + + matrix_free.initialize_dof_vector(dst); + matrix_free.initialize_dof_vector(dst2); + + matrix_free.template loop, + LinearAlgebra::distributed::Vector>( + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto &range) { + FEEvaluation fe_eval(matrix_free); + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.gather_evaluate(src, EvaluationFlags::gradients); + for (const unsigned int q : fe_eval.quadrature_point_indices()) + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + fe_eval.integrate_scatter(EvaluationFlags::gradients, dst); + } + }, + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto &range) { + FEFaceEvaluation fe_eval_m(matrix_free, true); + FEFaceEvaluation fe_eval_p(matrix_free, false); + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval_m.reinit(face); + fe_eval_p.reinit(face); + fe_eval_m.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_p.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + for (unsigned int q = 0; q < fe_eval_m.n_q_points; ++q) + do_flux_term(fe_eval_m, fe_eval_p, 1.0, q); + fe_eval_m.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + fe_eval_p.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }, + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto &range) { + FEFaceEvaluation fe_eval_m(matrix_free, true); + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval_m.reinit(face); + fe_eval_m.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + for (unsigned int q = 0; q < fe_eval_m.n_q_points; ++q) + { + const auto value = fe_eval_m.get_value(q); + const auto gradient = fe_eval_m.get_gradient(q); + + fe_eval_m.submit_value(gradient * fe_eval_m.normal_vector(q) + + value, + q); + fe_eval_m.submit_gradient(value * fe_eval_m.normal_vector(q), q); + } + fe_eval_m.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }, + dst, + src, + true); + + QGauss quad_cell(n_q_points); + std::vector> quad_vec_cells; + quad_vec_cells.reserve( + (matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) * + n_lanes); + + + std::vector::cell_iterator> vector_accessors; + vector_accessors.reserve( + (matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) * + n_lanes); + for (unsigned int cell_batch = 0; + cell_batch < + matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches(); + ++cell_batch) + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (v < matrix_free.n_active_entries_per_cell_batch(cell_batch)) + vector_accessors.push_back( + matrix_free.get_cell_iterator(cell_batch, v)); + else + vector_accessors.push_back( + matrix_free.get_cell_iterator(cell_batch, 0)); + + quad_vec_cells.push_back(quad_cell); + } + + QGauss quad_face(n_q_points); + std::vector> quad_vec_faces; + quad_vec_faces.reserve((matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()) * + n_lanes); + std::vector::cell_iterator, unsigned int>> + vector_face_accessors_m; + vector_face_accessors_m.reserve((matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()) * + n_lanes); + // fill container for inner face batches + unsigned int face_batch = 0; + for (; face_batch < matrix_free.n_inner_face_batches(); ++face_batch) + { + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (v < matrix_free.n_active_entries_per_face_batch(face_batch)) + vector_face_accessors_m.push_back( + matrix_free.get_face_iterator(face_batch, v)); + else + vector_face_accessors_m.push_back( + matrix_free.get_face_iterator(face_batch, 0)); + + quad_vec_faces.push_back(quad_face); + } + } + // and boundary face batches + for (; face_batch < (matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + ++face_batch) + { + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (v < matrix_free.n_active_entries_per_face_batch(face_batch)) + vector_face_accessors_m.push_back( + matrix_free.get_face_iterator(face_batch, v)); + else + vector_face_accessors_m.push_back( + matrix_free.get_face_iterator(face_batch, 0)); + + quad_vec_faces.push_back(quad_face); + } + } + + NonMatching::MappingInfo mapping_info_cells(mapping, + update_gradients | + update_JxW_values); + NonMatching::MappingInfo mapping_info_faces(mapping, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors); + + mapping_info_cells.reinit_cells(vector_accessors, quad_vec_cells); + mapping_info_faces.reinit_faces(vector_face_accessors_m, quad_vec_faces); + + FEPointEvaluation<1, dim, dim, double> fe_peval(mapping_info_cells, fe); + FEFacePointEvaluation<1, dim, dim, double> fe_peval_m(mapping_info_faces, + fe, + true); + FEFacePointEvaluation<1, dim, dim, double> fe_peval_p(mapping_info_faces, + fe, + false); + + matrix_free.template loop, + LinearAlgebra::distributed::Vector>( + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto &range) { + FEEvaluation fe_eval(matrix_free); + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + for (unsigned int v = 0; v < n_lanes; ++v) + { + fe_peval.reinit(cell * n_lanes + v); + fe_peval.evaluate(StridedArrayView( + &fe_eval.begin_dof_values()[0][v], + fe.dofs_per_cell), + EvaluationFlags::gradients); + for (const unsigned int q : fe_peval.quadrature_point_indices()) + fe_peval.submit_gradient(fe_peval.get_gradient(q), q); + fe_peval.integrate(StridedArrayView( + &fe_eval.begin_dof_values()[0][v], + fe.dofs_per_cell), + EvaluationFlags::gradients); + } + fe_eval.distribute_local_to_global(dst); + } + }, + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto &range) { + FEFaceEvaluation fe_eval_m(matrix_free, true); + FEFaceEvaluation fe_eval_p(matrix_free, false); + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval_m.reinit(face); + fe_eval_p.reinit(face); + + fe_eval_m.read_dof_values(src); + fe_eval_p.read_dof_values(src); + + fe_eval_m.project_to_face(EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_p.project_to_face(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int v = 0; v < n_lanes; ++v) + { + fe_peval_m.reinit(face * n_lanes + v); + fe_peval_p.reinit(face * n_lanes + v); + fe_peval_m.evaluate_in_face( + &fe_eval_m.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + fe_peval_p.evaluate_in_face( + &fe_eval_p.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + for (const unsigned int q : fe_peval_m.quadrature_point_indices()) + do_flux_term(fe_peval_m, fe_peval_p, 1.0, q); + fe_peval_m.integrate_in_face( + &fe_eval_m.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + fe_peval_p.integrate_in_face( + &fe_eval_p.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + } + + fe_eval_m.collect_from_face(EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_p.collect_from_face(EvaluationFlags::values | + EvaluationFlags::gradients); + + fe_eval_m.distribute_local_to_global(dst); + fe_eval_p.distribute_local_to_global(dst); + } + }, + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto &range) { + FEFaceEvaluation fe_eval_m(matrix_free, true); + for (unsigned int face = range.first; face < range.second; ++face) + { + fe_eval_m.reinit(face); + + fe_eval_m.read_dof_values(src); + + fe_eval_m.project_to_face(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int v = 0; v < n_lanes; ++v) + { + fe_peval_m.reinit(face * n_lanes + v); + fe_peval_m.evaluate_in_face( + &fe_eval_m.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + for (const unsigned int q : fe_peval_m.quadrature_point_indices()) + { + const auto value = fe_peval_m.get_value(q); + const auto gradient = fe_peval_m.get_gradient(q); + + fe_peval_m.submit_value( + gradient * fe_peval_m.normal_vector(q) + value, q); + fe_peval_m.submit_gradient(value * + fe_peval_m.normal_vector(q), + q); + } + fe_peval_m.integrate_in_face( + &fe_eval_m.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + } + + fe_eval_m.collect_from_face(EvaluationFlags::values | + EvaluationFlags::gradients); + + fe_eval_m.distribute_local_to_global(dst); + } + }, + dst2, + src, + true); + + + dst2 -= dst; + const double error = dst2.l2_norm() / dst.l2_norm(); + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "FEPointEvaluation verification: " << error << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + initlog(); + + test_dg_fcl<2>(1, false); + deallog << std::endl; + test_dg_fcl<2>(2, false); + deallog << std::endl; + test_dg_fcl<3>(1, false); + deallog << std::endl; + test_dg_fcl<3>(2, false); + deallog << std::endl; + test_dg_fcl<2>(1, true); + deallog << std::endl; + test_dg_fcl<2>(2, true); + // TODO: fix face orientation + // deallog << std::endl; + // test_dg_fcl<3>(1, true); + // deallog << std::endl; + // test_dg_fcl<3>(2, true); +} diff --git a/tests/non_matching/mapping_info_04.output b/tests/non_matching/mapping_info_04.output new file mode 100644 index 0000000000..a3b3a1fc17 --- /dev/null +++ b/tests/non_matching/mapping_info_04.output @@ -0,0 +1,18 @@ + +DEAL::Working with FE_DGQ<2>(1) and 64 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<2>(2) and 144 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<3>(1) and 512 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<3>(2) and 1728 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<2>(1) and 96 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<2>(2) and 216 dofs +DEAL::FEPointEvaluation verification: 0.0 diff --git a/tests/non_matching/mapping_info_05.cc b/tests/non_matching/mapping_info_05.cc new file mode 100644 index 0000000000..f7dc976cda --- /dev/null +++ b/tests/non_matching/mapping_info_05.cc @@ -0,0 +1,390 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +/* + * Test the NonMatching::MappingInfo class together with FEPointEvaluation and + * compare to FEEvaluation + */ + +#include + +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +do_flux_term_ecl(Integrator &evaluator_m, + Integrator &evaluator_p, + const Number2 &tau, + const unsigned int q) +{ + const auto gradient_m = evaluator_m.get_gradient(q); + const auto gradient_p = evaluator_p.get_gradient(q); + + const auto value_m = evaluator_m.get_value(q); + const auto value_p = evaluator_p.get_value(q); + + const auto normal = evaluator_m.normal_vector(q); + + const auto jump_value = (value_m - value_p) * normal; + + const auto central_flux_gradient = 0.5 * (gradient_m + gradient_p); + + const auto value_terms = normal * (central_flux_gradient - tau * jump_value); + + evaluator_m.submit_value(-value_terms, q); + + evaluator_m.submit_gradient(-0.5 * jump_value, q); +} + +template +void +test_dg_ecl(const unsigned int degree, const bool curved_mesh) +{ + constexpr unsigned int n_lanes = VectorizedArray::size(); + + const unsigned int n_q_points = degree + 1; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + if (curved_mesh && dim > 1) + GridGenerator::hyper_shell(tria, Point(), 0.5, 1, 6); + else + GridGenerator::subdivided_hyper_cube(tria, 2, 0, 1); + + tria.refine_global(1); + + FE_DGQ fe(degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + MappingQGeneric mapping(degree); + + AffineConstraints constraints; + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, Functions::ZeroFunction(), constraints); + constraints.close(); + + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags_inner_faces = + update_values | update_gradients; + additional_data.mapping_update_flags_boundary_faces = + update_values | update_gradients; + additional_data.hold_all_faces_to_owned_cells = true; + additional_data.mapping_update_flags_faces_by_cells = + additional_data.mapping_update_flags_inner_faces | + additional_data.mapping_update_flags_boundary_faces; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, QGauss<1>(n_q_points), additional_data); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "Working with " << fe.get_name() << " and " + << dof_handler.n_dofs() << " dofs" << std::endl; + + LinearAlgebra::distributed::Vector src, dst, dst2; + matrix_free.initialize_dof_vector(src); + for (auto &v : src) + v = static_cast(rand()) / RAND_MAX; + + matrix_free.initialize_dof_vector(dst); + matrix_free.initialize_dof_vector(dst2); + + matrix_free + .template loop_cell_centric, + LinearAlgebra::distributed::Vector>( + [&](const auto &matrix_free, + auto &dst, + const auto &src, + const auto &range) { + FEEvaluation fe_eval(matrix_free); + FEFaceEvaluation fe_eval_m(matrix_free, true); + FEFaceEvaluation fe_eval_p(matrix_free, false); + AlignedVector> vec_solution_values_in_m( + fe_eval.dofs_per_cell); + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + for (unsigned int i = 0; i < fe_eval.dofs_per_cell; ++i) + vec_solution_values_in_m[i] = fe_eval.begin_dof_values()[i]; + fe_eval.evaluate(EvaluationFlags::gradients); + + for (const unsigned int q : fe_eval.quadrature_point_indices()) + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + + fe_eval.integrate(EvaluationFlags::gradients); + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + // ask for boundary ids of face + const auto boundary_ids = + matrix_free.get_faces_by_cells_boundary_id(cell, f); + + // only internal faces have a neighbor, setup a mask + std::bitset mask; + VectorizedArray fluxes = 0.; + for (unsigned int v = 0; v < n_lanes; ++v) + { + mask[v] = + boundary_ids[v] == numbers::internal_face_boundary_id; + fluxes[v] = mask[v] == true ? 1. : 0.; + } + + fe_eval_m.reinit(cell, f); + fe_eval_p.reinit(cell, f); + + fe_eval_p.read_dof_values(src, 0, mask); + + fe_eval_m.evaluate(vec_solution_values_in_m.data(), + EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_p.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (const auto q : fe_eval_m.quadrature_point_indices()) + { + do_flux_term_ecl(fe_eval_m, fe_eval_p, 1.0, q); + + // clear lanes where face at boundary + fe_eval_m.begin_values()[q] *= fluxes; + for (unsigned int d = 0; d < dim; ++d) + fe_eval_m.begin_gradients()[q * dim + d] *= fluxes; + } + + fe_eval_m.integrate(EvaluationFlags::values | + EvaluationFlags::gradients, + fe_eval.begin_dof_values(), + true); + } + + fe_eval.distribute_local_to_global(dst); + } + }, + dst, + src, + true); + + QGauss quad_cell(n_q_points); + QGauss quad_face(n_q_points); + std::vector> quad_vec_cells; + quad_vec_cells.reserve( + (matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) * + n_lanes); + std::vector>> quad_vec_faces( + (matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) * + n_lanes); + + std::vector::cell_iterator> vector_accessors; + vector_accessors.reserve( + (matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) * + n_lanes); + for (unsigned int cell_batch = 0; + cell_batch < + matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches(); + ++cell_batch) + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (v < matrix_free.n_active_entries_per_cell_batch(cell_batch)) + vector_accessors.push_back( + matrix_free.get_cell_iterator(cell_batch, v)); + else + vector_accessors.push_back( + matrix_free.get_cell_iterator(cell_batch, 0)); + + quad_vec_cells.push_back(quad_cell); + + for (const auto f : GeometryInfo::face_indices()) + { + (void)f; + quad_vec_faces[cell_batch * n_lanes + v].push_back(quad_face); + } + } + + NonMatching::MappingInfo mapping_info_cells(mapping, + update_gradients | + update_JxW_values); + NonMatching::MappingInfo mapping_info_faces(mapping, + update_values | + update_gradients | + update_JxW_values | + update_normal_vectors); + + mapping_info_cells.reinit_cells(vector_accessors, quad_vec_cells); + mapping_info_faces.reinit_faces(vector_accessors, quad_vec_faces); + + FEPointEvaluation<1, dim, dim, double> fe_peval(mapping_info_cells, fe); + FEFacePointEvaluation<1, dim, dim, double> fe_peval_m(mapping_info_faces, fe); + FEFacePointEvaluation<1, dim, dim, double> fe_peval_p(mapping_info_faces, fe); + + matrix_free + .template loop_cell_centric, + LinearAlgebra::distributed::Vector>( + [&](const auto &matrix_free, + auto &dst, + const auto &src, + const auto &range) { + FEEvaluation fe_eval(matrix_free); + FEFaceEvaluation fe_eval_m(matrix_free, true); + FEFaceEvaluation fe_eval_p(matrix_free, false); + AlignedVector> vec_solution_values_in_m( + fe_eval.dofs_per_cell); + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + + for (unsigned int i = 0; i < fe_eval.dofs_per_cell; ++i) + vec_solution_values_in_m[i] = fe_eval.begin_dof_values()[i]; + + for (unsigned int v = 0; v < n_lanes; ++v) + { + fe_peval.reinit(cell * n_lanes + v); + fe_peval.evaluate(StridedArrayView( + &vec_solution_values_in_m[0][v], + fe.dofs_per_cell), + EvaluationFlags::gradients); + for (const unsigned int q : fe_peval.quadrature_point_indices()) + fe_peval.submit_gradient(fe_peval.get_gradient(q), q); + fe_peval.integrate(StridedArrayView( + &fe_eval.begin_dof_values()[0][v], + fe.dofs_per_cell), + EvaluationFlags::gradients); + } + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + fe_eval_m.reinit(cell, f); + + fe_eval_m.project_to_face(&vec_solution_values_in_m.begin()[0], + EvaluationFlags::values | + EvaluationFlags::gradients); + + // ask for boundary ids of face + const auto boundary_ids = + matrix_free.get_faces_by_cells_boundary_id(cell, f); + + // only internal faces have a neighbor, setup a mask + std::bitset mask; + for (unsigned int v = 0; v < n_lanes; ++v) + mask[v] = + boundary_ids[v] == numbers::internal_face_boundary_id; + + fe_eval_p.reinit(cell, f); + fe_eval_p.read_dof_values(src, 0, mask); + + fe_eval_p.project_to_face(EvaluationFlags::values | + EvaluationFlags::gradients); + + const auto &cell_indices_p = fe_eval_p.get_cell_ids(); + + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (mask[v] == false) + { + for (unsigned int i = 0; + i < 2 * fe_eval_m.dofs_per_face; + ++i) + fe_eval_m.get_scratch_data().begin()[i][v] = 0.; + continue; + } + + fe_peval_m.reinit(cell * n_lanes + v, f); + + fe_peval_p.reinit(cell_indices_p[v], + fe_eval_p.get_face_no(v)); + + fe_peval_m.evaluate_in_face( + &fe_eval_m.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + fe_peval_p.evaluate_in_face( + &fe_eval_p.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + + for (const unsigned int q : + fe_peval_m.quadrature_point_indices()) + do_flux_term_ecl(fe_peval_m, fe_peval_p, 1.0, q); + + fe_peval_m.integrate_in_face( + &fe_eval_m.get_scratch_data().begin()[0][v], + EvaluationFlags::values | EvaluationFlags::gradients); + } + + fe_eval_m.collect_from_face(EvaluationFlags::values | + EvaluationFlags::gradients, + fe_eval.begin_dof_values(), + true); + } + + fe_eval.distribute_local_to_global(dst); + } + }, + dst2, + src, + true); + + dst2 -= dst; + const double error = dst2.l2_norm() / dst.l2_norm(); + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + deallog << "FEPointEvaluation verification: " << error << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + initlog(); + + test_dg_ecl<2>(1, false); + deallog << std::endl; + test_dg_ecl<2>(2, false); + deallog << std::endl; + test_dg_ecl<3>(1, false); + deallog << std::endl; + test_dg_ecl<3>(2, false); + deallog << std::endl; + test_dg_ecl<2>(1, true); + deallog << std::endl; + test_dg_ecl<2>(2, true); + // TODO: fix face orientation + // deallog << std::endl; + // test_dg_ecl<3>(1, true); + // deallog << std::endl; + // test_dg_ecl<3>(2, true); +} diff --git a/tests/non_matching/mapping_info_05.output b/tests/non_matching/mapping_info_05.output new file mode 100644 index 0000000000..a3b3a1fc17 --- /dev/null +++ b/tests/non_matching/mapping_info_05.output @@ -0,0 +1,18 @@ + +DEAL::Working with FE_DGQ<2>(1) and 64 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<2>(2) and 144 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<3>(1) and 512 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<3>(2) and 1728 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<2>(1) and 96 dofs +DEAL::FEPointEvaluation verification: 0.0 +DEAL:: +DEAL::Working with FE_DGQ<2>(2) and 216 dofs +DEAL::FEPointEvaluation verification: 0.0