From bcd0ab7b9feeeebed407de0bc9bb3f017d5c59f9 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 28 Nov 2021 19:34:07 +0100 Subject: [PATCH] MatrixFree: store refinement configuration once per cell --- include/deal.II/matrix_free/dof_info.h | 6 ++ .../deal.II/matrix_free/dof_info.templates.h | 39 +++++++++--- include/deal.II/matrix_free/fe_evaluation.h | 35 +++++------ .../matrix_free/matrix_free.templates.h | 3 +- include/deal.II/matrix_free/tools.h | 10 +-- source/matrix_free/dof_info.cc | 61 +++++++++++-------- 6 files changed, 96 insertions(+), 58 deletions(-) diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 7b0583ea9e..1dc6024ff1 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -526,6 +526,12 @@ namespace internal */ std::vector dof_indices; + /** + * Supported components of all entries of the hp::FECollection object of + * the given DoFHandler. + */ + std::vector> hanging_node_constraint_masks_comp; + /** * Masks indicating for each cell and component if the optimized * hanging-node constraint is applicable and if yes which type. diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 63ff5f22a3..9e3571a154 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -283,13 +283,38 @@ namespace internal const TriaIterator> &cell, std::vector & dof_indices) { - const ArrayView mask_view( - hanging_node_constraint_masks.data() + - cell_number * cell->get_fe().n_components(), - cell->get_fe().n_components()); - - return hanging_nodes.setup_constraints( - cell, {}, lexicographic_mapping, dof_indices, mask_view); + // 1) check if finite elements support fast hanging-node algorithm + this->hanging_node_constraint_masks_comp = + hanging_nodes.compute_supported_components( + cell->get_dof_handler().get_fe_collection()); + + if ([](const auto &supported_components) { + return std::none_of(supported_components.begin(), + supported_components.end(), + [](const auto &a) { + return *std::max_element(a.begin(), a.end()); + }); + }(hanging_node_constraint_masks_comp)) + return false; + + // 2) determine the refinement configuration of the cell + const auto refinement_configuration = + hanging_nodes.compute_refinement_configuration(cell); + + if (refinement_configuration == ConstraintKinds::unconstrained) + return false; + + // 3) update DoF indices of cell for specified components + hanging_nodes.update_dof_indices(cell, + {}, + lexicographic_mapping, + hanging_node_constraint_masks_comp, + refinement_configuration, + dof_indices); + + hanging_node_constraint_masks[cell_number] = refinement_configuration; + + return true; } diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 647edf9208..e5f58db727 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3325,10 +3325,10 @@ FEEvaluationBase:: for (unsigned int v = 0; v < n_lanes; ++v) if (cells[v] != numbers::invalid_unsigned_int && this->dof_info->hanging_node_constraint_masks.size() > 0 && - this->dof_info->hanging_node_constraint_masks - [cells[v] * this->n_fe_components + - this->first_selected_component] != - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) + this->dof_info->hanging_node_constraint_masks[cells[v]] != + internal::MatrixFreeFunctions::ConstraintKinds::unconstrained && + this->dof_info->hanging_node_constraint_masks_comp + [this->active_fe_index][this->first_selected_component]) has_hn_constraints = true; } @@ -3427,10 +3427,10 @@ FEEvaluationBase:: has_constraints = true; if (this->dof_info->hanging_node_constraint_masks.size() > 0 && - this->dof_info->hanging_node_constraint_masks - [cells[v] * this->n_fe_components + - this->first_selected_component] != - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) + this->dof_info->hanging_node_constraint_masks[cells[v]] != + internal::MatrixFreeFunctions::ConstraintKinds::unconstrained && + this->dof_info->hanging_node_constraint_masks_comp + [this->active_fe_index][this->first_selected_component]) has_hn_constraints = true; Assert(my_index_start[n_components_read].first == @@ -3514,9 +3514,11 @@ FEEvaluationBase:: (this->dof_info->row_starts[cell_dof_index].second != this->dof_info->row_starts[cell_dof_index + n_components_read] .second || - (this->dof_info->hanging_node_constraint_masks.size() > 0 && - this->dof_info->hanging_node_constraint_masks[cell_dof_index] != - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained))) + ((this->dof_info->hanging_node_constraint_masks.size() > 0 && + this->dof_info->hanging_node_constraint_masks[cell_index] != + internal::MatrixFreeFunctions::ConstraintKinds::unconstrained) && + this->dof_info->hanging_node_constraint_masks_comp + [this->active_fe_index][this->first_selected_component]))) { Assert(this->dof_info->row_starts_plain_indices[cell_index] != numbers::invalid_unsigned_int, @@ -4152,7 +4154,9 @@ FEEvaluationBase:: apply_hanging_node_constraints(const bool transpose) const { if (this->dof_info == nullptr || - this->dof_info->hanging_node_constraint_masks.size() == 0) + this->dof_info->hanging_node_constraint_masks.size() == 0 || + this->dof_info->hanging_node_constraint_masks_comp + [this->active_fe_index][this->first_selected_component] == false) return; // nothing to do with faces constexpr unsigned int n_lanes = VectorizedArrayType::size(); @@ -4174,11 +4178,8 @@ FEEvaluationBase:: } const unsigned int cell_index = cells[v]; - const unsigned int cell_dof_index = - cell_index * this->n_fe_components + this->first_selected_component; - - const auto mask = - this->dof_info->hanging_node_constraint_masks[cell_dof_index]; + const auto mask = + this->dof_info->hanging_node_constraint_masks[cell_index]; constraint_mask[v] = mask; hn_available |= diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 7a341ee05e..f254e9ed87 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1236,8 +1236,7 @@ namespace internal hanging_nodes = std::make_unique< dealii::internal::MatrixFreeFunctions::HangingNodes>(tria); for (unsigned int no = 0; no < n_dof_handlers; ++no) - dof_info[no].hanging_node_constraint_masks.resize( - n_active_cells * dof_handler[no]->get_fe().n_components()); + dof_info[no].hanging_node_constraint_masks.resize(n_active_cells); } for (unsigned int counter = 0; counter < n_active_cells; ++counter) diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index a2429989e6..39d52486a3 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -447,13 +447,13 @@ namespace MatrixFreeTools locally_relevant_constrains.end()); // STEP 2c: apply hanging-node constraints - if (dof_info.hanging_node_constraint_masks.size() > 0) + if (dof_info.hanging_node_constraint_masks.size() > 0 && + dof_info + .hanging_node_constraint_masks_comp[phi.get_active_fe_index()] + [first_selected_component]) { const auto mask = - dof_info - .hanging_node_constraint_masks[(cell * n_lanes + v) * - n_fe_components + - first_selected_component]; + dof_info.hanging_node_constraint_masks[cell * n_lanes + v]; // cell has hanging nodes if (mask != dealii::internal::MatrixFreeFunctions:: diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc index 6bd7669f11..c9e2717eab 100644 --- a/source/matrix_free/dof_info.cc +++ b/source/matrix_free/dof_info.cc @@ -104,8 +104,9 @@ namespace internal // one const bool has_constraints = (hanging_node_constraint_masks.size() != 0 && - hanging_node_constraint_masks[ib] != - ConstraintKinds::unconstrained) || + hanging_node_constraint_masks[cell * n_vectorization + v] != + ConstraintKinds::unconstrained && + hanging_node_constraint_masks_comp[fe_index][0 /*TODO*/]) || (row_starts[ib].second != row_starts[ib + n_fe_components].second); auto do_copy = [&](const unsigned int *begin, @@ -226,13 +227,19 @@ namespace internal { bool has_hanging_nodes = false; - if (hanging_node_constraint_masks.size() > 0) + const unsigned int fe_index = + (cell_active_fe_index.size() == 0 || + dofs_per_cell.size() == 1) ? + 0 : + cell_active_fe_index[boundary_cells[i]]; + AssertIndexRange(fe_index, dofs_per_cell.size()); + + if (hanging_node_constraint_masks.size() > 0 && + hanging_node_constraint_masks[boundary_cells[i]] != + ConstraintKinds::unconstrained) for (unsigned int comp = 0; comp < n_components; ++comp) has_hanging_nodes |= - hanging_node_constraint_masks[boundary_cells[i] * - n_components + - comp] != - ConstraintKinds::unconstrained; + hanging_node_constraint_masks_comp[fe_index][comp]; if (has_hanging_nodes || row_starts[boundary_cells[i] * n_components].second != @@ -242,12 +249,6 @@ namespace internal unsigned int *data_ptr = plain_dof_indices.data() + row_starts_plain_indices[boundary_cells[i]]; - const unsigned int fe_index = - (cell_active_fe_index.size() == 0 || - dofs_per_cell.size() == 1) ? - 0 : - cell_active_fe_index[boundary_cells[i]]; - AssertIndexRange(fe_index, dofs_per_cell.size()); const unsigned int *row_end = data_ptr + dofs_per_cell[fe_index]; for (; data_ptr != row_end; ++data_ptr) @@ -340,7 +341,7 @@ namespace internal std::vector new_hanging_node_constraint_masks; new_hanging_node_constraint_masks.reserve( - new_hanging_node_constraint_masks.size()); + hanging_node_constraint_masks.size()); if (store_plain_indices == true) { @@ -373,6 +374,19 @@ namespace internal bool has_hanging_nodes = false; + if (hanging_node_constraint_masks.size() > 0) + { + const auto mask = + hanging_node_constraint_masks[renumbering[position_cell + + j]]; + new_hanging_node_constraint_masks.push_back(mask); + + if (mask != ConstraintKinds::unconstrained) + 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]; + } + for (unsigned int comp = 0; comp < n_components; ++comp) { new_row_starts[(i * vectorization_length + j) * n_components + @@ -382,15 +396,6 @@ namespace internal comp] .second = new_constraint_indicator.size(); - if (hanging_node_constraint_masks.size() > 0) - { - const auto mask = - hanging_node_constraint_masks[cell_no + comp]; - new_hanging_node_constraint_masks.push_back(mask); - has_hanging_nodes |= - mask != ConstraintKinds::unconstrained; - } - new_dof_indices.insert( new_dof_indices.end(), dof_indices.data() + row_starts[cell_no + comp].first, @@ -426,11 +431,13 @@ namespace internal new_row_starts[(i * vectorization_length + j) * n_components + comp] .second = new_constraint_indicator.size(); - - if (hanging_node_constraint_masks.size() > 0) - new_hanging_node_constraint_masks.push_back( - ConstraintKinds::unconstrained); } + + 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); + position_cell += n_vect; } AssertDimension(position_cell * n_components + 1, row_starts.size()); -- 2.39.5