From: Marc Fehling Date: Wed, 3 Mar 2021 21:13:10 +0000 (-0700) Subject: Added DoFHandler::prepare_for_coarsening_and_refinement(). X-Git-Tag: v9.3.0-rc1~380^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=baea52ad47cc93c5334b926cfd78381efab54993;p=dealii.git Added DoFHandler::prepare_for_coarsening_and_refinement(). --- diff --git a/doc/news/changes/minor/20210223Fehling b/doc/news/changes/minor/20210223Fehling new file mode 100644 index 0000000000..5b57fa6d39 --- /dev/null +++ b/doc/news/changes/minor/20210223Fehling @@ -0,0 +1,5 @@ +New: Member function DoFHandler::prepare_coarsening_and_refinement() +restricts the maximal level difference of neighboring cells with respect +to the finite element hierarchy of the registered hp::FECollection. +
+(Marc Fehling, 2021/02/23) diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h index d3dd609226..df5a1ec44b 100644 --- a/include/deal.II/dofs/dof_handler.h +++ b/include/deal.II/dofs/dof_handler.h @@ -701,6 +701,38 @@ 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. + * + * @note Only implemented for serial applications. + */ + 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 a8ae7ebdec..5f763b888a 100644 --- a/include/deal.II/hp/refinement.h +++ b/include/deal.II/hp/refinement.h @@ -147,9 +147,10 @@ namespace hp * Each cell flagged for h-refinement will also be flagged for p-refinement. * The same applies to coarsening. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -164,9 +165,10 @@ namespace hp * Each entry of the parameter @p p_flags needs to correspond to an active * cell. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -194,9 +196,10 @@ namespace hp * Each entry of the parameter @p criteria needs to correspond to an active * cell. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -236,9 +239,10 @@ namespace hp * cell. Parameters @p p_refine_fraction and @p p_coarsen_fraction need to be * in the interval $[0,1]$. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -279,9 +283,10 @@ namespace hp * cell. Parameters @p p_refine_fraction and @p p_coarsen_fraction need to be * in the interval $[0,1]$. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -316,9 +321,10 @@ namespace hp * * For more theoretical details see @cite ainsworth1998hp . * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -341,9 +347,10 @@ namespace hp * Each entry of the parameters @p criteria and @p references needs to * correspond to an active cell. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -553,7 +560,8 @@ namespace hp * * @note We want to predict the error by how adaptation will actually happen. * Thus, this function needs to be called after - * Triangulation::prepare_coarsening_and_refinement(). + * Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement(). */ template void @@ -579,9 +587,10 @@ namespace hp * Removes all refine and coarsen flags on cells that have a * @p future_fe_index assigned. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void @@ -626,9 +635,10 @@ namespace hp * the decision that Triangulation::prepare_coarsening_and_refinement() * would have made later on. * - * @note Triangulation::prepare_coarsening_and_refinement() may change - * refine and coarsen flags. Avoid calling it before this particular - * function. + * @note Triangulation::prepare_coarsening_and_refinement() and + * DoFHandler::prepare_coarsening_and_refinement() may change + * refine and coarsen flags as well as future finite element indices. + * Avoid calling them before this particular function. */ template void diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 5e440212fe..7f9aa9839e 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -2652,6 +2652,119 @@ 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; + + // communication of future fe indices on ghost cells necessary. + // thus, currently only implemented for serial triangulations. + Assert((dynamic_cast *>( + &(*tria)) == nullptr), + ExcNotImplemented()); + + // + // establish hierarchy + // + // - create bimap between hierarchy levels and FE indices + // - classify cells based on their level + + // 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 are not in the map + std::map hierarchy_level_for_fe_index; + for (unsigned int i = 0; i < fe_index_for_hierarchy_level.size(); ++i) + hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[i]] = i; + + // map from level in hierarchy to cells on that particular level + // cells that are not covered in the hierarchy are not in the map + std::vector> cells_on_hierarchy_level( + fe_index_for_hierarchy_level.size()); + if (fe_index_for_hierarchy_level.size() == fe_collection.size()) + { + // all FE indices are part of hierarchy + for (const auto &cell : active_cell_iterators()) + cells_on_hierarchy_level + [hierarchy_level_for_fe_index[cell->future_fe_index()]] + .push_back(cell); + } + else + { + // only some FE indices are part of hierarchy + typename decltype(hierarchy_level_for_fe_index)::iterator iter; + for (const auto &cell : active_cell_iterators()) + { + iter = hierarchy_level_for_fe_index.find(cell->future_fe_index()); + if (iter != hierarchy_level_for_fe_index.end()) + cells_on_hierarchy_level[iter->second].push_back(cell); + } + } + + // + // limit level difference of neighboring cells + // + // - always raise levels to match criterion, never lower them + // (hence iterating from highest to lowest level) + // - update future FE indices + // - move cells in level representation correspondigly + + bool changed_fe_indices = false; + for (unsigned int level = fe_index_for_hierarchy_level.size() - 1; + level > max_difference; + --level) + for (const auto &cell : cells_on_hierarchy_level[level]) + for (unsigned int f = 0; f < cell->n_faces(); ++f) + if (cell->face(f)->at_boundary() == false) + { + const auto neighbor = cell->neighbor(f); + + const auto neighbor_fe_and_level = + hierarchy_level_for_fe_index.find(neighbor->future_fe_index()); + + // ignore neighbors that are not part of the hierarchy + if (neighbor_fe_and_level == hierarchy_level_for_fe_index.end()) + continue; + + const auto neighbor_level = neighbor_fe_and_level->second; + + if ((level - max_difference) > neighbor_level) + { + // remove neighbor from old level representation + auto & old_cells = cells_on_hierarchy_level[neighbor_level]; + const auto iter = + std::find(old_cells.begin(), old_cells.end(), neighbor); + old_cells.erase(iter); + + // move neighbor to new level representation + const unsigned int new_level = level - max_difference; + cells_on_hierarchy_level[new_level].push_back(neighbor); + + // update future FE index + neighbor->set_future_fe_index( + fe_index_for_hierarchy_level[new_level]); + changed_fe_indices = true; + } + } + + return changed_fe_indices; +} + + + template void DoFHandler::initialize_local_block_info() diff --git a/tests/hp/prepare_coarsening_and_refinement_01.cc b/tests/hp/prepare_coarsening_and_refinement_01.cc new file mode 100644 index 0000000000..a702998be5 --- /dev/null +++ b/tests/hp/prepare_coarsening_and_refinement_01.cc @@ -0,0 +1,111 @@ +// --------------------------------------------------------------------- +// +// 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() +// +// set the center cell to the highest p-level in a hyper_cross geometry +// and verify that all other cells comply to the level difference + + +#include +#include + +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + +template +void +test(const unsigned int fes_size, const unsigned int max_difference) +{ + Assert(fes_size > 0, ExcInternalError()); + Assert(max_difference > 0, ExcInternalError()); + + // setup FE collection + 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); + + // setup cross-shaped mesh + Triangulation tria; + { + std::vector sizes(Utilities::pow(2, dim), + static_cast( + (sequence.size() - 1) / max_difference)); + GridGenerator::hyper_cross(tria, sizes); + } + + deallog << "ncells:" << tria.n_cells() << ", nfes:" << fes.size() << std::endl + << "sequence:" << sequence << std::endl; + + DoFHandler dofh(tria); + dofh.distribute_dofs(fes); + + // set center to highest p-level + const auto center_cell = dofh.begin(); + Assert(center_cell->center() == Point(), ExcInternalError()); + + center_cell->set_active_fe_index(sequence.back()); + dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index); + tria.execute_coarsening_and_refinement(); + + // display number of cells for each FE index + std::vector count(fes.size(), 0); + for (const auto &cell : dofh.active_cell_iterators()) + count[cell->active_fe_index()]++; + deallog << "fe count:" << count << std::endl; + +#ifdef DEBUG + // check each cell's active FE index by its distance from the center + for (const auto &cell : dofh.active_cell_iterators()) + { + const double distance = cell->center().distance(center_cell->center()); + const unsigned int expected_level = + (sequence.size() - 1) - + max_difference * static_cast(std::round(distance)); + + Assert(cell->active_fe_index() == sequence[expected_level], + ExcInternalError()); + } +#endif + + deallog << "OK" << std::endl; +} + + +int +main() +{ + initlog(); + + test<2>(4, 1); + test<2>(8, 2); + test<2>(12, 3); +} diff --git a/tests/hp/prepare_coarsening_and_refinement_01.output b/tests/hp/prepare_coarsening_and_refinement_01.output new file mode 100644 index 0000000000..a7bf63f5df --- /dev/null +++ b/tests/hp/prepare_coarsening_and_refinement_01.output @@ -0,0 +1,13 @@ + +DEAL::ncells:13, nfes:4 +DEAL::sequence:0 1 2 3 +DEAL::fe count:4 4 4 1 +DEAL::OK +DEAL::ncells:13, nfes:8 +DEAL::sequence:0 1 2 3 4 5 6 7 +DEAL::fe count:0 4 0 4 0 4 0 1 +DEAL::OK +DEAL::ncells:13, nfes:12 +DEAL::sequence:0 1 2 3 4 5 6 7 8 9 10 11 +DEAL::fe count:0 0 4 0 0 4 0 0 4 0 0 1 +DEAL::OK diff --git a/tests/hp/prepare_coarsening_and_refinement_02.cc b/tests/hp/prepare_coarsening_and_refinement_02.cc new file mode 100644 index 0000000000..42057d8603 --- /dev/null +++ b/tests/hp/prepare_coarsening_and_refinement_02.cc @@ -0,0 +1,118 @@ +// --------------------------------------------------------------------- +// +// 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() +// +// sequentially increase the p-level of the center cell in a hyper_cross +// geometry and verify that all other comply to the level difference + + +#include +#include + +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + +template +void +test(const unsigned int fes_size, const unsigned int max_difference) +{ + Assert(fes_size > 0, ExcInternalError()); + Assert(max_difference > 0, ExcInternalError()); + + // setup FE collection + 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); + + // setup cross-shaped mesh + Triangulation tria; + { + std::vector sizes(Utilities::pow(2, dim), + static_cast( + (sequence.size() - 1) / max_difference)); + GridGenerator::hyper_cross(tria, sizes); + } + + deallog << "ncells:" << tria.n_cells() << ", nfes:" << fes.size() << std::endl + << "sequence:" << sequence << std::endl; + + DoFHandler dofh(tria); + dofh.distribute_dofs(fes); + + // increase p-level in center subsequently + const auto center_cell = dofh.begin(); + Assert(center_cell->center() == Point(), ExcInternalError()); + + for (unsigned int i = 0; i < sequence.size() - 1; ++i) + { + center_cell->set_future_fe_index( + fes.next_in_hierarchy(center_cell->active_fe_index())); + + dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index); + + tria.execute_coarsening_and_refinement(); + + // display number of cells for each FE index + std::vector count(fes.size(), 0); + for (const auto &cell : dofh.active_cell_iterators()) + count[cell->active_fe_index()]++; + deallog << "cycle:" << i << ", fe count:" << count << std::endl; + } + Assert(center_cell->active_fe_index() == sequence.back(), ExcInternalError()); + +#ifdef DEBUG + // check each cell's active FE index by its distance from the center + for (const auto &cell : dofh.active_cell_iterators()) + { + const double distance = cell->center().distance(center_cell->center()); + const unsigned int expected_level = + (sequence.size() - 1) - + max_difference * static_cast(std::round(distance)); + + Assert(cell->active_fe_index() == sequence[expected_level], + ExcInternalError()); + } +#endif + + deallog << "OK" << std::endl; +} + + +int +main() +{ + initlog(); + + test<2>(4, 1); + test<2>(8, 2); + test<2>(12, 3); +} diff --git a/tests/hp/prepare_coarsening_and_refinement_02.output b/tests/hp/prepare_coarsening_and_refinement_02.output new file mode 100644 index 0000000000..1c07c11cb5 --- /dev/null +++ b/tests/hp/prepare_coarsening_and_refinement_02.output @@ -0,0 +1,31 @@ + +DEAL::ncells:13, nfes:4 +DEAL::sequence:0 1 2 3 +DEAL::cycle:0, fe count:12 1 0 0 +DEAL::cycle:1, fe count:8 4 1 0 +DEAL::cycle:2, fe count:4 4 4 1 +DEAL::OK +DEAL::ncells:13, nfes:8 +DEAL::sequence:0 1 2 3 4 5 6 7 +DEAL::cycle:0, fe count:12 1 0 0 0 0 0 0 +DEAL::cycle:1, fe count:12 0 1 0 0 0 0 0 +DEAL::cycle:2, fe count:8 4 0 1 0 0 0 0 +DEAL::cycle:3, fe count:8 0 4 0 1 0 0 0 +DEAL::cycle:4, fe count:4 4 0 4 0 1 0 0 +DEAL::cycle:5, fe count:4 0 4 0 4 0 1 0 +DEAL::cycle:6, fe count:0 4 0 4 0 4 0 1 +DEAL::OK +DEAL::ncells:13, nfes:12 +DEAL::sequence:0 1 2 3 4 5 6 7 8 9 10 11 +DEAL::cycle:0, fe count:12 1 0 0 0 0 0 0 0 0 0 0 +DEAL::cycle:1, fe count:12 0 1 0 0 0 0 0 0 0 0 0 +DEAL::cycle:2, fe count:12 0 0 1 0 0 0 0 0 0 0 0 +DEAL::cycle:3, fe count:8 4 0 0 1 0 0 0 0 0 0 0 +DEAL::cycle:4, fe count:8 0 4 0 0 1 0 0 0 0 0 0 +DEAL::cycle:5, fe count:8 0 0 4 0 0 1 0 0 0 0 0 +DEAL::cycle:6, fe count:4 4 0 0 4 0 0 1 0 0 0 0 +DEAL::cycle:7, fe count:4 0 4 0 0 4 0 0 1 0 0 0 +DEAL::cycle:8, fe count:4 0 0 4 0 0 4 0 0 1 0 0 +DEAL::cycle:9, fe count:0 4 0 0 4 0 0 4 0 0 1 0 +DEAL::cycle:10, fe count:0 0 4 0 0 4 0 0 4 0 0 1 +DEAL::OK