From: Peter Munch Date: Wed, 23 Mar 2022 12:05:31 +0000 (+0100) Subject: Use compressed_constraint_kind in MatrixFree (CPU) X-Git-Tag: v9.4.0-rc1~341^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=274c5381df741b551a41fb6cd466195dcaf5b90f;p=dealii.git Use compressed_constraint_kind in MatrixFree (CPU) --- diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h index dbb635a2f6..dd5c97a483 100644 --- a/include/deal.II/matrix_free/constraint_info.h +++ b/include/deal.II/matrix_free/constraint_info.h @@ -100,9 +100,9 @@ namespace internal std::vector constraint_pool_data; std::vector constraint_pool_row_index; - std::vector> shape_infos; - std::vector hanging_node_constraint_masks; - std::vector active_fe_indices; + std::vector> shape_infos; + std::vector hanging_node_constraint_masks; + std::vector active_fe_indices; private: inline const typename Number::value_type * @@ -231,7 +231,7 @@ namespace internal hanging_nodes->setup_constraints( cell, {}, lexicographic_numbering, local_dof_indices_lex, mask); - hanging_node_constraint_masks[cell_no] = mask[0]; + hanging_node_constraint_masks[cell_no] = compress(mask[0], dim); active_fe_indices[cell_no] = cell->active_fe_index(); } @@ -353,12 +353,12 @@ namespace internal dof_indices_per_cell.clear(); constraint_indicator_per_cell.clear(); - if (hanging_nodes && std::all_of(hanging_node_constraint_masks.begin(), - hanging_node_constraint_masks.end(), - [](const auto i) { - return i == - ConstraintKinds::unconstrained; - })) + if (hanging_nodes && + std::all_of(hanging_node_constraint_masks.begin(), + hanging_node_constraint_masks.end(), + [](const auto i) { + return i == unconstrained_compressed_constraint_kind; + })) hanging_node_constraint_masks.clear(); } @@ -459,7 +459,9 @@ namespace internal if (hanging_node_constraint_masks.size() == 0) return; - std::array constraint_mask; + std::array + constraint_mask; bool hn_available = false; @@ -469,13 +471,13 @@ namespace internal constraint_mask[v] = mask; - hn_available |= (mask != ConstraintKinds::unconstrained); + hn_available |= (mask != unconstrained_compressed_constraint_kind); } if (hn_available == true) { for (unsigned int v = n_lanes_filled; v < Number::size(); ++v) - constraint_mask[v] = ConstraintKinds::unconstrained; + constraint_mask[v] = unconstrained_compressed_constraint_kind; for (unsigned int i = 1; i < n_lanes_filled; ++i) AssertDimension(active_fe_indices[first_cell], diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 1dc6024ff1..0f986441d3 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -536,7 +536,7 @@ namespace internal * Masks indicating for each cell and component if the optimized * hanging-node constraint is applicable and if yes which type. */ - std::vector hanging_node_constraint_masks; + std::vector hanging_node_constraint_masks; /** * This variable describes the position of constraints in terms of the diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 2b4d859e16..f884701828 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -301,7 +301,8 @@ namespace internal refinement_configuration, dof_indices); - hanging_node_constraint_masks[cell_number] = refinement_configuration; + hanging_node_constraint_masks[cell_number] = + compress(refinement_configuration, dim); return true; } diff --git a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h index b5137069cf..12a97f2bf9 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h +++ b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h @@ -181,13 +181,13 @@ namespace internal public: template static void - run_internal(const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const std::array & constraint_mask, - Number * values) + run_internal( + const unsigned int n_components, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const std::array & constraint_mask, + Number * values) { - using Kinds = MatrixFreeFunctions::ConstraintKinds; const unsigned int given_degree = fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; @@ -205,7 +205,7 @@ namespace internal for (unsigned int v = 0; v < Number::size(); ++v) { - const auto kind = static_cast(constraint_mask[v]); + const auto kind = constraint_mask[v]; const bool subcell_x = (kind >> 0) & 1; const bool subcell_y = (kind >> 1) & 1; const bool face_x = (kind >> 3) & 1; @@ -278,16 +278,17 @@ namespace internal for (unsigned int v = 0; v < Number::size(); ++v) { - const auto kind = static_cast(constraint_mask[v]); + const auto kind = constraint_mask[v]; + const bool subcell_x = (kind >> 0) & 1; const bool subcell_y = (kind >> 1) & 1; const bool subcell_z = (kind >> 2) & 1; - const bool face_x = (kind >> 3) & 1; - const bool face_y = (kind >> 4) & 1; - const bool face_z = (kind >> 5) & 1; - const bool edge_x = (kind >> 6) & 1; - const bool edge_y = (kind >> 7) & 1; - const bool edge_z = (kind >> 8) & 1; + const bool face_x = ((kind >> 3) & 1) ? (kind >> 5) & 1 : 0; + const bool face_y = ((kind >> 3) & 1) ? (kind >> 6) & 1 : 0; + const bool face_z = ((kind >> 3) & 1) ? (kind >> 7) & 1 : 0; + const bool edge_x = ((kind >> 4) & 1) ? (kind >> 5) & 1 : 0; + const bool edge_y = ((kind >> 4) & 1) ? (kind >> 6) & 1 : 0; + const bool edge_z = ((kind >> 4) & 1) ? (kind >> 7) & 1 : 0; if (subcell_x) mask_weights[0][v] = 1; @@ -471,11 +472,11 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE unsigned int - create( - const std::array mask, - const std::array - mask_new, - const unsigned int v) + create(const std::array mask, + const std::array mask_new, + const unsigned int v) { (void)mask; (void)mask_new; @@ -483,7 +484,8 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + do_break(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)v; (void)kind; @@ -491,17 +493,19 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_continue(unsigned int v, - const MatrixFreeFunctions::ConstraintKinds &kind) + do_continue(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)v; - return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + return kind == + MatrixFreeFunctions::unconstrained_compressed_constraint_kind; } static inline DEAL_II_ALWAYS_INLINE_RELEASE - std::array + std::array create_mask( - const std::array mask) + const std::array mask) { return mask; } @@ -543,7 +547,8 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + do_break(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)v; (void)kind; @@ -551,19 +556,20 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_continue(unsigned int v, - const MatrixFreeFunctions::ConstraintKinds &kind) + do_continue(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)v; - return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + return kind == + MatrixFreeFunctions::unconstrained_compressed_constraint_kind; } static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type - create( - const std::array mask, - const std::array - mask_new, - const unsigned int v) + create(const std::array mask, + const std::array mask_new, + const unsigned int v) { (void)mask; (void)mask_new; @@ -573,9 +579,10 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE - std::array + std::array create_mask( - const std::array mask) + const std::array mask) { return mask; } @@ -608,26 +615,29 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + do_break(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)v; - return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + return kind == + MatrixFreeFunctions::unconstrained_compressed_constraint_kind; } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_continue(unsigned int v, - const MatrixFreeFunctions::ConstraintKinds &kind) + do_continue(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)v; - return kind == MatrixFreeFunctions::ConstraintKinds::unconstrained; + return kind == + MatrixFreeFunctions::unconstrained_compressed_constraint_kind; } static inline DEAL_II_ALWAYS_INLINE_RELEASE index_type - create( - const std::array mask, - const std::array - mask_new, - const unsigned int v) + create(const std::array mask, + const std::array mask_new, + const unsigned int v) { T1 result; @@ -638,16 +648,17 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE - std::array + std::array create_mask( - const std::array mask) + const std::array mask) { auto new_mask = mask; std::sort(new_mask.begin(), new_mask.end()); std::fill(std::unique(new_mask.begin(), new_mask.end()), new_mask.end(), - MatrixFreeFunctions::ConstraintKinds::unconstrained); + MatrixFreeFunctions::unconstrained_compressed_constraint_kind); return new_mask; } @@ -680,15 +691,16 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_break(unsigned int v, const MatrixFreeFunctions::ConstraintKinds &kind) + do_break(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)kind; return v > 0; } static inline DEAL_II_ALWAYS_INLINE_RELEASE bool - do_continue(unsigned int v, - const MatrixFreeFunctions::ConstraintKinds &kind) + do_continue(unsigned int v, + const MatrixFreeFunctions::compressed_constraint_kind &kind) { (void)kind; @@ -698,11 +710,11 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE T1 - create( - const std::array mask, - const std::array - mask_new, - const unsigned int v) + create(const std::array mask, + const std::array mask_new, + const unsigned int v) { (void)mask; (void)mask_new; @@ -711,9 +723,10 @@ namespace internal } static inline DEAL_II_ALWAYS_INLINE_RELEASE - std::array + std::array create_mask( - const std::array mask) + const std::array mask) { return mask; } @@ -1514,11 +1527,12 @@ namespace internal public: template static void - run_internal(const unsigned int n_desired_components, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const std::array & constraint_mask, - Number * values) + run_internal( + const unsigned int n_desired_components, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const std::array & constraint_mask, + Number * values) { const unsigned int given_degree = fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; @@ -1548,25 +1562,18 @@ namespace internal if (dim == 2) // 2D: only faces { - const auto is_set = [](const auto a, const auto b) -> bool { - return (a & b) == b; - }; + const bool subcell_x = (mask >> 0) & 1; + const bool subcell_y = (mask >> 1) & 1; + const bool face_x = (mask >> 3) & 1; + const bool face_y = (mask >> 4) & 1; // direction 0: - if ((mask & MatrixFreeFunctions::ConstraintKinds::face_y) != - MatrixFreeFunctions::ConstraintKinds::unconstrained) + if (face_y) { - const bool is_subface_0 = - (mask & - MatrixFreeFunctions::ConstraintKinds::subcell_x) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - const auto *weights = - interpolation_matrices[is_subface_0].data(); + interpolation_matrices[!subcell_x].data(); - if (is_set( - mask, - MatrixFreeFunctions::ConstraintKinds::subcell_y)) + if (subcell_y) interpolate_2D<2, transpose>(given_degree, vv, weights, @@ -1579,20 +1586,12 @@ namespace internal } // direction 1: - if ((mask & MatrixFreeFunctions::ConstraintKinds::face_x) != - MatrixFreeFunctions::ConstraintKinds::unconstrained) + if (face_x) { - const bool is_subface_0 = - (mask & - MatrixFreeFunctions::ConstraintKinds::subcell_y) == - MatrixFreeFunctions::ConstraintKinds::unconstrained; - const auto *weights = - interpolation_matrices[is_subface_0].data(); + interpolation_matrices[!subcell_y].data(); - if (is_set( - mask, - MatrixFreeFunctions::ConstraintKinds::subcell_x)) + if (subcell_x) interpolate_2D<0, transpose>(given_degree, vv, weights, @@ -1606,13 +1605,14 @@ namespace internal } else if (dim == 3) // 3D faces and edges { - const auto m = static_cast(mask); + const bool type_x = (mask >> 0) & 1; + const bool type_y = (mask >> 1) & 1; + const bool type_z = (mask >> 2) & 1; - const bool type_x = (m >> 0) & 1; - const bool type_y = (m >> 1) & 1; - const bool type_z = (m >> 2) & 1; - const auto faces = (m >> 3) & 7; - const auto edges = (m >> 6); + const auto flag_0 = (mask >> 3) & 3; + const auto flag_1 = (mask >> 5) & 7; + const auto faces = (flag_0 & 0b01) ? flag_1 : 0; + const auto edges = (flag_0 & 0b10) ? flag_1 : 0; Helper &shape_info, const bool transpose, - const std::array - & c_mask, - Number *values) + const std::array & c_mask, + Number * values) { using RunnerType = FEEvaluationImplHangingNodesRunner(), diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index 6ea92bf884..af6d8e29d4 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -131,7 +131,7 @@ namespace internal const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo &shape_info, const bool transpose, - const std::array & c_mask, VectorizedArrayType * values); }; diff --git a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h index c1622f2859..234e187961 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h @@ -36,7 +36,7 @@ namespace internal const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo &shape_info, const bool transpose, - const std::array & c_mask, VectorizedArrayType * values) { diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index cbb87aa3f5..159f5e5d29 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3327,7 +3327,8 @@ FEEvaluationBase:: this->dof_info->hanging_node_constraint_masks.size() > 0 && this->dof_info->hanging_node_constraint_masks_comp.size() > 0 && this->dof_info->hanging_node_constraint_masks[cells[v]] != - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained && + internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind && this->dof_info->hanging_node_constraint_masks_comp [this->active_fe_index][this->first_selected_component]) has_hn_constraints = true; @@ -3430,7 +3431,8 @@ FEEvaluationBase:: if (this->dof_info->hanging_node_constraint_masks.size() > 0 && this->dof_info->hanging_node_constraint_masks_comp.size() > 0 && this->dof_info->hanging_node_constraint_masks[cells[v]] != - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained && + internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind && this->dof_info->hanging_node_constraint_masks_comp [this->active_fe_index][this->first_selected_component]) has_hn_constraints = true; @@ -3519,7 +3521,8 @@ FEEvaluationBase:: ((this->dof_info->hanging_node_constraint_masks.size() > 0 && this->dof_info->hanging_node_constraint_masks_comp.size() > 0 && this->dof_info->hanging_node_constraint_masks[cell_index] != - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) && + internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind) && this->dof_info->hanging_node_constraint_masks_comp [this->active_fe_index][this->first_selected_component]))) { @@ -4164,7 +4167,7 @@ FEEvaluationBase:: return; // nothing to do with faces constexpr unsigned int n_lanes = VectorizedArrayType::size(); - std::array + std::array constraint_mask; bool hn_available = false; @@ -4176,8 +4179,8 @@ FEEvaluationBase:: { if (cells[v] == numbers::invalid_unsigned_int) { - constraint_mask[v] = - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained; + constraint_mask[v] = internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind; continue; } @@ -4186,8 +4189,8 @@ FEEvaluationBase:: this->dof_info->hanging_node_constraint_masks[cell_index]; constraint_mask[v] = mask; - hn_available |= - (mask != internal::MatrixFreeFunctions::ConstraintKinds::unconstrained); + hn_available |= (mask != internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind); } if (hn_available == false) diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index a2f1951528..f406690258 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1379,8 +1379,8 @@ namespace internal { auto &vec = dof_info[no].hanging_node_constraint_masks; if (std::all_of(vec.begin(), vec.end(), [](const auto i) { - return i == internal::MatrixFreeFunctions::ConstraintKinds:: - unconstrained; + return i == internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind; })) vec.clear(); } diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index fc62fe3a3e..9ea4677e80 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -458,7 +458,7 @@ namespace MatrixFreeTools // cell has hanging nodes if (mask != dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) + unconstrained_compressed_constraint_kind) { // check if hanging node internpolation matrix has been set // up @@ -471,9 +471,13 @@ namespace MatrixFreeTools dofs_per_component); std::array constraint_mask; + constraint_mask.fill( + dealii::internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind); + constraint_mask[0] = mask; std::vector< @@ -796,8 +800,9 @@ namespace MatrixFreeTools // constraints! std::array, n_lanes> diagonals_local_constrained; - std::map>> + std::map< + dealii::internal::MatrixFreeFunctions::compressed_constraint_kind, + std::vector>> locally_relevant_constrains_hn_map; }; diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc index 38cfa797ef..c1098500c2 100644 --- a/source/matrix_free/dof_info.cc +++ b/source/matrix_free/dof_info.cc @@ -106,7 +106,7 @@ namespace internal (hanging_node_constraint_masks.size() != 0 && hanging_node_constraint_masks_comp.size() != 0 && hanging_node_constraint_masks[cell * n_vectorization + v] != - ConstraintKinds::unconstrained && + unconstrained_compressed_constraint_kind && hanging_node_constraint_masks_comp[fe_index][0 /*TODO*/]) || (row_starts[ib].second != row_starts[ib + n_fe_components].second); @@ -238,7 +238,7 @@ namespace internal if (hanging_node_constraint_masks.size() > 0 && hanging_node_constraint_masks_comp.size() > 0 && hanging_node_constraint_masks[boundary_cells[i]] != - ConstraintKinds::unconstrained) + unconstrained_compressed_constraint_kind) for (unsigned int comp = 0; comp < n_components; ++comp) has_hanging_nodes |= hanging_node_constraint_masks_comp[fe_index][comp]; @@ -341,7 +341,7 @@ namespace internal new_dof_indices.reserve(dof_indices.size()); new_constraint_indicator.reserve(constraint_indicator.size()); - std::vector new_hanging_node_constraint_masks; + std::vector new_hanging_node_constraint_masks; new_hanging_node_constraint_masks.reserve( hanging_node_constraint_masks.size()); @@ -384,7 +384,7 @@ namespace internal j]]; new_hanging_node_constraint_masks.push_back(mask); - if (mask != ConstraintKinds::unconstrained) + if (mask != unconstrained_compressed_constraint_kind) for (unsigned int comp = 0; comp < n_components; ++comp) has_hanging_nodes |= hanging_node_constraint_masks_comp [have_hp ? cell_active_fe_index[i] : 0][comp]; @@ -439,7 +439,7 @@ namespace internal for (unsigned int j = n_vect; j < vectorization_length; ++j) if (hanging_node_constraint_masks.size() > 0) new_hanging_node_constraint_masks.push_back( - ConstraintKinds::unconstrained); + unconstrained_compressed_constraint_kind); position_cell += n_vect; } diff --git a/tests/matrix_free/hanging_node_kernels_01.cc b/tests/matrix_free/hanging_node_kernels_01.cc index da65765ed6..12feedb4d7 100644 --- a/tests/matrix_free/hanging_node_kernels_01.cc +++ b/tests/matrix_free/hanging_node_kernels_01.cc @@ -450,14 +450,23 @@ test(const unsigned int degree, FEEvaluation eval(matrix_free); eval.reinit(0); - std::array::size()> cmask; + std::fill(cmask.begin(), + cmask.end(), + dealii::internal::MatrixFreeFunctions:: + unconstrained_compressed_constraint_kind); + cmask[0] = dealii::internal::MatrixFreeFunctions::compress(mask_value, dim); + + std::array::size()> + cmask_; std::fill( - cmask.begin(), - cmask.end(), + cmask_.begin(), + cmask_.end(), dealii::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained); - cmask[0] = mask_value; + cmask_[0] = mask_value; for (unsigned int b = 0; b < 2; ++b) { @@ -477,7 +486,7 @@ test(const unsigned int degree, internal::FEEvaluationImplHangingNodesReference< dim, VectorizedArray, - false>::template run<-1, -1>(eval, b == 1, cmask, values1.data()); + false>::template run<-1, -1>(eval, b == 1, cmask_, values1.data()); internal::FEEvaluationImplHangingNodes>:: template run<-1, -1>( 1, eval.get_shape_info(), b == 1, cmask, values2.data());