From: Peter Munch Date: Wed, 12 Aug 2020 08:14:21 +0000 (+0200) Subject: ECL: enable contiguous data access X-Git-Tag: v9.3.0-rc1~1183^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=eb3aabd4b7447b0979201ef307a5aa906d0570a9;p=dealii.git ECL: enable contiguous data access --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 8c3ad32cac..1a06ec4a52 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2069,11 +2069,11 @@ namespace internal const bool evaluate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, - const unsigned int cell, - const std::array face_no, + const std::array cells, + const std::array face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array face_orientation, + const std::array face_orientations, const Table<2, unsigned int> & orientation_map) { return process_and_io( // @@ -2088,11 +2088,11 @@ namespace internal evaluate_gradients, active_fe_index, first_selected_component, - cell, - face_no, + cells, + face_nos, subface_index, dof_access_index, - face_orientation, + face_orientations, orientation_map, [](auto & temp_1, auto & temp_2, @@ -2194,11 +2194,11 @@ namespace internal const bool integrate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, - const unsigned int cell, - const std::array face_no, + const std::array cells, + const std::array face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array face_orientation, + const std::array face_orientations, const Table<2, unsigned int> & orientation_map) { return process_and_io( // @@ -2213,11 +2213,11 @@ namespace internal integrate_gradients, active_fe_index, first_selected_component, - cell, - face_no, + cells, + face_nos, subface_index, dof_access_index, - face_orientation, + face_orientations, orientation_map, [](const auto &temp_1, const auto &temp_2, @@ -2270,14 +2270,14 @@ namespace internal // case 5: default vector access, must be handled separately, just do // the face-normal interpolation - AssertDimension(face_no.size(), 1); + AssertDimension(face_nos.size(), 1); FEFaceNormalEvaluationImpl:: template interpolate( - data, temp1, values_array, integrate_gradients, face_no[0]); + data, temp1, values_array, integrate_gradients, face_nos[0]); }, [&](auto &temp1, const auto &dofs_per_face) { if (fe_degree > -1 && @@ -2344,11 +2344,11 @@ namespace internal const bool do_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, - const unsigned int cell, - const std::array face_no, + const std::array cells, + const std::array face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array face_orientation, + const std::array face_orientations, const Table<2, unsigned int> & orientation_map, const Function1a & function_1a, const Function1b & function_1b, @@ -2359,14 +2359,14 @@ namespace internal const Function5 & function_5, const Function0 & function_0) { - (void)subface_index; + const unsigned int cell = cells[0]; if (integrate && - (face_orientation[0] > 0 && + (face_orientations[0] > 0 && subface_index < GeometryInfo::max_children_per_cell)) { - AssertDimension(face_orientation.size(), 1); - adjust_for_face_orientation(face_orientation[0], + AssertDimension(face_orientations.size(), 1); + adjust_for_face_orientation(face_orientations[0], orientation_map, true, do_values, @@ -2385,8 +2385,8 @@ namespace internal data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ? data.data.front() .shape_data_on_face[0][fe_degree + (integrate ? - (2 - (face_no[0] % 2)) : - (1 + (face_no[0] % 2)))] : + (2 - (face_nos[0] % 2)) : + (1 + (face_nos[0] % 2)))] : VectorizedArrayType(0.0 /*dummy*/); constexpr unsigned int static_dofs_per_component = @@ -2419,388 +2419,340 @@ namespace internal // re-orientation std::array orientation; - orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ? - &data.face_orientations[face_orientation[0]][0] : - &dummy; + if (n_face_orientations == 1) + orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ? + &data.face_orientations[face_orientations[0]][0] : + &dummy; + else + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + // the loop breaks once an invalid_unsigned_int is hit for + // all cases except the exterior faces in the ECL loop (where + // some faces might be at the boundaries but others not) + if (cells[v] == numbers::invalid_unsigned_int) + continue; + + orientation[v] = + (data.data.front().nodal_at_cell_boundaries == true) ? + &data.face_orientations[face_orientations[v]][0] : + &dummy; + } + } // face_to_cell_index_hermite std::array index_array_hermite; - index_array_hermite[0] = - (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 && - data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ? - &data.face_to_cell_index_hermite(face_no[0], 0) : - &dummy; + if (n_face_orientations == 1) + index_array_hermite[0] = + (data.data.front().nodal_at_cell_boundaries == true && + fe_degree > 1 && + data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ? + &data.face_to_cell_index_hermite(face_nos[0], 0) : + &dummy; if (n_face_orientations > 1 && data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 && data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) { - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - for (unsigned int v = 1; v < n_filled_lanes; ++v) + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + grad_weight[v] = data.data.front().shape_data_on_face - [0][fe_degree + (integrate ? (2 - (face_no[v] % 2)) : - (1 + (face_no[v] % 2)))][v]; + [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) : + (1 + (face_nos[v] % 2)))][v]; index_array_hermite[v] = - &data.face_to_cell_index_hermite(face_no[v], 0); + &data.face_to_cell_index_hermite(face_nos[v], 0); } } // face_to_cell_index_nodal std::array index_array_nodal; - index_array_nodal[0] = - (data.data.front().nodal_at_cell_boundaries == true) ? - &data.face_to_cell_index_nodal(face_no[0], 0) : - &dummy; + if (n_face_orientations == 1) + index_array_nodal[0] = + (data.data.front().nodal_at_cell_boundaries == true) ? + &data.face_to_cell_index_nodal(face_nos[0], 0) : + &dummy; if (n_face_orientations > 1 && (data.data.front().nodal_at_cell_boundaries == true)) { - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + if (cells[v] == numbers::invalid_unsigned_int) + continue; - for (unsigned int v = 1; v < n_filled_lanes; ++v) - index_array_nodal[v] = - &data.face_to_cell_index_nodal(face_no[v], 0); + index_array_nodal[v] = + &data.face_to_cell_index_nodal(face_nos[v], 0); + } } const auto reorientate = [&](const unsigned int v, const unsigned int i) { - return (dim < 3 || face_orientation[0] == 0 || + return (dim < 3 || + face_orientations[n_face_orientations == 1 ? 0 : v] == 0 || subface_index < GeometryInfo::max_children_per_cell) ? i : orientation[v][i]; }; - // case 1: contiguous and interleaved indices - if (((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous) + if ((do_gradients == false && + data.data.front().nodal_at_cell_boundaries == true) || + (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) { - AssertDimension( - dof_info.n_vectorization_lanes_filled[dof_access_index][cell], - VectorizedArrayType::size()); - const unsigned int dof_index = - dof_info - .dof_indices_contiguous[dof_access_index] - [cell * VectorizedArrayType::size()] + - dof_info.component_dof_indices_offset[active_fe_index] + // case 1: contiguous and interleaved indices + if (n_face_orientations == 1 && + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous) + { + AssertDimension(n_face_orientations, 1); + + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::size()); + const unsigned int dof_index = + dof_info + .dof_indices_contiguous[dof_access_index] + [cell * VectorizedArrayType::size()] + + dof_info + .component_dof_indices_offset[active_fe_index] [first_selected_component] * - VectorizedArrayType::size(); + VectorizedArrayType::size(); - if (fe_degree > 1 && do_gradients == true) - { - for (unsigned int i = 0; i < dofs_per_face; ++i) + if (fe_degree > 1 && do_gradients == true) { - if (n_face_orientations == 1) - { - const unsigned int ind1 = index_array_hermite[0][2 * i]; - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1]; - AssertIndexRange(ind1, data.dofs_per_component_on_cell); - AssertIndexRange(ind2, data.dofs_per_component_on_cell); - const unsigned int i_ = reorientate(0, i); - for (unsigned int comp = 0; comp < n_components; ++comp) - function_1a( - temp1[i_ + 2 * comp * dofs_per_face], - temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], - global_vector_ptr + dof_index + - (ind1 + comp * static_dofs_per_component) * - VectorizedArrayType::size(), - global_vector_ptr + dof_index + - (ind2 + comp * static_dofs_per_component) * - VectorizedArrayType::size(), - grad_weight); - } - else + for (unsigned int i = 0; i < dofs_per_face; ++i) { - Assert(false, ExcNotImplemented()); - - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index][cell]; - - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3a( - temp1[reorientate(v, i) + 2 * comp * dofs_per_face] - [v], - temp1[reorientate(v, i) + dofs_per_face + - 2 * comp * dofs_per_face][v], - global_vector_ptr[dof_index + - (index_array_hermite[v][2 * i] + - comp * - static_dofs_per_component) * - VectorizedArrayType::size() + - v], - global_vector_ptr - [dof_index + - (index_array_hermite[v][2 * i + 1] + - comp * static_dofs_per_component) * - VectorizedArrayType::size() + - v], - grad_weight[v]); + if (n_face_orientations == 1) + { + const unsigned int ind1 = + index_array_hermite[0][2 * i]; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1]; + AssertIndexRange(ind1, + data.dofs_per_component_on_cell); + AssertIndexRange(ind2, + data.dofs_per_component_on_cell); + const unsigned int i_ = reorientate(0, i); + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_1a( + temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + + 2 * comp * dofs_per_face], + global_vector_ptr + dof_index + + (ind1 + comp * static_dofs_per_component) * + VectorizedArrayType::size(), + global_vector_ptr + dof_index + + (ind2 + comp * static_dofs_per_component) * + VectorizedArrayType::size(), + grad_weight); + } + else + { + Assert(false, ExcNotImplemented()); + } } } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = index_array_nodal[0][i]; - for (unsigned int comp = 0; comp < n_components; ++comp) - function_1b(temp1[i_ + 2 * comp * dofs_per_face], - global_vector_ptr + dof_index + - (ind + comp * static_dofs_per_component) * - VectorizedArrayType::size()); - } - else + for (unsigned int i = 0; i < dofs_per_face; ++i) { - Assert(false, ExcNotImplemented()); - - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index][cell]; - - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3b( - temp1[reorientate(v, i) + 2 * comp * dofs_per_face] - [v], - global_vector_ptr[dof_index + - (index_array_nodal[v][i] + - comp * - static_dofs_per_component) * - VectorizedArrayType::size() + - v]); + if (n_face_orientations == 1) + { + const unsigned int i_ = reorientate(0, i); + const unsigned int ind = index_array_nodal[0][i]; + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_1b(temp1[i_ + 2 * comp * dofs_per_face], + global_vector_ptr + dof_index + + (ind + + comp * static_dofs_per_component) * + VectorizedArrayType::size()); + } + else + { + Assert(false, ExcNotImplemented()); + } } } } - } - // case 2: contiguous and interleaved indices with fixed stride - else if (((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_strided) - { - AssertDimension( - dof_info.n_vectorization_lanes_filled[dof_access_index][cell], - VectorizedArrayType::size()); - const unsigned int *indices = - &dof_info - .dof_indices_contiguous[dof_access_index] - [cell * VectorizedArrayType::size()]; - if (fe_degree > 1 && do_gradients == true) + // case 2: contiguous and interleaved indices with fixed stride + else if (n_face_orientations == 1 && + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided) { - for (unsigned int i = 0; i < dofs_per_face; ++i) + AssertDimension(n_face_orientations, 1); + + AssertDimension( + dof_info.n_vectorization_lanes_filled[dof_access_index][cell], + VectorizedArrayType::size()); + const unsigned int *indices = + &dof_info + .dof_indices_contiguous[dof_access_index] + [cell * VectorizedArrayType::size()]; + if (fe_degree > 1 && do_gradients == true) { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - - const unsigned int ind1 = index_array_hermite[0][2 * i] * - VectorizedArrayType::size(); - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1] * - VectorizedArrayType::size(); - for (unsigned int comp = 0; comp < n_components; ++comp) - function_2a( - temp1[i_ + 2 * comp * dofs_per_face], - temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], - global_vector_ptr + ind1 + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size(), - global_vector_ptr + ind2 + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size(), - grad_weight, - indices, - indices); - } - else + for (unsigned int i = 0; i < dofs_per_face; ++i) { - Assert(false, ExcNotImplemented()); - - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index][cell]; + if (n_face_orientations == 1) + { + const unsigned int i_ = reorientate(0, i); - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int comp = 0; comp < n_components; ++comp) - { - const unsigned int i_ = reorientate(0, i); - function_3a(temp1[i_ + 2 * comp * dofs_per_face][v], + const unsigned int ind1 = + index_array_hermite[0][2 * i] * + VectorizedArrayType::size(); + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1] * + VectorizedArrayType::size(); + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_2a(temp1[i_ + 2 * comp * dofs_per_face], temp1[i_ + dofs_per_face + - 2 * comp * dofs_per_face][v], - global_vector_ptr - [index_array_hermite[v][2 * i] * - VectorizedArrayType::size() + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index] - [first_selected_component] * - VectorizedArrayType::size() + - indices[v] + v], - global_vector_ptr - [index_array_hermite[v][2 * i + 1] * - VectorizedArrayType::size() + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index] - [first_selected_component] * - VectorizedArrayType::size() + - indices[v] + v], - grad_weight[v]); - } + 2 * comp * dofs_per_face], + global_vector_ptr + ind1 + + comp * static_dofs_per_component * + VectorizedArrayType::size() + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] * + VectorizedArrayType::size(), + global_vector_ptr + ind2 + + comp * static_dofs_per_component * + VectorizedArrayType::size() + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] * + VectorizedArrayType::size(), + grad_weight, + indices, + indices); + } + else + { + Assert(false, ExcNotImplemented()); + } } } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) + else { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - - const unsigned int ind = - index_array_nodal[0][i] * VectorizedArrayType::size(); - for (unsigned int comp = 0; comp < n_components; ++comp) - function_2b( - temp1[i_ + 2 * comp * dofs_per_face], - global_vector_ptr + ind + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size(), - indices); - } - else - { - Assert(false, ExcNotImplemented()); - - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index][cell]; - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int comp = 0; comp < n_components; ++comp) - function_3b( - temp1[reorientate(v, i) + 2 * comp * dofs_per_face] - [v], - global_vector_ptr - [index_array_nodal[v][i] + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size() + - indices[v] + v]); - } - } - } - } - - // case 3: contiguous and interleaved indices with mixed stride - else if (((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_mixed_strides) - { - const unsigned int *strides = - &dof_info.dof_indices_interleave_strides - [dof_access_index][cell * VectorizedArrayType::size()]; - unsigned int indices[VectorizedArrayType::size()]; - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - indices[v] = - dof_info.dof_indices_contiguous - [dof_access_index][cell * VectorizedArrayType::size() + v] + - dof_info.component_dof_indices_offset[active_fe_index] - [first_selected_component] * - strides[v]; - const unsigned int n_filled_lanes = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - if (fe_degree > 1 && do_gradients == true) - { - if (n_filled_lanes == VectorizedArrayType::size()) - for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_face; ++i) { if (n_face_orientations == 1) { const unsigned int i_ = reorientate(0, i); - unsigned int ind1[VectorizedArrayType::size()]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++v) - ind1[v] = indices[v] + - (comp * static_dofs_per_component + - index_array_hermite[0 /*TODO*/][2 * i]) * - strides[v]; - unsigned int ind2[VectorizedArrayType::size()]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++v) - ind2[v] = - indices[v] + - (comp * static_dofs_per_component + - index_array_hermite[0 /*TODO*/][2 * i + 1]) * - strides[v]; - function_2a(temp1[i_ + 2 * comp * dofs_per_face], - temp1[i_ + dofs_per_face + - 2 * comp * dofs_per_face], - global_vector_ptr, - global_vector_ptr, - grad_weight, - ind1, - ind2); + + const unsigned int ind = index_array_nodal[0][i] * + VectorizedArrayType::size(); + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_2b(temp1[i_ + 2 * comp * dofs_per_face], + global_vector_ptr + ind + + comp * static_dofs_per_component * + VectorizedArrayType::size() + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] * + VectorizedArrayType::size(), + indices); } else { Assert(false, ExcNotImplemented()); + } + } + } + } - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index] - [cell]; + // case 3: contiguous and interleaved indices with mixed stride + else if (n_face_orientations == 1 && + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_mixed_strides) + { + AssertDimension(n_face_orientations, 1); + + const unsigned int *strides = + &dof_info.dof_indices_interleave_strides + [dof_access_index][cell * VectorizedArrayType::size()]; + unsigned int indices[VectorizedArrayType::size()]; + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + indices[v] = + dof_info.dof_indices_contiguous + [dof_access_index][cell * VectorizedArrayType::size() + v] + + dof_info + .component_dof_indices_offset[active_fe_index] + [first_selected_component] * + strides[v]; + const unsigned int n_filled_lanes = + dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (fe_degree > 1 && do_gradients == true) + { + if (n_filled_lanes == VectorizedArrayType::size()) + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + if (n_face_orientations == 1) + { + const unsigned int i_ = reorientate(0, i); + unsigned int ind1[VectorizedArrayType::size()]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++v) + ind1[v] = + indices[v] + + (comp * static_dofs_per_component + + index_array_hermite[0 /*TODO*/][2 * i]) * + strides[v]; + unsigned int ind2[VectorizedArrayType::size()]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++v) + ind2[v] = + indices[v] + + (comp * static_dofs_per_component + + index_array_hermite[0 /*TODO*/][2 * i + 1]) * + strides[v]; + function_2a(temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + + 2 * comp * dofs_per_face], + global_vector_ptr, + global_vector_ptr, + grad_weight, + ind1, + ind2); + } + else { - const unsigned int i_ = reorientate(v, i); + Assert(false, ExcNotImplemented()); + } + } + else + { + if (integrate == false) + for (unsigned int i = 0; + i < n_components * 2 * dofs_per_face; + ++i) + temp1[i] = VectorizedArrayType(); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int i_ = + reorientate(n_face_orientations == 1 ? 0 : v, + i); function_3a( temp1[i_ + 2 * comp * dofs_per_face][v], temp1[i_ + dofs_per_face + @@ -2808,241 +2760,280 @@ namespace internal global_vector_ptr [indices[v] + (comp * static_dofs_per_component + - index_array_hermite[v][2 * i]) * + index_array_hermite + [n_face_orientations == 1 ? 0 : v] + [2 * i]) * strides[v]], global_vector_ptr [indices[v] + (comp * static_dofs_per_component + - index_array_hermite[v][2 * i + 1]) * + index_array_hermite + [n_face_orientations == 1 ? 0 : v] + [2 * i + 1]) * strides[v]], - grad_weight[v]); + grad_weight[n_face_orientations == 1 ? 0 : v]); } - } } + } else { - if (integrate == false) - for (unsigned int i = 0; - i < n_components * 2 * dofs_per_face; - ++i) - temp1[i] = VectorizedArrayType(); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) + if (n_filled_lanes == VectorizedArrayType::size()) for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_face; ++i) { - const unsigned int i_ = - reorientate(n_face_orientations == 1 ? 0 : v, i); - function_3a( - temp1[i_ + 2 * comp * dofs_per_face][v], - temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face] - [v], - global_vector_ptr - [indices[v] + - (comp * static_dofs_per_component + - index_array_hermite - [n_face_orientations == 1 ? 0 : v][2 * i]) * - strides[v]], - global_vector_ptr - [indices[v] + - (comp * static_dofs_per_component + - index_array_hermite[n_face_orientations == 1 ? - 0 : - v][2 * i + 1]) * - strides[v]], - grad_weight[n_face_orientations == 1 ? 0 : v]); + if (n_face_orientations == 1) + { + unsigned int ind[VectorizedArrayType::size()]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++v) + ind[v] = indices[v] + + (comp * static_dofs_per_component + + index_array_nodal[0][i]) * + strides[v]; + const unsigned int i_ = reorientate(0, i); + function_2b(temp1[i_ + 2 * comp * dofs_per_face], + global_vector_ptr, + ind); + } + else + { + Assert(false, ExcNotImplemented()); + } } + else + { + if (integrate == false) + for (unsigned int i = 0; + i < n_components * dofs_per_face; + ++i) + temp1[i] = VectorizedArrayType(); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + function_3b( + temp1[reorientate( + n_face_orientations == 1 ? 0 : v, i) + + 2 * comp * dofs_per_face][v], + global_vector_ptr + [indices[v] + + (comp * static_dofs_per_component + + index_array_nodal + [n_face_orientations == 1 ? 0 : v][i]) * + strides[v]]); + } } } - else + + // case 4: contiguous indices without interleaving + else if (n_face_orientations > 1 || + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous) { - if (n_filled_lanes == VectorizedArrayType::size()) - for (unsigned int comp = 0; comp < n_components; ++comp) + const unsigned int *indices = + &dof_info + .dof_indices_contiguous[dof_access_index] + [cell * VectorizedArrayType::size()]; + + if (do_gradients == true && + data.element_type == + MatrixFreeFunctions::tensor_symmetric_hermite) + { for (unsigned int i = 0; i < dofs_per_face; ++i) { - if (n_face_orientations == 1) + if (n_face_orientations == 1 && + dof_info + .n_vectorization_lanes_filled[dof_access_index] + [cell] == + VectorizedArrayType::size()) { - unsigned int ind[VectorizedArrayType::size()]; - DEAL_II_OPENMP_SIMD_PRAGMA - for (unsigned int v = 0; - v < VectorizedArrayType::size(); - ++v) - ind[v] = - indices[v] + (comp * static_dofs_per_component + - index_array_nodal[0 /*TODO*/][i]) * - strides[v]; + const unsigned int ind1 = + index_array_hermite[0][2 * i]; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1]; const unsigned int i_ = reorientate(0, i); - function_2b(temp1[i_ + 2 * comp * dofs_per_face], - global_vector_ptr, - ind); + + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_2a( + temp1[i_ + 2 * comp * dofs_per_face], + temp1[i_ + dofs_per_face + + 2 * comp * dofs_per_face], + global_vector_ptr + + comp * static_dofs_per_component + ind1 + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component], + global_vector_ptr + + comp * static_dofs_per_component + ind2 + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component], + grad_weight, + indices, + indices); + } + else if (n_face_orientations == 1) + { + const unsigned int ind1 = + index_array_hermite[0][2 * i]; + const unsigned int ind2 = + index_array_hermite[0][2 * i + 1]; + const unsigned int i_ = reorientate(0, i); + + const unsigned int n_filled_lanes = + dof_info + .n_vectorization_lanes_filled[dof_access_index] + [cell]; + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_3a( + temp1[i_ + 2 * comp * dofs_per_face][v], + temp1[i_ + dofs_per_face + + 2 * comp * dofs_per_face][v], + global_vector_ptr + [comp * static_dofs_per_component + ind1 + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] + + indices[v]], + global_vector_ptr + [comp * static_dofs_per_component + ind2 + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] + + indices[v]], + grad_weight[v]); + + if (integrate == false) + for (unsigned int v = n_filled_lanes; + v < VectorizedArrayType::size(); + ++v) + for (unsigned int comp = 0; comp < n_components; + ++comp) + { + temp1[i_ + 2 * comp * dofs_per_face][v] = 0.0; + temp1[i_ + dofs_per_face + + 2 * comp * dofs_per_face][v] = 0.0; + } } else { Assert(false, ExcNotImplemented()); - const unsigned int n_filled_lanes = dof_info .n_vectorization_lanes_filled[dof_access_index] [cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3b( - temp1[reorientate(v, i) + - 2 * comp * dofs_per_face][v], - global_vector_ptr[indices[v] + - (comp * - static_dofs_per_component + - index_array_nodal[v][i]) * - strides[v]]); + for (unsigned int comp = 0; comp < n_components; + ++comp) + function_3a( + temp1[reorientate(v, i) + + 2 * comp * dofs_per_face][v], + temp1[reorientate(v, i) + dofs_per_face + + 2 * comp * dofs_per_face][v], + global_vector_ptr + [comp * static_dofs_per_component + + index_array_hermite[v][2 * i] + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] + + indices[v]], + global_vector_ptr + [comp * static_dofs_per_component + + index_array_hermite[v][2 * i + 1] + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] + + indices[v]], + grad_weight[v]); } } + } else { - if (integrate == false) - for (unsigned int i = 0; i < n_components * dofs_per_face; - ++i) - temp1[i] = VectorizedArrayType(); - - for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int i = 0; i < dofs_per_face; ++i) for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_face; ++i) - function_3b( - temp1[reorientate(n_face_orientations == 1 ? 0 : v, - i) + - 2 * comp * dofs_per_face][v], - global_vector_ptr - [indices[v] + - (comp * static_dofs_per_component + - index_array_nodal - [n_face_orientations == 1 ? 0 : v][i]) * - strides[v]]); - } - } - } - - // case 4: contiguous indices without interleaving - else if (((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) && - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - contiguous && - dof_info.n_vectorization_lanes_filled[dof_access_index][cell] == - VectorizedArrayType::size()) - { - const unsigned int *indices = - &dof_info - .dof_indices_contiguous[dof_access_index] - [cell * VectorizedArrayType::size()]; - - if (do_gradients == true && - data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite) - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - if (n_face_orientations == 1) - { - const unsigned int ind1 = index_array_hermite[0][2 * i]; - const unsigned int ind2 = - index_array_hermite[0][2 * i + 1]; - const unsigned int i_ = reorientate(0, i); - - for (unsigned int comp = 0; comp < n_components; ++comp) - function_2a( - temp1[i_ + 2 * comp * dofs_per_face], - temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face], - global_vector_ptr + comp * static_dofs_per_component + - ind1 + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component], - global_vector_ptr + comp * static_dofs_per_component + - ind2 + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component], - grad_weight, - indices, - indices); - } - else - { - Assert(false, ExcNotImplemented()); + { + if (n_face_orientations == 1 && + dof_info + .n_vectorization_lanes_filled[dof_access_index] + [cell] == + VectorizedArrayType::size()) + { + const unsigned int ind = index_array_nodal[0][i]; + const unsigned int i_ = reorientate(0, i); + + function_2b( + temp1[i_ + 2 * comp * dofs_per_face], + global_vector_ptr + + comp * static_dofs_per_component + ind + + dof_info.component_dof_indices_offset + [active_fe_index][first_selected_component], + indices); + } + else if (n_face_orientations == 1) + { + const unsigned int ind = index_array_nodal[0][i]; + const unsigned int i_ = reorientate(0, i); - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index][cell]; + const unsigned int n_filled_lanes = + dof_info + .n_vectorization_lanes_filled[dof_access_index] + [cell]; - for (unsigned int v = 0; v < n_filled_lanes; ++v) - for (unsigned int comp = 0; comp < n_components; ++comp) - function_3a( - temp1[reorientate(v, i) + 2 * comp * dofs_per_face] - [v], - temp1[reorientate(v, i) + dofs_per_face + - 2 * comp * dofs_per_face][v], - global_vector_ptr - [comp * static_dofs_per_component + - index_array_hermite[v][2 * i] + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] + - indices[v]], - global_vector_ptr - [comp * static_dofs_per_component + - index_array_hermite[v][2 * i + 1] + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] + - indices[v]], - grad_weight[v]); - } + for (unsigned int v = 0; v < n_filled_lanes; ++v) + function_3b( + temp1[i_ + 2 * comp * dofs_per_face][v], + global_vector_ptr + [comp * static_dofs_per_component + ind + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] + + indices[v]]); + + if (integrate == false) + for (unsigned int v = n_filled_lanes; + v < VectorizedArrayType::size(); + ++v) + temp1[i_ + 2 * comp * dofs_per_face][v] = 0.0; + } + else + { + for (unsigned int v = 0; + v < VectorizedArrayType::size(); + ++v) + if (cells[v] != numbers::invalid_unsigned_int) + function_3b( + temp1[reorientate(v, i) + + 2 * comp * dofs_per_face][v], + global_vector_ptr + [index_array_nodal[v][i] + + comp * static_dofs_per_component + + dof_info.component_dof_indices_offset + [active_fe_index] + [first_selected_component] + + dof_info.dof_indices_contiguous + [dof_access_index][cells[v]]]); + } + } } } else { - for (unsigned int i = 0; i < dofs_per_face; ++i) - for (unsigned int comp = 0; comp < n_components; ++comp) - { - if (n_face_orientations == 1) - { - const unsigned int ind = index_array_nodal[0][i]; - const unsigned int i_ = reorientate(0, i); - - function_2b( - temp1[i_ + 2 * comp * dofs_per_face], - global_vector_ptr + comp * static_dofs_per_component + - ind + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component], - indices); - } - else - { - Assert(false, ExcNotImplemented()); - - const unsigned int n_filled_lanes = - dof_info - .n_vectorization_lanes_filled[dof_access_index] - [cell]; - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3b( - temp1[reorientate(v, i) + 2 * comp * dofs_per_face] - [v], - global_vector_ptr - [comp * static_dofs_per_component + - index_array_nodal[v][i] + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] + - indices[v]]); - } - } + // case 5: default vector access + function_5(temp1); + return false; } } - - // case 5: default vector access else { + // case 5: default vector access function_5(temp1); return false; } @@ -3051,11 +3042,11 @@ namespace internal function_0(temp1, dofs_per_face); if (!integrate && - (face_orientation[0] > 0 && + (face_orientations[0] > 0 && subface_index < GeometryInfo::max_children_per_cell)) { - AssertDimension(face_orientation.size(), 1); - adjust_for_face_orientation(face_orientation[0], + AssertDimension(face_orientations.size(), 1); + adjust_for_face_orientation(face_orientations[0], orientation_map, false, do_values, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 488637848f..21136c51f8 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -754,6 +754,14 @@ public: std::array get_cell_ids() const; + + /** + * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is + * associated with. + */ + std::array + get_cell_or_face_ids() const; + //@} /** @@ -3142,6 +3150,20 @@ public: * points is inaccurate and this value must be used instead. */ const unsigned int n_q_points; + + +private: + /** + * Return face number of each face of the current face batch. + */ + std::array + compute_face_no_data(); + + /** + * Determine the orientation of each face of the current face batch. + */ + std::array + compute_face_orientations(); }; @@ -3849,6 +3871,65 @@ namespace internal +template +std::array +FEEvaluationBase:: + get_cell_or_face_ids() const +{ + const unsigned int v_len = VectorizedArrayType::size(); + std::array cells; + + // initialize array + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = numbers::invalid_unsigned_int; + + if (is_face && + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false) + { + // cell-based face-loop: plus face + for (unsigned int i = 0; i < v_len; i++) + { + // compute actual (non vectorized) cell ID + const unsigned int cell_this = this->cell * v_len + i; + // compute face ID + unsigned int fn = + this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, + this->face_no, + i); + + if (fn == numbers::invalid_unsigned_int) + continue; // invalid face ID: no neighbor on boundary + + // get cell ID on both sides of face + auto cell_m = this->matrix_info->get_face_info(fn / v_len) + .cells_interior[fn % v_len]; + auto cell_p = this->matrix_info->get_face_info(fn / v_len) + .cells_exterior[fn % v_len]; + + // compare the IDs with the given cell ID + if (cell_m == cell_this) + cells[i] = cell_p; // neighbor has the other ID + else if (cell_p == cell_this) + cells[i] = cell_m; + } + } + else + { + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = cell * v_len + i; + } + + return cells; +} + + + template :: - gather_evaluate(input_vector.begin(), - *this->data, - *this->dof_info, - this->begin_values(), - this->begin_gradients(), - this->scratch_data, - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - this->active_fe_index, - this->first_selected_component, - this->cell, - std::array{{this->face_no}}, - this->subface_index, - this->dof_access_index, - std::array{{this->face_orientation}}, - this->mapping_data->descriptor[this->active_fe_index] - .face_orientations)) + const auto fu = [&]() { + if (this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false) + { + const auto cells = this->get_cell_or_face_ids(); + const auto face_nos = this->compute_face_no_data(); + const auto face_orientations = this->compute_face_orientations(); + + return internal::FEFaceEvaluationSelector:: + gather_evaluate(input_vector.begin(), + *this->data, + *this->dof_info, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluation_flag & EvaluationFlags::values, + evaluation_flag & EvaluationFlags::gradients, + this->active_fe_index, + this->first_selected_component, + cells, + face_nos, + this->subface_index, + this->dof_access_index, + face_orientations, + this->mapping_data->descriptor[this->active_fe_index] + .face_orientations); + } + else + { + return internal::FEFaceEvaluationSelector:: + gather_evaluate(input_vector.begin(), + *this->data, + *this->dof_info, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluation_flag & EvaluationFlags::values, + evaluation_flag & EvaluationFlags::gradients, + this->active_fe_index, + this->first_selected_component, + std::array{{this->cell}}, + std::array{{this->face_no}}, + this->subface_index, + this->dof_access_index, + std::array{{this->face_orientation}}, + this->mapping_data->descriptor[this->active_fe_index] + .face_orientations); + } + }; + + if (!fu()) { this->read_dof_values(input_vector); this->evaluate(evaluation_flag); @@ -8235,6 +8355,11 @@ FEFaceEvaluationdof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false) == false, + ExcNotImplemented()); + if (!internal::FEFaceEvaluationSelectoractive_fe_index, this->first_selected_component, - this->cell, + std::array{{this->cell}}, std::array{{this->face_no}}, this->subface_index, this->dof_access_index, @@ -8318,6 +8443,161 @@ FEFaceEvaluation +std::array +FEFaceEvaluation::compute_face_no_data() +{ + std::array face_no_data; + + if (this->dof_access_index != + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell || + this->is_interior_face == true) + { + std::fill(face_no_data.begin(), + face_no_data.begin() + + this->dof_info->n_vectorization_lanes_filled + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell], + this->face_no); + } + else + { + std::fill(face_no_data.begin(), + face_no_data.end(), + numbers::invalid_unsigned_int); + + for (unsigned int i = 0; + i < this->dof_info->n_vectorization_lanes_filled + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell]; + i++) + { + // compute actual (non vectorized) cell ID + const unsigned int cell_this = + this->cell * VectorizedArrayType::size() + i; + // compute face ID + const unsigned int face_index = + this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, + this->face_no, + i); + + Assert(face_index != numbers::invalid_unsigned_int, + ExcNotInitialized()); + + // get cell ID on both sides of face + auto cell_m = + this->matrix_info + ->get_face_info(face_index / VectorizedArrayType::size()) + .cells_interior[face_index % VectorizedArrayType::size()]; + + // compare the IDs with the given cell ID + face_no_data[i] = + (cell_m == cell_this) ? + this->matrix_info + ->get_face_info(face_index / VectorizedArrayType::size()) + .exterior_face_no : + this->matrix_info + ->get_face_info(face_index / VectorizedArrayType::size()) + .interior_face_no; + } + } + + return face_no_data; +} + + + +template +std::array +FEFaceEvaluation::compute_face_orientations() +{ + std::array face_no_data; + + if (this->dof_access_index != + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell || + this->is_interior_face == true) + { + Assert(false, ExcNotImplemented()); + } + else + { + std::fill(face_no_data.begin(), face_no_data.end(), 0); + + if (dim == 3) + { + for (unsigned int i = 0; + i < this->dof_info->n_vectorization_lanes_filled + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell]; + i++) + { + // compute actual (non vectorized) cell ID + const unsigned int cell_this = + this->cell * VectorizedArrayType::size() + i; + // compute face ID + const unsigned int face_index = + this->matrix_info->get_cell_and_face_to_plain_faces()( + this->cell, this->face_no, i); + + Assert(face_index != numbers::invalid_unsigned_int, + ExcNotInitialized()); + + const unsigned int macro = + face_index / VectorizedArrayType::size(); + const unsigned int lane = + face_index % VectorizedArrayType::size(); + + const auto &faces = this->matrix_info->get_face_info(macro); + + // get cell ID on both sides of face + auto cell_m = faces.cells_interior[lane]; + + const bool is_interior_face = cell_m != cell_this; + const bool fo_interior_face = faces.face_orientation >= 8; + + unsigned int face_orientation = faces.face_orientation % 8; + + if (is_interior_face != fo_interior_face) + { + // invert (see also: + // Triangulation::update_periodic_face_map()) + static const std::array table{ + {0, 1, 0, 3, 6, 5, 4, 7}}; + + face_orientation = table[face_orientation]; + } + + // compare the IDs with the given cell ID + face_no_data[i] = face_orientation; + } + } + } + + return face_no_data; +} + + + /*------------------------- end FEFaceEvaluation ------------------------- */ diff --git a/tests/matrix_free/ecl_02.cc b/tests/matrix_free/ecl_02.cc new file mode 100644 index 0000000000..8a0df9a52b --- /dev/null +++ b/tests/matrix_free/ecl_02.cc @@ -0,0 +1,153 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + +// tests matrix-free read_cell_data for locally refined mesh + +template > +void +test(const unsigned int n_refinements = 1) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + GridGenerator::subdivided_hyper_rectangle( + tria, {2, 1}, Point(0.0, 0.0), Point(2.0, 1.0), true); + + std::vector::cell_iterator>> + periodic_faces; + + if (dim >= 1) + dealii::GridTools::collect_periodic_faces(tria, 0, 1, 0, periodic_faces); + + if (dim >= 2) + dealii::GridTools::collect_periodic_faces(tria, 2, 3, 1, periodic_faces); + + tria.add_periodicity(periodic_faces); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + QGauss<1> quad(n_points); + + AffineConstraints constraint; + + using MF = MatrixFree; + + typename MF::AdditionalData additional_data; + additional_data.mapping_update_flags = update_values; + additional_data.mapping_update_flags_inner_faces = update_values; + additional_data.mapping_update_flags_boundary_faces = update_values; + additional_data.mapping_update_flags_faces_by_cells = update_values; + additional_data.hold_all_faces_to_owned_cells = true; + + MF matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + FEEvaluation phi( + matrix_free); + FEFaceEvaluation + phi_m(matrix_free, true); + FEFaceEvaluation + phi_p(matrix_free, false); + + + for (unsigned int i = 0; i < src.size() / 2; i++) + src[i] = 1; + + for (unsigned int i = src.size() / 2; i < src.size(); i++) + src[i] = 2; + + dst = 0.0; + + /** + * Element-centric loop + */ + matrix_free.template loop_cell_centric( + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + face++) + { + phi_m.reinit(cell, face); + phi_p.reinit(cell, face); + + phi_m.read_dof_values(src); + + for (unsigned int i = 0; i < phi_m.static_dofs_per_component; i++) + deallog << static_cast(phi_m.begin_dof_values()[i][0]) + << " "; + deallog << std::endl; + phi_m.gather_evaluate(src, true, false); + for (unsigned int i = 0; i < phi_m.static_n_q_points; i++) + deallog << static_cast(phi_m.begin_values()[i][0]) << " "; + deallog << std::endl; + + phi_p.read_dof_values(src); + for (unsigned int i = 0; i < phi_p.static_dofs_per_component; i++) + deallog << static_cast(phi_p.begin_dof_values()[i][0]) + << " "; + deallog << std::endl; + + phi_p.gather_evaluate(src, true, false); + for (unsigned int i = 0; i < phi_p.static_n_q_points; i++) + deallog << static_cast(phi_p.begin_values()[i][0]) << " "; + deallog << std::endl << std::endl; + } + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::gradients); +} + +int +main() +{ + initlog(); + test<2, 1, 2, double, VectorizedArray>(); +} diff --git a/tests/matrix_free/ecl_02.output b/tests/matrix_free/ecl_02.output new file mode 100644 index 0000000000..d584d53cc8 --- /dev/null +++ b/tests/matrix_free/ecl_02.output @@ -0,0 +1,21 @@ + +DEAL::1 1 1 1 +DEAL::1 1 +DEAL::1 1 1 1 +DEAL::2 2 +DEAL:: +DEAL::1 1 1 1 +DEAL::1 1 +DEAL::1 1 1 1 +DEAL::2 2 +DEAL:: +DEAL::1 1 1 1 +DEAL::1 1 +DEAL::1 1 1 1 +DEAL::1 1 +DEAL:: +DEAL::1 1 1 1 +DEAL::1 1 +DEAL::1 1 1 1 +DEAL::1 1 +DEAL::