From bad4c4d9f25daca2e9f255718c7d47940747b69d Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Wed, 17 Mar 2021 19:59:12 -0600 Subject: [PATCH] Move DoFhandler::prepare_c_and_r() to hp::Refinement::limit_p_level_difference(). --- doc/news/changes/minor/20210223Fehling | 2 +- include/deal.II/dofs/dof_handler.h | 30 --- include/deal.II/hp/refinement.h | 44 +++- source/dofs/dof_handler.cc | 204 ---------------- source/hp/refinement.cc | 219 +++++++++++++++++- source/hp/refinement.inst.in | 7 + .../prepare_coarsening_and_refinement_01.cc | 5 +- .../prepare_coarsening_and_refinement_02.cc | 6 +- .../prepare_coarsening_and_refinement_03.cc | 5 +- .../prepare_coarsening_and_refinement_01.cc | 5 +- .../prepare_coarsening_and_refinement_02.cc | 6 +- .../prepare_coarsening_and_refinement_03.cc | 5 +- 12 files changed, 292 insertions(+), 246 deletions(-) diff --git a/doc/news/changes/minor/20210223Fehling b/doc/news/changes/minor/20210223Fehling index 5b57fa6d39..79d740b4ea 100644 --- a/doc/news/changes/minor/20210223Fehling +++ b/doc/news/changes/minor/20210223Fehling @@ -1,4 +1,4 @@ -New: Member function DoFHandler::prepare_coarsening_and_refinement() +New: Function hp::Refinement::limit_p_level_difference() restricts the maximal level difference of neighboring cells with respect to the finite element hierarchy of the registered hp::FECollection.
diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h index 256ba91808..c4c15ab2e8 100644 --- a/include/deal.II/dofs/dof_handler.h +++ b/include/deal.II/dofs/dof_handler.h @@ -701,36 +701,6 @@ public: void distribute_mg_dofs(); - /** - * Prepare all cells for p-adaptation in hp-mode. Does nothing in non-hp-mode. - * - * Essentially does to future FE indices what - * Triangulation::prepare_coarsening_and_refinement() does to refinement - * flags. - * - * In detail, this function limits the level difference of neighboring cells - * and thus smoothes the overall function space. Future FE indices will be - * raised (and never lowered) so that the level difference to neighboring - * cells is never larger than @p max_difference. - * - * Multiple FE hierarchies might have been registered via - * hp::FECollection::set_hierarchy(). This function operates on only one - * hierarchy, namely the one that contains the FE index @p contains_fe_index. - * Cells with future FE indices that are not part of the corresponding - * hierarchy will be ignored. - * - * The function can optionally be called before performing adaptation with - * Triangulation::execute_coarsening_and_refinement(). It is not necessary to - * call this function, nor will it be automatically invoked in any part of the - * library (contrary to its Triangulation counterpart). - * - * Returns whether any future FE indices have been changed by this function. - */ - bool - prepare_coarsening_and_refinement( - const unsigned int max_difference = 1, - const unsigned int contains_fe_index = 0) const; - /** * Returns whether this DoFHandler has hp-capabilities. */ diff --git a/include/deal.II/hp/refinement.h b/include/deal.II/hp/refinement.h index 5f763b888a..b674183cd0 100644 --- a/include/deal.II/hp/refinement.h +++ b/include/deal.II/hp/refinement.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2019 - 2020 by the deal.II authors +// Copyright (C) 2019 - 2021 by the deal.II authors // // This file is part of the deal.II library. // @@ -361,6 +361,7 @@ namespace hp const ComparisonFunction::type> &compare_refine, const ComparisonFunction::type> &compare_coarsen); + /** * @} */ @@ -647,6 +648,47 @@ namespace hp /** * @} */ + + /** + * @name Optimiize p-level distribution + * @{ + */ + + /** + * Limit p-level differences between neighboring cells. + * + * Essentially does to future FE indices what + * Triangulation::prepare_coarsening_and_refinement() does to refinement + * flags. + * + * In detail, this function limits the level difference of neighboring cells + * and thus smoothes the overall function space. Future FE indices will be + * raised (and never lowered) so that the level difference to neighboring + * cells is never larger than @p max_difference. + * + * Multiple FE hierarchies might have been registered via + * hp::FECollection::set_hierarchy(). This function operates on only one + * hierarchy, namely the one that contains the FE index @p contains_fe_index. + * Cells with future FE indices that are not part of the corresponding + * hierarchy will be ignored. + * + * The function can optionally be called before performing adaptation with + * Triangulation::execute_coarsening_and_refinement(). It is not necessary + * to call this function, nor will it be automatically invoked in any part + * of the library (contrary to its Triangulation counterpart). + * + * Returns whether any future FE indices have been changed by this function. + */ + template + bool + limit_p_level_difference( + const dealii::DoFHandler &dof_handler, + const unsigned int max_difference = 1, + const unsigned int contains_fe_index = 0); + + /** + * @} + */ } // namespace Refinement } // namespace hp diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index e7ba6ba8a2..a2f32baae5 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -32,8 +32,6 @@ #include #include -#include - #include #include #include @@ -2657,208 +2655,6 @@ DoFHandler::distribute_mg_dofs() -template -bool -DoFHandler::prepare_coarsening_and_refinement( - const unsigned int max_difference, - const unsigned int contains_fe_index) const -{ - Assert(max_difference > 0, - ExcMessage( - "This function does not serve any purpose for max_difference = 0.")); - AssertIndexRange(contains_fe_index, fe_collection.size()); - - if (!hp_capability_enabled) - // nothing to do - return false; - - - // - // establish hierarchy - // - // - create bimap between hierarchy levels and FE indices - - // there can be as many levels in the hierarchy as active FE indices are - // possible - using level_type = active_fe_index_type; - static const level_type invalid_level = invalid_active_fe_index; - - // map from FE index to level in hierarchy - // FE indices that are not covered in the hierarchy are not in the map - const std::vector fe_index_for_hierarchy_level = - fe_collection.get_hierarchy_sequence(contains_fe_index); - - // map from level in hierarchy to FE index - // FE indices that are not covered in the hierarchy will be mapped to - // invalid_level - std::vector hierarchy_level_for_fe_index(fe_collection.size(), - invalid_level); - for (unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l) - hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l; - - - // - // parallelization - // - // - create distributed vector of level indices - // - update ghost values in each iteration (see later) - // - no need to compress, since the owning processor will have the correct - // level index - - // HOTFIX: dealii::Vector does not accept integral types - LinearAlgebra::distributed::Vector future_levels; - if (const auto parallel_tria = - dynamic_cast *>( - &(*tria))) - { - future_levels.reinit( - parallel_tria->global_active_cell_index_partitioner().lock()); - } - else - { - future_levels.reinit(tria->n_active_cells()); - } - - for (const auto &cell : active_cell_iterators()) - if (cell->is_locally_owned()) - future_levels[cell->global_active_cell_index()] = - hierarchy_level_for_fe_index[cell->future_fe_index()]; - - - // - // limit level difference of neighboring cells - // - // - go over all locally relevant cells, and adjust the level indices of - // locally owned neighbors to match the level difference (as a consequence, - // indices on ghost cells will be updated only on the owning processor) - // - always raise levels to match criterion, never lower them - // - exchange level indices on ghost cells - - bool levels_changed = false; - bool levels_changed_in_cycle; - do - { - levels_changed_in_cycle = false; - - future_levels.update_ghost_values(); - - for (const auto &cell : active_cell_iterators()) - if (!cell->is_artificial()) - { - const level_type cell_level = static_cast( - future_levels[cell->global_active_cell_index()]); - - // ignore cells that are not part of the hierarchy - if (cell_level == invalid_level) - continue; - - // ignore lowest levels of the hierarchy that always fulfill the - // max_difference criterion - if (cell_level <= max_difference) - continue; - - for (unsigned int f = 0; f < cell->n_faces(); ++f) - if (cell->face(f)->at_boundary() == false) - { - if (cell->face(f)->has_children()) - { - for (unsigned int sf = 0; - sf < cell->face(f)->n_children(); - ++sf) - { - const auto neighbor = - cell->neighbor_child_on_subface(f, sf); - - // We only care about locally owned neighbors. If - // neighbor is a ghost cell, its future FE index will - // be updated on the owning process and communicated - // at the next loop iteration. - if (neighbor->is_locally_owned()) - { - const level_type neighbor_level = - static_cast( - future_levels - [neighbor->global_active_cell_index()]); - - // ignore neighbors that are not part of the - // hierarchy - if (neighbor_level == invalid_level) - continue; - - if ((cell_level - max_difference) > - neighbor_level) - { - future_levels - [neighbor->global_active_cell_index()] = - cell_level - max_difference; - - levels_changed_in_cycle = true; - } - } - } - } - else - { - const auto neighbor = cell->neighbor(f); - - // We only care about locally owned neighbors. If neighbor - // is a ghost cell, its future FE index will be updated on - // the owning process and communicated at the next loop - // iteration. - if (neighbor->is_locally_owned()) - { - const level_type neighbor_level = - static_cast( - future_levels[neighbor - ->global_active_cell_index()]); - - // ignore neighbors that are not part of the hierarchy - if (neighbor_level == invalid_level) - continue; - - if ((cell_level - max_difference) > neighbor_level) - { - future_levels[neighbor - ->global_active_cell_index()] = - cell_level - max_difference; - - levels_changed_in_cycle = true; - } - } - } - } - } - - levels_changed_in_cycle = - Utilities::MPI::logical_or(levels_changed_in_cycle, - tria->get_communicator()); - levels_changed |= levels_changed_in_cycle; - } - while (levels_changed_in_cycle); - - // update future FE indices on locally owned cells - for (const auto &cell : active_cell_iterators()) - if (cell->is_locally_owned()) - { - const level_type cell_level = static_cast( - future_levels[cell->global_active_cell_index()]); - - if (cell_level != invalid_level) - { - const active_fe_index_type fe_index = - fe_index_for_hierarchy_level[cell_level]; - - // only update if necessary - if (fe_index != cell->active_fe_index()) - cell->set_future_fe_index(fe_index); - } - } - - return levels_changed; -} - - - template void DoFHandler::initialize_local_block_info() diff --git a/source/hp/refinement.cc b/source/hp/refinement.cc index ac0549bd48..edc5cb1794 100644 --- a/source/hp/refinement.cc +++ b/source/hp/refinement.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2019 - 2020 by the deal.II authors +// Copyright (C) 2019 - 2021 by the deal.II authors // // This file is part of the deal.II library. // @@ -22,12 +22,14 @@ #include #include +#include + #include #include -#include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -752,6 +754,219 @@ namespace hp } } } + + + + /** + * Optimize p-level distribution + */ + template + bool + limit_p_level_difference( + const dealii::DoFHandler &dof_handler, + const unsigned int max_difference, + const unsigned int contains_fe_index) + { + Assert( + dof_handler.has_hp_capabilities(), + (typename dealii::DoFHandler::ExcOnlyAvailableWithHP())); + Assert( + max_difference > 0, + ExcMessage( + "This function does not serve any purpose for max_difference = 0.")); + AssertIndexRange(contains_fe_index, + dof_handler.get_fe_collection().size()); + + // + // establish hierarchy + // + // - create bimap between hierarchy levels and FE indices + + // there can be as many levels in the hierarchy as active FE indices are + // possible + using level_type = + typename dealii::DoFHandler::active_fe_index_type; + static const level_type invalid_level = + dealii::DoFHandler::invalid_active_fe_index; + + // map from FE index to level in hierarchy + // FE indices that are not covered in the hierarchy are not in the map + const std::vector fe_index_for_hierarchy_level = + dof_handler.get_fe_collection().get_hierarchy_sequence( + contains_fe_index); + + // map from level in hierarchy to FE index + // FE indices that are not covered in the hierarchy will be mapped to + // invalid_level + std::vector hierarchy_level_for_fe_index( + dof_handler.get_fe_collection().size(), invalid_level); + for (unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l) + hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l; + + + // + // parallelization + // + // - create distributed vector of level indices + // - update ghost values in each iteration (see later) + // - no need to compress, since the owning processor will have the correct + // level index + + // HOTFIX: dealii::Vector does not accept integral types + LinearAlgebra::distributed::Vector future_levels; + if (const auto parallel_tria = + dynamic_cast *>( + &(dof_handler.get_triangulation()))) + { + future_levels.reinit( + parallel_tria->global_active_cell_index_partitioner().lock()); + } + else + { + future_levels.reinit( + dof_handler.get_triangulation().n_active_cells()); + } + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + future_levels[cell->global_active_cell_index()] = + hierarchy_level_for_fe_index[cell->future_fe_index()]; + + + // + // limit level difference of neighboring cells + // + // - go over all locally relevant cells, and adjust the level indices of + // locally owned neighbors to match the level difference (as a + // consequence, indices on ghost cells will be updated only on the + // owning processor) + // - always raise levels to match criterion, never lower them + // - exchange level indices on ghost cells + + bool levels_changed = false; + bool levels_changed_in_cycle; + do + { + levels_changed_in_cycle = false; + + future_levels.update_ghost_values(); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (!cell->is_artificial()) + { + const level_type cell_level = static_cast( + future_levels[cell->global_active_cell_index()]); + + // ignore cells that are not part of the hierarchy + if (cell_level == invalid_level) + continue; + + // ignore lowest levels of the hierarchy that always fulfill the + // max_difference criterion + if (cell_level <= max_difference) + continue; + + for (unsigned int f = 0; f < cell->n_faces(); ++f) + if (cell->face(f)->at_boundary() == false) + { + if (cell->face(f)->has_children()) + { + for (unsigned int sf = 0; + sf < cell->face(f)->n_children(); + ++sf) + { + const auto neighbor = + cell->neighbor_child_on_subface(f, sf); + + // We only care about locally owned neighbors. If + // neighbor is a ghost cell, its future FE index + // will be updated on the owning process and + // communicated at the next loop iteration. + if (neighbor->is_locally_owned()) + { + const level_type neighbor_level = + static_cast( + future_levels + [neighbor->global_active_cell_index()]); + + // ignore neighbors that are not part of the + // hierarchy + if (neighbor_level == invalid_level) + continue; + + if ((cell_level - max_difference) > + neighbor_level) + { + future_levels + [neighbor->global_active_cell_index()] = + cell_level - max_difference; + + levels_changed_in_cycle = true; + } + } + } + } + else + { + const auto neighbor = cell->neighbor(f); + + // We only care about locally owned neighbors. If + // neighbor is a ghost cell, its future FE index will + // be updated on the owning process and communicated + // at the next loop iteration. + if (neighbor->is_locally_owned()) + { + const level_type neighbor_level = + static_cast( + future_levels + [neighbor->global_active_cell_index()]); + + // ignore neighbors that are not part of the + // hierarchy + if (neighbor_level == invalid_level) + continue; + + if ((cell_level - max_difference) > + neighbor_level) + { + future_levels + [neighbor->global_active_cell_index()] = + cell_level - max_difference; + + levels_changed_in_cycle = true; + } + } + } + } + } + + levels_changed_in_cycle = + Utilities::MPI::logical_or(levels_changed_in_cycle, + dof_handler.get_communicator()); + levels_changed |= levels_changed_in_cycle; + } + while (levels_changed_in_cycle); + + // update future FE indices on locally owned cells + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + const level_type cell_level = static_cast( + future_levels[cell->global_active_cell_index()]); + + if (cell_level != invalid_level) + { + const unsigned int fe_index = + fe_index_for_hierarchy_level[cell_level]; + + // only update if necessary + if (fe_index != cell->active_fe_index()) + cell->set_future_fe_index(fe_index); + } + } + + return levels_changed; + } } // namespace Refinement } // namespace hp diff --git a/source/hp/refinement.inst.in b/source/hp/refinement.inst.in index a9042f0c59..b2ec253d4e 100644 --- a/source/hp/refinement.inst.in +++ b/source/hp/refinement.inst.in @@ -42,6 +42,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) choose_p_over_h( const dealii::DoFHandler &); + + template bool + limit_p_level_difference( + const dealii::DoFHandler + &, + const unsigned int, + const unsigned int); \} \} #endif diff --git a/tests/hp/prepare_coarsening_and_refinement_01.cc b/tests/hp/prepare_coarsening_and_refinement_01.cc index 10d685e8ad..e0494dff5a 100644 --- a/tests/hp/prepare_coarsening_and_refinement_01.cc +++ b/tests/hp/prepare_coarsening_and_refinement_01.cc @@ -32,6 +32,7 @@ #include #include +#include #include @@ -74,7 +75,9 @@ test(const unsigned int fes_size, const unsigned int max_difference) center_cell->set_active_fe_index(sequence.back()); const bool fe_indices_changed = - dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index); + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); tria.execute_coarsening_and_refinement(); (void)fe_indices_changed; diff --git a/tests/hp/prepare_coarsening_and_refinement_02.cc b/tests/hp/prepare_coarsening_and_refinement_02.cc index f46d847f00..7038552e2c 100644 --- a/tests/hp/prepare_coarsening_and_refinement_02.cc +++ b/tests/hp/prepare_coarsening_and_refinement_02.cc @@ -32,6 +32,7 @@ #include #include +#include #include @@ -78,8 +79,9 @@ test(const unsigned int fes_size, const unsigned int max_difference) fes.next_in_hierarchy(center_cell->active_fe_index())); const bool fe_indices_changed = - dofh.prepare_coarsening_and_refinement(max_difference, - contains_fe_index); + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); tria.execute_coarsening_and_refinement(); (void)fe_indices_changed; diff --git a/tests/hp/prepare_coarsening_and_refinement_03.cc b/tests/hp/prepare_coarsening_and_refinement_03.cc index 806e656307..3d34890e66 100644 --- a/tests/hp/prepare_coarsening_and_refinement_03.cc +++ b/tests/hp/prepare_coarsening_and_refinement_03.cc @@ -28,6 +28,7 @@ #include #include +#include #include "../tests.h" @@ -90,7 +91,9 @@ test(const unsigned int fes_size, const unsigned int max_difference) dofh.distribute_dofs(fes); const bool fe_indices_changed = - dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index); + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); tria.execute_coarsening_and_refinement(); (void)fe_indices_changed; diff --git a/tests/mpi/prepare_coarsening_and_refinement_01.cc b/tests/mpi/prepare_coarsening_and_refinement_01.cc index 26851fc263..f9e36450ee 100644 --- a/tests/mpi/prepare_coarsening_and_refinement_01.cc +++ b/tests/mpi/prepare_coarsening_and_refinement_01.cc @@ -35,6 +35,7 @@ #include #include +#include #include @@ -79,7 +80,9 @@ test(parallel::TriangulationBase &tria, cell->set_active_fe_index(sequence.back()); const bool fe_indices_changed = - dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index); + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); tria.execute_coarsening_and_refinement(); (void)fe_indices_changed; diff --git a/tests/mpi/prepare_coarsening_and_refinement_02.cc b/tests/mpi/prepare_coarsening_and_refinement_02.cc index e324e09652..e660e42fe8 100644 --- a/tests/mpi/prepare_coarsening_and_refinement_02.cc +++ b/tests/mpi/prepare_coarsening_and_refinement_02.cc @@ -35,6 +35,7 @@ #include #include +#include #include @@ -83,8 +84,9 @@ test(parallel::TriangulationBase &tria, fes.next_in_hierarchy(cell->active_fe_index())); const bool fe_indices_changed = - dofh.prepare_coarsening_and_refinement(max_difference, - contains_fe_index); + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); tria.execute_coarsening_and_refinement(); (void)fe_indices_changed; diff --git a/tests/mpi/prepare_coarsening_and_refinement_03.cc b/tests/mpi/prepare_coarsening_and_refinement_03.cc index 0ab5d6f59f..ece1b8681d 100644 --- a/tests/mpi/prepare_coarsening_and_refinement_03.cc +++ b/tests/mpi/prepare_coarsening_and_refinement_03.cc @@ -30,6 +30,7 @@ #include #include +#include #include "../tests.h" @@ -94,7 +95,9 @@ test(parallel::TriangulationBase &tria, dofh.distribute_dofs(fes); const bool fe_indices_changed = - dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index); + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); tria.execute_coarsening_and_refinement(); (void)fe_indices_changed; -- 2.39.5