From 3248660525d07b78410a8de86594f46480de18d8 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 22 Sep 2020 08:33:15 +0200 Subject: [PATCH] Remove n_face_orientations template argument within FEEvaluation Conflicts: include/deal.II/matrix_free/evaluation_kernels.h --- .../deal.II/matrix_free/evaluation_kernels.h | 942 +++++++++++------- .../matrix_free/evaluation_template_factory.h | 40 +- .../evaluation_template_factory.templates.h | 40 +- include/deal.II/matrix_free/fe_evaluation.h | 211 ++-- .../mg_transfer_global_coarsening.templates.h | 4 +- .../evaluation_template_factory.inst.in | 47 - 6 files changed, 736 insertions(+), 548 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 8b26c3f956..7877ae06b4 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -181,11 +181,10 @@ namespace internal (Eval::n_rows_of_product > Eval::n_columns_of_product ? Eval::n_rows_of_product : Eval::n_columns_of_product); - Number *temp1; + Number *temp1 = scratch_data; Number *temp2; if (temp_size == 0) { - temp1 = scratch_data; temp2 = temp1 + std::max(Utilities::fixed_power( shape_info.data.front().fe_degree + 1), Utilities::fixed_power( @@ -193,7 +192,6 @@ namespace internal } else { - temp1 = scratch_data; temp2 = temp1 + temp_size; } @@ -455,11 +453,10 @@ namespace internal (Eval::n_rows_of_product > Eval::n_columns_of_product ? Eval::n_rows_of_product : Eval::n_columns_of_product); - Number *temp1; + Number *temp1 = scratch_data; Number *temp2; if (temp_size == 0) { - temp1 = scratch_data; temp2 = temp1 + std::max(Utilities::fixed_power( shape_info.data.front().fe_degree + 1), Utilities::fixed_power( @@ -467,7 +464,6 @@ namespace internal } else { - temp1 = scratch_data; temp2 = temp1 + temp_size; } @@ -2490,50 +2486,35 @@ namespace internal - template + 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) + fe_face_evaluation_process_and_io(Processor &proc) { + auto n_components = proc.n_components; + auto integrate = proc.integrate; + auto global_vector_ptr = proc.global_vector_ptr; + auto &data = proc.data; + auto &dof_info = proc.dof_info; + auto values_quad = proc.values_quad; + auto gradients_quad = proc.gradients_quad; + auto scratch_data = proc.scratch_data; + auto do_values = proc.do_values; + auto do_gradients = proc.do_gradients; + auto active_fe_index = proc.active_fe_index; + auto first_selected_component = proc.first_selected_component; + auto cells = proc.cells; + auto face_nos = proc.face_nos; + auto subface_index = proc.subface_index; + auto dof_access_index = proc.dof_access_index; + auto face_orientations = proc.face_orientations; + auto &orientation_map = proc.orientation_map; + + static const int dim = Processor::dim_; + static const int fe_degree = Processor::fe_degree_; + using VectorizedArrayType = typename Processor::VectorizedArrayType_; + + using Number2_ = typename Processor::Number2_; + const unsigned int cell = cells[0]; // In the case of integration, we do not need to reshuffle the @@ -2600,28 +2581,14 @@ namespace internal 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 *DEAL_II_RESTRICT temp1; - VectorizedArrayType - temp_data[static_dofs_per_face < stack_array_size_threshold ? - 2 * dofs_per_face : - 1]; - if (static_dofs_per_face < stack_array_size_threshold) - temp1 = &temp_data[0]; - else - temp1 = scratch_data; + VectorizedArrayType *temp1 = scratch_data; const unsigned int dummy = 0; // re-orientation std::array orientation; + std::fill(orientation.begin(), orientation.end(), nullptr); + if (n_face_orientations == 1) orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ? &data.face_orientations[face_orientations[0]][0] : @@ -2645,6 +2612,7 @@ namespace internal // face_to_cell_index_hermite std::array index_array_hermite; + std::fill(index_array_hermite.begin(), index_array_hermite.end(), nullptr); if (n_face_orientations == 1) index_array_hermite[0] = @@ -2674,6 +2642,7 @@ namespace internal // face_to_cell_index_nodal std::array index_array_nodal; + std::fill(index_array_nodal.begin(), index_array_nodal.end(), nullptr); if (n_face_orientations == 1) index_array_nodal[0] = @@ -2710,7 +2679,7 @@ namespace internal for (unsigned int comp = 0; comp < n_components; ++comp) { if (integrate) - function_0(temp1, comp); + proc.function_0(temp1, comp); if ((do_gradients == false && data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == @@ -2754,13 +2723,12 @@ namespace internal 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); + proc.function_1a( + temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1 * VectorizedArrayType::size(), + vector_ptr + ind2 * VectorizedArrayType::size(), + grad_weight); } else { @@ -2776,9 +2744,10 @@ namespace internal { 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()); + proc.function_1b(temp1[i_], + vector_ptr + + ind * + VectorizedArrayType::size()); } else { @@ -2823,13 +2792,13 @@ namespace internal 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); + proc.function_2a(temp1[i_], + temp1[i_ + dofs_per_face], + vector_ptr + ind1, + vector_ptr + ind2, + grad_weight, + indices, + indices); } else { @@ -2847,7 +2816,9 @@ namespace internal const unsigned int ind = index_array_nodal[0][i] * VectorizedArrayType::size(); - function_2b(temp1[i_], vector_ptr + ind, indices); + proc.function_2b(temp1[i_], + vector_ptr + ind, + indices); } else { @@ -2908,13 +2879,13 @@ namespace internal 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); + proc.function_2a(temp1[i_], + temp1[i_ + dofs_per_face], + global_vector_ptr, + global_vector_ptr, + grad_weight, + ind1, + ind2); } else { @@ -2933,7 +2904,7 @@ namespace internal const unsigned int i_ = reorientate(n_face_orientations == 1 ? 0 : v, i); - function_3a( + proc.function_3a( temp1[i_][v], temp1[i_ + dofs_per_face][v], global_vector_ptr @@ -2967,7 +2938,9 @@ namespace internal 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); + proc.function_2b(temp1[i_], + global_vector_ptr, + ind); } else { @@ -2982,7 +2955,7 @@ namespace internal for (unsigned int v = 0; v < n_filled_lanes; ++v) for (unsigned int i = 0; i < dofs_per_face; ++i) - function_3b( + proc.function_3b( temp1[reorientate( n_face_orientations == 1 ? 0 : v, i)][v], global_vector_ptr @@ -3028,13 +3001,13 @@ namespace internal 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); + proc.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) { @@ -3050,11 +3023,11 @@ namespace internal [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]); + proc.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; @@ -3075,7 +3048,7 @@ namespace internal [cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3a( + proc.function_3a( temp1[reorientate(v, i)][v], temp1[reorientate(v, i) + dofs_per_face][v], vector_ptr[index_array_hermite[v][2 * i] + @@ -3099,7 +3072,9 @@ namespace internal const unsigned int ind = index_array_nodal[0][i]; const unsigned int i_ = reorientate(0, i); - function_2b(temp1[i_], vector_ptr + ind, indices); + proc.function_2b(temp1[i_], + vector_ptr + ind, + indices); } else if (n_face_orientations == 1) { @@ -3112,8 +3087,8 @@ namespace internal [cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) - function_3b(temp1[i_][v], - vector_ptr[ind + indices[v]]); + proc.function_3b(temp1[i_][v], + vector_ptr[ind + indices[v]]); if (integrate == false) for (unsigned int v = n_filled_lanes; @@ -3127,7 +3102,7 @@ namespace internal v < VectorizedArrayType::size(); ++v) if (cells[v] != numbers::invalid_unsigned_int) - function_3b( + proc.function_3b( temp1[reorientate(v, i)][v], vector_ptr[index_array_nodal[v][i] + dof_info.dof_indices_contiguous @@ -3139,6 +3114,7 @@ namespace internal else { // case 5: default vector access + AssertDimension(n_face_orientations, 1); // for the integrate_scatter path (integrate == true), we // need to only prepare the data in this function for all @@ -3146,7 +3122,7 @@ namespace internal // 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); + proc.function_5(temp1, comp); if (integrate) accesses_global_vector = false; else @@ -3156,7 +3132,9 @@ namespace internal else { // case 5: default vector access - function_5(temp1, comp); + AssertDimension(n_face_orientations, 1); + + proc.function_5(temp1, comp); if (integrate) accesses_global_vector = false; else @@ -3164,14 +3142,14 @@ namespace internal } if (!integrate) - function_0(temp1, comp); + proc.function_0(temp1, comp); } if (!integrate && (face_orientations[0] > 0 && subface_index < GeometryInfo::max_children_per_cell)) { - AssertDimension(face_orientations.size(), 1); + AssertDimension(n_face_orientations, 1); adjust_for_face_orientation(dim, n_components, face_orientations[0], @@ -3195,10 +3173,11 @@ namespace internal typename Number2 = Number> struct FEFaceEvaluationImplGatherEvaluateSelector { - template + template static bool - run(const unsigned int n_components, - const Number2 * src_ptr, + run(const unsigned int n_components, + const unsigned int n_face_orientations, + const Number2 * src_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -3208,130 +3187,234 @@ namespace internal 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) + 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 fe_face_evaluation_process_and_io( // - n_components, - false /*=evaluate*/, - src_ptr, - data, - dof_info, - values_quad, - gradients_quad, - scratch_data, - evaluate_values, - evaluate_gradients, - active_fe_index, - first_selected_component, - cells, - face_nos, - subface_index, - dof_access_index, - face_orientations, - orientation_map, - [](VectorizedArrayType & temp_1, - VectorizedArrayType & temp_2, - const Number * src_ptr_1, - const Number * src_ptr_2, - const VectorizedArrayType &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); - }, - [](VectorizedArrayType &temp, const Number *src_ptr) { - // case 1b) - do_vectorized_read(src_ptr, temp); - }, - [](VectorizedArrayType & temp_1, - VectorizedArrayType & temp_2, - const Number * src_ptr_1, - const Number * src_ptr_2, - const VectorizedArrayType &grad_weight, - const unsigned int * indices_1, - const unsigned int * 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); - }, - [](VectorizedArrayType &temp, - const Number * src_ptr, - const unsigned int * indices) { - // case 2b) - do_vectorized_gather(src_ptr, indices, temp); - }, - [](Number & temp_1, - Number & temp_2, - const Number src_ptr_1, - const Number src_ptr_2, - Number grad_weight) { - // case 3a) - temp_1 = src_ptr_1; - temp_2 = grad_weight * (temp_1 - src_ptr_2); - }, - [](Number &temp, const Number src_ptr) { - // case 3b) - temp = src_ptr; - }, - [&](const VectorizedArrayType *, const unsigned int) { - // case 5) - }, - [&](auto &temp1, const unsigned int comp) { - const unsigned int dofs_per_face = - fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); - const unsigned int n_q_points = - fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : - data.n_q_points_face; - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - data.element_type <= MatrixFreeFunctions::tensor_symmetric) - FEFaceEvaluationImpl:: - evaluate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + comp * dim * n_q_points, - scratch_data + 2 * dofs_per_face, - evaluate_values, - evaluate_gradients, - subface_index); - else - FEFaceEvaluationImpl:: - evaluate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + comp * dim * n_q_points, - scratch_data + 2 * dofs_per_face, - evaluate_values, - evaluate_gradients, - subface_index); - }); + Processor p(n_components, + false, + src_ptr, + data, + dof_info, + values_quad, + gradients_quad, + scratch_data, + evaluate_values, + evaluate_gradients, + active_fe_index, + first_selected_component, + cells, + face_nos, + subface_index, + dof_access_index, + face_orientations, + orientation_map); + + if (n_face_orientations == VectorizedArrayType::size()) + return fe_face_evaluation_process_and_io( + p); + else + return fe_face_evaluation_process_and_io<1>(p); } + + private: + template + struct Processor + { + static const int dim_ = dim; + static const int fe_degree_ = fe_degree; + static const int n_q_points_1d_ = n_q_points_1d; + using VectorizedArrayType_ = VectorizedArrayType; + using Number2_ = const Number2; + + Processor( + const unsigned int n_components, + const bool integrate, + const 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) + : n_components(n_components) + , integrate(integrate) + , global_vector_ptr(global_vector_ptr) + , data(data) + , dof_info(dof_info) + , values_quad(values_quad) + , gradients_quad(gradients_quad) + , scratch_data(scratch_data) + , do_values(do_values) + , do_gradients(do_gradients) + , active_fe_index(active_fe_index) + , first_selected_component(first_selected_component) + , cells(cells) + , face_nos(face_nos) + , subface_index(subface_index) + , dof_access_index(dof_access_index) + , face_orientations(face_orientations) + , orientation_map(orientation_map) + {} + + template + void + function_1a(T0 & temp_1, + T0 & temp_2, + const T1 src_ptr_1, + const T1 src_ptr_2, + const T2 &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); + } + + template + void + function_1b(T1 &temp, const T2 src_ptr) + { + // case 1b) + do_vectorized_read(src_ptr, temp); + } + + template + void + function_2a(T0 & temp_1, + T0 & temp_2, + const T1 src_ptr_1, + const T1 src_ptr_2, + const T2 &grad_weight, + const T3 &indices_1, + const T3 &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); + } + + template + void + function_2b(T0 &temp, const T1 src_ptr, const T2 &indices) + { + // case 2b) + do_vectorized_gather(src_ptr, indices, temp); + } + + template + void + function_3a(T0 & temp_1, + T0 & temp_2, + const T1 &src_ptr_1, + const T2 &src_ptr_2, + const T2 &grad_weight) + { + // case 3a) + temp_1 = src_ptr_1; + temp_2 = grad_weight * (temp_1 - src_ptr_2); + } + + template + void + function_3b(T1 &temp, const T2 &src_ptr) + { + // case 3b) + temp = src_ptr; + } + + template + void + function_5(const T1 &, const unsigned int) + { + // case 5) + } + + template + void + function_0(T1 &temp1, const unsigned int comp) + { + const unsigned int dofs_per_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(data.data.front().fe_degree + 1, dim - 1); + const unsigned int n_q_points = + fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : + data.n_q_points_face; + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= MatrixFreeFunctions::tensor_symmetric) + FEFaceEvaluationImpl:: + evaluate_in_face(/* n_components */ 1, + data, + temp1, + values_quad + comp * n_q_points, + gradients_quad + comp * dim * n_q_points, + scratch_data + 2 * dofs_per_face, + do_values, + do_gradients, + subface_index); + else + FEFaceEvaluationImpl:: + evaluate_in_face(/* n_components */ 1, + data, + temp1, + values_quad + comp * n_q_points, + gradients_quad + comp * dim * n_q_points, + scratch_data + 2 * dofs_per_face, + do_values, + do_gradients, + subface_index); + } + + const unsigned int n_components; + const bool integrate; + const 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; + }; }; template struct FEFaceEvaluationImplIntegrateScatterSelector { - template + template static bool - run(const unsigned int n_components, - Number2 * dst_ptr, + run(const unsigned int n_components, + const unsigned int n_face_orientations, + Number2 * dst_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, @@ -3354,17 +3438,18 @@ namespace internal 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) + 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) { - AssertDimension(face_nos.size(), 1); - AssertDimension(face_orientations.size(), 1); + AssertDimension(n_face_orientations, 1); + AssertDimension(n_face_orientations, 1); // for block vectors simply integrate FEFaceEvaluationImplIntegrateSelector:: @@ -3385,135 +3470,244 @@ namespace internal return false; } - return fe_face_evaluation_process_and_io( // - n_components, - true /*=integrate*/, - dst_ptr, - data, - dof_info, - values_quad, - gradients_quad, - scratch_data, - integrate_values, - integrate_gradients, - active_fe_index, - first_selected_component, - cells, - face_nos, - subface_index, - dof_access_index, - face_orientations, - orientation_map, - [](const VectorizedArrayType &temp_1, - const VectorizedArrayType &temp_2, - Number * dst_ptr_1, - Number * dst_ptr_2, - const VectorizedArrayType &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); - }, - [](VectorizedArrayType & temp_1, - VectorizedArrayType & temp_2, - Number * dst_ptr_1, - Number * dst_ptr_2, - const VectorizedArrayType &grad_weight, - const unsigned int * indices_1, - const unsigned int * 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 VectorizedArrayType &temp, - Number * dst_ptr, - const unsigned int * indices) { - // case 2b) - do_vectorized_scatter_add(temp, indices, dst_ptr); - }, - [](const Number temp_1, - const Number temp_2, - Number & dst_ptr_1, - Number & dst_ptr_2, - const Number 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 Number temp, Number &dst_ptr) { - // case 3b) - dst_ptr += temp; - }, - [&](const VectorizedArrayType *temp1, const unsigned int comp) { - // case 5: default vector access, must be handled separately, just do - // the face-normal interpolation - - AssertDimension(face_nos.size(), 1); - - FEFaceNormalEvaluationImpl:: - template interpolate( - /* n_components */ 1, - data, - temp1, - values_array + comp * data.dofs_per_component_on_cell, - integrate_gradients, - face_nos[0]); - }, - [&](auto &temp1, const unsigned int comp) { - const unsigned int dofs_per_face = - fe_degree > -1 ? - Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(data.data.front().fe_degree + 1, dim - 1); - const unsigned int n_q_points = - fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : - data.n_q_points_face; - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - data.element_type <= - internal::MatrixFreeFunctions::tensor_symmetric) - internal::FEFaceEvaluationImpl:: - integrate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + dim * comp * n_q_points, - scratch_data + 2 * dofs_per_face, - integrate_values, - integrate_gradients, - subface_index); - else - internal::FEFaceEvaluationImpl:: - integrate_in_face(/* n_components */ 1, - data, - temp1, - values_quad + comp * n_q_points, - gradients_quad + dim * comp * n_q_points, - scratch_data + 2 * dofs_per_face, - integrate_values, - integrate_gradients, - subface_index); - }); + + Processor p(values_array, + n_components, + true, + dst_ptr, + data, + dof_info, + values_quad, + gradients_quad, + scratch_data, + integrate_values, + integrate_gradients, + active_fe_index, + first_selected_component, + cells, + face_nos, + subface_index, + dof_access_index, + face_orientations, + orientation_map); + + if (n_face_orientations == VectorizedArrayType::size()) + return fe_face_evaluation_process_and_io( + p); + else + return fe_face_evaluation_process_and_io<1>(p); } + + private: + template + struct Processor + { + static const int dim_ = dim; + static const int fe_degree_ = fe_degree; + static const int n_q_points_1d_ = n_q_points_1d; + using VectorizedArrayType_ = VectorizedArrayType; + using Number2_ = Number2; + + + Processor( + VectorizedArrayType *values_array, + 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) + : values_array(values_array) + , n_components(n_components) + , integrate(integrate) + , global_vector_ptr(global_vector_ptr) + , data(data) + , dof_info(dof_info) + , values_quad(values_quad) + , gradients_quad(gradients_quad) + , scratch_data(scratch_data) + , do_values(do_values) + , do_gradients(do_gradients) + , active_fe_index(active_fe_index) + , first_selected_component(first_selected_component) + , cells(cells) + , face_nos(face_nos) + , subface_index(subface_index) + , dof_access_index(dof_access_index) + , face_orientations(face_orientations) + , orientation_map(orientation_map) + {} + + template + void + function_1a(const T0 &temp_1, + const T1 &temp_2, + T2 dst_ptr_1, + T3 dst_ptr_2, + const T4 &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); + } + + template + void + function_1b(const T0 &temp, T1 dst_ptr) + { + // case 1b) + do_vectorized_add(temp, dst_ptr); + } + + template + void + function_2a(const T0 &temp_1, + const T0 &temp_2, + T1 dst_ptr_1, + T1 dst_ptr_2, + const T2 &grad_weight, + const T3 &indices_1, + const T3 &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); + } + + template + void + function_2b(const T0 &temp, T1 dst_ptr, const T2 &indices) + { + // case 2b) + do_vectorized_scatter_add(temp, indices, dst_ptr); + } + + template + void + function_3a(const T0 &temp_1, + const T0 &temp_2, + T1 & dst_ptr_1, + T1 & dst_ptr_2, + const T2 &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; + } + + template + void + function_3b(const T0 &temp, T1 &dst_ptr) + { + // case 3b) + dst_ptr += temp; + } + + template + void + function_5(const T0 &temp1, const unsigned int comp) + { + // case 5: default vector access, must be handled separately, just do + // the face-normal interpolation + + FEFaceNormalEvaluationImpl:: + template interpolate( + /* n_components */ 1, + data, + temp1, + values_array + comp * data.dofs_per_component_on_cell, + do_gradients, + face_nos[0]); + } + + template + void + function_0(T0 &temp1, const unsigned int comp) + { + const unsigned int dofs_per_face = + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(data.data.front().fe_degree + 1, dim - 1); + const unsigned int n_q_points = + fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : + data.n_q_points_face; + if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + data.element_type <= + internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl:: + integrate_in_face(/* n_components */ 1, + data, + temp1, + values_quad + comp * n_q_points, + gradients_quad + dim * comp * n_q_points, + scratch_data + 2 * dofs_per_face, + do_values, + do_gradients, + subface_index); + else + internal::FEFaceEvaluationImpl:: + integrate_in_face(/* n_components */ 1, + data, + temp1, + values_quad + comp * n_q_points, + gradients_quad + dim * comp * n_q_points, + scratch_data + 2 * dofs_per_face, + do_values, + do_gradients, + subface_index); + } + + VectorizedArrayType *values_array; + + + 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; + }; }; diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index a7484b2f25..dbb4e57a42 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -97,11 +97,11 @@ namespace internal const unsigned int face_orientation, const Table<2, unsigned int> &orientation_map); - template static bool gather_evaluate( - const unsigned int n_components, - const Number * src_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + const Number * src_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -111,18 +111,19 @@ namespace internal 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); - - template + 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); + static bool integrate_scatter( - const unsigned int n_components, - Number * dst_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + Number * dst_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, @@ -133,12 +134,13 @@ namespace internal 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); + 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); }; diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index d730f89e1c..168cf1c82c 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -198,11 +198,11 @@ namespace internal template - template bool FEFaceEvaluationFactory::gather_evaluate( - const unsigned int n_components, - const Number * src_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + const Number * src_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -212,12 +212,13 @@ namespace internal 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) + 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) { return instantiation_helper_run< 1, @@ -227,6 +228,7 @@ namespace internal data.data[0].fe_degree, data.data[0].n_q_points_1d, n_components, + n_face_orientations, src_ptr, data, dof_info, @@ -248,11 +250,11 @@ namespace internal template - template bool FEFaceEvaluationFactory::integrate_scatter( - const unsigned int n_components, - Number * dst_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + Number * dst_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, @@ -263,12 +265,13 @@ namespace internal 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) + 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) { return instantiation_helper_run< 1, @@ -278,6 +281,7 @@ namespace internal data.data[0].fe_degree, data.data[0].n_q_points_1d, n_components, + n_face_orientations, dst_ptr, data, dof_info, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index c6e6e6027e..049c8c9cf9 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -8470,36 +8470,53 @@ FEFaceEvaluation:: - 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); + VectorizedArrayType>::template run( + n_components, + VectorizedArrayType::size(), + 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 { + // TODO: this copying should not be necessary once we have introduced + // an internal-data structure + std::array cells_; + std::array face_no_; + std::array face_orientation_; + + std::fill(cells_.begin(), cells_.end(), -1); + std::fill(face_no_.begin(), face_no_.end(), -1); + std::fill(face_orientation_.begin(), face_orientation_.end(), -1); + + cells_[0] = this->cell; + face_no_[0] = this->face_no; + face_orientation_[0] = this->face_orientation; + if (fe_degree > -1) - return internal::FEFaceEvaluationImplGatherEvaluateSelector< - dim, - Number, - VectorizedArrayType>:: - template run( + { + return internal::FEFaceEvaluationImplGatherEvaluateSelector< + dim, + Number, + VectorizedArrayType>::template run( n_components, + 1, internal::get_beginning(input_vector), *this->data, *this->dof_info, @@ -8510,35 +8527,37 @@ FEFaceEvaluationactive_fe_index, this->first_selected_component, - std::array{{this->cell}}, - std::array{{this->face_no}}, + cells_, + face_no_, this->subface_index, this->dof_access_index, - std::array{{this->face_orientation}}, + face_orientation_, this->mapping_data->descriptor[this->active_fe_index] .face_orientations); + } else return internal:: FEFaceEvaluationFactory:: - 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); + gather_evaluate(n_components, + 1, + 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_no_, + this->subface_index, + this->dof_access_index, + face_orientation_, + this->mapping_data + ->descriptor[this->active_fe_index] + .face_orientations); } }; @@ -8608,32 +8627,47 @@ FEFaceEvaluationis_interior_face == false) == false, ExcNotImplemented()); + // TODO: this copying should not be necessary once we have introduced + // an internal-data structure + std::array cells_; + std::array face_no_; + std::array face_orientation_; + + std::fill(cells_.begin(), cells_.end(), -1); + std::fill(face_no_.begin(), face_no_.end(), -1); + std::fill(face_orientation_.begin(), face_orientation_.end(), -1); + + cells_[0] = this->cell; + face_no_[0] = this->face_no; + face_orientation_[0] = this->face_orientation; + if (fe_degree > -1) { 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)) + VectorizedArrayType>::template run( + n_components, + 1, + 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, + cells_, + face_no_, + this->subface_index, + this->dof_access_index, + 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 @@ -8645,26 +8679,27 @@ FEFaceEvaluation:: - 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)) + integrate_scatter(n_components, + 1, + 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, + cells_, + face_no_, + this->subface_index, + this->dof_access_index, + face_orientation_, + this->mapping_data + ->descriptor[this->active_fe_index] + .face_orientations)) { this->distribute_local_to_global(destination); } diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index 3814c31af6..a3414772be 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -182,8 +182,8 @@ namespace internal get_mpi_comm(const MeshType &mesh) { const auto *tria_parallel = dynamic_cast< - const parallel::TriangulationBase *>( + const dealii::parallel::TriangulationBase *>( &(mesh.get_triangulation())); return tria_parallel != nullptr ? tria_parallel->get_communicator() : diff --git a/source/matrix_free/evaluation_template_factory.inst.in b/source/matrix_free/evaluation_template_factory.inst.in index a0019c9c57..468f6f3724 100644 --- a/source/matrix_free/evaluation_template_factory.inst.in +++ b/source/matrix_free/evaluation_template_factory.inst.in @@ -32,51 +32,4 @@ for (deal_II_dimension : DIMENSIONS; deal_II_dimension, deal_II_scalar_vectorized::value_type, deal_II_scalar_vectorized>; - - template bool dealii::internal::FEFaceEvaluationFactory< - deal_II_dimension, - deal_II_scalar_vectorized::value_type, - deal_II_scalar_vectorized>:: - gather_evaluate<1>( - const unsigned int, - const deal_II_scalar_vectorized::value_type *, - const MatrixFreeFunctions::ShapeInfo &, - const MatrixFreeFunctions::DoFInfo &, - deal_II_scalar_vectorized *, - deal_II_scalar_vectorized *, - deal_II_scalar_vectorized *, - const bool, - const bool, - const unsigned int, - const unsigned int, - const std::array, - const std::array, - const unsigned int, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex, - const std::array, - const Table<2, unsigned int> &); - - template bool dealii::internal::FEFaceEvaluationFactory< - deal_II_dimension, - deal_II_scalar_vectorized::value_type, - deal_II_scalar_vectorized>:: - integrate_scatter<1>( - const unsigned int, - deal_II_scalar_vectorized::value_type *, - const MatrixFreeFunctions::ShapeInfo &, - const MatrixFreeFunctions::DoFInfo &, - deal_II_scalar_vectorized *, - deal_II_scalar_vectorized *, - deal_II_scalar_vectorized *, - deal_II_scalar_vectorized *, - const bool, - const bool, - const unsigned int, - const unsigned int, - const std::array, - const std::array, - const unsigned int, - const MatrixFreeFunctions::DoFInfo::DoFAccessIndex, - const std::array, - const Table<2, unsigned int> &); } -- 2.39.5