From 1e85c0e2a88f22602c42192b9d9fbffbc5909378 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 16 Sep 2020 16:17:28 +0200 Subject: [PATCH] Refactor access to FEEvaluationImpl and friends --- .../deal.II/matrix_free/evaluation_kernels.h | 1977 ++++++++++------- .../deal.II/matrix_free/evaluation_selector.h | 220 +- include/deal.II/matrix_free/fe_evaluation.h | 181 +- .../matrix_free/mapping_info.templates.h | 65 +- 4 files changed, 1289 insertions(+), 1154 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index b630b89f2f..da6a399cc8 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -1377,6 +1377,300 @@ namespace internal + /** + * This class chooses an appropriate evaluation strategy based on the + * template parameters and the shape_info variable which contains runtime + * parameters for the strategy underlying FEEvaluation::evaluate(), i.e. + * this calls internal::FEEvaluationImpl::evaluate(), + * internal::FEEvaluationImplCollocation::evaluate() or + * internal::FEEvaluationImplTransformToCollocation::evaluate() with + * appropriate template parameters. In case the template parameters + * fe_degree and n_q_points_1d contain valid information (i.e. fe_degree>-1 + * and n_q_points_1d>0), we simply pass these values to the respective + * template specializations. Otherwise, we perform a runtime matching of + * the runtime parameters to find the correct specialization. This matching + * currently supports $0\leq fe\_degree \leq 9$ and $degree+1\leq + * n\_q\_points\_1d\leq fe\_degree+2$. + */ + template + struct FEEvaluationImplEvaluateSelector + { + template + static void + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const internal::MatrixFreeFunctions::ShapeInfo &shape_info, + Number *values_dofs_actual, + Number *values_quad, + Number *gradients_quad, + Number *hessians_quad, + Number *scratch_data) + { + // We enable a transformation to collocation for derivatives if it gives + // correct results (first condition), if it is the most efficient choice + // in terms of operation counts (second condition) and if we were able to + // initialize the fields in shape_info.templates.h from the polynomials + // (third condition). + static constexpr bool use_collocation = + n_q_points_1d > fe_degree && n_q_points_1d <= 3 * fe_degree / 2 + 1 && + n_q_points_1d < 200; + + if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && + shape_info.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_collocation) + { + internal::FEEvaluationImplCollocation:: + evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } + // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see + // shape_info.h for more details + else if (fe_degree >= 0 && use_collocation && + shape_info.element_type <= + internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImplTransformToCollocation< + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } + else if (fe_degree >= 0 && + shape_info.element_type <= + internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::tensor_symmetric, + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0, + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::truncated_tensor) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::truncated_tensor, + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::tensor_general) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::tensor_general, + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); + } + else + AssertThrow(false, ExcNotImplemented()); + } + }; + + + + /** + * This class chooses an appropriate evaluation strategy based on the + * template parameters and the shape_info variable which contains runtime + * parameters for the strategy underlying FEEvaluation::integrate(), i.e. + * this calls internal::FEEvaluationImpl::integrate(), + * internal::FEEvaluationImplCollocation::integrate() or + * internal::FEEvaluationImplTransformToCollocation::integrate() with + * appropriate template parameters. In case the template parameters + * fe_degree and n_q_points_1d contain valid information (i.e. fe_degree>-1 + * and n_q_points_1d>0), we simply pass these values to the respective + * template specializations. Otherwise, we perform a runtime matching of + * the runtime parameters to find the correct specialization. This matching + * currently supports $0\leq fe\_degree \leq 9$ and $degree+1\leq + * n\_q\_points\_1d\leq fe\_degree+2$. + */ + template + struct FEEvaluationImplIntegrateSelector + { + template + static void + run(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const internal::MatrixFreeFunctions::ShapeInfo &shape_info, + Number * values_dofs_actual, + Number * values_quad, + Number * gradients_quad, + Number * scratch_data, + const bool sum_into_values_array) + { + // We enable a transformation to collocation for derivatives if it gives + // correct results (first condition), if it is the most efficient choice + // in terms of operation counts (second condition) and if we were able to + // initialize the fields in shape_info.templates.h from the polynomials + // (third condition). + constexpr bool use_collocation = n_q_points_1d > fe_degree && + n_q_points_1d <= 3 * fe_degree / 2 + 1 && + n_q_points_1d < 200; + + if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && + shape_info.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_collocation) + { + internal::FEEvaluationImplCollocation:: + integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } + // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see + // shape_info.h for more details + else if (fe_degree >= 0 && use_collocation && + shape_info.element_type <= + internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImplTransformToCollocation< + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } + else if (fe_degree >= 0 && + shape_info.element_type <= + internal::MatrixFreeFunctions::tensor_symmetric) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::tensor_symmetric, + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0, + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::truncated_tensor) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::truncated_tensor, + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } + else if (shape_info.element_type == + internal::MatrixFreeFunctions::tensor_general) + { + internal::FEEvaluationImpl< + internal::MatrixFreeFunctions::tensor_general, + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); + } + else + AssertThrow(false, ExcNotImplemented()); + } + }; + + + template - struct FEFaceEvaluationSelector + template + void + adjust_for_face_orientation(const unsigned int dim, + const unsigned int n_components, + const unsigned int face_orientation, + const Table<2, unsigned int> &orientation_map, + const bool integrate, + const bool values, + const bool gradients, + const unsigned int n_q_points, + Number * tmp_values, + Number * values_quad, + Number * gradients_quad) { + Assert(face_orientation, ExcInternalError()); + const unsigned int *orientation = &orientation_map[face_orientation][0]; + for (unsigned int c = 0; c < n_components; ++c) + { + if (values == true) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = values_quad[c * n_q_points + orientation[q]]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = values_quad[c * n_q_points + q]; + for (unsigned int q = 0; q < n_q_points; ++q) + values_quad[c * n_q_points + q] = tmp_values[q]; + } + if (gradients == true) + for (unsigned int d = 0; d < dim; ++d) + { + if (integrate) + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[q] = + gradients_quad[(c * dim + d) * n_q_points + orientation[q]]; + else + for (unsigned int q = 0; q < n_q_points; ++q) + tmp_values[orientation[q]] = + gradients_quad[(c * dim + d) * n_q_points + q]; + for (unsigned int q = 0; q < n_q_points; ++q) + gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q]; + } + } + } + + + + template + struct FEFaceEvaluationImplEvaluateSelector + { + template static void - evaluate(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - const VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map) + run(const unsigned int n_components, + const MatrixFreeFunctions::ShapeInfo &data, + const VectorizedArrayType * values_array, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool evaluate_values, + const bool evaluate_gradients, + const unsigned int face_no, + const unsigned int subface_index, + const unsigned int face_orientation, + const Table<2, unsigned int> &orientation_map) { constexpr unsigned int static_dofs_per_face = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : @@ -2043,7 +2381,8 @@ namespace internal subface_index); if (face_orientation) - adjust_for_face_orientation(n_components, + adjust_for_face_orientation(dim, + n_components, face_orientation, orientation_map, false, @@ -2054,23 +2393,31 @@ namespace internal values_quad, gradients_quad); } + }; + + + template + struct FEFaceEvaluationImplIntegrateSelector + { + template static void - integrate(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &data, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const unsigned int face_no, - const unsigned int subface_index, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map) + run(const unsigned int n_components, + const MatrixFreeFunctions::ShapeInfo &data, + VectorizedArrayType * values_array, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const bool integrate_values, + const bool integrate_gradients, + const unsigned int face_no, + const unsigned int subface_index, + const unsigned int face_orientation, + const Table<2, unsigned int> &orientation_map) { if (face_orientation) - adjust_for_face_orientation(n_components, + adjust_for_face_orientation(dim, + n_components, face_orientation, orientation_map, true, @@ -2138,34 +2485,744 @@ namespace internal integrate_gradients, face_no); } + }; + + + + template + static bool + fe_face_evaluation_process_and_io( + const unsigned int n_components, + 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 std::array cells, + const std::array face_nos, + const unsigned int subface_index, + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, + const std::array face_orientations, + 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) + { + const unsigned int cell = cells[0]; + + // In the case of integration, we do not need to reshuffle the + // data at the quadrature points to adjust for the face + // orientation if the shape functions are nodal at the cell + // boundaries (and we only requested the integration of the + // values) or Hermite shape functions are used. These cases are + // handled later when the values are written back into the + // glrobal vector. + if (integrate && + (face_orientations[0] > 0 && + (subface_index < GeometryInfo::max_children_per_cell || + !(((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 || + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided || + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_mixed_strides || + dof_info.index_storage_variants[dof_access_index][cell] == + MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous))))) + { + AssertDimension(face_orientations.size(), 1); + adjust_for_face_orientation(dim, + n_components, + face_orientations[0], + orientation_map, + true, + do_values, + do_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + } + + // 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. + VectorizedArrayType grad_weight = + (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 && + data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ? + data.data.front() + .shape_data_on_face[0][fe_degree + (integrate ? + (2 - (face_nos[0] % 2)) : + (1 + (face_nos[0] % 2)))] : + VectorizedArrayType(0.0 /*dummy*/); + + 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); + + // we allocate small amounts of data on the stack to signal the compiler + // that this temporary data is only needed for the calculations but the + // final results can be discarded and need not be written back to + // memory. For large sizes or when the dofs per face is not a + // compile-time constant, however, we want to go to the heap in the + // `scratch_data` variable to not risk a stack overflow. + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArrayType + temp_data[static_dofs_per_face < stack_array_size_threshold ? + 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; + + const unsigned int dummy = 0; + + // re-orientation + std::array orientation; + 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; + + 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) + { + 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_nos[v] % 2)) : + (1 + (face_nos[v] % 2)))][v]; + + index_array_hermite[v] = + &data.face_to_cell_index_hermite(face_nos[v], 0); + } + } - template + // face_to_cell_index_nodal + std::array index_array_nodal; + + 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)) + { + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + + 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_orientations[n_face_orientations == 1 ? 0 : v] == 0 || + subface_index < GeometryInfo::max_children_per_cell) ? + i : + orientation[v][i]; + }; + + // this variable keeps track of whether we are able to directly write + // the results into the result (function returns true) or not, requiring + // an additional call to another function + bool accesses_global_vector = true; + + for (unsigned int comp = 0; comp < n_components; ++comp) + { + if (integrate) + function_0(temp1, comp); + if ((do_gradients == false && + data.data.front().nodal_at_cell_boundaries == true) || + (data.element_type == + MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) + { + // 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()); + Number2_ *vector_ptr = + global_vector_ptr + + dof_info.dof_indices_contiguous[dof_access_index] + [cell * + VectorizedArrayType::size()] + + (dof_info + .component_dof_indices_offset[active_fe_index] + [first_selected_component] + + comp * static_dofs_per_component) * + VectorizedArrayType::size(); + + if (fe_degree > 1 && do_gradients == true) + { + 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]; + AssertIndexRange(ind1, + data.dofs_per_component_on_cell); + AssertIndexRange(ind2, + data.dofs_per_component_on_cell); + const unsigned int i_ = reorientate(0, i); + function_1a(temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + + ind1 * VectorizedArrayType::size(), + vector_ptr + + ind2 * VectorizedArrayType::size(), + grad_weight); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + else + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + if (n_face_orientations == 1) + { + const unsigned int i_ = reorientate(0, i); + const unsigned int ind = index_array_nodal[0][i]; + function_1b(temp1[i_], + vector_ptr + + ind * VectorizedArrayType::size()); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + } + + // 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) + { + 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()]; + Number2_ *vector_ptr = + global_vector_ptr + + (comp * static_dofs_per_component + + dof_info + .component_dof_indices_offset[active_fe_index] + [first_selected_component]) * + VectorizedArrayType::size(); + if (fe_degree > 1 && do_gradients == true) + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + 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(); + function_2a(temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1, + vector_ptr + ind2, + grad_weight, + indices, + indices); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + else + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + if (n_face_orientations == 1) + { + const unsigned int i_ = reorientate(0, i); + const unsigned int ind = + index_array_nodal[0][i] * + VectorizedArrayType::size(); + function_2b(temp1[i_], vector_ptr + ind, indices); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + } + } + + // 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] + + comp * static_dofs_per_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 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] + + 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] + + index_array_hermite[0 /*TODO*/][2 * i + 1] * + strides[v]; + function_2a(temp1[i_], + temp1[i_ + dofs_per_face], + global_vector_ptr, + global_vector_ptr, + grad_weight, + ind1, + ind2); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + else + { + if (integrate == false) + for (unsigned int i = 0; i < 2 * 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) + { + const unsigned int i_ = + reorientate(n_face_orientations == 1 ? 0 : v, + i); + function_3a( + temp1[i_][v], + temp1[i_ + dofs_per_face][v], + global_vector_ptr + [indices[v] + + index_array_hermite + [n_face_orientations == 1 ? 0 : v] + [2 * i] * + strides[v]], + global_vector_ptr + [indices[v] + + index_array_hermite + [n_face_orientations == 1 ? 0 : v] + [2 * i + 1] * + strides[v]], + grad_weight[n_face_orientations == 1 ? 0 : v]); + } + } + } + else + { + if (n_filled_lanes == VectorizedArrayType::size()) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + 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] + + index_array_nodal[0][i] * strides[v]; + const unsigned int i_ = reorientate(0, i); + function_2b(temp1[i_], global_vector_ptr, ind); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + else + { + if (integrate == false) + for (unsigned int i = 0; i < 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) + function_3b( + temp1[reorientate( + n_face_orientations == 1 ? 0 : v, i)][v], + global_vector_ptr + [indices[v] + + index_array_nodal + [n_face_orientations == 1 ? 0 : v][i] * + strides[v]]); + } + } + } + + // 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) + { + const unsigned int *indices = + &dof_info.dof_indices_contiguous[dof_access_index] + [cell * + VectorizedArrayType::size()]; + Number2_ *vector_ptr = + global_vector_ptr + comp * static_dofs_per_component + + dof_info + .component_dof_indices_offset[active_fe_index] + [first_selected_component]; + + 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 && + dof_info + .n_vectorization_lanes_filled[dof_access_index] + [cell] == + VectorizedArrayType::size()) + { + 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_2a(temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1, + vector_ptr + ind2, + 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) + function_3a(temp1[i_][v], + temp1[i_ + dofs_per_face][v], + vector_ptr[ind1 + indices[v]], + vector_ptr[ind2 + indices[v]], + grad_weight[v]); + + if (integrate == false) + for (unsigned int v = n_filled_lanes; + v < VectorizedArrayType::size(); + ++v) + { + temp1[i_][v] = 0.0; + temp1[i_ + 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_3a( + temp1[reorientate(v, i)][v], + temp1[reorientate(v, i) + dofs_per_face][v], + vector_ptr[index_array_hermite[v][2 * i] + + indices[v]], + vector_ptr[index_array_hermite[v][2 * i + 1] + + indices[v]], + grad_weight[v]); + } + } + } + else + { + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + 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_], vector_ptr + ind, 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]; + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + function_3b(temp1[i_][v], + vector_ptr[ind + indices[v]]); + + if (integrate == false) + for (unsigned int v = n_filled_lanes; + v < VectorizedArrayType::size(); + ++v) + temp1[i_][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)][v], + vector_ptr[index_array_nodal[v][i] + + dof_info.dof_indices_contiguous + [dof_access_index][cells[v]]]); + } + } + } + } + else + { + // case 5: default vector access + + // for the integrate_scatter path (integrate == true), we + // need to only prepare the data in this function for all + // components to later call distribute_local_to_global(); + // for the gather_evaluate path (integrate == false), we + // instead want to leave early because we need to get the + // vector data from somewhere else + function_5(temp1, comp); + if (integrate) + accesses_global_vector = false; + else + return false; + } + } + else + { + // case 5: default vector access + function_5(temp1, comp); + if (integrate) + accesses_global_vector = false; + else + return false; + } + + if (!integrate) + function_0(temp1, comp); + } + + if (!integrate && + (face_orientations[0] > 0 && + subface_index < GeometryInfo::max_children_per_cell)) + { + AssertDimension(face_orientations.size(), 1); + adjust_for_face_orientation(dim, + n_components, + face_orientations[0], + orientation_map, + false, + do_values, + do_gradients, + data.n_q_points_face, + scratch_data, + values_quad, + gradients_quad); + } + + return accesses_global_vector; + } + + + template + struct FEFaceEvaluationImplGatherEvaluateSelector + { + template static bool - gather_evaluate( - const unsigned int n_components, - const Number2 * src_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - 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_orientations, - const Table<2, unsigned int> & orientation_map) + run(const unsigned int n_components, + const Number2 * src_ptr, + const MatrixFreeFunctions::ShapeInfo &data, + const MatrixFreeFunctions::DoFInfo & dof_info, + VectorizedArrayType * values_quad, + VectorizedArrayType *gradients_quad, + VectorizedArrayType *scratch_data, + const bool evaluate_values, + const bool evaluate_gradients, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + 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_orientations, + const Table<2, unsigned int> & orientation_map) { if (src_ptr == nullptr) { return false; } - return process_and_io( // + return fe_face_evaluation_process_and_io( // n_components, false /*=evaluate*/, src_ptr, @@ -2272,28 +3329,34 @@ namespace internal subface_index); }); } + }; - template + template + struct FEFaceEvaluationImplIntegrateScatterSelector + { + template static bool - integrate_scatter( - const unsigned int n_components, - Number2 * dst_ptr, - const MatrixFreeFunctions::ShapeInfo &data, - const MatrixFreeFunctions::DoFInfo & dof_info, - VectorizedArrayType * values_array, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad, - VectorizedArrayType * scratch_data, - const bool integrate_values, - const bool integrate_gradients, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - 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_orientations, - const Table<2, unsigned int> & orientation_map) + run(const unsigned int n_components, + Number2 * dst_ptr, + const MatrixFreeFunctions::ShapeInfo &data, + const MatrixFreeFunctions::DoFInfo & dof_info, + VectorizedArrayType * values_array, + VectorizedArrayType * values_quad, + VectorizedArrayType *gradients_quad, + VectorizedArrayType *scratch_data, + const bool integrate_values, + const bool integrate_gradients, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + 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_orientations, + const Table<2, unsigned int> & orientation_map) { if (dst_ptr == nullptr) { @@ -2301,24 +3364,28 @@ namespace internal AssertDimension(face_orientations.size(), 1); // for block vectors simply integrate - integrate(n_components, - data, - values_array, - values_quad, - gradients_quad, - scratch_data, - integrate_values, - integrate_gradients, - face_nos[0], - subface_index, - face_orientations[0], - orientation_map); + FEFaceEvaluationImplIntegrateSelector:: + template run(n_components, + data, + values_array, + values_quad, + gradients_quad, + scratch_data, + integrate_values, + integrate_gradients, + face_nos[0], + subface_index, + face_orientations[0], + orientation_map); // default vector access return false; } - return process_and_io( // + return fe_face_evaluation_process_and_io( // n_components, true /*=integrate*/, dst_ptr, @@ -2442,743 +3509,7 @@ namespace internal subface_index); }); } - - private: - template - static bool - process_and_io( - const unsigned int n_components, - 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 std::array cells, - const std::array face_nos, - const unsigned int subface_index, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, - const std::array face_orientations, - 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) - { - const unsigned int cell = cells[0]; - - // In the case of integration, we do not need to reshuffle the - // data at the quadrature points to adjust for the face - // orientation if the shape functions are nodal at the cell - // boundaries (and we only requested the integration of the - // values) or Hermite shape functions are used. These cases are - // handled later when the values are written back into the - // glrobal vector. - if (integrate && - (face_orientations[0] > 0 && - (subface_index < GeometryInfo::max_children_per_cell || - !(((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 || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_strided || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - interleaved_contiguous_mixed_strides || - dof_info.index_storage_variants[dof_access_index][cell] == - MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: - contiguous))))) - { - AssertDimension(face_orientations.size(), 1); - adjust_for_face_orientation(n_components, - face_orientations[0], - orientation_map, - true, - do_values, - do_gradients, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad); - } - - // 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. - VectorizedArrayType grad_weight = - (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 && - data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ? - data.data.front() - .shape_data_on_face[0][fe_degree + (integrate ? - (2 - (face_nos[0] % 2)) : - (1 + (face_nos[0] % 2)))] : - VectorizedArrayType(0.0 /*dummy*/); - - 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); - - // we allocate small amounts of data on the stack to signal the compiler - // that this temporary data is only needed for the calculations but the - // final results can be discarded and need not be written back to - // memory. For large sizes or when the dofs per face is not a - // compile-time constant, however, we want to go to the heap in the - // `scratch_data` variable to not risk a stack overflow. - constexpr unsigned int stack_array_size_threshold = 100; - - VectorizedArrayType - temp_data[static_dofs_per_face < stack_array_size_threshold ? - 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; - - const unsigned int dummy = 0; - - // re-orientation - std::array orientation; - 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; - - 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) - { - 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_nos[v] % 2)) : - (1 + (face_nos[v] % 2)))][v]; - - index_array_hermite[v] = - &data.face_to_cell_index_hermite(face_nos[v], 0); - } - } - - // face_to_cell_index_nodal - std::array index_array_nodal; - - 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)) - { - for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) - { - if (cells[v] == numbers::invalid_unsigned_int) - continue; - - 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_orientations[n_face_orientations == 1 ? 0 : v] == 0 || - subface_index < GeometryInfo::max_children_per_cell) ? - i : - orientation[v][i]; - }; - - // this variable keeps track of whether we are able to directly write - // the results into the result (function returns true) or not, requiring - // an additional call to another function - bool accesses_global_vector = true; - - for (unsigned int comp = 0; comp < n_components; ++comp) - { - if (integrate) - function_0(temp1, comp); - if ((do_gradients == false && - data.data.front().nodal_at_cell_boundaries == true) || - (data.element_type == - MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1)) - { - // 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()); - Number2_ *vector_ptr = - global_vector_ptr + - dof_info.dof_indices_contiguous - [dof_access_index][cell * VectorizedArrayType::size()] + - (dof_info - .component_dof_indices_offset[active_fe_index] - [first_selected_component] + - comp * static_dofs_per_component) * - VectorizedArrayType::size(); - - if (fe_degree > 1 && do_gradients == true) - { - 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]; - AssertIndexRange(ind1, - data.dofs_per_component_on_cell); - AssertIndexRange(ind2, - data.dofs_per_component_on_cell); - const unsigned int i_ = reorientate(0, i); - function_1a(temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + - ind1 * VectorizedArrayType::size(), - vector_ptr + - ind2 * VectorizedArrayType::size(), - grad_weight); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = index_array_nodal[0][i]; - function_1b(temp1[i_], - vector_ptr + - ind * VectorizedArrayType::size()); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - } - } - - // 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) - { - 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()]; - Number2_ *vector_ptr = - global_vector_ptr + - (comp * static_dofs_per_component + - dof_info.component_dof_indices_offset - [active_fe_index][first_selected_component]) * - VectorizedArrayType::size(); - if (fe_degree > 1 && do_gradients == true) - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - 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(); - function_2a(temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1, - vector_ptr + ind2, - grad_weight, - indices, - indices); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - if (n_face_orientations == 1) - { - const unsigned int i_ = reorientate(0, i); - const unsigned int ind = - index_array_nodal[0][i] * - VectorizedArrayType::size(); - function_2b(temp1[i_], vector_ptr + ind, indices); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - } - } - - // 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] + - comp * static_dofs_per_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 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] + - 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] + - index_array_hermite[0 /*TODO*/][2 * i + 1] * - strides[v]; - function_2a(temp1[i_], - temp1[i_ + dofs_per_face], - global_vector_ptr, - global_vector_ptr, - grad_weight, - ind1, - ind2); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else - { - if (integrate == false) - for (unsigned int i = 0; i < 2 * 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) - { - const unsigned int i_ = - reorientate(n_face_orientations == 1 ? 0 : v, - i); - function_3a( - temp1[i_][v], - temp1[i_ + dofs_per_face][v], - global_vector_ptr - [indices[v] + - index_array_hermite - [n_face_orientations == 1 ? 0 : v] - [2 * i] * - strides[v]], - global_vector_ptr - [indices[v] + - index_array_hermite - [n_face_orientations == 1 ? 0 : v] - [2 * i + 1] * - strides[v]], - grad_weight[n_face_orientations == 1 ? 0 : - v]); - } - } - } - else - { - if (n_filled_lanes == VectorizedArrayType::size()) - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - 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] + - index_array_nodal[0][i] * strides[v]; - const unsigned int i_ = reorientate(0, i); - function_2b(temp1[i_], global_vector_ptr, ind); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else - { - if (integrate == false) - for (unsigned int i = 0; i < 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) - function_3b( - temp1[reorientate( - n_face_orientations == 1 ? 0 : v, i)][v], - global_vector_ptr - [indices[v] + - index_array_nodal - [n_face_orientations == 1 ? 0 : v][i] * - strides[v]]); - } - } - } - - // 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) - { - const unsigned int *indices = - &dof_info.dof_indices_contiguous - [dof_access_index][cell * VectorizedArrayType::size()]; - Number2_ *vector_ptr = - global_vector_ptr + comp * static_dofs_per_component + - dof_info - .component_dof_indices_offset[active_fe_index] - [first_selected_component]; - - 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 && - dof_info.n_vectorization_lanes_filled - [dof_access_index][cell] == - VectorizedArrayType::size()) - { - 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_2a(temp1[i_], - temp1[i_ + dofs_per_face], - vector_ptr + ind1, - vector_ptr + ind2, - 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) - function_3a(temp1[i_][v], - temp1[i_ + dofs_per_face][v], - vector_ptr[ind1 + indices[v]], - vector_ptr[ind2 + indices[v]], - grad_weight[v]); - - if (integrate == false) - for (unsigned int v = n_filled_lanes; - v < VectorizedArrayType::size(); - ++v) - { - temp1[i_][v] = 0.0; - temp1[i_ + 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_3a( - temp1[reorientate(v, i)][v], - temp1[reorientate(v, i) + dofs_per_face][v], - vector_ptr[index_array_hermite[v][2 * i] + - indices[v]], - vector_ptr[index_array_hermite[v][2 * i + 1] + - indices[v]], - grad_weight[v]); - } - } - } - else - { - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - 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_], vector_ptr + ind, 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]; - - for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3b(temp1[i_][v], - vector_ptr[ind + indices[v]]); - - if (integrate == false) - for (unsigned int v = n_filled_lanes; - v < VectorizedArrayType::size(); - ++v) - temp1[i_][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)][v], - vector_ptr[index_array_nodal[v][i] + - dof_info.dof_indices_contiguous - [dof_access_index][cells[v]]]); - } - } - } - } - else - { - // case 5: default vector access - - // for the integrate_scatter path (integrate == true), we - // need to only prepare the data in this function for all - // components to later call distribute_local_to_global(); - // for the gather_evaluate path (integrate == false), we - // instead want to leave early because we need to get the - // vector data from somewhere else - function_5(temp1, comp); - if (integrate) - accesses_global_vector = false; - else - return false; - } - } - else - { - // case 5: default vector access - function_5(temp1, comp); - if (integrate) - accesses_global_vector = false; - else - return false; - } - - if (!integrate) - function_0(temp1, comp); - } - - if (!integrate && - (face_orientations[0] > 0 && - subface_index < GeometryInfo::max_children_per_cell)) - { - AssertDimension(face_orientations.size(), 1); - adjust_for_face_orientation(n_components, - face_orientations[0], - orientation_map, - false, - do_values, - do_gradients, - data.n_q_points_face, - scratch_data, - values_quad, - gradients_quad); - } - - return accesses_global_vector; - } - - static void - adjust_for_face_orientation(const unsigned int n_components, - const unsigned int face_orientation, - const Table<2, unsigned int> &orientation_map, - const bool integrate, - const bool values, - const bool gradients, - const unsigned int n_q_points, - VectorizedArrayType * tmp_values, - VectorizedArrayType * values_quad, - VectorizedArrayType * gradients_quad) - { - Assert(face_orientation, ExcInternalError()); - const unsigned int *orientation = &orientation_map[face_orientation][0]; - for (unsigned int c = 0; c < n_components; ++c) - { - if (values == true) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = values_quad[c * n_q_points + orientation[q]]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = values_quad[c * n_q_points + q]; - for (unsigned int q = 0; q < n_q_points; ++q) - values_quad[c * n_q_points + q] = tmp_values[q]; - } - if (gradients == true) - for (unsigned int d = 0; d < dim; ++d) - { - if (integrate) - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points + - orientation[q]]; - else - for (unsigned int q = 0; q < n_q_points; ++q) - tmp_values[orientation[q]] = - gradients_quad[(c * dim + d) * n_q_points + q]; - for (unsigned int q = 0; q < n_q_points; ++q) - gradients_quad[(c * dim + d) * n_q_points + q] = - tmp_values[q]; - } - } - } - }; // namespace internal + }; diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h index 08089ceabf..4aac6af5d2 100644 --- a/include/deal.II/matrix_free/evaluation_selector.h +++ b/include/deal.II/matrix_free/evaluation_selector.h @@ -567,107 +567,15 @@ SelectEvaluator::evaluate( { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); - if (fe_degree + 1 == n_q_points_1d && - shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_collocation) - { - internal::FEEvaluationImplCollocation::evaluate( - n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); - } - // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see - // shape_info.h for more details - else if (use_collocation && shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - { - internal::FEEvaluationImplTransformToCollocation< - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); - } - else if (shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); - } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) - { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); - } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::truncated_tensor) - { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::truncated_tensor, - dim, - fe_degree, - n_q_points_1d, - Number>::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); - } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_general) - { - internal::FEEvaluationImpl::evaluate(n_components, - evaluation_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - hessians_quad, - scratch_data); - } - else - AssertThrow(false, ExcNotImplemented()); + internal::FEEvaluationImplEvaluateSelector:: + template run(n_components, + evaluation_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + hessians_quad, + scratch_data); } @@ -686,107 +594,15 @@ SelectEvaluator::integrate( { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); - if (fe_degree + 1 == n_q_points_1d && - shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_collocation) - { - internal::FEEvaluationImplCollocation::integrate( - n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - scratch_data, - sum_into_values_array); - } - // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see - // shape_info.h for more details - else if (use_collocation && shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - { - internal::FEEvaluationImplTransformToCollocation< - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - scratch_data, - sum_into_values_array); - } - else if (shape_info.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - scratch_data, - sum_into_values_array); - } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) - { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - scratch_data, - sum_into_values_array); - } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::truncated_tensor) - { - internal::FEEvaluationImpl< - internal::MatrixFreeFunctions::truncated_tensor, - dim, - fe_degree, - n_q_points_1d, - Number>::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - scratch_data, - sum_into_values_array); - } - else if (shape_info.element_type == - internal::MatrixFreeFunctions::tensor_general) - { - internal::FEEvaluationImpl::integrate(n_components, - integration_flag, - shape_info, - values_dofs_actual, - values_quad, - gradients_quad, - scratch_data, - sum_into_values_array); - } - else - AssertThrow(false, ExcNotImplemented()); + internal::FEEvaluationImplIntegrateSelector:: + template run(n_components, + integration_flag, + shape_info, + values_dofs_actual, + values_quad, + gradients_quad, + scratch_data, + sum_into_values_array); } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 38bb982967..e7fc1e075f 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -8105,25 +8105,20 @@ FEFaceEvaluation::evaluate(n_components, - *this->data, - values_array, - this->begin_values(), - this->begin_gradients(), - this->scratch_data, - evaluation_flag & EvaluationFlags::values, - evaluation_flag & EvaluationFlags::gradients, - this->face_no, - this->subface_index, - this->face_orientation, - this->mapping_data - ->descriptor[this->active_fe_index] - .face_orientations); + internal::FEFaceEvaluationImplEvaluateSelector:: + template run( + n_components, + *this->data, + values_array, + this->begin_values(), + this->begin_gradients(), + this->scratch_data, + evaluation_flag & EvaluationFlags::values, + evaluation_flag & EvaluationFlags::gradients, + this->face_no, + this->subface_index, + this->face_orientation, + this->mapping_data->descriptor[this->active_fe_index].face_orientations); # ifdef DEBUG if (evaluation_flag & EvaluationFlags::values) @@ -8236,12 +8231,8 @@ FEFaceEvaluation:: - integrate( + internal::FEFaceEvaluationImplIntegrateSelector:: + template run( n_components, *this->data, values_array, @@ -8318,55 +8309,55 @@ FEFaceEvaluationcompute_face_no_data(); const auto face_orientations = this->compute_face_orientations(); - return internal::FEFaceEvaluationSelector:: - gather_evaluate(n_components, - internal::get_beginning(input_vector), - *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); + return internal::FEFaceEvaluationImplGatherEvaluateSelector< + dim, + Number, + VectorizedArrayType>:: + template run( + n_components, + internal::get_beginning(input_vector), + *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(n_components, - internal::get_beginning(input_vector), - *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); + return internal::FEFaceEvaluationImplGatherEvaluateSelector< + dim, + Number, + VectorizedArrayType>:: + template run( + n_components, + internal::get_beginning(input_vector), + *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); } }; @@ -8436,30 +8427,30 @@ FEFaceEvaluationis_interior_face == false) == false, ExcNotImplemented()); - if (!internal::FEFaceEvaluationSelector:: - integrate_scatter(n_components, - internal::get_beginning(destination), - *this->data, - *this->dof_info, - this->begin_dof_values(), - 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 (!internal::FEFaceEvaluationImplIntegrateScatterSelector< + dim, + Number, + VectorizedArrayType>:: + template run( + n_components, + internal::get_beginning(destination), + *this->data, + *this->dof_info, + this->begin_dof_values(), + 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 we arrive here, writing into the destination vector did not succeed // because some of the assumptions in integrate_scatter were not diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index a73568fea2..99b0349a47 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -2126,21 +2126,21 @@ namespace internal // now let the matrix-free evaluators provide us with the // data on faces - FEFaceEvaluationSelector:: - evaluate(dim, - shape_info, - cell_points.data(), - face_quads.data(), - face_grads.data(), - scratch_data.data(), - true, - true, - face_no, - GeometryInfo::max_children_per_cell, - faces[face].face_orientation > 8 ? - faces[face].face_orientation - 8 : - 0, - my_data.descriptor[0].face_orientations); + FEFaceEvaluationImplEvaluateSelector:: + template run<-1, 0>(dim, + shape_info, + cell_points.data(), + face_quads.data(), + face_grads.data(), + scratch_data.data(), + true, + true, + face_no, + GeometryInfo::max_children_per_cell, + faces[face].face_orientation > 8 ? + faces[face].face_orientation - 8 : + 0, + my_data.descriptor[0].face_orientations); if (update_flags_faces & update_quadrature_points) @@ -2243,25 +2243,22 @@ namespace internal start_indices, cell_points.data()); - FEFaceEvaluationSelector:: - evaluate(dim, - shape_info, - cell_points.data(), - face_quads.data(), - face_grads.data(), - scratch_data.data(), - false, - true, - faces[face].exterior_face_no, - faces[face].subface_index, - faces[face].face_orientation < 8 ? - faces[face].face_orientation : - 0, - my_data.descriptor[0].face_orientations); + FEFaceEvaluationImplEvaluateSelector:: + template run<-1, 0>( + dim, + shape_info, + cell_points.data(), + face_quads.data(), + face_grads.data(), + scratch_data.data(), + false, + true, + faces[face].exterior_face_no, + faces[face].subface_index, + faces[face].face_orientation < 8 ? + faces[face].face_orientation : + 0, + my_data.descriptor[0].face_orientations); for (unsigned int q = 0; q < n_points_compute; ++q) { -- 2.39.5