From: Peter Munch Date: Sun, 6 Mar 2022 16:52:56 +0000 (+0100) Subject: Use fast HN algorithm for polynomial GC X-Git-Tag: v9.4.0-rc1~401^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=53016842c904a81a43c133be9449b826efd0d784;p=dealii.git Use fast HN algorithm for polynomial GC --- 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 4e4e88db04..b1ffe7c3d4 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -395,6 +395,74 @@ namespace internal return matrix; } + template + bool + use_fast_hanging_node_algorithm(const DoFHandler &dof_handler_coarse, + const unsigned int mg_level_coarse) + { + // algorithm is only needed on active levels + bool use_fast_hanging_node_algorithm = + mg_level_coarse == numbers::invalid_unsigned_int; + + // algorithm can be only used on meshes consisting of hypercube and + // simplices + if (use_fast_hanging_node_algorithm) + { + const auto &reference_cells = + dof_handler_coarse.get_triangulation().get_reference_cells(); + use_fast_hanging_node_algorithm = + std::all_of(reference_cells.begin(), + reference_cells.end(), + [](const auto &r) { + return r.is_hyper_cube() || r.is_simplex(); + }); + } + + // local p-refinement is not supported + if (use_fast_hanging_node_algorithm) + { + const auto &fes = dof_handler_coarse.get_fe_collection(); + + use_fast_hanging_node_algorithm &= + std::all_of(fes.begin(), fes.end(), [&fes](const auto &fe) { + return fes[0].compare_for_domination(fe) == + FiniteElementDomination::Domination:: + either_element_can_dominate; + }); + } + + // check that all components are either supported or not + if (use_fast_hanging_node_algorithm) + { + const std::vector> supported_components = + internal::MatrixFreeFunctions::HangingNodes< + dim>::compute_supported_components(dof_handler_coarse + .get_fe_collection()); + + use_fast_hanging_node_algorithm &= std::any_of( + supported_components.begin(), + supported_components.end(), + [](const auto &supported_components_per_fe) { + return std::all_of(supported_components_per_fe.begin(), + supported_components_per_fe.end(), + [](const auto &a) { return a == true; }); + }); + + use_fast_hanging_node_algorithm &= std::all_of( + supported_components.begin(), + supported_components.end(), + [](const auto &supported_components_per_fe) { + return std::all_of(supported_components_per_fe.begin(), + supported_components_per_fe.end(), + [&supported_components_per_fe](const auto &a) { + return a == supported_components_per_fe[0]; + }); + }); + } + + return use_fast_hanging_node_algorithm; + } + } // namespace class FineDoFHandlerViewCell @@ -1288,105 +1356,32 @@ namespace internal hanging_nodes; std::vector> lexicographic_mappings; - { - // algorithm is only needed on active levels - bool use_fast_hanging_node_algorithm = - mg_level_coarse == numbers::invalid_unsigned_int; - - const auto &tria = dof_handler_coarse.get_triangulation(); - - // algorithm can be only used on meshes consisting of hypercube and - // simplices - if (use_fast_hanging_node_algorithm) - { - const auto &reference_cells = tria.get_reference_cells(); - use_fast_hanging_node_algorithm = - std::all_of(reference_cells.begin(), - reference_cells.end(), - [](const auto &r) { - return r.is_hyper_cube() || r.is_simplex(); - }); - } - - // local p-refinement is not supported - if (use_fast_hanging_node_algorithm) - { - const auto &fes = dof_handler_coarse.get_fe_collection(); - - use_fast_hanging_node_algorithm &= - std::all_of(fes.begin(), fes.end(), [&fes](const auto &fe) { - return fes[0].compare_for_domination(fe) == - FiniteElementDomination::Domination:: - either_element_can_dominate; - }); - } - - // check that all components are either supported or not - if (use_fast_hanging_node_algorithm) - { - const std::vector> supported_components = - internal::MatrixFreeFunctions::HangingNodes< - dim>::compute_supported_components(dof_handler_coarse - .get_fe_collection()); - - use_fast_hanging_node_algorithm &= std::any_of( - supported_components.begin(), - supported_components.end(), - [](const auto &supported_components_per_fe) { - return std::all_of(supported_components_per_fe.begin(), - supported_components_per_fe.end(), - [](const auto &a) { return a == true; }); - }); - - use_fast_hanging_node_algorithm &= - std::all_of(supported_components.begin(), - supported_components.end(), - [](const auto &supported_components_per_fe) { - return std::all_of( - supported_components_per_fe.begin(), - supported_components_per_fe.end(), - [&supported_components_per_fe](const auto &a) { - return a == supported_components_per_fe[0]; - }); - }); - } - - use_fast_hanging_node_algorithm = false; - - // create helper class - if (use_fast_hanging_node_algorithm) - { - hanging_nodes = std::make_unique< - internal::MatrixFreeFunctions::HangingNodes>(tria); - - const auto &fes = dof_handler_coarse.get_fe_collection(); - lexicographic_mappings.resize(fes.size()); + // create helper class + if (use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse)) + { + hanging_nodes = + std::make_unique>( + dof_handler_coarse.get_triangulation()); - for (unsigned int i = 0; i < fes.size(); ++i) - if (fes[i].reference_cell().is_hyper_cube()) - { - const Quadrature<1> dummy_quadrature( - std::vector>(1, Point<1>())); - internal::MatrixFreeFunctions::ShapeInfo shape_info; - shape_info.reinit(dummy_quadrature, fes[i], 0); - lexicographic_mappings[i] = - shape_info.lexicographic_numbering; - } + const auto &fes = dof_handler_coarse.get_fe_collection(); + lexicographic_mappings.resize(fes.size()); - if (fes[fe_index_coarse].reference_cell().is_hyper_cube()) + for (unsigned int i = 0; i < fes.size(); ++i) + if (fes[i].reference_cell().is_hyper_cube()) { const Quadrature<1> dummy_quadrature( std::vector>(1, Point<1>())); internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray> shape_info; - shape_info.reinit(dummy_quadrature, fes[fe_index_coarse], 0); + shape_info.reinit(dummy_quadrature, fes[i], 0); + lexicographic_mappings[i] = shape_info.lexicographic_numbering; - transfer.schemes[0].shape_info_coarse = - transfer.schemes[1].shape_info_coarse = shape_info; + if (i == fe_index_coarse) + transfer.schemes[0].shape_info_coarse = + transfer.schemes[1].shape_info_coarse = shape_info; } - } - } + } // create partitioners and vectors for internal purposes { @@ -2002,9 +1997,35 @@ namespace internal return fe.n_base_elements() == 1 && (dynamic_cast *>( &fe.base_element(0)) != nullptr); - ; }); + std::unique_ptr> + hanging_nodes; + std::vector> lexicographic_mappings; + + // create helper class + if (use_fast_hanging_node_algorithm(dof_handler_coarse, mg_level_coarse)) + { + hanging_nodes = + std::make_unique>( + dof_handler_coarse.get_triangulation()); + + const auto &fes = dof_handler_coarse.get_fe_collection(); + lexicographic_mappings.resize(fes.size()); + + for (unsigned int i = 0; i < fes.size(); ++i) + if (fes[i].reference_cell().is_hyper_cube()) + { + const Quadrature<1> dummy_quadrature( + std::vector>(1, Point<1>())); + internal::MatrixFreeFunctions::ShapeInfo< + VectorizedArray> + shape_info; + shape_info.reinit(dummy_quadrature, fes[i], 0); + lexicographic_mappings[i] = shape_info.lexicographic_numbering; + } + } + const auto process_cells = [&](const auto &fu) { loop_over_active_or_level_cells( dof_handler_coarse, mg_level_coarse, [&](const auto &cell_coarse) { @@ -2099,6 +2120,8 @@ namespace internal fe_index_pairs.size()); std::vector> local_dof_indices_coarse(fe_index_pairs.size()); + std::vector> + local_dof_indices_coarse_lex(fe_index_pairs.size()); std::vector> local_dof_indices_fine(fe_index_pairs.size()); @@ -2106,6 +2129,8 @@ namespace internal { local_dof_indices_coarse[fe_index_pair.second].resize( transfer.schemes[fe_index_pair.second].n_dofs_per_cell_coarse); + local_dof_indices_coarse_lex[fe_index_pair.second].resize( + transfer.schemes[fe_index_pair.second].n_dofs_per_cell_coarse); local_dof_indices_fine[fe_index_pair.second].resize( transfer.schemes[fe_index_pair.second].n_dofs_per_cell_fine); @@ -2130,7 +2155,9 @@ namespace internal { const Quadrature<1> dummy_quadrature( std::vector>(1, Point<1>())); - internal::MatrixFreeFunctions::ShapeInfo shape_info; + internal::MatrixFreeFunctions::ShapeInfo< + VectorizedArray> + shape_info; shape_info.reinit(dummy_quadrature, dof_handler_fine.get_fe( fe_index_pair.first.second), @@ -2144,6 +2171,10 @@ namespace internal 0); lexicographic_numbering_coarse[fe_index_pair.second] = shape_info.lexicographic_numbering; + + if (hanging_nodes) + transfer.schemes[fe_index_pair.second].shape_info_coarse = + shape_info; } else { @@ -2174,24 +2205,55 @@ namespace internal } transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back()); + if (hanging_nodes) + transfer.level_dof_indices_coarse.resize(n_dof_indices_coarse.back()); transfer.level_dof_indices_coarse_plain.resize( n_dof_indices_coarse.back()); // ------------------------------ indices ----------------------------- - std::vector level_dof_indices_coarse_( + std::vector level_dof_indices_coarse( + fe_index_pairs.size()); + std::vector level_dof_indices_coarse_plain( fe_index_pairs.size()); - std::vector level_dof_indices_fine_( + std::vector level_dof_indices_fine( fe_index_pairs.size()); for (unsigned int i = 0; i < fe_index_pairs.size(); ++i) { - level_dof_indices_fine_[i] = + level_dof_indices_fine[i] = transfer.level_dof_indices_fine.data() + n_dof_indices_fine[i]; - level_dof_indices_coarse_[i] = + level_dof_indices_coarse[i] = + transfer.level_dof_indices_coarse.data() + + n_dof_indices_coarse[i]; + level_dof_indices_coarse_plain[i] = transfer.level_dof_indices_coarse_plain.data() + n_dof_indices_coarse[i]; } + std::vector + coarse_cell_refinement_configurations; + + if (hanging_nodes) + { + unsigned int n_coarse_cells = 0; + for (unsigned int i = 0; i < transfer.schemes.size(); ++i) + n_coarse_cells += transfer.schemes[i].n_coarse_cells; + transfer.coarse_cell_refinement_configurations.resize( + n_coarse_cells); + + coarse_cell_refinement_configurations.resize(fe_index_pairs.size()); + + if (fe_index_pairs.size() > 0) + { + coarse_cell_refinement_configurations[0] = + transfer.coarse_cell_refinement_configurations.data(); + for (unsigned int i = 1; i < transfer.schemes.size(); ++i) + coarse_cell_refinement_configurations[i] = + coarse_cell_refinement_configurations[i - 1] + + transfer.schemes[i - 1].n_coarse_cells; + } + } + bool fine_indices_touch_remote_dofs = false; const IndexSet locally_owned_dofs = mg_level_fine == numbers::invalid_unsigned_int ? @@ -2203,40 +2265,86 @@ namespace internal fe_index_pairs[std::pair( cell_coarse->active_fe_index(), cell_fine->active_fe_index())]; - if (mg_level_coarse == numbers::invalid_unsigned_int) - cell_coarse->get_dof_indices(local_dof_indices_coarse[fe_pair_no]); - else - cell_coarse->get_mg_dof_indices( - local_dof_indices_coarse[fe_pair_no]); - - for (unsigned int i = 0; - i < transfer.schemes[fe_pair_no].n_dofs_per_cell_coarse; - i++) - level_dof_indices_coarse_[fe_pair_no][i] = - transfer.partitioner_coarse->global_to_local( + // parent + { + if (mg_level_coarse == numbers::invalid_unsigned_int) + cell_coarse->get_dof_indices( + local_dof_indices_coarse[fe_pair_no]); + else + cell_coarse->get_mg_dof_indices( + local_dof_indices_coarse[fe_pair_no]); + + for (unsigned int i = 0; + i < local_dof_indices_coarse[fe_pair_no].size(); + ++i) + local_dof_indices_coarse_lex[fe_pair_no][i] = local_dof_indices_coarse - [fe_pair_no][lexicographic_numbering_coarse[fe_pair_no][i]]); + [fe_pair_no][lexicographic_numbering_coarse[fe_pair_no][i]]; - cell_fine->get_dof_indices(local_dof_indices_fine[fe_pair_no]); - for (unsigned int i = 0; - i < transfer.schemes[fe_pair_no].n_dofs_per_cell_fine; - i++) - { - level_dof_indices_fine_[fe_pair_no][i] = - transfer.partitioner_fine->global_to_local( - local_dof_indices_fine - [fe_pair_no][lexicographic_numbering_fine[fe_pair_no][i]]); - if (!locally_owned_dofs.is_element( + for (auto &i : local_dof_indices_coarse_lex[fe_pair_no]) + i = transfer.partitioner_coarse->global_to_local(i); + + if (hanging_nodes) + { + std::vector + mask(transfer.n_components); + + local_dof_indices_coarse[fe_pair_no] = + local_dof_indices_coarse_lex[fe_pair_no]; + + hanging_nodes->setup_constraints( + cell_coarse, + transfer.partitioner_coarse, + lexicographic_mappings, + local_dof_indices_coarse[fe_pair_no], + mask); + + for (unsigned int i = 0; + i < transfer.schemes[fe_pair_no].n_dofs_per_cell_coarse; + i++) + level_dof_indices_coarse[fe_pair_no][i] = + local_dof_indices_coarse[fe_pair_no][i]; + + coarse_cell_refinement_configurations[fe_pair_no][0] = mask[0]; + coarse_cell_refinement_configurations[fe_pair_no] += 1; + } + + for (unsigned int i = 0; + i < transfer.schemes[fe_pair_no].n_dofs_per_cell_coarse; + i++) + level_dof_indices_coarse_plain[fe_pair_no][i] = + local_dof_indices_coarse_lex[fe_pair_no][i]; + } + + // child + { + cell_fine->get_dof_indices(local_dof_indices_fine[fe_pair_no]); + for (unsigned int i = 0; + i < transfer.schemes[fe_pair_no].n_dofs_per_cell_fine; + i++) + { + level_dof_indices_fine[fe_pair_no][i] = + transfer.partitioner_fine->global_to_local( local_dof_indices_fine [fe_pair_no] - [lexicographic_numbering_fine[fe_pair_no][i]])) - fine_indices_touch_remote_dofs = true; - } + [lexicographic_numbering_fine[fe_pair_no][i]]); + if (!locally_owned_dofs.is_element( + local_dof_indices_fine + [fe_pair_no] + [lexicographic_numbering_fine[fe_pair_no][i]])) + fine_indices_touch_remote_dofs = true; + } + } - level_dof_indices_coarse_[fe_pair_no] += - transfer.schemes[fe_pair_no].n_dofs_per_cell_coarse; - level_dof_indices_fine_[fe_pair_no] += - transfer.schemes[fe_pair_no].n_dofs_per_cell_fine; + // move pointers + { + level_dof_indices_coarse[fe_pair_no] += + transfer.schemes[fe_pair_no].n_dofs_per_cell_coarse; + level_dof_indices_coarse_plain[fe_pair_no] += + transfer.schemes[fe_pair_no].n_dofs_per_cell_coarse; + level_dof_indices_fine[fe_pair_no] += + transfer.schemes[fe_pair_no].n_dofs_per_cell_fine; + } }); // if all access goes to the locally owned dofs on all ranks, we do @@ -2245,6 +2353,15 @@ namespace internal fine_indices_touch_remote_dofs), comm) == 0) transfer.vec_fine.reinit(0); + + if (hanging_nodes && + std::all_of(transfer.coarse_cell_refinement_configurations.begin(), + transfer.coarse_cell_refinement_configurations.end(), + [](const auto i) { + return i == internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained; + })) + transfer.coarse_cell_refinement_configurations.clear(); } // ------------------------- prolongation matrix ------------------------- @@ -2700,10 +2817,8 @@ MGTwoLevelTransfer>:: const unsigned int *indices_coarse = level_dof_indices_coarse.size() > 0 ? level_dof_indices_coarse.data() : level_dof_indices_coarse_plain.data(); - const unsigned int *indices_coarse_plain = - level_dof_indices_coarse_plain.data(); - const unsigned int *indices_fine = level_dof_indices_fine.data(); - const Number * weights = nullptr; + const unsigned int *indices_fine = level_dof_indices_fine.data(); + const Number * weights = nullptr; const std::array, Utilities::pow(3, dim)> *weights_compressed = nullptr; @@ -2748,7 +2863,6 @@ MGTwoLevelTransfer>:: evaluation_data_coarse[i][v] = read_dof_values(indices_coarse[i], this->vec_coarse); indices_coarse += scheme.n_dofs_per_cell_coarse; - indices_coarse_plain += scheme.n_dofs_per_cell_coarse; } // ... fast hanging-node-constraints algorithm @@ -2869,10 +2983,8 @@ MGTwoLevelTransfer>:: const unsigned int *indices_coarse = level_dof_indices_coarse.size() > 0 ? level_dof_indices_coarse.data() : level_dof_indices_coarse_plain.data(); - const unsigned int *indices_coarse_plain = - level_dof_indices_coarse_plain.data(); - const unsigned int *indices_fine = level_dof_indices_fine.data(); - const Number * weights = nullptr; + const unsigned int *indices_fine = level_dof_indices_fine.data(); + const Number * weights = nullptr; const std::array, Utilities::pow(3, dim)> *weights_compressed = nullptr; @@ -2979,7 +3091,6 @@ MGTwoLevelTransfer>:: evaluation_data_coarse[i][v], this->vec_coarse); indices_coarse += scheme.n_dofs_per_cell_coarse; - indices_coarse_plain += scheme.n_dofs_per_cell_coarse; } } }