--- /dev/null
+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.
+<br>
+(Marc Fehling, 2021/02/23)
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.
*/
* 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 <int dim, int spacedim>
void
* 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 <int dim, int spacedim>
void
* 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 <int dim, typename Number, int spacedim>
void
* 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 <int dim, typename Number, int spacedim>
void
* 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 <int dim, typename Number, int spacedim>
void
*
* 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 <int dim, typename Number, int spacedim>
void
* 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 <int dim, typename Number, int spacedim>
void
*
* @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 <int dim, typename Number, int spacedim>
void
* 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 <int dim, int spacedim>
void
* 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 <int dim, int spacedim>
void
+template <int dim, int spacedim>
+bool
+DoFHandler<dim, spacedim>::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<const parallel::TriangulationBase<dim, spacedim> *>(
+ &(*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<unsigned int> 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<unsigned int, unsigned int> 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<std::list<active_cell_iterator>> 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 <int dim, int spacedim>
void
DoFHandler<dim, spacedim>::initialize_local_block_info()
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+template <int dim>
+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<dim> fes;
+ while (fes.size() < fes_size)
+ fes.push_back(FE_Q<dim>(1));
+
+ const unsigned int contains_fe_index = 0;
+ const auto sequence = fes.get_hierarchy_sequence(contains_fe_index);
+
+ // setup cross-shaped mesh
+ Triangulation<dim> tria;
+ {
+ std::vector<unsigned int> sizes(Utilities::pow(2, dim),
+ static_cast<unsigned int>(
+ (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<dim> dofh(tria);
+ dofh.distribute_dofs(fes);
+
+ // set center to highest p-level
+ const auto center_cell = dofh.begin();
+ Assert(center_cell->center() == Point<dim>(), 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<unsigned int> 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<unsigned int>(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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+template <int dim>
+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<dim> fes;
+ while (fes.size() < fes_size)
+ fes.push_back(FE_Q<dim>(1));
+
+ const unsigned int contains_fe_index = 0;
+ const auto sequence = fes.get_hierarchy_sequence(contains_fe_index);
+
+ // setup cross-shaped mesh
+ Triangulation<dim> tria;
+ {
+ std::vector<unsigned int> sizes(Utilities::pow(2, dim),
+ static_cast<unsigned int>(
+ (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<dim> dofh(tria);
+ dofh.distribute_dofs(fes);
+
+ // increase p-level in center subsequently
+ const auto center_cell = dofh.begin();
+ Assert(center_cell->center() == Point<dim>(), 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<unsigned int> 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<unsigned int>(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);
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff