From 6422d921d195a9034fa3a868b67d77e7b48ccc62 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 10 Dec 2021 11:46:12 +0100 Subject: [PATCH] Simplify extract of face information, preparing for use in evaluation_kernels --- examples/step-76/step-76.cc | 6 +- include/deal.II/fe/mapping_q_internal.h | 2 +- .../deal.II/matrix_free/evaluation_kernels.h | 16 +- include/deal.II/matrix_free/fe_evaluation.h | 627 +++++------------- .../matrix_free/fe_evaluation_base_data.h | 298 +++++++-- .../matrix_free/mapping_info.templates.h | 6 +- .../matrix_free/interpolate_quadrature_01.cc | 4 +- 7 files changed, 418 insertions(+), 541 deletions(-) diff --git a/examples/step-76/step-76.cc b/examples/step-76/step-76.cc index c7b2a8f987..dd437c1867 100644 --- a/examples/step-76/step-76.cc +++ b/examples/step-76/step-76.cc @@ -790,11 +790,10 @@ namespace Euler_DG VectorizedArrayType>:: template interpolate_quadrature( dim + 2, + EvaluationFlags::values, data.get_shape_info(), buffer.data(), phi_m.begin_values(), - false, - false, face); // Check if the face is an internal or a boundary face and @@ -893,11 +892,10 @@ namespace Euler_DG VectorizedArrayType>:: template interpolate_quadrature( dim + 2, + EvaluationFlags::values, data.get_shape_info(), phi_m.begin_values(), phi.begin_values(), - false, - false, face); } diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 3d40e5154d..1fed26230c 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -1136,7 +1136,7 @@ namespace internal temp_data.descriptor.resize(1); temp_data.descriptor[0].n_q_points = n_q_points; FEEvaluationBaseData eval( - {&data.shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0); + {&data.shape_info, nullptr, &temp_data, 0, 0}); // prepare arrays if (evaluation_flag != EvaluationFlags::nothing) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 121f5a83d5..576187616b 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2553,7 +2553,7 @@ namespace internal private: template - static void + static DEAL_II_ALWAYS_INLINE void interpolate_generic(const unsigned int n_components, const Number * input, Number * output, @@ -3137,6 +3137,20 @@ namespace internal // re-orientation std::array orientation = {}; +# ifdef DEBUG + // currently on structured meshes are supported -> face numbers and + // orientations have to be the same for all filled lanes + for (unsigned int v = 1; v < n_lanes; ++v) + { + if (this->all_face_numbers[v] != static_cast(-1)) + AssertDimension(this->all_face_numbers[0], + this->all_face_numbers[v]); + if (this->all_face_orientations[v] != static_cast(-1)) + AssertDimension(this->all_face_orientations[0], + this->all_face_orientations[v]); + } +# endif + if (n_face_orientations == 1) orientation[0] = &eval.get_orientation_map()(eval.get_face_orientation(), 0); diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 81f6aaa3d4..a3192566af 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -41,6 +41,8 @@ #include #include +#include + DEAL_II_NAMESPACE_OPEN @@ -615,26 +617,6 @@ public: AlignedVector> &array, const std::array & value) const; - /** - * Return the id of the cells this FEEvaluation or FEFaceEvaluation is - * associated with. - */ - std::array - get_cell_ids() const; - - /** - * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is - * associated with. - */ - std::array - get_cell_or_face_ids() const; - - /** - * Return the first selected component. - */ - unsigned int - get_first_selected_component() const; - /** * Return the underlying MatrixFree object. */ @@ -792,18 +774,6 @@ protected: */ const MatrixFree *matrix_info; - /** - * Stores the number of components in the finite element as detected in the - * MatrixFree storage class for comparison with the template argument. - */ - const unsigned int n_fe_components; - - /** - * For a FiniteElement with more than one base element, select at which - * component this data structure should start. - */ - const unsigned int first_selected_component; - /** * A temporary data structure necessary to read degrees of freedom when no * MatrixFree object was given at initialization. @@ -2802,20 +2772,6 @@ public: * points is inaccurate and this value must be used instead. */ const unsigned int n_q_points; - - -private: - /** - * Return face number of each face of the current face batch. - */ - std::array - compute_face_no_data(); - - /** - * Determine the orientation of each face of the current face batch. - */ - std::array - compute_face_orientations(); }; @@ -2935,13 +2891,12 @@ inline FEEvaluationBasedof_info->start_components.back()) - , first_selected_component(first_selected_component) { this->set_data_pointers(scratch_data_array, n_components_); Assert( @@ -3002,13 +2957,11 @@ inline FEEvaluationBase>( mapping, quadrature, - update_flags)) + update_flags), + n_components_, + first_selected_component) , scratch_data_array(new AlignedVector()) , matrix_info(nullptr) - , n_fe_components(n_components_) - // keep the number of the selected component within the current base element - // for reading dof values - , first_selected_component(first_selected_component) { const unsigned int base_element_number = fe.component_to_base_index(first_selected_component).first; @@ -3052,8 +3005,6 @@ inline FEEvaluationBase() : other.matrix_info->acquire_scratch_data()) , matrix_info(other.matrix_info) - , n_fe_components(other.n_fe_components) - , first_selected_component(other.first_selected_component) { if (other.matrix_info == nullptr) { @@ -3143,9 +3094,6 @@ operator=(const FEEvaluationBaseacquire_scratch_data(); } - AssertDimension(n_fe_components, other.n_fe_components); - AssertDimension(first_selected_component, other.first_selected_component); - this->set_data_pointers(scratch_data_array, n_components_); return *this; @@ -3198,86 +3146,6 @@ FEEvaluationBase:: } - -template -inline std::array -FEEvaluationBase:: - get_cell_ids() const -{ - Assert(this->matrix_info != nullptr, ExcNotInitialized()); - - const unsigned int n_lanes = VectorizedArrayType::size(); - std::array cells; - - // initialize array - for (unsigned int i = 0; i < n_lanes; ++i) - cells[i] = numbers::invalid_unsigned_int; - - if ((is_face == false) || - (is_face && - this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face)) - { - // cell or interior face face (element-centric loop) - for (unsigned int i = 0; i < n_lanes; ++i) - cells[i] = this->cell * n_lanes + i; - } - else if (is_face && - this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face == false) - { - // exterior face (element-centric loop): for this case, we need to - // look into the FaceInfo field that collects information from both - // sides of a face once for the global mesh, and pick the face id that - // is not the local one (cell_this). - for (unsigned int i = 0; i < n_lanes; ++i) - { - // compute actual (non vectorized) cell ID - const unsigned int cell_this = this->cell * n_lanes + i; - // compute face ID - unsigned int face_index = - this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, - this->face_no, - i); - - if (face_index == numbers::invalid_unsigned_int) - continue; // invalid face ID: no neighbor on boundary - - // get cell ID on both sides of face - auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes) - .cells_interior[face_index % n_lanes]; - auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes) - .cells_exterior[face_index % n_lanes]; - - // compare the IDs with the given cell ID - if (cell_m == cell_this) - cells[i] = cell_p; // neighbor has the other ID - else if (cell_p == cell_this) - cells[i] = cell_m; - } - } - else if (is_face) - { - // face-centric faces - const unsigned int *cells_ = - this->is_interior_face ? - &this->matrix_info->get_face_info(this->cell).cells_interior[0] : - &this->matrix_info->get_face_info(this->cell).cells_exterior[0]; - for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) - if (cells_[i] != numbers::invalid_unsigned_int) - cells[i] = cells_[i]; - } - - return cells; -} - - namespace internal { template -std::array -FEEvaluationBase:: - get_cell_or_face_ids() const -{ - const unsigned int v_len = VectorizedArrayType::size(); - std::array cells; - - // initialize array - for (unsigned int i = 0; i < v_len; ++i) - cells[i] = numbers::invalid_unsigned_int; - - if ((is_face && - this->dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && - this->is_interior_face == false) || - (!is_face && !this->is_interior_face)) - { - // cell-based face-loop: plus face - for (unsigned int i = 0; i < v_len; ++i) - { - // compute actual (non vectorized) cell ID - const unsigned int cell_this = this->cell * v_len + i; - // compute face ID - unsigned int fn = - this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, - this->face_no, - i); - - if (fn == numbers::invalid_unsigned_int) - continue; // invalid face ID: no neighbor on boundary - - // get cell ID on both sides of face - auto cell_m = this->matrix_info->get_face_info(fn / v_len) - .cells_interior[fn % v_len]; - auto cell_p = this->matrix_info->get_face_info(fn / v_len) - .cells_exterior[fn % v_len]; - - // compare the IDs with the given cell ID - if (cell_m == cell_this) - cells[i] = cell_p; // neighbor has the other ID - else if (cell_p == cell_this) - cells[i] = cell_m; - } - } - else - { - for (unsigned int i = 0; i < v_len; ++i) - cells[i] = this->cell * v_len + i; - } - - return cells; -} - - - template :: Assert(this->dof_info != nullptr, ExcNotInitialized()); Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized()); - if (n_fe_components == 1) + if (this->n_fe_components == 1) for (unsigned int comp = 0; comp < n_components; ++comp) { Assert(src[comp] != nullptr, @@ -3613,10 +3421,9 @@ FEEvaluationBase:: { for (unsigned int v = 0; v < n_vectorization_actual; ++v) if (this->dof_info->hanging_node_constraint_masks.size() > 0 && - this->dof_info - ->hanging_node_constraint_masks[(this->cell * n_lanes + v) * - n_fe_components + - first_selected_component] != + this->dof_info->hanging_node_constraint_masks + [(this->cell * n_lanes + v) * this->n_fe_components + + this->first_selected_component] != internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) has_hn_constraints = true; } @@ -3640,13 +3447,13 @@ FEEvaluationBase:: { const unsigned int *dof_indices = this->dof_info->dof_indices_interleaved.data() + - this->dof_info->row_starts[this->cell * n_fe_components * n_lanes] + this->dof_info->row_starts[this->cell * this->n_fe_components * n_lanes] .first + this->dof_info ->component_dof_indices_offset[this->active_fe_index] [this->first_selected_component] * n_lanes; - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int i = 0; i < dofs_per_component; ++i, dof_indices += n_lanes) for (unsigned int comp = 0; comp < n_components; ++comp) @@ -3670,8 +3477,9 @@ FEEvaluationBase:: // to the dof indices of the cells on all lanes unsigned int cells_copied[n_lanes]; const unsigned int *cells; - bool has_constraints = false; - const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1; + bool has_constraints = false; + const unsigned int n_components_read = + this->n_fe_components > 1 ? n_components : 1; if (is_face) { @@ -3691,8 +3499,8 @@ FEEvaluationBase:: Assert(cells[v] < this->dof_info->row_starts.size() - 1, ExcInternalError()); const std::pair *my_index_start = - &this->dof_info->row_starts[cells[v] * n_fe_components + - first_selected_component]; + &this->dof_info->row_starts[cells[v] * this->n_fe_components + + this->first_selected_component]; // check whether any of the SIMD lanes has constraints, i.e., the // constraint indicator which is the second entry of row_starts @@ -3709,23 +3517,22 @@ FEEvaluationBase:: } else { - AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components, + AssertIndexRange((this->cell + 1) * n_lanes * this->n_fe_components, this->dof_info->row_starts.size()); for (unsigned int v = 0; v < n_vectorization_actual; ++v) { const std::pair *my_index_start = &this->dof_info - ->row_starts[(this->cell * n_lanes + v) * n_fe_components + - first_selected_component]; + ->row_starts[(this->cell * n_lanes + v) * this->n_fe_components + + this->first_selected_component]; if (my_index_start[n_components_read].second != my_index_start[0].second) has_constraints = true; if (this->dof_info->hanging_node_constraint_masks.size() > 0 && - this->dof_info - ->hanging_node_constraint_masks[(this->cell * n_lanes + v) * - n_fe_components + - first_selected_component] != + this->dof_info->hanging_node_constraint_masks + [(this->cell * n_lanes + v) * this->n_fe_components + + this->first_selected_component] != internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) has_hn_constraints = true; @@ -3750,7 +3557,7 @@ FEEvaluationBase:: for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_component; ++i) operation.process_empty(values_dofs[comp][i]); - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) { for (unsigned int v = 0; v < n_vectorization_actual; ++v) for (unsigned int i = 0; i < dofs_per_component; ++i) @@ -3786,9 +3593,9 @@ FEEvaluationBase:: const unsigned int cell_index = is_face ? cells[v] : this->cell * n_lanes + v; const unsigned int cell_dof_index = - cell_index * n_fe_components + first_selected_component; + cell_index * this->n_fe_components + this->first_selected_component; const unsigned int n_components_read = - n_fe_components > 1 ? n_components : 1; + this->n_fe_components > 1 ? n_components : 1; unsigned int index_indicators = this->dof_info->row_starts[cell_dof_index].second; unsigned int next_index_indicators = @@ -3816,7 +3623,7 @@ FEEvaluationBase:: next_index_indicators = index_indicators; } - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) { unsigned int ind_local = 0; for (; index_indicators != next_index_indicators; ++index_indicators) @@ -3950,7 +3757,7 @@ FEEvaluationBase:: Assert(!local_dof_indices.empty(), ExcNotInitialized()); const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell; - unsigned int index = first_selected_component * dofs_per_component; + unsigned int index = this->first_selected_component * dofs_per_component; for (unsigned int comp = 0; comp < n_components; ++comp) { for (unsigned int i = 0; i < dofs_per_component; ++i, ++index) @@ -4021,10 +3828,11 @@ FEEvaluationBase:: { const unsigned int dof_index = dof_indices_cont[this->cell * VectorizedArrayType::size()] + - this->dof_info->component_dof_indices_offset[this->active_fe_index] - [first_selected_component] * + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * VectorizedArrayType::size(); - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int comp = 0; comp < n_components; ++comp) operation.process_dofs_vectorized(dofs_per_component, dof_index, @@ -4040,7 +3848,7 @@ FEEvaluationBase:: return; } - std::array cells = + const std::array &cells = this->get_cell_or_face_ids(); // More general case: Must go through the components one by one and apply @@ -4102,7 +3910,7 @@ FEEvaluationBase:: if (n_filled_lanes == VectorizedArrayType::size() && n_lanes == VectorizedArrayType::size() && !is_ecl) { - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) { for (unsigned int comp = 0; comp < n_components; ++comp) { @@ -4128,12 +3936,12 @@ FEEvaluationBase:: for (unsigned int comp = 0; comp < n_components; ++comp) { auto vector_ptrs = compute_vector_ptrs( - (n_components == 1 || n_fe_components == 1) ? comp : 0); + (n_components == 1 || this->n_fe_components == 1) ? comp : 0); for (unsigned int i = 0; i < dofs_per_component; ++i) operation.process_empty(values_dofs[comp][i]); - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) { for (unsigned int v = 0; v < n_filled_lanes; ++v) if (mask[v] == true) @@ -4179,7 +3987,7 @@ FEEvaluationBase:: internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: contiguous) { - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int comp = 0; comp < n_components; ++comp) operation.process_dofs_vectorized_transpose(dofs_per_component, dof_indices, @@ -4198,7 +4006,7 @@ FEEvaluationBase:: internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous_strided) { - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int i = 0; i < dofs_per_component; ++i) { for (unsigned int comp = 0; comp < n_components; ++comp) @@ -4229,7 +4037,7 @@ FEEvaluationBase:: const unsigned int *offsets = &this->dof_info->dof_indices_interleave_strides [ind][VectorizedArrayType::size() * this->cell]; - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int i = 0; i < dofs_per_component; ++i) { for (unsigned int comp = 0; comp < n_components; ++comp) @@ -4266,7 +4074,7 @@ FEEvaluationBase:: internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: contiguous) { - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) { for (unsigned int v = 0; v < n_filled_lanes; ++v) if (mask[v] == true) @@ -4293,7 +4101,7 @@ FEEvaluationBase:: [ind][VectorizedArrayType::size() * this->cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1); - if (n_components == 1 || n_fe_components == 1) + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int v = 0; v < n_filled_lanes; ++v) { if (mask[v] == true) @@ -4472,7 +4280,7 @@ FEEvaluationBase:: const unsigned int cell_index = is_face ? cells[v] : this->cell * n_lanes + v; const unsigned int cell_dof_index = - cell_index * n_fe_components + first_selected_component; + cell_index * this->n_fe_components + this->first_selected_component; const auto mask = this->dof_info->hanging_node_constraint_masks[cell_dof_index]; @@ -4662,18 +4470,6 @@ FEEvaluationBase:: /*------------------------------ access to data fields ----------------------*/ -template -inline unsigned int -FEEvaluationBase:: - get_first_selected_component() const -{ - return first_selected_component; -} - template jacobian = &this->mapping_data->jacobians[0][offsets]; this->J_value = &this->mapping_data->JxW_values[offsets]; + for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) + this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i; + # ifdef DEBUG this->dof_values_initialized = false; this->values_quad_initialized = false; @@ -7353,7 +7152,6 @@ namespace internal * Implementation for standard vectors (that have the begin() methods). */ template ::type * = nullptr> bool try_gather_evaluate_inplace( - T phi, - const VectorType & input_vector, - const unsigned int cell, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const internal::MatrixFreeFunctions::DoFInfo *dof_info, - const EvaluationFlags::EvaluationFlags evaluation_flag) + T & phi, + const VectorType & input_vector, + const EvaluationFlags::EvaluationFlags evaluation_flag) { + const unsigned int cell = phi.get_current_cell_index(); + const auto & dof_info = phi.get_dof_info(); + using VectorizedArrayType = + typename std::remove_reference::type; + // If the index storage is interleaved and contiguous and the vector storage // has the correct alignment, we can directly pass the pointer into the // vector to the evaluate() call, without reading the vector entries into a // separate data field. This saves some operations. if (std::is_same::value && - dof_info->index_storage_variants + dof_info.index_storage_variants [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous && reinterpret_cast( input_vector.begin() + - dof_info->dof_indices_contiguous + dof_info.dof_indices_contiguous [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] [cell * VectorizedArrayType::size()]) % sizeof(VectorizedArrayType) == @@ -7391,14 +7190,15 @@ namespace internal const VectorizedArrayType *vec_values = reinterpret_cast( input_vector.begin() + - dof_info->dof_indices_contiguous + dof_info.dof_indices_contiguous [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] [cell * VectorizedArrayType::size()] + - dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component] * + dof_info.component_dof_indices_offset + [phi.get_active_fe_index()] + [phi.get_first_selected_component()] * VectorizedArrayType::size()); - phi->evaluate(vec_values, evaluation_flag); + phi.evaluate(vec_values, evaluation_flag); return true; } @@ -7410,7 +7210,6 @@ namespace internal * Implementation for block vectors. */ template ::value, VectorType>::type * = nullptr> bool - try_gather_evaluate_inplace(T, + try_gather_evaluate_inplace(T &, const VectorType &, - const unsigned int, - const unsigned int, - const unsigned int, - const internal::MatrixFreeFunctions::DoFInfo *, const EvaluationFlags::EvaluationFlags) { return false; @@ -7434,7 +7229,6 @@ namespace internal * Implementation for vectors that have the begin() methods. */ template ::type * = nullptr> bool try_integrate_scatter_inplace( - T phi, - VectorType & destination, - const unsigned int cell, - const unsigned int active_fe_index, - const unsigned int first_selected_component, - const internal::MatrixFreeFunctions::DoFInfo *dof_info, - const EvaluationFlags::EvaluationFlags evaluation_flag) + T & phi, + VectorType & destination, + const EvaluationFlags::EvaluationFlags evaluation_flag) { + const unsigned int cell = phi.get_current_cell_index(); + const auto & dof_info = phi.get_dof_info(); + using VectorizedArrayType = + typename std::remove_reference::type; + // If the index storage is interleaved and contiguous and the vector storage // has the correct alignment, we can directly pass the pointer into the // vector to the evaluate() call, without reading the vector entries into a // separate data field. This saves some operations. if (std::is_same::value && - dof_info->index_storage_variants + dof_info.index_storage_variants [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous && reinterpret_cast( destination.begin() + - dof_info->dof_indices_contiguous + dof_info.dof_indices_contiguous [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] [cell * VectorizedArrayType::size()]) % sizeof(VectorizedArrayType) == @@ -7475,11 +7270,12 @@ namespace internal dof_info->dof_indices_contiguous [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] [cell * VectorizedArrayType::size()] + - dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component] * + dof_info->component_dof_indices_offset + [phi.get_active_fe_index()] + [phi.get_first_selected_component()] * VectorizedArrayType::size()); - phi->integrate(evaluation_flag, vec_values, true); + phi.integrate(evaluation_flag, vec_values, true); return true; } @@ -7491,7 +7287,6 @@ namespace internal * Implementation for all other vectors like block vectors. */ template ( - this, - input_vector, - this->cell, - this->active_fe_index, - this->first_selected_component, - this->dof_info, - evaluation_flag) == false) + if (internal::try_gather_evaluate_inplace(*this, + input_vector, + evaluation_flag) == false) { this->read_dof_values(input_vector); evaluate(this->begin_dof_values(), evaluation_flag); @@ -7743,14 +7529,8 @@ FEEvaluation( - this, - destination, - this->cell, - this->active_fe_index, - this->first_selected_component, - this->dof_info, - integration_flag) == false) + if (internal::try_integrate_scatter_inplace( + *this, destination, integration_flag) == false) { integrate(integration_flag, this->begin_dof_values()); this->distribute_local_to_global(destination); @@ -7871,6 +7651,8 @@ FEFaceEvaluationreinit_face(this->matrix_info->get_face_info(face_index)); + for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) + this->cell_or_face_ids[i] = face_index * VectorizedArrayType::size() + i; this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index]; const unsigned int offsets = @@ -7925,12 +7707,80 @@ FEFaceEvaluationcell_type = this->matrix_info->get_mapping_info().cell_type[cell_index]; this->cell = cell_index; - this->face_orientation = 0; - this->subface_index = GeometryInfo::max_children_per_cell; - this->face_no = face_number; + this->subface_index = GeometryInfo::max_children_per_cell; this->dof_access_index = internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; + constexpr unsigned int n_lanes = VectorizedArrayType::size(); + + if (this->is_interior_face == false) + { + // for this case, we need to look into the FaceInfo field that collects + // information from both sides of a face once for the global mesh, and + // pick the face id that is not the local one (cell_this). + for (unsigned int i = 0; i < n_lanes; ++i) + { + // compute actual (non vectorized) cell ID + const unsigned int cell_this = cell_index * n_lanes + i; + // compute face ID + unsigned int face_index = + this->matrix_info->get_cell_and_face_to_plain_faces()(cell_index, + face_number, + i); + + if (face_index == numbers::invalid_unsigned_int) + { + this->cell_ids[i] = numbers::invalid_unsigned_int; + this->all_face_numbers[i] = numbers::invalid_unsigned_int; + this->all_face_orientations[i] = numbers::invalid_unsigned_int; + continue; // invalid face ID: no neighbor on boundary + } + + const auto &faces = + this->matrix_info->get_face_info(face_index / n_lanes); + // get cell ID on both sides of face + auto cell_m = faces.cells_interior[face_index % n_lanes]; + auto cell_p = faces.cells_exterior[face_index % n_lanes]; + + const bool is_interior_face = cell_m != cell_this; + + Assert(cell_m == cell_this || cell_p == cell_this, + ExcInternalError()); + + // compare the IDs with the given cell ID + if (is_interior_face) + { + this->cell_ids[i] = cell_m; // neighbor has the other ID + this->all_face_numbers[i] = faces.interior_face_no; + } + else + { + this->cell_ids[i] = cell_p; + this->all_face_numbers[i] = faces.exterior_face_no; + } + + const bool orientation_interior_face = faces.face_orientation >= 8; + unsigned int face_orientation = faces.face_orientation % 8; + if (is_interior_face != orientation_interior_face) + { + constexpr std::array table{ + {0, 1, 0, 3, 6, 5, 4, 7}}; + face_orientation = table[face_orientation]; + } + this->all_face_orientations[i] = face_orientation; + } + + this->face_no = this->all_face_numbers[0]; + this->face_orientation = this->all_face_orientations[0]; + } + else + { + this->face_orientation = 0; + this->face_no = face_number; + for (unsigned int i = 0; i < n_lanes; ++i) + this->cell_ids[i] = cell_index * n_lanes + i; + } + const unsigned int offsets = this->matrix_info->get_mapping_info() .face_data_by_cells[this->quad_no] @@ -8470,173 +8320,6 @@ FEFaceEvaluation -std::array -FEFaceEvaluation::compute_face_no_data() -{ - std::array face_no_data; - - if (this->dof_access_index != - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell || - this->is_interior_face == true) - { - std::fill(face_no_data.begin(), - face_no_data.begin() + - this->dof_info->n_vectorization_lanes_filled - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell], - this->face_no); - } - else - { - std::fill(face_no_data.begin(), - face_no_data.end(), - numbers::invalid_unsigned_int); - - for (unsigned int i = 0; - i < this->dof_info->n_vectorization_lanes_filled - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell]; - ++i) - { - // compute actual (non vectorized) cell ID - const unsigned int cell_this = - this->cell * VectorizedArrayType::size() + i; - // compute face ID - const unsigned int face_index = - this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, - this->face_no, - i); - - Assert(face_index != numbers::invalid_unsigned_int, - ExcNotInitialized()); - - // get cell ID on both sides of face - auto cell_m = - this->matrix_info - ->get_face_info(face_index / VectorizedArrayType::size()) - .cells_interior[face_index % VectorizedArrayType::size()]; - - // compare the IDs with the given cell ID - face_no_data[i] = - (cell_m == cell_this) ? - this->matrix_info - ->get_face_info(face_index / VectorizedArrayType::size()) - .exterior_face_no : - this->matrix_info - ->get_face_info(face_index / VectorizedArrayType::size()) - .interior_face_no; - } - } - - return face_no_data; -} - - - -template -std::array -FEFaceEvaluation::compute_face_orientations() -{ - std::array face_no_data; - - if (this->dof_access_index != - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell || - this->is_interior_face == true) - { - Assert(false, ExcNotImplemented()); - } - else - { - std::fill(face_no_data.begin(), - face_no_data.end(), - numbers::invalid_unsigned_int); - - if (dim == 3) - { - for (unsigned int i = 0; - i < this->dof_info->n_vectorization_lanes_filled - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell]; - ++i) - { - // compute actual (non vectorized) cell ID - const unsigned int cell_this = - this->cell * VectorizedArrayType::size() + i; - // compute face ID - const unsigned int face_index = - this->matrix_info->get_cell_and_face_to_plain_faces()( - this->cell, this->face_no, i); - - Assert(face_index != numbers::invalid_unsigned_int, - ExcNotInitialized()); - - const unsigned int macro = - face_index / VectorizedArrayType::size(); - const unsigned int lane = - face_index % VectorizedArrayType::size(); - - const auto &faces = this->matrix_info->get_face_info(macro); - - // get cell ID on both sides of face - auto cell_m = faces.cells_interior[lane]; - - const bool is_interior_face = cell_m != cell_this; - const bool fo_interior_face = faces.face_orientation >= 8; - - unsigned int face_orientation = faces.face_orientation % 8; - - if (is_interior_face != fo_interior_face) - { - // invert (see also: - // Triangulation::update_periodic_face_map()) - static const std::array table{ - {0, 1, 2, 3, 6, 5, 4, 7}}; - - face_orientation = table[face_orientation]; - } - - // compare the IDs with the given cell ID - face_no_data[i] = face_orientation; - } - } - else - { - std::fill( - face_no_data.begin(), - face_no_data.begin() + - this->dof_info->n_vectorization_lanes_filled - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell], - 0); - } - } - - return face_no_data; -} - - - template &info, - const unsigned int quad_no, - const bool is_interior_face, - const unsigned int face_type); + const bool is_interior_face = true, + const unsigned int quad_no = 0, + const unsigned int face_type = 0, + const unsigned int first_selected_component = 0); /** * Copy constructor. @@ -404,6 +407,12 @@ public: unsigned int get_active_quadrature_index() const; + /** + * Return the first selected component in a multi-component system. + */ + unsigned int + get_first_selected_component() const; + /** * If is_face is true, this function returns the face number within a cell for * the face this object was initialized to. On cells where the face makes no @@ -471,6 +480,46 @@ public: bool get_is_interior_face() const; + /** + * Return the id of the cells this FEEvaluation or FEFaceEvaluation is + * associated with. + */ + const std::array & + get_cell_ids() const; + + /** + * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is + * associated with. + */ + const std::array & + get_cell_or_face_ids() const; + + /** + * Return the (non-vectorized) number of faces within cells in case of ECL + * for the exterior cells where the single number accessible via + * get_face_no() is not enough to cover all cases. + * + * @note Only available for `dof_access_index == dof_access_cell` and + * `is_interior_face == false`. + * + * @note This function depends on the internal representation of data, which + * is not stable between releases of deal.II, and is hence mostly for + * internal use. + */ + const std::array & + get_all_face_numbers() const; + + /** + * Store the orientation of the neighbor's faces with respect to the current + * cell for the case of exterior faces on ECL with possibly different + * orientations behind different cells. + * + * @note Only available for `dof_access_index == dof_access_cell` and + * `is_interior_face == false`. + */ + const std::array & + get_all_face_orientations() const; + //@} protected: @@ -481,7 +530,9 @@ protected: FEEvaluationBaseData( const std::shared_ptr< internal::MatrixFreeFunctions:: - MappingDataOnTheFly> &mapping_data); + MappingDataOnTheFly> &mapping_data, + const unsigned int n_fe_components, + const unsigned int first_selected_component); /** * A pointer to the unit cell shape data, i.e., values, gradients and @@ -513,6 +564,18 @@ protected: */ const unsigned int quad_no; + /** + * Stores the number of components in the finite element as detected in the + * MatrixFree storage class for comparison with the template argument. + */ + const unsigned int n_fe_components; + + /** + * For a FiniteElement with more than one base element, select at which + * component this data structure should start. + */ + const unsigned int first_selected_component; + /** * The active FE index for this class for efficient indexing in the hp-case. */ @@ -731,6 +794,38 @@ protected: */ internal::MatrixFreeFunctions::GeometryType cell_type; + /** + * Stores the (non-vectorized) id of the cells or faces this object is + * initialized to. Relevant for ECL. + */ + std::array cell_ids; + + /** + * Stores the (non-vectorized) id of the cells or faces this object is + * initialized to. Relevant for ECL. + */ + std::array cell_or_face_ids; + + /** + * Stores the (non-vectorized) number of faces within cells in case of ECL + * for the exterior cells where the single number `face_no` is not enough to + * cover all cases. + * + * @note Only available for `dof_access_index == dof_access_cell` and + * `is_interior_face == false`. + */ + std::array all_face_numbers; + + /** + * Store the orientation of the neighbor's faces with respect to the current + * cell for the case of exterior faces on ECL with possibly different + * orientations behind different cells. + * + * @note Only available for `dof_access_index == dof_access_cell` and + * `is_interior_face == false`. + */ + std::array all_face_orientations; + /** * Geometry data that can be generated FEValues on the fly with the * respective constructor, as an alternative to the entry point with @@ -762,13 +857,16 @@ inline FEEvaluationBaseData:: const MappingInfoStorageType *, unsigned int, unsigned int> &shape_dof_info, - const unsigned int quad_no, const bool is_interior_face, - const unsigned int face_type) + const unsigned int quad_no, + const unsigned int face_type, + const unsigned int first_selected_component) : data(std::get<0>(shape_dof_info)) , dof_info(std::get<1>(shape_dof_info)) , mapping_data(std::get<2>(shape_dof_info)) , quad_no(quad_no) + , n_fe_components(dof_info != nullptr ? dof_info->start_components.back() : 0) + , first_selected_component(first_selected_component) , active_fe_index(std::get<3>(shape_dof_info)) , active_quad_index(std::get<4>(shape_dof_info)) , descriptor( @@ -810,11 +908,15 @@ inline FEEvaluationBaseData:: FEEvaluationBaseData( const std::shared_ptr< internal::MatrixFreeFunctions:: - MappingDataOnTheFly> &mapped_geometry) + MappingDataOnTheFly> &mapped_geometry, + const unsigned int n_fe_components, + const unsigned int first_selected_component) : data(nullptr) , dof_info(nullptr) , mapping_data(nullptr) , quad_no(numbers::invalid_unsigned_int) + , n_fe_components(n_fe_components) + , first_selected_component(first_selected_component) , active_fe_index(numbers::invalid_unsigned_int) , active_quad_index(numbers::invalid_unsigned_int) , descriptor(nullptr) @@ -844,6 +946,8 @@ FEEvaluationBaseData::operator=( const FEEvaluationBaseData &other) { AssertDimension(quad_no, other.quad_no); + AssertDimension(n_fe_components, other.n_fe_components); + AssertDimension(first_selected_component, other.first_selected_component); AssertDimension(active_fe_index, other.active_fe_index); AssertDimension(active_quad_index, other.active_quad_index); AssertDimension(n_quadrature_points, descriptor->n_q_points); @@ -883,6 +987,78 @@ FEEvaluationBaseData::operator=( +template +inline void +FEEvaluationBaseData:: + set_data_pointers(AlignedVector *scratch_data_array, + const unsigned int n_components) +{ + Assert(scratch_data_array != nullptr, ExcInternalError()); + + const unsigned int tensor_dofs_per_component = + Utilities::fixed_power(data->data.front().fe_degree + 1); + const unsigned int dofs_per_component = data->dofs_per_component_on_cell; + + const unsigned int size_scratch_data = + std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components * + 3 + + 2 * n_quadrature_points; + const unsigned int size_data_arrays = + n_components * dofs_per_component + + (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 1) * + n_quadrature_points); + + const unsigned int allocated_size = size_scratch_data + size_data_arrays; + scratch_data_array->resize_fast(allocated_size); + scratch_data.reinit(scratch_data_array->begin() + size_data_arrays, + size_scratch_data); + + // set the pointers to the correct position in the data array + values_dofs = scratch_data_array->begin(); + values_quad = scratch_data_array->begin() + n_components * dofs_per_component; + gradients_quad = scratch_data_array->begin() + + n_components * (dofs_per_component + n_quadrature_points); + gradients_from_hessians_quad = + scratch_data_array->begin() + + n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); + hessians_quad = + scratch_data_array->begin() + + n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points); +} + + + +template +inline void +FEEvaluationBaseData::reinit_face( + const internal::MatrixFreeFunctions::FaceToCellTopology< + VectorizedArrayType::size()> &face) +{ + Assert(is_face == true, + ExcMessage("Faces can only be set if the is_face template parameter " + "is true")); + face_no = (is_interior_face ? face.interior_face_no : face.exterior_face_no); + subface_index = is_interior_face == true ? + GeometryInfo::max_children_per_cell : + face.subface_index; + + // First check if interior or exterior cell has non-standard orientation + // (i.e. the third bit is one or not). Then set zero if this cell has + // standard-orientation else copy the first three bits + // (which is equivalent to modulo 8). See also the documentation of + // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation. + face_orientation = (is_interior_face == (face.face_orientation >= 8)) ? + (face.face_orientation % 8) : + 0; + + if (is_interior_face) + cell_ids = face.cells_interior; + else + cell_ids = face.cells_exterior; +} + + + template inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> FEEvaluationBaseData:: @@ -1147,6 +1323,16 @@ FEEvaluationBaseData:: +template +inline unsigned int +FEEvaluationBaseData:: + get_first_selected_component() const +{ + return first_selected_component; +} + + + template inline ArrayView FEEvaluationBaseData:: @@ -1219,65 +1405,63 @@ FEEvaluationBaseData:: template -inline void -FEEvaluationBaseData:: - set_data_pointers(AlignedVector *scratch_data_array, - const unsigned int n_components) +inline const std::array & +FEEvaluationBaseData::get_cell_ids() + const { - Assert(scratch_data_array != nullptr, ExcInternalError()); - - const unsigned int tensor_dofs_per_component = - Utilities::fixed_power(data->data.front().fe_degree + 1); - const unsigned int dofs_per_component = data->dofs_per_component_on_cell; + Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + return cell_ids; +} - const unsigned int size_scratch_data = - std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components * - 3 + - 2 * n_quadrature_points; - const unsigned int size_data_arrays = - n_components * dofs_per_component + - (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 1) * - n_quadrature_points); - const unsigned int allocated_size = size_scratch_data + size_data_arrays; - scratch_data_array->resize_fast(allocated_size); - scratch_data.reinit(scratch_data_array->begin() + size_data_arrays, - size_scratch_data); - // set the pointers to the correct position in the data array - values_dofs = scratch_data_array->begin(); - values_quad = scratch_data_array->begin() + n_components * dofs_per_component; - gradients_quad = scratch_data_array->begin() + - n_components * (dofs_per_component + n_quadrature_points); - gradients_from_hessians_quad = - scratch_data_array->begin() + - n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); - hessians_quad = - scratch_data_array->begin() + - n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points); +template +inline const std::array & +FEEvaluationBaseData:: + get_cell_or_face_ids() const +{ + Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + if (!is_face || dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) + return cell_ids; + else + return cell_or_face_ids; } template -inline void -FEEvaluationBaseData::reinit_face( - const internal::MatrixFreeFunctions::FaceToCellTopology< - VectorizedArrayType::size()> &face) +inline const std::array & +FEEvaluationBaseData:: + get_all_face_numbers() const { - face_no = (is_interior_face ? face.interior_face_no : face.exterior_face_no); - subface_index = is_interior_face == true ? - GeometryInfo::max_children_per_cell : - face.subface_index; + Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + Assert(is_face && + dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + is_interior_face == false, + ExcMessage("All face numbers can only be queried for ECL at exterior " + "faces. Use get_face_no() in other cases.")); + + return all_face_numbers; +} - // First check if interior or exterior cell has non-standard orientation - // (i.e. the third bit is one or not). Then set zero if this cell has - // standard-orientation else copy the first three bits - // (which is equivalent to modulo 8). See also the documentation of - // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation. - face_orientation = (is_interior_face == (face.face_orientation >= 8)) ? - (face.face_orientation % 8) : - 0; + + +template +inline const std::array & +FEEvaluationBaseData:: + get_all_face_orientations() const +{ + Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + Assert(is_face && + dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + is_interior_face == false, + ExcMessage("All face numbers can only be queried for ECL at exterior " + "faces. Use get_face_no() in other cases.")); + + return all_face_orientations; } diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index cf562d5740..a434fe9b01 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -1179,7 +1179,7 @@ namespace internal MappingInfoStorage temp_data; temp_data.copy_descriptor(my_data); FEEvaluationBaseData eval( - {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0); + {&shape_info, nullptr, &temp_data, 0, 0}); AlignedVector evaluation_data; eval.set_data_pointers(&evaluation_data, dim); @@ -2091,9 +2091,9 @@ namespace internal MappingInfoStorage temp_data; temp_data.copy_descriptor(my_data); FEEvaluationBaseData eval_int( - {&shape_info, nullptr, &temp_data, 0, 0}, 0, true, 0); + {&shape_info, nullptr, &temp_data, 0, 0}, true); FEEvaluationBaseData eval_ext( - {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0); + {&shape_info, nullptr, &temp_data, 0, 0}, false); AlignedVector evaluation_data, evaluation_data_ext; eval_int.set_data_pointers(&evaluation_data, dim); diff --git a/tests/matrix_free/interpolate_quadrature_01.cc b/tests/matrix_free/interpolate_quadrature_01.cc index 6706195df7..4ce2186a4e 100644 --- a/tests/matrix_free/interpolate_quadrature_01.cc +++ b/tests/matrix_free/interpolate_quadrature_01.cc @@ -120,14 +120,12 @@ test(const unsigned int n_refinements = 1) VectorizedArrayType>:: template interpolate_quadrature( 1, + EvaluationFlags::values, matrix_free.get_shape_info(), phi.begin_values(), temp.data(), - false, - false, face); - for (unsigned int q = 0; q < phi_m.n_q_points; ++q) { const auto u_cell = temp[q]; -- 2.39.5