From: Marc Fehling <mafehling.git@gmail.com>
Date: Sat, 6 Mar 2021 04:30:16 +0000 (-0700)
Subject: Parallelize DoFHandler::prepare_coarsening_and_refinement().
X-Git-Tag: v9.3.0-rc1~349^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f5eaf27774143a2d6ae7914b49440164f44a6ba7;p=dealii.git

Parallelize DoFHandler::prepare_coarsening_and_refinement().
---

diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h
index b1494ef54b..1043a8ea05 100644
--- a/include/deal.II/dofs/dof_handler.h
+++ b/include/deal.II/dofs/dof_handler.h
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1998 - 2020 by the deal.II authors
+// Copyright (C) 1998 - 2021 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -725,8 +725,6 @@ public:
    * 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(
diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc
index dc48a04d56..e16d10eb70 100644
--- a/source/dofs/dof_handler.cc
+++ b/source/dofs/dof_handler.cc
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1998 - 2020 by the deal.II authors
+// Copyright (C) 1998 - 2021 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -32,6 +32,8 @@
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_levels.h>
 
+#include <deal.II/lac/la_parallel_vector.templates.h>
+
 #include <algorithm>
 #include <memory>
 #include <set>
@@ -2670,12 +2672,6 @@ DoFHandler<dim, spacedim>::prepare_coarsening_and_refinement(
     // 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
   //
@@ -2692,67 +2688,136 @@ DoFHandler<dim, spacedim>::prepare_coarsening_and_refinement(
   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;
 
+  //
+  // 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
+
+  // 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;
+
+  LinearAlgebra::distributed::Vector<level_type> future_levels;
+  if (const auto parallel_tria =
+        dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
+          &(*tria)))
+    {
+      const auto &partitioner =
+        *parallel_tria->global_active_cell_index_partitioner().lock();
+      future_levels.reinit(partitioner.locally_owned_range(),
+                           partitioner.ghost_indices(),
+                           partitioner.get_mpi_communicator());
+    }
+  else
+    {
+      future_levels.reinit(tria->n_active_cells());
+    }
+
+  for (const auto &cell : active_cell_iterators())
+    if (cell->is_locally_owned())
+      {
+        // set level if FE is part of hierarchy
+        const auto cell_fe_and_level =
+          hierarchy_level_for_fe_index.find(cell->future_fe_index());
+
+        future_levels[cell->global_active_cell_index()] =
+          (cell_fe_and_level != hierarchy_level_for_fe_index.end()) ?
+            cell_fe_and_level->second :
+            invalid_level;
+      }
+
+
   //
   // 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
-  //   (hence iterating from highest to lowest level)
-  // - update future FE indices
+  // - exchange level indices on ghost cells
 
-  bool fe_indices_changed = false;
-  bool fe_indices_changed_in_cycle;
+  bool levels_changed = false;
+  bool levels_changed_in_cycle;
   do
     {
-      fe_indices_changed_in_cycle = false;
+      levels_changed_in_cycle = false;
 
-      for (const auto &cell : active_cell_iterators())
-        {
-          const auto cell_fe_and_level =
-            hierarchy_level_for_fe_index.find(cell->future_fe_index());
+      future_levels.update_ghost_values();
 
-          // ignore cells that are not part of the hierarchy
-          if (cell_fe_and_level == hierarchy_level_for_fe_index.end())
-            continue;
+      for (const auto &cell : active_cell_iterators())
+        if (!cell->is_artificial())
+          {
+            const level_type cell_level =
+              future_levels[cell->global_active_cell_index()];
 
-          const unsigned int cell_level = cell_fe_and_level->second;
+            // 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;
+            // 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)
-              {
-                const auto neighbor = cell->neighbor(f);
+            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());
+                  // 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 =
+                        future_levels[neighbor->global_active_cell_index()];
 
-                // ignore neighbors that are not part of the hierarchy
-                if (neighbor_fe_and_level == hierarchy_level_for_fe_index.end())
-                  continue;
+                      // ignore neighbors that are not part of the hierarchy
+                      if (neighbor_level == invalid_level)
+                        continue;
 
-                const unsigned int neighbor_level =
-                  neighbor_fe_and_level->second;
+                      if ((cell_level - max_difference) > neighbor_level)
+                        {
+                          // update future level
+                          future_levels[neighbor->global_active_cell_index()] =
+                            cell_level - max_difference;
 
-                if ((cell_level - max_difference) > neighbor_level)
-                  {
-                    // update future FE index
-                    const unsigned int new_level = cell_level - max_difference;
-                    neighbor->set_future_fe_index(
-                      fe_index_for_hierarchy_level[new_level]);
+                          levels_changed_in_cycle = true;
+                        }
+                    }
+                }
+          }
 
-                    fe_indices_changed_in_cycle = true;
-                    fe_indices_changed          = true;
-                  }
-              }
-        }
+      levels_changed_in_cycle =
+        Utilities::MPI::logical_or(levels_changed_in_cycle,
+                                   tria->get_communicator());
+      levels_changed |= levels_changed_in_cycle;
     }
