From: Peter Munch Date: Sun, 9 Aug 2020 17:49:39 +0000 (+0200) Subject: Merge implementation of gather_evaluate and integrate_scatter X-Git-Tag: v9.3.0-rc1~1201^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10811%2Fhead;p=dealii.git Merge implementation of gather_evaluate and integrate_scatter --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 88cc9bfef2..1a101a63d0 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2075,457 +2075,108 @@ namespace internal const unsigned int face_orientation, const Table<2, unsigned int> & orientation_map) { - const unsigned int side = face_no % 2; - - constexpr unsigned int static_dofs_per_component = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) : - numbers::invalid_unsigned_int; - constexpr unsigned int static_dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - numbers::invalid_unsigned_int; - const unsigned int dofs_per_face = - fe_degree > -1 ? - static_dofs_per_face : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); - - constexpr unsigned int stack_array_size_threshold = 100; - - VectorizedArrayType - temp_data[static_dofs_per_face < stack_array_size_threshold ? - n_components * 2 * dofs_per_face : - 1]; - VectorizedArrayType *__restrict temp1; - if (static_dofs_per_face < stack_array_size_threshold) - temp1 = &temp_data[0]; - else - temp1 = scratch_data; - - // case 1: contiguous and interleaved indices - if (((evaluate_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) - { - 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(); - - if (fe_degree > 1 && evaluate_gradients == true) - { - // we know that the gradient weights for the Hermite case on the - // right (side==1) are the negative from the value at the left - // (side==0), so we only read out one of them. - const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; - AssertDimension(data.face_to_cell_index_hermite.size(1), - 2 * dofs_per_face); - const unsigned int *index_array = - &data.face_to_cell_index_hermite(face_no, 0); - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = index_array[2 * i]; - const unsigned int ind2 = index_array[2 * i + 1]; - AssertIndexRange(ind1, data.dofs_per_component_on_cell); - AssertIndexRange(ind2, data.dofs_per_component_on_cell); - for (unsigned int comp = 0; comp < n_components; ++comp) - { - do_vectorized_read(src_ptr + dof_index + - (ind1 + - comp * static_dofs_per_component) * - VectorizedArrayType::size(), - temp1[i + 2 * comp * dofs_per_face]); - do_vectorized_read( - src_ptr + dof_index + - (ind2 + comp * static_dofs_per_component) * - VectorizedArrayType::size(), - temp1[dofs_per_face + i + 2 * comp * dofs_per_face]); - temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = - grad_weight * - (temp1[i + 2 * comp * dofs_per_face] - - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]); - } - } - } - else - { - AssertDimension(data.face_to_cell_index_nodal.size(1), - dofs_per_face); - const unsigned int *index_array = - &data.face_to_cell_index_nodal(face_no, 0); - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind = index_array[i]; - for (unsigned int comp = 0; comp < n_components; ++comp) - do_vectorized_read(src_ptr + dof_index + - (ind + - comp * static_dofs_per_component) * - VectorizedArrayType::size(), - temp1[i + 2 * comp * dofs_per_face]); - } - } - } - - // case 2: contiguous and interleaved indices with fixed stride - else if (((evaluate_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 && evaluate_gradients == true) - { - // we know that the gradient weights for the Hermite case on the - // right (side==1) are the negative from the value at the left - // (side==0), so we only read out one of them. - const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; - AssertDimension(data.face_to_cell_index_hermite.size(1), - 2 * dofs_per_face); - - const unsigned int *index_array = - &data.face_to_cell_index_hermite(face_no, 0); - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = - index_array[2 * i] * VectorizedArrayType::size(); - const unsigned int ind2 = - index_array[2 * i + 1] * VectorizedArrayType::size(); - for (unsigned int comp = 0; comp < n_components; ++comp) - { - do_vectorized_gather( - src_ptr + ind1 + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size(), - indices, - temp1[i + 2 * comp * dofs_per_face]); - VectorizedArrayType grad; - do_vectorized_gather( - src_ptr + ind2 + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size(), - indices, - grad); - temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = - grad_weight * - (temp1[i + 2 * comp * dofs_per_face] - grad); - } - } - } - else - { - AssertDimension(data.face_to_cell_index_nodal.size(1), - dofs_per_face); - const unsigned int *index_array = - &data.face_to_cell_index_nodal(face_no, 0); - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind = - index_array[i] * VectorizedArrayType::size(); - for (unsigned int comp = 0; comp < n_components; ++comp) - do_vectorized_gather( - src_ptr + ind + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size(), - indices, - temp1[i + 2 * comp * dofs_per_face]); - } - } - } - - // case 3: contiguous and interleaved indices with mixed stride - else if (((evaluate_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 nvec = - dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - - if (fe_degree > 1 && evaluate_gradients == true) - { - // we know that the gradient weights for the Hermite case on the - // right (side==1) are the negative from the value at the left - // (side==0), so we only read out one of them. - const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; - AssertDimension(data.face_to_cell_index_hermite.size(1), - 2 * dofs_per_face); - - const unsigned int *index_array = - &data.face_to_cell_index_hermite(face_no, 0); - if (nvec == VectorizedArrayType::size()) - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_face; ++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[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[2 * i + 1]) * - strides[v]; - do_vectorized_gather(src_ptr, - ind1, - temp1[i + 2 * comp * dofs_per_face]); - VectorizedArrayType grad; - do_vectorized_gather(src_ptr, ind2, grad); - temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = - grad_weight * - (temp1[i + 2 * comp * dofs_per_face] - grad); - } - else - { - for (unsigned int i = 0; i < n_components * 2 * dofs_per_face; - ++i) - temp1[i] = VectorizedArrayType(); - for (unsigned int v = 0; v < nvec; ++v) - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = - indices[v] + (comp * static_dofs_per_component + - index_array[2 * i]) * - strides[v]; - const unsigned int ind2 = - indices[v] + (comp * static_dofs_per_component + - index_array[2 * i + 1]) * - strides[v]; - temp1[i + 2 * comp * dofs_per_face][v] = - src_ptr[ind1]; - const Number grad = src_ptr[ind2]; - temp1[i + dofs_per_face + - 2 * comp * dofs_per_face][v] = - grad_weight[0] * - (temp1[i + 2 * comp * dofs_per_face][v] - grad); - } - } - } - else - { - AssertDimension(data.face_to_cell_index_nodal.size(1), - dofs_per_face); - const unsigned int *index_array = - &data.face_to_cell_index_nodal(face_no, 0); - if (nvec == VectorizedArrayType::size()) - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - 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[i]) * - strides[v]; - do_vectorized_gather(src_ptr, - ind, - temp1[i + 2 * comp * dofs_per_face]); - } - else - { - for (unsigned int i = 0; i < n_components * dofs_per_face; - ++i) - temp1[i] = VectorizedArrayType(); - for (unsigned int v = 0; v < nvec; ++v) - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = - indices[v] + (comp * static_dofs_per_component + - index_array[i]) * - strides[v]; - temp1[i + 2 * comp * dofs_per_face][v] = - src_ptr[ind1]; - } - } - } - } - - // case 4: contiguous indices without interleaving - else if (((evaluate_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 (evaluate_gradients == true && - data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite) - { - // we know that the gradient weights for the Hermite case on the - // right (side==1) are the negative from the value at the left - // (side==0), so we only read out one of them. - const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; - AssertDimension(data.face_to_cell_index_hermite.size(1), - 2 * dofs_per_face); - - const unsigned int *index_array = - &data.face_to_cell_index_hermite(face_no, 0); - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind1 = index_array[2 * i]; - const unsigned int ind2 = index_array[2 * i + 1]; - for (unsigned int comp = 0; comp < n_components; ++comp) - { - do_vectorized_gather( - src_ptr + ind1 + comp * static_dofs_per_component + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component], - indices, - temp1[i + 2 * comp * dofs_per_face]); - VectorizedArrayType grad; - do_vectorized_gather( - src_ptr + ind2 + comp * static_dofs_per_component + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component], - indices, - grad); - temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = - grad_weight * - (temp1[i + 2 * comp * dofs_per_face] - grad); - } - } - } + return process_and_io( // + false /*=evaluate*/, + src_ptr, + data, + dof_info, + values_quad, + gradients_quad, + scratch_data, + evaluate_values, + evaluate_gradients, + active_fe_index, + first_selected_component, + cell, + face_no, + subface_index, + dof_access_index, + face_orientation, + orientation_map, + [](auto & temp_1, + auto & temp_2, + const auto src_ptr_1, + const auto src_ptr_2, + const auto &grad_weight) { + // case 1a) + do_vectorized_read(src_ptr_1, temp_1); + do_vectorized_read(src_ptr_2, temp_2); + temp_2 = grad_weight * (temp_1 - temp_2); + }, + [](auto &temp, const auto src_ptr) { + // case 1b) + do_vectorized_read(src_ptr, temp); + }, + [](auto & temp_1, + auto & temp_2, + const auto src_ptr_1, + const auto src_ptr_2, + const auto &grad_weight, + const auto &indices_1, + const auto &indices_2) { + // case 2a) + do_vectorized_gather(src_ptr_1, indices_1, temp_1); + do_vectorized_gather(src_ptr_2, indices_2, temp_2); + temp_2 = grad_weight * (temp_1 - temp_2); + }, + [](auto &temp, const auto src_ptr, const auto &indices) { + // case 2b) + do_vectorized_gather(src_ptr, indices, temp); + }, + [](auto & temp_1, + auto & temp_2, + const auto &src_ptr_1, + const auto &src_ptr_2, + const auto &grad_weight) { + // case 3a) + temp_1 = src_ptr_1; + temp_2 = grad_weight * (temp_1 - src_ptr_2); + }, + [](auto &temp, const auto &src_ptr) { + // case 3b) + temp = src_ptr; + }, + [&](const auto &) { + // case 5) + }, + [&](auto &temp1, const auto &dofs_per_face) { + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::evaluate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + evaluate_values, + evaluate_gradients, + subface_index); else - { - AssertDimension(data.face_to_cell_index_nodal.size(1), - dofs_per_face); - const unsigned int *index_array = - &data.face_to_cell_index_nodal(face_no, 0); - for (unsigned int i = 0; i < dofs_per_face; ++i) - for (unsigned int comp = 0; comp < n_components; ++comp) - { - const unsigned int ind = index_array[i]; - do_vectorized_gather( - src_ptr + ind + comp * static_dofs_per_component + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component], - indices, - temp1[i + comp * 2 * dofs_per_face]); - } - } - } - - // case 5: default vector access - else - { - return false; - } - - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - data.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d, - n_components, - VectorizedArrayType>::evaluate_in_face(data, - temp1, - values_quad, - gradients_quad, - scratch_data + 2 * - n_components * - dofs_per_face, - evaluate_values, - evaluate_gradients, - subface_index); - else - FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d, - n_components, - VectorizedArrayType>::evaluate_in_face(data, - temp1, - values_quad, - gradients_quad, - scratch_data + 2 * - n_components * - dofs_per_face, - evaluate_values, - evaluate_gradients, - subface_index); - - if (face_orientation) - adjust_for_face_orientation(face_orientation, - orientation_map, - false, - evaluate_values, - evaluate_gradients, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad); - - return true; + FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::evaluate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + evaluate_values, + evaluate_gradients, + subface_index); + }); } static bool @@ -2548,82 +2199,211 @@ namespace internal const unsigned int face_orientation, const Table<2, unsigned int> & orientation_map) { - if (face_orientation) + return process_and_io( // + true /*=integrate*/, + dst_ptr, + data, + dof_info, + values_quad, + gradients_quad, + scratch_data, + integrate_values, + integrate_gradients, + active_fe_index, + first_selected_component, + cell, + face_no, + subface_index, + dof_access_index, + face_orientation, + orientation_map, + [](const auto &temp_1, + const auto &temp_2, + auto dst_ptr_1, + auto dst_ptr_2, + const auto &grad_weight) { + // case 1a) + const VectorizedArrayType val = temp_1 - grad_weight * temp_2; + const VectorizedArrayType grad = grad_weight * temp_2; + do_vectorized_add(val, dst_ptr_1); + do_vectorized_add(grad, dst_ptr_2); + }, + [](const auto &temp, auto dst_ptr) { + // case 1b) + do_vectorized_add(temp, dst_ptr); + }, + [](const auto &temp_1, + const auto &temp_2, + auto dst_ptr_1, + auto dst_ptr_2, + const auto &grad_weight, + const auto &indices_1, + const auto &indices_2) { + // case 2a) + const VectorizedArrayType val = temp_1 - grad_weight * temp_2; + const VectorizedArrayType grad = grad_weight * temp_2; + do_vectorized_scatter_add(val, indices_1, dst_ptr_1); + do_vectorized_scatter_add(grad, indices_2, dst_ptr_2); + }, + [](const auto &temp, auto dst_ptr, const auto &indices) { + // case 2b) + do_vectorized_scatter_add(temp, indices, dst_ptr); + }, + [](const auto &temp_1, + const auto &temp_2, + auto & dst_ptr_1, + auto & dst_ptr_2, + const auto &grad_weight) { + // case 3a) + const Number val = temp_1 - grad_weight * temp_2; + const Number grad = grad_weight * temp_2; + dst_ptr_1 += val; + dst_ptr_2 += grad; + }, + [](const auto &temp, auto &dst_ptr) { + // case 3b) + dst_ptr += temp; + }, + [&](const auto &temp1) { + // case 5: default vector access, must be handled separately, just do + // the face-normal interpolation + FEFaceNormalEvaluationImpl:: + template interpolate( + data, temp1, values_array, integrate_gradients, face_no); + }, + [&](auto &temp1, const auto &dofs_per_face) { + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= + internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl< + true, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::integrate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + integrate_values, + integrate_gradients, + subface_index); + else + internal::FEFaceEvaluationImpl< + false, + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArrayType>::integrate_in_face(data, + temp1, + values_quad, + gradients_quad, + scratch_data + + 2 * n_components * + dofs_per_face, + integrate_values, + integrate_gradients, + subface_index); + }); + } + + private: + template + static bool + process_and_io( + const bool integrate, + Number2_ * global_vector_ptr, + const MatrixFreeFunctions::ShapeInfo &data, + const MatrixFreeFunctions::DoFInfo & dof_info, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool do_values, + const bool do_gradients, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + const unsigned int cell, + const unsigned int face_no, + const unsigned int subface_index, + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, + const unsigned int face_orientation, + const Table<2, unsigned int> & orientation_map, + const Function1a & function_1a, + const Function1b & function_1b, + const Function2a & function_2a, + const Function2b & function_2b, + const Function3a & function_3a, + const Function3b & function_3b, + const Function5 & function_5, + const Function0 & function_0) + { + (void)subface_index; + + if (integrate && (face_orientation > 0)) adjust_for_face_orientation(face_orientation, orientation_map, true, - integrate_values, - integrate_gradients, + do_values, + do_gradients, data.n_q_points_face, scratch_data, values_quad, gradients_quad); - const unsigned int side = face_no % 2; + const unsigned int side = face_no % 2; + + const unsigned int side_ = integrate ? (2 - side) : (1 + side); + constexpr unsigned int static_dofs_per_component = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) : numbers::invalid_unsigned_int; + constexpr unsigned int static_dofs_per_face = + fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : + numbers::invalid_unsigned_int; const unsigned int dofs_per_face = fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : + static_dofs_per_face : Utilities::pow(data.data.front().fe_degree + 1, dim - 1); constexpr unsigned int stack_array_size_threshold = 100; - VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ? - n_components * 2 * dofs_per_face : - 1]; + VectorizedArrayType + temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components * 2 * dofs_per_face : + 1]; VectorizedArrayType *__restrict temp1; - if (dofs_per_face < stack_array_size_threshold) + if (static_dofs_per_face < stack_array_size_threshold) temp1 = &temp_data[0]; else temp1 = scratch_data; - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - data.element_type <= internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl< - true, - dim, - fe_degree, - n_q_points_1d, - n_components, - VectorizedArrayType>::integrate_in_face(data, - temp1, - values_quad, - gradients_quad, - scratch_data + - 2 * n_components * - dofs_per_face, - integrate_values, - integrate_gradients, - subface_index); - else - internal::FEFaceEvaluationImpl< - false, - dim, - fe_degree, - n_q_points_1d, - n_components, - VectorizedArrayType>::integrate_in_face(data, - temp1, - values_quad, - gradients_quad, - scratch_data + - 2 * n_components * - dofs_per_face, - integrate_values, - integrate_gradients, - subface_index); + if (integrate) + function_0(temp1, dofs_per_face); // case 1: contiguous and interleaved indices - if (((integrate_gradients == false && + if (((do_gradients == false && data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite && + MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && dof_info.index_storage_variants[dof_access_index][cell] == - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous) { AssertDimension( @@ -2637,13 +2417,13 @@ namespace internal [first_selected_component] * VectorizedArrayType::size(); - if (fe_degree > 1 && integrate_gradients == true) + if (fe_degree > 1 && do_gradients == true) { // we know that the gradient weights for the Hermite case on the // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + side_]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); const unsigned int *index_array = @@ -2655,25 +2435,16 @@ namespace internal AssertIndexRange(ind1, data.dofs_per_component_on_cell); AssertIndexRange(ind2, data.dofs_per_component_on_cell); for (unsigned int comp = 0; comp < n_components; ++comp) - { - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType grad = - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - do_vectorized_add(val, - dst_ptr + dof_index + - (ind1 + - comp * static_dofs_per_component) * - VectorizedArrayType::size()); - do_vectorized_add(grad, - dst_ptr + dof_index + - (ind2 + - comp * static_dofs_per_component) * - VectorizedArrayType::size()); - } + 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 @@ -2686,23 +2457,22 @@ namespace internal { const unsigned int ind = index_array[i]; for (unsigned int comp = 0; comp < n_components; ++comp) - do_vectorized_add(temp1[i + 2 * comp * dofs_per_face], - dst_ptr + dof_index + - (ind + - comp * static_dofs_per_component) * - VectorizedArrayType::size()); + function_1b(temp1[i + 2 * comp * dofs_per_face], + global_vector_ptr + dof_index + + (ind + comp * static_dofs_per_component) * + VectorizedArrayType::size()); } } } // case 2: contiguous and interleaved indices with fixed stride - else if (((integrate_gradients == false && + else if (((do_gradients == false && data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite && + MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && dof_info.index_storage_variants[dof_access_index][cell] == - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous_strided) { AssertDimension( @@ -2712,13 +2482,13 @@ namespace internal &dof_info .dof_indices_contiguous[dof_access_index] [cell * VectorizedArrayType::size()]; - if (fe_degree > 1 && integrate_gradients == true) + if (fe_degree > 1 && do_gradients == true) { // we know that the gradient weights for the Hermite case on the // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + side_]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2731,33 +2501,24 @@ namespace internal const unsigned int ind2 = index_array[2 * i + 1] * VectorizedArrayType::size(); for (unsigned int comp = 0; comp < n_components; ++comp) - { - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType grad = - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - do_vectorized_scatter_add( - val, - indices, - dst_ptr + ind1 + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size()); - do_vectorized_scatter_add( - grad, - indices, - dst_ptr + ind2 + - comp * static_dofs_per_component * - VectorizedArrayType::size() + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component] * - VectorizedArrayType::size()); - } + 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 @@ -2771,27 +2532,27 @@ namespace internal const unsigned int ind = index_array[i] * VectorizedArrayType::size(); for (unsigned int comp = 0; comp < n_components; ++comp) - do_vectorized_scatter_add( + function_2b( temp1[i + 2 * comp * dofs_per_face], - indices, - dst_ptr + ind + + 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()); + VectorizedArrayType::size(), + indices); } } } // case 3: contiguous and interleaved indices with mixed stride - else if (((integrate_gradients == false && + else if (((do_gradients == false && data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite && + MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && dof_info.index_storage_variants[dof_access_index][cell] == - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous_mixed_strides) { const unsigned int *strides = @@ -2808,13 +2569,13 @@ namespace internal const unsigned int n_filled_lanes = dof_info.n_vectorization_lanes_filled[dof_access_index][cell]; - if (fe_degree > 1 && integrate_gradients == true) + if (fe_degree > 1 && do_gradients == true) { // we know that the gradient weights for the Hermite case on the // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + side_]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2840,18 +2601,23 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[2 * i + 1]) * strides[v]; - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType grad = - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - do_vectorized_scatter_add(val, ind1, dst_ptr); - do_vectorized_scatter_add(grad, ind2, dst_ptr); + 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 { + 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) @@ -2864,15 +2630,12 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[2 * i + 1]) * strides[v]; - Number val = - temp1[i + 2 * comp * dofs_per_face][v] - - grad_weight[0] * temp1[i + dofs_per_face + - 2 * comp * dofs_per_face][v]; - Number grad = - grad_weight[0] * temp1[i + dofs_per_face + - 2 * comp * dofs_per_face][v]; - dst_ptr[ind1] += val; - dst_ptr[ind2] += grad; + function_3a(temp1[i + 2 * comp * dofs_per_face][v], + temp1[i + dofs_per_face + + 2 * comp * dofs_per_face][v], + global_vector_ptr[ind1], + global_vector_ptr[ind2], + grad_weight[0]); } } } @@ -2894,11 +2657,17 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[i]) * strides[v]; - do_vectorized_scatter_add( - temp1[i + 2 * comp * dofs_per_face], ind, dst_ptr); + function_2b(temp1[i + 2 * comp * dofs_per_face], + global_vector_ptr, + ind); } 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) @@ -2907,21 +2676,21 @@ namespace internal indices[v] + (comp * static_dofs_per_component + index_array[i]) * strides[v]; - dst_ptr[ind1] += - temp1[i + 2 * comp * dofs_per_face][v]; + function_3b(temp1[i + 2 * comp * dofs_per_face][v], + global_vector_ptr[ind1]); } } } } // case 4: contiguous indices without interleaving - else if (((integrate_gradients == false && + else if (((do_gradients == false && data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite && + MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && dof_info.index_storage_variants[dof_access_index][cell] == - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: contiguous && dof_info.n_vectorization_lanes_filled[dof_access_index][cell] == VectorizedArrayType::size()) @@ -2931,17 +2700,18 @@ namespace internal .dof_indices_contiguous[dof_access_index] [cell * VectorizedArrayType::size()]; - if (integrate_gradients == true && + if (do_gradients == true && data.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_hermite) + MatrixFreeFunctions::tensor_symmetric_hermite) { // we know that the gradient weights for the Hermite case on the // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + side_]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); + const unsigned int *index_array = &data.face_to_cell_index_hermite(face_no, 0); for (unsigned int i = 0; i < dofs_per_face; ++i) @@ -2949,27 +2719,20 @@ namespace internal const unsigned int ind1 = index_array[2 * i]; const unsigned int ind2 = index_array[2 * i + 1]; for (unsigned int comp = 0; comp < n_components; ++comp) - { - VectorizedArrayType val = - temp1[i + 2 * comp * dofs_per_face] - - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - VectorizedArrayType grad = - grad_weight * - temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; - do_vectorized_scatter_add( - val, - indices, - dst_ptr + comp * static_dofs_per_component + ind1 + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component]); - do_vectorized_scatter_add( - grad, - indices, - dst_ptr + comp * static_dofs_per_component + ind2 + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component]); - } + 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 @@ -2979,32 +2742,40 @@ namespace internal const unsigned int *index_array = &data.face_to_cell_index_nodal(face_no, 0); for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int ind = index_array[i]; - for (unsigned int comp = 0; comp < n_components; ++comp) - do_vectorized_scatter_add( - temp1[i + 2 * comp * dofs_per_face], - indices, - dst_ptr + comp * static_dofs_per_component + ind + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component]); - } + for (unsigned int comp = 0; comp < n_components; ++comp) + { + const unsigned int ind = index_array[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); + } } } - // case 5: default vector access, must be handled separately, just do - // the face-normal interpolation + // case 5: default vector access else { - FEFaceNormalEvaluationImpl:: - template interpolate( - data, temp1, values_array, integrate_gradients, face_no); + function_5(temp1); return false; } + if (!integrate) + function_0(temp1, dofs_per_face); + + if (!integrate && (face_orientation > 0)) + adjust_for_face_orientation(face_orientation, + orientation_map, + false, + do_values, + do_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + return true; } @@ -3051,7 +2822,7 @@ namespace internal } } } - }; + }; // namespace internal