From 9c59501d722b01ad19bb2ddc96e68ab48fa3332d Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Mon, 29 Mar 2021 14:09:55 -0600 Subject: [PATCH] Fix hp-coarsening on p:s:Triangulation. --- include/deal.II/dofs/dof_handler.h | 21 +++ source/dofs/dof_handler.cc | 160 +++++++++++++++--- source/dofs/dof_handler.inst.in | 34 +++- .../sharedtria/limit_p_level_difference_04.cc | 126 ++++++++++++++ ...ifference_04.with_mpi=true.mpirun=1.output | 75 ++++++++ ...ifference_04.with_mpi=true.mpirun=4.output | 141 +++++++++++++++ 6 files changed, 527 insertions(+), 30 deletions(-) create mode 100644 tests/sharedtria/limit_p_level_difference_04.cc create mode 100644 tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=1.output create mode 100644 tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=4.output diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h index 35a7c578a1..0678b01006 100644 --- a/include/deal.II/dofs/dof_handler.h +++ b/include/deal.II/dofs/dof_handler.h @@ -1820,6 +1820,21 @@ namespace internal { namespace DoFHandlerImplementation { + /** + * Given a DoFHandler object in hp-mode, make sure that the + * future FE indices that a user has set for locally owned cells are + * communicated to all other relevant cells as well. + * + * For parallel::shared::Triangulation objects, + * this information is distributed on both ghost and artificial cells. + * + * In case a parallel::distributed::Triangulation is used, + * indices are communicated only to ghost cells. + */ + template + void + communicate_future_fe_indices(DoFHandler &dof_handler); + /** * Return the index of the finite element from the entire hp::FECollection * that is dominated by those assigned as future finite elements to the @@ -1837,6 +1852,12 @@ namespace internal * * @note This function can only be called on direct parent cells, i.e., * non-active cells whose children are all active. + * + * @note On parallel::shared::Triangulation objects where sibling cells + * can be ghost cells, make sure that future FE indices have been properly + * communicated with communicate_future_fe_indices() first. Otherwise, + * results might differ on different processors. There is no check for + * consistency of future FE indices. */ template unsigned int diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 9cffc4ccc9..b22490e765 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -1809,6 +1809,94 @@ namespace internal + /** + * Same as above, but for future FE indices. + * + * Given a DoFHandler object in hp-mode, make sure that the + * future FE indices that a user has set for locally owned cells are + * communicated to all other relevant cells as well. + * + * For parallel::shared::Triangulation objects, + * this information is distributed on both ghost and artificial cells. + * + * In case a parallel::distributed::Triangulation is used, + * indices are communicated only to ghost cells. + */ + template + static void + communicate_future_fe_indices(DoFHandler &dof_handler) + { + Assert( + dof_handler.hp_capability_enabled == true, + (typename DoFHandler::ExcOnlyAvailableWithHP())); + + using active_fe_index_type = + typename dealii::DoFHandler::active_fe_index_type; + + if (const dealii::parallel::shared::Triangulation *tr = + dynamic_cast< + const dealii::parallel::shared::Triangulation + *>(&dof_handler.get_triangulation())) + { + std::vector future_fe_indices( + tr->n_active_cells(), 0u); + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + future_fe_indices[cell->active_cell_index()] = + dof_handler + .hp_cell_future_fe_indices[cell->level()][cell->index()]; + + Utilities::MPI::sum(future_fe_indices, + tr->get_communicator(), + future_fe_indices); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (!cell->is_locally_owned()) + dof_handler + .hp_cell_future_fe_indices[cell->level()][cell->index()] = + future_fe_indices[cell->active_cell_index()]; + } + else if (const dealii::parallel:: + DistributedTriangulationBase *tr = + dynamic_cast< + const dealii::parallel:: + DistributedTriangulationBase *>( + &dof_handler.get_triangulation())) + { + auto pack = + [&dof_handler]( + const typename dealii::DoFHandler:: + active_cell_iterator &cell) -> active_fe_index_type { + return dof_handler + .hp_cell_future_fe_indices[cell->level()][cell->index()]; + }; + + auto unpack = + [&dof_handler]( + const typename dealii::DoFHandler:: + active_cell_iterator & cell, + const active_fe_index_type future_fe_index) -> void { + dof_handler + .hp_cell_future_fe_indices[cell->level()][cell->index()] = + future_fe_index; + }; + + GridTools::exchange_cell_data_to_ghosts< + active_fe_index_type, + dealii::DoFHandler>(dof_handler, pack, unpack); + } + else + { + Assert( + (dynamic_cast< + const dealii::parallel::TriangulationBase *>( + &dof_handler.get_triangulation()) == nullptr), + ExcInternalError()); + } + } + + + /** * Collect all finite element indices on cells that will be affected by * future refinement and coarsening. Further, prepare those indices to @@ -1873,8 +1961,8 @@ namespace internal dim>::ExcInconsistentCoarseningFlags()); #endif - const unsigned int fe_index = - dealii::internal::hp::DoFHandlerImplementation:: + const unsigned int fe_index = dealii::internal::hp:: + DoFHandlerImplementation::Implementation:: dominated_future_fe_on_children( parent); @@ -1926,11 +2014,11 @@ namespace internal const auto &parent = refine.first; for (const auto &child : parent->child_iterators()) - { - Assert(child->is_locally_owned() && child->is_active(), - ExcInternalError()); - child->set_active_fe_index(refine.second); - } + if (child->is_locally_owned()) + { + Assert(child->is_active(), ExcInternalError()); + child->set_active_fe_index(refine.second); + } } // Set active FE indices on coarsened cells that have been determined @@ -1938,9 +2026,12 @@ namespace internal for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index) { const auto &cell = coarsen.first; - Assert(cell->is_locally_owned() && cell->is_active(), - ExcInternalError()); - cell->set_active_fe_index(coarsen.second); + + if (cell->is_locally_owned()) + { + Assert(cell->is_active(), ExcInternalError()); + cell->set_active_fe_index(coarsen.second); + } } } @@ -1982,8 +2073,7 @@ namespace internal /** * Return the index of the finite element from the entire * hp::FECollection that is dominated by those assigned as future finite - * elements to the - * children of @p parent. + * elements to the children of @p parent. * * See documentation in the header file for more information. */ @@ -1998,6 +2088,11 @@ namespace internal "You ask for information on children of this cell which is only " "available for active cells. This cell has no children.")); + const auto &dof_handler = parent->get_dof_handler(); + Assert( + dof_handler.has_hp_capabilities(), + (typename DoFHandler::ExcOnlyAvailableWithHP())); + std::set future_fe_indices_children; for (const auto &child : parent->child_iterators()) { @@ -2006,18 +2101,23 @@ namespace internal ExcMessage( "You ask for information on children of this cell which is only " "available for active cells. One of its children is not active.")); - Assert( - child->is_locally_owned(), - ExcMessage( - "You ask for information on children of this cell which is only " - "available for locally owned cells. One of its children is not " - "locally owned.")); - future_fe_indices_children.insert(child->future_fe_index()); + + // Ghost siblings might occur on parallel::shared::Triangulation + // objects. The public interface does not allow to access future + // FE indices on ghost cells. However, we need this information + // here and thus call the internal function that does not check + // for cell ownership. This requires that future FE indices have + // been communicated prior to calling this function. + const unsigned int future_fe_index_child = + dealii::internal::DoFCellAccessorImplementation:: + Implementation::future_fe_index(*child); + + future_fe_indices_children.insert(future_fe_index_child); } Assert(!future_fe_indices_children.empty(), ExcInternalError()); const unsigned int future_fe_index = - parent->get_dof_handler().fe_collection.find_dominated_fe_extended( + dof_handler.fe_collection.find_dominated_fe_extended( future_fe_indices_children, /*codim=*/0); @@ -2031,7 +2131,20 @@ namespace internal /** - * Public wrapper for the above function. + * Public wrapper of Implementation::communicate_future_fe_indices(). + */ + template + void + communicate_future_fe_indices(DoFHandler &dof_handler) + { + Implementation::communicate_future_fe_indices( + dof_handler); + } + + + + /** + * Public wrapper of Implementation::dominated_future_fe_on_children(). */ template unsigned int @@ -3159,6 +3272,11 @@ DoFHandler::connect_to_triangulation_signals() })); // refinement signals + this->tria_listeners_for_transfer.push_back( + this->tria->signals.pre_refinement.connect([this]() { + internal::hp::DoFHandlerImplementation::Implementation:: + communicate_future_fe_indices(*this); + })); this->tria_listeners_for_transfer.push_back( this->tria->signals.pre_refinement.connect( [this] { this->pre_transfer_action(); })); diff --git a/source/dofs/dof_handler.inst.in b/source/dofs/dof_handler.inst.in index b6fb6cb840..e78524dd31 100644 --- a/source/dofs/dof_handler.inst.in +++ b/source/dofs/dof_handler.inst.in @@ -19,15 +19,31 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) #if deal_II_dimension <= deal_II_space_dimension template class DoFHandler; - template std::string internal::policy_to_string( - const dealii::internal::DoFHandlerImplementation::Policy:: - PolicyBase &policy); - - template unsigned int internal::hp::DoFHandlerImplementation:: - dominated_future_fe_on_children( - const typename DoFHandler::cell_iterator &); + namespace internal + \{ + template std::string + policy_to_string(const DoFHandlerImplementation::Policy::PolicyBase< + deal_II_dimension, + deal_II_space_dimension> &); + + namespace hp + \{ + namespace DoFHandlerImplementation + \{ + template void + communicate_future_fe_indices( + DoFHandler &); + + template unsigned int + dominated_future_fe_on_children( + const typename DoFHandler::cell_iterator + &); + \} + \} + \} #endif } diff --git a/tests/sharedtria/limit_p_level_difference_04.cc b/tests/sharedtria/limit_p_level_difference_04.cc new file mode 100644 index 0000000000..d2fc146234 --- /dev/null +++ b/tests/sharedtria/limit_p_level_difference_04.cc @@ -0,0 +1,126 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// verify restrictions on level differences imposed by +// DoFHandler::prepare_coarsening_and_refinement() +// on h-coarsened grids + + +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + + +template +void +test(const unsigned int max_difference, const bool allow_artificial_cells) +{ + Assert(max_difference > 0, ExcInternalError()); + + // setup FE collection + const unsigned int fes_size = 2 * max_difference + 2; + + hp::FECollection fes; + while (fes.size() < fes_size) + fes.push_back(FE_Q(1)); + + const unsigned int contains_fe_index = 0; + const auto sequence = fes.get_hierarchy_sequence(contains_fe_index); + + // set up line grid + // - refine central cell and flag them for coarsening + // - assign highest p-level to leftmost cell + // + // +-------+---+---+-------+ + // | | c | c | | + // | +---+---+ | + // | | c | c | | + // +-------+---+---+-------+ + // + // after prepare_coarsening_and_refinement(), the p-levels should propagate + // through the central cells as if they were already coarsened + + parallel::shared::Triangulation tria(MPI_COMM_WORLD, + Triangulation::none, + allow_artificial_cells); + TestGrids::hyper_line(tria, 3); + + // refine the central cell + for (const auto &cell : tria.active_cell_iterators()) + if (cell->center()[0] > 1. && cell->center()[0] < 2.) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + // now flag these cells for coarsening + for (const auto &cell : tria.active_cell_iterators()) + if (cell->center()[0] > 1. && cell->center()[0] < 2.) + cell->set_coarsen_flag(); + + DoFHandler dofh(tria); + for (const auto &cell : dofh.active_cell_iterators()) + if (cell->is_locally_owned() && cell->center()[0] < 1.) + cell->set_active_fe_index(sequence.back()); + dofh.distribute_dofs(fes); + + const bool fe_indices_changed = + hp::Refinement::limit_p_level_difference(dofh, + max_difference, + contains_fe_index); + (void)fe_indices_changed; + Assert(fe_indices_changed, ExcInternalError()); + + deallog << "future FE indices before adaptation:" << std::endl; + for (const auto &cell : dofh.active_cell_iterators()) + if (cell->is_locally_owned()) + deallog << " " << cell->id().to_string() << " " << cell->future_fe_index() + << std::endl; + + tria.execute_coarsening_and_refinement(); + + deallog << "active FE indices after adaptation:" << std::endl; + for (const auto &cell : dofh.active_cell_iterators()) + if (cell->is_locally_owned()) + deallog << " " << cell->id().to_string() << " " << cell->active_fe_index() + << std::endl; + + deallog << "OK" << std::endl; +} + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + deallog << std::boolalpha; + for (const bool allow_artificial_cells : {false, true}) + { + deallog << "artificial cells: " << allow_artificial_cells << std::endl; + test<2>(1, allow_artificial_cells); + test<2>(2, allow_artificial_cells); + test<2>(3, allow_artificial_cells); + } +} diff --git a/tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=1.output b/tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=1.output new file mode 100644 index 0000000000..a44dd666d8 --- /dev/null +++ b/tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=1.output @@ -0,0 +1,75 @@ + +DEAL:0::artificial cells: false +DEAL:0::future FE indices before adaptation: +DEAL:0:: 0_0: 3 +DEAL:0:: 2_0: 1 +DEAL:0:: 1_1:0 2 +DEAL:0:: 1_1:1 2 +DEAL:0:: 1_1:2 2 +DEAL:0:: 1_1:3 2 +DEAL:0::active FE indices after adaptation: +DEAL:0:: 0_0: 3 +DEAL:0:: 1_0: 2 +DEAL:0:: 2_0: 1 +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 0_0: 5 +DEAL:0:: 2_0: 1 +DEAL:0:: 1_1:0 3 +DEAL:0:: 1_1:1 3 +DEAL:0:: 1_1:2 3 +DEAL:0:: 1_1:3 3 +DEAL:0::active FE indices after adaptation: +DEAL:0:: 0_0: 5 +DEAL:0:: 1_0: 3 +DEAL:0:: 2_0: 1 +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 0_0: 7 +DEAL:0:: 2_0: 1 +DEAL:0:: 1_1:0 4 +DEAL:0:: 1_1:1 4 +DEAL:0:: 1_1:2 4 +DEAL:0:: 1_1:3 4 +DEAL:0::active FE indices after adaptation: +DEAL:0:: 0_0: 7 +DEAL:0:: 1_0: 4 +DEAL:0:: 2_0: 1 +DEAL:0::OK +DEAL:0::artificial cells: true +DEAL:0::future FE indices before adaptation: +DEAL:0:: 0_0: 3 +DEAL:0:: 2_0: 1 +DEAL:0:: 1_1:0 2 +DEAL:0:: 1_1:1 2 +DEAL:0:: 1_1:2 2 +DEAL:0:: 1_1:3 2 +DEAL:0::active FE indices after adaptation: +DEAL:0:: 0_0: 3 +DEAL:0:: 1_0: 2 +DEAL:0:: 2_0: 1 +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 0_0: 5 +DEAL:0:: 2_0: 1 +DEAL:0:: 1_1:0 3 +DEAL:0:: 1_1:1 3 +DEAL:0:: 1_1:2 3 +DEAL:0:: 1_1:3 3 +DEAL:0::active FE indices after adaptation: +DEAL:0:: 0_0: 5 +DEAL:0:: 1_0: 3 +DEAL:0:: 2_0: 1 +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 0_0: 7 +DEAL:0:: 2_0: 1 +DEAL:0:: 1_1:0 4 +DEAL:0:: 1_1:1 4 +DEAL:0:: 1_1:2 4 +DEAL:0:: 1_1:3 4 +DEAL:0::active FE indices after adaptation: +DEAL:0:: 0_0: 7 +DEAL:0:: 1_0: 4 +DEAL:0:: 2_0: 1 +DEAL:0::OK diff --git a/tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=4.output b/tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=4.output new file mode 100644 index 0000000000..e571804fa0 --- /dev/null +++ b/tests/sharedtria/limit_p_level_difference_04.with_mpi=true.mpirun=4.output @@ -0,0 +1,141 @@ + +DEAL:0::artificial cells: false +DEAL:0::future FE indices before adaptation: +DEAL:0:: 1_1:3 2 +DEAL:0::active FE indices after adaptation: +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 1_1:3 3 +DEAL:0::active FE indices after adaptation: +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 1_1:3 4 +DEAL:0::active FE indices after adaptation: +DEAL:0::OK +DEAL:0::artificial cells: true +DEAL:0::future FE indices before adaptation: +DEAL:0:: 1_1:3 2 +DEAL:0::active FE indices after adaptation: +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 1_1:3 3 +DEAL:0::active FE indices after adaptation: +DEAL:0::OK +DEAL:0::future FE indices before adaptation: +DEAL:0:: 1_1:3 4 +DEAL:0::active FE indices after adaptation: +DEAL:0::OK + +DEAL:1::artificial cells: false +DEAL:1::future FE indices before adaptation: +DEAL:1:: 2_0: 1 +DEAL:1:: 1_1:1 2 +DEAL:1::active FE indices after adaptation: +DEAL:1:: 1_0: 2 +DEAL:1::OK +DEAL:1::future FE indices before adaptation: +DEAL:1:: 2_0: 1 +DEAL:1:: 1_1:1 3 +DEAL:1::active FE indices after adaptation: +DEAL:1:: 1_0: 3 +DEAL:1::OK +DEAL:1::future FE indices before adaptation: +DEAL:1:: 2_0: 1 +DEAL:1:: 1_1:1 4 +DEAL:1::active FE indices after adaptation: +DEAL:1:: 1_0: 4 +DEAL:1::OK +DEAL:1::artificial cells: true +DEAL:1::future FE indices before adaptation: +DEAL:1:: 2_0: 1 +DEAL:1:: 1_1:1 2 +DEAL:1::active FE indices after adaptation: +DEAL:1:: 1_0: 2 +DEAL:1::OK +DEAL:1::future FE indices before adaptation: +DEAL:1:: 2_0: 1 +DEAL:1:: 1_1:1 3 +DEAL:1::active FE indices after adaptation: +DEAL:1:: 1_0: 3 +DEAL:1::OK +DEAL:1::future FE indices before adaptation: +DEAL:1:: 2_0: 1 +DEAL:1:: 1_1:1 4 +DEAL:1::active FE indices after adaptation: +DEAL:1:: 1_0: 4 +DEAL:1::OK + + +DEAL:2::artificial cells: false +DEAL:2::future FE indices before adaptation: +DEAL:2:: 1_1:2 2 +DEAL:2::active FE indices after adaptation: +DEAL:2:: 2_0: 1 +DEAL:2::OK +DEAL:2::future FE indices before adaptation: +DEAL:2:: 1_1:2 3 +DEAL:2::active FE indices after adaptation: +DEAL:2:: 2_0: 1 +DEAL:2::OK +DEAL:2::future FE indices before adaptation: +DEAL:2:: 1_1:2 4 +DEAL:2::active FE indices after adaptation: +DEAL:2:: 2_0: 1 +DEAL:2::OK +DEAL:2::artificial cells: true +DEAL:2::future FE indices before adaptation: +DEAL:2:: 1_1:2 2 +DEAL:2::active FE indices after adaptation: +DEAL:2:: 2_0: 1 +DEAL:2::OK +DEAL:2::future FE indices before adaptation: +DEAL:2:: 1_1:2 3 +DEAL:2::active FE indices after adaptation: +DEAL:2:: 2_0: 1 +DEAL:2::OK +DEAL:2::future FE indices before adaptation: +DEAL:2:: 1_1:2 4 +DEAL:2::active FE indices after adaptation: +DEAL:2:: 2_0: 1 +DEAL:2::OK + + +DEAL:3::artificial cells: false +DEAL:3::future FE indices before adaptation: +DEAL:3:: 0_0: 3 +DEAL:3:: 1_1:0 2 +DEAL:3::active FE indices after adaptation: +DEAL:3:: 0_0: 3 +DEAL:3::OK +DEAL:3::future FE indices before adaptation: +DEAL:3:: 0_0: 5 +DEAL:3:: 1_1:0 3 +DEAL:3::active FE indices after adaptation: +DEAL:3:: 0_0: 5 +DEAL:3::OK +DEAL:3::future FE indices before adaptation: +DEAL:3:: 0_0: 7 +DEAL:3:: 1_1:0 4 +DEAL:3::active FE indices after adaptation: +DEAL:3:: 0_0: 7 +DEAL:3::OK +DEAL:3::artificial cells: true +DEAL:3::future FE indices before adaptation: +DEAL:3:: 0_0: 3 +DEAL:3:: 1_1:0 2 +DEAL:3::active FE indices after adaptation: +DEAL:3:: 0_0: 3 +DEAL:3::OK +DEAL:3::future FE indices before adaptation: +DEAL:3:: 0_0: 5 +DEAL:3:: 1_1:0 3 +DEAL:3::active FE indices after adaptation: +DEAL:3:: 0_0: 5 +DEAL:3::OK +DEAL:3::future FE indices before adaptation: +DEAL:3:: 0_0: 7 +DEAL:3:: 1_1:0 4 +DEAL:3::active FE indices after adaptation: +DEAL:3:: 0_0: 7 +DEAL:3::OK + -- 2.39.5