-  while (fe_indices_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 =
+          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 fe_indices_changed;
+  return levels_changed;
 }
 
 
diff --git a/tests/hp/prepare_coarsening_and_refinement_01.cc b/tests/hp/prepare_coarsening_and_refinement_01.cc
index a702998be5..10d685e8ad 100644
--- a/tests/hp/prepare_coarsening_and_refinement_01.cc
+++ b/tests/hp/prepare_coarsening_and_refinement_01.cc
@@ -73,9 +73,13 @@ test(const unsigned int fes_size, const unsigned int max_difference)
   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);
+  const bool fe_indices_changed =
+    dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index);
   tria.execute_coarsening_and_refinement();
 
+  (void)fe_indices_changed;
+  Assert(fe_indices_changed, ExcInternalError());
+
   // display number of cells for each FE index
   std::vector<unsigned int> count(fes.size(), 0);
   for (const auto &cell : dofh.active_cell_iterators())
diff --git a/tests/hp/prepare_coarsening_and_refinement_02.cc b/tests/hp/prepare_coarsening_and_refinement_02.cc
index 42057d8603..f46d847f00 100644
--- a/tests/hp/prepare_coarsening_and_refinement_02.cc
+++ b/tests/hp/prepare_coarsening_and_refinement_02.cc
@@ -77,10 +77,15 @@ test(const unsigned int fes_size, const unsigned int max_difference)
       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);
-
+      const bool fe_indices_changed =
+        dofh.prepare_coarsening_and_refinement(max_difference,
+                                               contains_fe_index);
       tria.execute_coarsening_and_refinement();
 
+      (void)fe_indices_changed;
+      if (i >= max_difference)
+        Assert(fe_indices_changed, ExcInternalError());
+
       // display number of cells for each FE index
       std::vector<unsigned int> count(fes.size(), 0);
       for (const auto &cell : dofh.active_cell_iterators())
@@ -93,7 +98,7 @@ test(const unsigned int fes_size, const unsigned int max_difference)
   // 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 double       distance = cell->center().distance(Point<dim>());
       const unsigned int expected_level =
         (sequence.size() - 1) -
         max_difference * static_cast<unsigned int>(std::round(distance));
diff --git a/tests/mpi/prepare_coarsening_and_refinement_01.cc b/tests/mpi/prepare_coarsening_and_refinement_01.cc
new file mode 100644
index 0000000000..26851fc263
--- /dev/null
+++ b/tests/mpi/prepare_coarsening_and_refinement_01.cc
@@ -0,0 +1,160 @@
+// ---------------------------------------------------------------------
+//
+// 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/distributed/shared_tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/tria_base.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/hp/fe_collection.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(parallel::TriangulationBase<dim> &tria,
+     const unsigned int                fes_size,
+     const unsigned int                max_difference)
+{
+  Assert(tria.n_levels() == 0, ExcInternalError());
+  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
+  {
+    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 cell to highest p-level
+  for (const auto &cell : dofh.active_cell_iterators())
+    if (cell->is_locally_owned() && cell->center() == Point<dim>())
+      cell->set_active_fe_index(sequence.back());
+
+  const bool fe_indices_changed =
+    dofh.prepare_coarsening_and_refinement(max_difference, contains_fe_index);
+  tria.execute_coarsening_and_refinement();
+
+  (void)fe_indices_changed;
+  Assert(fe_indices_changed, ExcInternalError());
+
+  // display number of cells for each FE index
+  std::vector<unsigned int> count(fes.size(), 0);
+  for (const auto &cell : dofh.active_cell_iterators())
+    if (cell->is_locally_owned())
+      count[cell->active_fe_index()]++;
+  Utilities::MPI::sum(count, tria.get_communicator(), count);
+  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())
+    if (cell->is_locally_owned())
+      {
+        const double       distance = cell->center().distance(Point<dim>());
+        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(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    initlog();
+
+  constexpr const unsigned int dim = 2;
+
+  deallog << "parallel::shared::Triangulation" << std::endl;
+  {
+    parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+    test<dim>(tria, 4, 1);
+    tria.clear();
+    test<dim>(tria, 8, 2);
+    tria.clear();
+    test<dim>(tria, 12, 3);
+  }
+
+  deallog << "parallel::shared::Triangulation with artificial cells"
+          << std::endl;
+  {
+    parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD,
+                                              Triangulation<dim>::none,
+                                              /*allow_artificial_cells=*/true);
+
+    test<dim>(tria, 4, 1);
+    tria.clear();
+    test<dim>(tria, 8, 2);
+    tria.clear();
+    test<dim>(tria, 12, 3);
+  }
+
+  deallog << "parallel::distributed::Triangulation" << std::endl;
+  {
+    parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+    test<dim>(tria, 4, 1);
+    tria.clear();
+    test<dim>(tria, 8, 2);
+    tria.clear();
+    test<dim>(tria, 12, 3);
+  }
+}
diff --git a/tests/mpi/prepare_coarsening_and_refinement_01.with_p4est=true.mpirun=1.output b/tests/mpi/prepare_coarsening_and_refinement_01.with_p4est=true.mpirun=1.output
new file mode 100644
index 0000000000..065ac3fadb
--- /dev/null
+++ b/tests/mpi/prepare_coarsening_and_refinement_01.with_p4est=true.mpirun=1.output
@@ -0,0 +1,40 @@
+
+DEAL::parallel::shared::Triangulation
+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
+DEAL::parallel::shared::Triangulation with artificial cells
+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
+DEAL::parallel::distributed::Triangulation
+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/mpi/prepare_coarsening_and_refinement_01.with_p4est=true.mpirun=4.output b/tests/mpi/prepare_coarsening_and_refinement_01.with_p4est=true.mpirun=4.output
new file mode 100644
index 0000000000..065ac3fadb
--- /dev/null
+++ b/tests/mpi/prepare_coarsening_and_refinement_01.with_p4est=true.mpirun=4.output
@@ -0,0 +1,40 @@
+
+DEAL::parallel::shared::Triangulation
+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
+DEAL::parallel::shared::Triangulation with artificial cells
+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
+DEAL::parallel::distributed::Triangulation
+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/mpi/prepare_coarsening_and_refinement_02.cc b/tests/mpi/prepare_coarsening_and_refinement_02.cc
new file mode 100644
index 0000000000..e324e09652
--- /dev/null
+++ b/tests/mpi/prepare_coarsening_and_refinement_02.cc
@@ -0,0 +1,167 @@
+// ---------------------------------------------------------------------
+//
+// 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/distributed/shared_tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/tria_base.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/hp/fe_collection.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(parallel::TriangulationBase<dim> &tria,
+     const unsigned int                fes_size,
+     const unsigned int                max_difference)
+{
+  Assert(tria.n_levels() == 0, ExcInternalError());
+  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
+  {
+    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
+  for (unsigned int i = 0; i < sequence.size() - 1; ++i)
+    {
+      // find center cell
+      for (const auto &cell : dofh.active_cell_iterators())
+        if (cell->is_locally_owned() && cell->center() == Point<dim>())
+          cell->set_future_fe_index(
+            fes.next_in_hierarchy(cell->active_fe_index()));
+
+      const bool fe_indices_changed =
+        dofh.prepare_coarsening_and_refinement(max_difference,
+                                               contains_fe_index);
+      tria.execute_coarsening_and_refinement();
+
+      (void)fe_indices_changed;
+      if (i >= max_difference)
+        Assert(fe_indices_changed, ExcInternalError());
+
+      // display number of cells for each FE index
+      std::vector<unsigned int> count(fes.size(), 0);
+      for (const auto &cell : dofh.active_cell_iterators())
+        if (cell->is_locally_owned())
+          count[cell->active_fe_index()]++;
+      Utilities::MPI::sum(count, tria.get_communicator(), count);
+      deallog << "cycle:" << i << ", 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())
+    if (cell->is_locally_owned())
+      {
+        const double       distance = cell->center().distance(Point<dim>());
+        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(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    initlog();
+
+  constexpr const unsigned int dim = 2;
+
+  deallog << "parallel::shared::Triangulation" << std::endl;
+  {
+    parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+    test<dim>(tria, 4, 1);
+    tria.clear();
+    test<dim>(tria, 8, 2);
+    tria.clear();
+    test<dim>(tria, 12, 3);
+  }
+
+  deallog << "parallel::shared::Triangulation with artificial cells"
+          << std::endl;
+  {
+    parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD,
+                                              Triangulation<dim>::none,
+                                              /*allow_artificial_cells=*/true);
+
+    test<dim>(tria, 4, 1);
+    tria.clear();
+    test<dim>(tria, 8, 2);
+    tria.clear();
+    test<dim>(tria, 12, 3);
+  }
+
+  deallog << "parallel::distributed::Triangulation" << std::endl;
+  {
+    parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+    test<dim>(tria, 4, 1);
+    tria.clear();
+    test<dim>(tria, 8, 2);
+    tria.clear();
+    test<dim>(tria, 12, 3);
+  }
+}
diff --git a/tests/mpi/prepare_coarsening_and_refinement_02.with_p4est=true.mpirun=1.output b/tests/mpi/prepare_coarsening_and_refinement_02.with_p4est=true.mpirun=1.output
new file mode 100644
index 0000000000..562cc53c58
--- /dev/null
+++ b/tests/mpi/prepare_coarsening_and_refinement_02.with_p4est=true.mpirun=1.output
@@ -0,0 +1,94 @@
+
+DEAL::parallel::shared::Triangulation
+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
+DEAL::parallel::shared::Triangulation with artificial cells
+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
+DEAL::parallel::distributed::Triangulation
+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
diff --git a/tests/mpi/prepare_coarsening_and_refinement_02.with_p4est=true.mpirun=4.output b/tests/mpi/prepare_coarsening_and_refinement_02.with_p4est=true.mpirun=4.output
new file mode 100644
index 0000000000..562cc53c58
--- /dev/null
+++ b/tests/mpi/prepare_coarsening_and_refinement_02.with_p4est=true.mpirun=4.output
@@ -0,0 +1,94 @@
+
+DEAL::parallel::shared::Triangulation
+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
+DEAL::parallel::shared::Triangulation with artificial cells
+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
+DEAL::parallel::distributed::Triangulation
+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