Restructure TriangulationDescription::Utilities::create_description_from_triangulatio...

diff --git a/source/grid/tria_description.cc b/source/grid/tria_description.cc
index fa6df293558e62b1b4634334e1217d980de76b6e..37210bbb8cc5854f567f673c7a9ae07d5007d11c 100644 (file)
--- a/source/grid/tria_description.cc
+++ b/source/grid/tria_description.cc
@@ -82,6 +82,244 @@ namespace TriangulationDescription
   {
     namespace
     {
+      /**
+       * A temporal class similar to TriangulationDescription::Description.
+       */
+      template <int dim, int spacedim>
+      struct DescriptionTemp
+      {
+        /**
+         * Serialization function for packing and unpacking the content of this
+         * class.
+         */
+        template <class Archive>
+        void
+        serialize(Archive &ar, const unsigned int /*version*/)
+        {
+          ar &coarse_cells;
+          ar &coarse_cell_vertices;
+          ar &coarse_cell_index_to_coarse_cell_id;
+          ar &cell_infos;
+        }
+
+        /**
+         * Send the DescriptionTemp objects to the right process.
+         */
+        void
+        collect(
+          const std::vector<unsigned int> &                  relevant_processes,
+          const std::vector<DescriptionTemp<dim, spacedim>> &description_temp,
+          const MPI_Comm &                                   comm)
+        {
+          dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess<char,
+                                                                        char>
+            process(
+              [&]() { return relevant_processes; },
+              [&](const unsigned int other_rank,
+                  std::vector<char> &send_buffer) {
+                const auto ptr = std::find(relevant_processes.begin(),
+                                           relevant_processes.end(),
+                                           other_rank);
+
+                Assert(ptr != relevant_processes.end(), ExcInternalError());
+
+                const auto other_rank_index =
+                  std::distance(relevant_processes.begin(), ptr);
+
+                send_buffer =
+                  dealii::Utilities::pack(description_temp[other_rank_index],
+                                          false);
+              },
+              [&](const unsigned int &,
+                  const std::vector<char> &recv_buffer,
+                  std::vector<char> &) {
+                this->merge(
+                  dealii::Utilities::unpack<DescriptionTemp<dim, spacedim>>(
+                    recv_buffer, false));
+              });
+
+          dealii::Utilities::MPI::ConsensusAlgorithms::Selector<char, char>(
+            process, comm)
+            .run();
+        }
+
+        /**
+         * Merge the given DescriptionTemp object into the current one. Note
+         * that no duplicate cells, vertices, ... are removed at this stage.
+         * This done by the reduce() function.
+         */
+        void
+        merge(const DescriptionTemp<dim, spacedim> &other)
+        {
+          this->cell_infos.resize(
+            std::max(other.cell_infos.size(), this->cell_infos.size()));
+
+          this->coarse_cells.insert(this->coarse_cells.end(),
+                                    other.coarse_cells.begin(),
+                                    other.coarse_cells.end());
+          this->coarse_cell_vertices.insert(this->coarse_cell_vertices.end(),
+                                            other.coarse_cell_vertices.begin(),
+                                            other.coarse_cell_vertices.end());
+          this->coarse_cell_index_to_coarse_cell_id.insert(
+            this->coarse_cell_index_to_coarse_cell_id.end(),
+            other.coarse_cell_index_to_coarse_cell_id.begin(),
+            other.coarse_cell_index_to_coarse_cell_id.end());
+
+          for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
+            this->cell_infos[i].insert(this->cell_infos[i].end(),
+                                       other.cell_infos[i].begin(),
+                                       other.cell_infos[i].end());
+        }
+
+        /**
+         * Remove all duplicate information obtained during merge().
+         */
+        void
+        reduce()
+        {
+          // make coarse cells unique
+          {
+            std::vector<std::tuple<types::coarse_cell_id,
+                                   dealii::CellData<dim>,
+                                   unsigned int>>
+              temp;
+
+            for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
+              temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
+                                this->coarse_cells[i],
+                                i);
+
+            std::sort(temp.begin(),
+                      temp.end(),
+                      [](const auto &a, const auto &b) {
+                        return std::get<0>(a) < std::get<0>(b);
+                      });
+            temp.erase(std::unique(temp.begin(),
+                                   temp.end(),
+                                   [](const auto &a, const auto &b) {
+                                     return std::get<0>(a) == std::get<0>(b);
+                                   }),
+                       temp.end());
+            std::sort(temp.begin(),
+                      temp.end(),
+                      [](const auto &a, const auto &b) {
+                        return std::get<2>(a) < std::get<2>(b);
+                      });
+
+            this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
+            this->coarse_cells.resize(temp.size());
+
+            for (unsigned int i = 0; i < temp.size(); ++i)
+              {
+                this->coarse_cell_index_to_coarse_cell_id[i] =
+                  std::get<0>(temp[i]);
+                this->coarse_cells[i] = std::get<1>(temp[i]);
+              }
+          }
+
+          // make coarse cell vertices unique
+          {
+            std::sort(this->coarse_cell_vertices.begin(),
+                      this->coarse_cell_vertices.end(),
+                      [](const auto &a, const auto &b) {
+                        return a.first < b.first;
+                      });
+            this->coarse_cell_vertices.erase(
+              std::unique(this->coarse_cell_vertices.begin(),
+                          this->coarse_cell_vertices.end(),
+                          [](const auto &a, const auto &b) {
+                            return a.first == b.first;
+                          }),
+              this->coarse_cell_vertices.end());
+          }
+
+          // make cells unique
+          for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
+            {
+              if (this->cell_infos[i].size() == 0)
+                continue;
+
+              std::sort(this->cell_infos[i].begin(),
+                        this->cell_infos[i].end(),
+                        [](const auto &a, const auto &b) {
+                          return a.id < b.id;
+                        });
+
+              std::vector<CellData<dim>> temp;
+              temp.push_back(this->cell_infos[i][0]);
+
+              for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j)
+                if (temp.back().id == cell_infos[i][j].id)
+                  {
+                    temp.back().subdomain_id =
+                      std::min(temp.back().subdomain_id,
+                               this->cell_infos[i][j].subdomain_id);
+                    temp.back().level_subdomain_id =
+                      std::min(temp.back().level_subdomain_id,
+                               this->cell_infos[i][j].level_subdomain_id);
+                  }
+                else
+                  {
+                    temp.push_back(this->cell_infos[i][j]);
+                  }
+
+              this->cell_infos[i] = temp;
+            }
+        }
+
+        /**
+         * Convert the content of this class to a
+         * TriangulationDescription::Description, which can be used during
+         * Triangulation::create_triangulation().
+         */
+        Description<dim, spacedim>
+        convert(const MPI_Comm comm,
+                const typename Triangulation<dim, spacedim>::MeshSmoothing
+                                                         mesh_smoothing,
+                const TriangulationDescription::Settings settings)
+        {
+          Description<dim, spacedim> description;
+
+          // copy communicator
+          description.comm = comm;
+
+          description.settings = settings;
+
+          // use mesh smoothing from base triangulation
+          description.smoothing = mesh_smoothing;
+
+          std::map<unsigned int, unsigned int> map;
+
+          for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
+            {
+              description.coarse_cell_vertices.push_back(
+                this->coarse_cell_vertices[i].second);
+              map[this->coarse_cell_vertices[i].first] = i;
+            }
+
+          description.coarse_cells = this->coarse_cells;
+
+          for (auto &cell : description.coarse_cells)
+            for (unsigned int v = 0; v < cell.vertices.size(); ++v)
+              cell.vertices[v] = map[cell.vertices[v]];
+
+          description.coarse_cell_index_to_coarse_cell_id =
+            this->coarse_cell_index_to_coarse_cell_id;
+          description.cell_infos = this->cell_infos;
+
+          return description;
+        }
+
+        std::vector<dealii::CellData<dim>> coarse_cells;
+
+        std::vector<std::pair<unsigned int, Point<spacedim>>>
+          coarse_cell_vertices;
+
+        std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
+
+        std::vector<std::vector<CellData<dim>>> cell_infos;
+      };
+
       /**
        * Set the user_flag of a cell and of all its parent cells.
        */
@@ -94,6 +332,345 @@ namespace TriangulationDescription
         if (cell->level() != 0)
           set_user_flag_and_of_its_parents(cell->parent());
       }
+
+      /**
+       * Helper calls for the
+       * TriangulationDescription::Utilities::create_description_from_triangulation()
+       * functions.
+       */
+      template <int dim, int spacedim>
+      class CreateDescriptionFromTriangulationHelper
+      {
+      public:
+        CreateDescriptionFromTriangulationHelper(
+          const dealii::Triangulation<dim, spacedim> &tria,
+          const std::function<types::subdomain_id(
+            const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+              &)> &                                   subdomain_id_function,
+          const std::function<types::subdomain_id(
+            const typename dealii::Triangulation<dim, spacedim>::cell_iterator
+              &)> &                                level_subdomain_id_function,
+          const MPI_Comm &                         comm,
+          const TriangulationDescription::Settings settings)
+          : tria(tria)
+          , subdomain_id_function(subdomain_id_function)
+          , level_subdomain_id_function(level_subdomain_id_function)
+          , comm(comm)
+          , settings(settings)
+          , construct_multigrid(
+              settings &
+              TriangulationDescription::Settings::construct_multigrid_hierarchy)
+        {
+          Assert(
+            !(settings & TriangulationDescription::Settings::
+                           construct_multigrid_hierarchy) ||
+              (tria.get_mesh_smoothing() &
+               Triangulation<dim,
+                             spacedim>::limit_level_difference_at_vertices),
+            ExcMessage(
+              "Source triangulation has to be set up with "
+              "limit_level_difference_at_vertices if the construction of the "
+              "multigrid hierarchy is requested!"));
+
+          GridTools::collect_coinciding_vertices(
+            tria, coinciding_vertex_groups, vertex_to_coinciding_vertex_group);
+        }
+
+        template <typename T>
+        T
+        create_description_for_rank(const unsigned int my_rank) const
+        {
+          T construction_data;
+
+          set_additional_data(construction_data);
+
+          // helper function, which collects all vertices belonging to a cell
+          // (also taking into account periodicity)
+          const auto
+            add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
+              [this](const auto &       cell,
+                     std::vector<bool> &vertices_owned_by_locally_owned_cells) {
+                // add local vertices
+                for (const auto v : cell->vertex_indices())
+                  {
+                    vertices_owned_by_locally_owned_cells[cell->vertex_index(
+                      v)] = true;
+                    const auto coinciding_vertex_group =
+                      vertex_to_coinciding_vertex_group.find(
+                        cell->vertex_index(v));
+                    if (coinciding_vertex_group !=
+                        vertex_to_coinciding_vertex_group.end())
+                      for (const auto &co_vertex : coinciding_vertex_groups.at(
+                             coinciding_vertex_group->second))
+                        vertices_owned_by_locally_owned_cells[co_vertex] = true;
+                  }
+              };
+
+          // 1) collect locally relevant cells (set user_flag)
+          std::vector<bool> old_user_flags;
+          tria.save_user_flags(old_user_flags);
+
+          // 1a) clear user_flags
+          const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
+            .clear_user_flags();
+
+          // 1b) loop over levels (from fine to coarse) and mark on each level
+          //     the locally relevant cells
+          for (int level = tria.get_triangulation().n_global_levels() - 1;
+               level >= 0;
+               --level)
+            {
+              // collect vertices connected to a (on any level) locally owned
+              // cell
+              std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
+                tria.n_vertices());
+              for (const auto &cell : tria.cell_iterators_on_level(level))
+                if (construct_multigrid &&
+                    (level_subdomain_id_function(cell) == my_rank))
+                  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
+                    cell, vertices_owned_by_locally_owned_cells_on_level);
+
+              for (const auto &cell : tria.active_cell_iterators())
+                if (subdomain_id_function(cell) == my_rank)
+                  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
+                    cell, vertices_owned_by_locally_owned_cells_on_level);
+
+              // helper function to determine if cell is locally relevant
+              // (i.e. a cell which is connected to a vertex via a locally
+              // owned cell)
+              const auto is_locally_relevant_on_level = [&](const auto &cell) {
+                for (const auto v : cell->vertex_indices())
+                  if (vertices_owned_by_locally_owned_cells_on_level
+                        [cell->vertex_index(v)])
+                    return true;
+                return false;
+              };
+
+              // mark all locally relevant cells
+              for (const auto &cell : tria.cell_iterators_on_level(level))
+                if (is_locally_relevant_on_level(cell))
+                  set_user_flag_and_of_its_parents(cell);
+            }
+
+          // 2) set_up coarse-grid triangulation
+          {
+            std::vector<bool> vertices_locally_relevant(tria.n_vertices(),
+                                                        false);
+
+            // a) loop over all cells
+            for (const auto &cell : tria.cell_iterators_on_level(0))
+              {
+                if (!cell->user_flag_set())
+                  continue;
+
+                // extract cell definition (with old numbering of vertices)
+                dealii::CellData<dim> cell_data(cell->n_vertices());
+                cell_data.material_id = cell->material_id();
+                cell_data.manifold_id = cell->manifold_id();
+                for (const auto v : cell->vertex_indices())
+                  cell_data.vertices[v] = cell->vertex_index(v);
+                construction_data.coarse_cells.push_back(cell_data);
+
+                // save indices of each vertex of this cell
+                for (const auto v : cell->vertex_indices())
+                  vertices_locally_relevant[cell->vertex_index(v)] = true;
+
+                // save translation for corase grid: lid -> gid
+                construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
+                  cell->id().get_coarse_cell_id());
+              }
+
+            add_vertices(construction_data, vertices_locally_relevant);
+          }
+
+
+          // 3) collect info of each cell
+          construction_data.cell_infos.resize(
+            tria.get_triangulation().n_global_levels());
+
+          // collect local vertices on active level
+          std::vector<bool> vertices_owned_by_locally_owned_active_cells(
+            tria.n_vertices());
+          for (const auto &cell : tria.active_cell_iterators())
+            if (subdomain_id_function(cell) == my_rank)
+              add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
+                cell, vertices_owned_by_locally_owned_active_cells);
+
+          // helper function to determine if cell is locally relevant
+          // on active level
+          const auto is_locally_relevant_on_active_level =
+            [&](const auto &cell) {
+              if (cell->is_active())
+                for (const auto v : cell->vertex_indices())
+                  if (vertices_owned_by_locally_owned_active_cells
+                        [cell->vertex_index(v)])
+                    return true;
+              return false;
+            };
+
+          for (unsigned int level = 0;
+               level < tria.get_triangulation().n_global_levels();
+               ++level)
+            {
+              // collect local vertices on level
+              std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
+                tria.n_vertices());
+              for (const auto &cell : tria.cell_iterators_on_level(level))
+                if ((construct_multigrid &&
+                     (level_subdomain_id_function(cell) == my_rank)) ||
+                    (cell->is_active() &&
+                     subdomain_id_function(cell) == my_rank))
+                  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
+                    cell, vertices_owned_by_locally_owned_cells_on_level);
+
+              // helper function to determine if cell is locally relevant
+              // on level
+              const auto is_locally_relevant_on_level = [&](const auto &cell) {
+                for (const auto v : cell->vertex_indices())
+                  if (vertices_owned_by_locally_owned_cells_on_level
+                        [cell->vertex_index(v)])
+                    return true;
+                return false;
+              };
+
+              auto &level_cell_infos = construction_data.cell_infos[level];
+              for (const auto &cell : tria.cell_iterators_on_level(level))
+                {
+                  // check if cell is locally relevant
+                  if (!(cell->user_flag_set()))
+                    continue;
+
+                  CellData<dim> cell_info;
+
+                  // save coarse-cell id
+                  cell_info.id = cell->id().template to_binary<dim>();
+
+                  // save boundary_ids of each face of this cell
+                  for (const auto f : cell->face_indices())
+                    {
+                      types::boundary_id boundary_ind =
+                        cell->face(f)->boundary_id();
+                      if (boundary_ind != numbers::internal_face_boundary_id)
+                        cell_info.boundary_ids.emplace_back(f, boundary_ind);
+                    }
+
+                  // save manifold id
+                  {
+                    // ... of cell
+                    cell_info.manifold_id = cell->manifold_id();
+
+                    // ... of lines
+                    if (dim >= 2)
+                      for (const auto line : cell->line_indices())
+                        cell_info.manifold_line_ids[line] =
+                          cell->line(line)->manifold_id();
+
+                    // ... of quads
+                    if (dim == 3)
+                      for (const auto f : cell->face_indices())
+                        cell_info.manifold_quad_ids[f] =
+                          cell->quad(f)->manifold_id();
+                  }
+
+                  // subdomain and level subdomain id
+                  cell_info.subdomain_id = numbers::artificial_subdomain_id;
+                  cell_info.level_subdomain_id =
+                    numbers::artificial_subdomain_id;
+
+                  if (is_locally_relevant_on_active_level(cell))
+                    {
+                      cell_info.subdomain_id = subdomain_id_function(cell);
+
+                      cell_info.level_subdomain_id =
+                        level_subdomain_id_function(cell);
+                    }
+                  else if (is_locally_relevant_on_level(cell))
+                    {
+                      cell_info.level_subdomain_id =
+                        level_subdomain_id_function(cell);
+                    }
+                  else
+                    {
+                      // cell is locally relevant but an artificial cell
+                    }
+
+                  level_cell_infos.emplace_back(cell_info);
+                }
+            }
+
+          const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
+            .load_user_flags(old_user_flags);
+
+          return construction_data;
+        }
+
+      private:
+        void
+        set_additional_data(Description<dim, spacedim> &construction_data) const
+        {
+          construction_data.comm      = comm;
+          construction_data.smoothing = tria.get_mesh_smoothing();
+          construction_data.settings  = settings;
+        }
+
+        void
+        set_additional_data(DescriptionTemp<dim, spacedim> &) const
+        {
+          // nothing to do
+        }
+
+        void
+        add_vertices(Description<dim, spacedim> &construction_data,
+                     const std::vector<bool> &vertices_locally_relevant) const
+        {
+          std::vector<unsigned int> vertices_locally_relevant_indices(
+            vertices_locally_relevant.size());
+
+          // enumerate locally relevant vertices
+          unsigned int vertex_counter = 0;
+          for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
+            if (vertices_locally_relevant[i])
+              {
+                construction_data.coarse_cell_vertices.push_back(
+                  tria.get_vertices()[i]);
+                vertices_locally_relevant_indices[i] = vertex_counter++;
+              }
+
+          // correct vertices of cells (make them local)
+          for (auto &cell : construction_data.coarse_cells)
+            for (unsigned int v = 0; v < cell.vertices.size(); ++v)
+              cell.vertices[v] =
+                vertices_locally_relevant_indices[cell.vertices[v]];
+        }
+
+        void
+        add_vertices(DescriptionTemp<dim, spacedim> &construction_data,
+                     const std::vector<bool> &vertices_locally_relevant) const
+        {
+          for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
+            if (vertices_locally_relevant[i])
+              construction_data.coarse_cell_vertices.emplace_back(
+                i, tria.get_vertices()[i]);
+        }
+
+
+        const dealii::Triangulation<dim, spacedim> &tria;
+        const std::function<types::subdomain_id(
+          const typename dealii::Triangulation<dim, spacedim>::cell_iterator &)>
+          subdomain_id_function;
+        const std::function<types::subdomain_id(
+          const typename dealii::Triangulation<dim, spacedim>::cell_iterator &)>
+          level_subdomain_id_function;
+
+        const MPI_Comm &                         comm;
+        const TriangulationDescription::Settings settings;
+        const bool                               construct_multigrid;
+
+        std::map<unsigned int, std::vector<unsigned int>>
+                                             coinciding_vertex_groups;
+        std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
+      };
+
     } // namespace
 
 
@@ -105,7 +682,7 @@ namespace TriangulationDescription
       const TriangulationDescription::Settings    settings,
       const unsigned int                          my_rank_in)
     {
-      if (auto tria_pdt = dynamic_cast<
+      if (const auto tria_pdt = dynamic_cast<
             const parallel::distributed::Triangulation<dim, spacedim> *>(&tria))
         Assert(comm == tria_pdt->get_communicator(),
                ExcMessage("MPI communicators do not match."));
@@ -141,255 +718,22 @@ namespace TriangulationDescription
                  ExcMessage("This type of triangulation is not supported!"));
         }
 
-      Description<dim, spacedim> construction_data;
-
-      // store the communicator
-      construction_data.comm = comm;
-
-      std::map<unsigned int, std::vector<unsigned int>>
-                                           coinciding_vertex_groups;
-      std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
-
-      GridTools::collect_coinciding_vertices(tria,
-                                             coinciding_vertex_groups,
-                                             vertex_to_coinciding_vertex_group);
-
-      // helper function, which collects all vertices belonging to a cell
-      // (also taking into account periodicity)
-      const auto add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
-        [&coinciding_vertex_groups, &vertex_to_coinciding_vertex_group](
-          const auto &       cell,
-          std::vector<bool> &vertices_owned_by_locally_owned_cells) {
-          // add local vertices
-          for (const auto v : cell->vertex_indices())
-            {
-              vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] =
-                true;
-              const auto coinciding_vertex_group =
-                vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
-              if (coinciding_vertex_group !=
-                  vertex_to_coinciding_vertex_group.end())
-                for (const auto &co_vertex : coinciding_vertex_groups.at(
-                       coinciding_vertex_group->second))
-                  vertices_owned_by_locally_owned_cells[co_vertex] = true;
-            }
-        };
-
-      construction_data.smoothing = tria.get_mesh_smoothing();
-      construction_data.settings  = settings;
-
-      const bool construct_multigrid =
-        settings &
-        TriangulationDescription::Settings::construct_multigrid_hierarchy;
-
-      Assert(
-        !(settings &
-          TriangulationDescription::Settings::construct_multigrid_hierarchy) ||
-          (tria.get_mesh_smoothing() &
-           Triangulation<dim, spacedim>::limit_level_difference_at_vertices),
-        ExcMessage(
-          "Source triangulation has to be set up with "
-          "limit_level_difference_at_vertices if the construction of the "
-          "multigrid hierarchy is requested!"));
-
-      // 1) collect locally relevant cells (set user_flag)
-      std::vector<bool> old_user_flags;
-      tria.save_user_flags(old_user_flags);
-
-      // 1a) clear user_flags
-      const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
-        .clear_user_flags();
-
-      // 1b) loop over levels (from fine to coarse) and mark on each level
-      //     the locally relevant cells
-      for (int level = tria.get_triangulation().n_global_levels() - 1;
-           level >= 0;
-           --level)
-        {
-          // collect vertices connected to a (on any level) locally owned
-          // cell
-          std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
-            tria.n_vertices());
-          for (const auto &cell : tria.cell_iterators_on_level(level))
-            if (construct_multigrid && (cell->level_subdomain_id() == my_rank))
-              add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
-                cell, vertices_owned_by_locally_owned_cells_on_level);
-
-          for (const auto &cell : tria.active_cell_iterators())
-            if (cell->subdomain_id() == my_rank)
-              add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
-                cell, vertices_owned_by_locally_owned_cells_on_level);
-
-          // helper function to determine if cell is locally relevant
-          // (i.e. a cell which is connected to a vertex via a locally
-          // owned cell)
-          const auto is_locally_relevant_on_level = [&](const auto &cell) {
-            for (const auto v : cell->vertex_indices())
-              if (vertices_owned_by_locally_owned_cells_on_level
-                    [cell->vertex_index(v)])
-                return true;
-            return false;
-          };
-
-          // mark all locally relevant cells
-          for (const auto &cell : tria.cell_iterators_on_level(level))
-            if (is_locally_relevant_on_level(cell))
-              set_user_flag_and_of_its_parents(cell);
-        }
-
-      // 2) set_up coarse-grid triangulation
-      {
-        std::map<unsigned int, unsigned int> vertices_locally_relevant;
-
-        // a) loop over all cells
-        for (const auto &cell : tria.cell_iterators_on_level(0))
-          {
-            if (!cell->user_flag_set())
-              continue;
-
-            // extract cell definition (with old numbering of vertices)
-            dealii::CellData<dim> cell_data(cell->n_vertices());
-            cell_data.material_id = cell->material_id();
-            cell_data.manifold_id = cell->manifold_id();
-            for (const auto v : cell->vertex_indices())
-              cell_data.vertices[v] = cell->vertex_index(v);
-            construction_data.coarse_cells.push_back(cell_data);
-
-            // save indices of each vertex of this cell
-            for (const auto v : cell->vertex_indices())
-              vertices_locally_relevant[cell->vertex_index(v)] =
-                numbers::invalid_unsigned_int;
-
-            // save translation for corase grid: lid -> gid
-            construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
-              cell->id().get_coarse_cell_id());
-          }
-
-        // b) enumerate locally relevant vertices
-        unsigned int vertex_counter = 0;
-        for (auto &vertex : vertices_locally_relevant)
-          {
-            construction_data.coarse_cell_vertices.push_back(
-              tria.get_vertices()[vertex.first]);
-            vertex.second = vertex_counter++;
-          }
-
-        // c) correct vertices of cells (make them local)
-        for (auto &cell : construction_data.coarse_cells)
-          for (unsigned int v = 0; v < cell.vertices.size(); ++v)
-            cell.vertices[v] = vertices_locally_relevant[cell.vertices[v]];
-      }
-
-
-      // 3) collect info of each cell
-      construction_data.cell_infos.resize(
-        tria.get_triangulation().n_global_levels());
-
-      // collect local vertices on active level
-      std::vector<bool> vertices_owned_by_locally_owned_active_cells(
-        tria.n_vertices());
-      for (const auto &cell : tria.active_cell_iterators())
-        if (cell->subdomain_id() == my_rank)
-          add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
-            cell, vertices_owned_by_locally_owned_active_cells);
-
-      // helper function to determine if cell is locally relevant
-      // on active level
-      const auto is_locally_relevant_on_active_level = [&](const auto &cell) {
-        if (cell->is_active())
-          for (const auto v : cell->vertex_indices())
-            if (vertices_owned_by_locally_owned_active_cells[cell->vertex_index(
-                  v)])
-              return true;
-        return false;
+      const auto subdomain_id_function = [](const auto &cell) {
+        return cell->subdomain_id();
       };
 
-      for (unsigned int level = 0;
-           level < tria.get_triangulation().n_global_levels();
-           ++level)
-        {
-          // collect local vertices on level
-          std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
-            tria.n_vertices());
-          for (const auto &cell : tria.cell_iterators_on_level(level))
-            if ((construct_multigrid &&
-                 (cell->level_subdomain_id() == my_rank)) ||
-                (cell->is_active() && cell->subdomain_id() == my_rank))
-              add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
-                cell, vertices_owned_by_locally_owned_cells_on_level);
-
-          // helper function to determine if cell is locally relevant
-          // on level
-          const auto is_locally_relevant_on_level = [&](const auto &cell) {
-            for (const auto v : cell->vertex_indices())
-              if (vertices_owned_by_locally_owned_cells_on_level
-                    [cell->vertex_index(v)])
-                return true;
-            return false;
-          };
-
-          auto &level_cell_infos = construction_data.cell_infos[level];
-          for (const auto &cell : tria.cell_iterators_on_level(level))
-            {
-              // check if cell is locally relevant
-              if (!(cell->user_flag_set()))
-                continue;
-
-              CellData<dim> cell_info;
-
-              // save coarse-cell id
-              cell_info.id = cell->id().template to_binary<dim>();
-
-              // save boundary_ids of each face of this cell
-              for (const auto f : cell->face_indices())
-                {
-                  types::boundary_id boundary_ind =
-                    cell->face(f)->boundary_id();
-                  if (boundary_ind != numbers::internal_face_boundary_id)
-                    cell_info.boundary_ids.emplace_back(f, boundary_ind);
-                }
-
-              // save manifold id
-              {
-                // ... of cell
-                cell_info.manifold_id = cell->manifold_id();
-
-                // ... of lines
-                if (dim >= 2)
-                  for (const auto line : cell->line_indices())
-                    cell_info.manifold_line_ids[line] =
-                      cell->line(line)->manifold_id();
-
-                // ... of quads
-                if (dim == 3)
-                  for (const auto f : cell->face_indices())
-                    cell_info.manifold_quad_ids[f] =
-                      cell->quad(f)->manifold_id();
-              }
-
-              // subdomain and level subdomain id
-              cell_info.subdomain_id       = numbers::artificial_subdomain_id;
-              cell_info.level_subdomain_id = numbers::artificial_subdomain_id;
-
-              if (is_locally_relevant_on_active_level(cell))
-                {
-                  cell_info.level_subdomain_id = cell->level_subdomain_id();
-                  cell_info.subdomain_id       = cell->subdomain_id();
-                }
-              else if (is_locally_relevant_on_level(cell))
-                {
-                  cell_info.level_subdomain_id = cell->level_subdomain_id();
-                }
-              // else: cell is locally relevant but an artificial cell
-
-              level_cell_infos.emplace_back(cell_info);
-            }
-        }
-
-      const_cast<dealii::Triangulation<dim, spacedim> &>(tria).load_user_flags(
-        old_user_flags);
+      const auto level_subdomain_id_function = [](const auto &cell) {
+        return cell->level_subdomain_id();
+      };
 
-      return construction_data;
+      return CreateDescriptionFromTriangulationHelper<dim, spacedim>(
+               tria,
+               subdomain_id_function,
+               level_subdomain_id_function,
+               comm,
+               settings)
+        .template create_description_for_rank<Description<dim, spacedim>>(
+          my_rank);
     }
 
 
@@ -522,211 +866,6 @@ namespace TriangulationDescription
 #endif
     }
 
-    template <int dim, int spacedim>
-    struct DescriptionTemp
-    {
-      template <class Archive>
-      void
-      serialize(Archive &ar, const unsigned int /*version*/)
-      {
-        ar &coarse_cells;
-        ar &coarse_cell_vertices;
-        ar &coarse_cell_index_to_coarse_cell_id;
-        ar &cell_infos;
-      }
-
-      void
-      merge(const DescriptionTemp<dim, spacedim> &other)
-      {
-        this->cell_infos.resize(
-          std::max(other.cell_infos.size(), this->cell_infos.size()));
-
-        this->coarse_cells.insert(this->coarse_cells.end(),
-                                  other.coarse_cells.begin(),
-                                  other.coarse_cells.end());
-        this->coarse_cell_vertices.insert(this->coarse_cell_vertices.end(),
-                                          other.coarse_cell_vertices.begin(),
-                                          other.coarse_cell_vertices.end());
-        this->coarse_cell_index_to_coarse_cell_id.insert(
-          this->coarse_cell_index_to_coarse_cell_id.end(),
-          other.coarse_cell_index_to_coarse_cell_id.begin(),
-          other.coarse_cell_index_to_coarse_cell_id.end());
-
-        for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
-          this->cell_infos[i].insert(this->cell_infos[i].end(),
-                                     other.cell_infos[i].begin(),
-                                     other.cell_infos[i].end());
-      }
-
-      void
-      reduce()
-      {
-        // make coarse cells unique
-        {
-          std::vector<std::tuple<types::coarse_cell_id,
-                                 dealii::CellData<dim>,
-                                 unsigned int>>
-            temp;
-
-          for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
-            temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
-                              this->coarse_cells[i],
-                              i);
-
-          std::sort(temp.begin(), temp.end(), [](const auto &a, const auto &b) {
-            return std::get<0>(a) < std::get<0>(b);
-          });
-          temp.erase(std::unique(temp.begin(),
-                                 temp.end(),
-                                 [](const auto &a, const auto &b) {
-                                   return std::get<0>(a) == std::get<0>(b);
-                                 }),
-                     temp.end());
-          std::sort(temp.begin(), temp.end(), [](const auto &a, const auto &b) {
-            return std::get<2>(a) < std::get<2>(b);
-          });
-
-          this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
-          this->coarse_cells.resize(temp.size());
-
-          for (unsigned int i = 0; i < temp.size(); ++i)
-            {
-              this->coarse_cell_index_to_coarse_cell_id[i] =
-                std::get<0>(temp[i]);
-              this->coarse_cells[i] = std::get<1>(temp[i]);
-            }
-        }
-
-        // make coarse cell vertices unique
-        {
-          std::sort(this->coarse_cell_vertices.begin(),
-                    this->coarse_cell_vertices.end(),
-                    [](const auto &a, const auto &b) {
-                      return a.first < b.first;
-                    });
-          this->coarse_cell_vertices.erase(
-            std::unique(this->coarse_cell_vertices.begin(),
-                        this->coarse_cell_vertices.end(),
-                        [](const auto &a, const auto &b) {
-                          return a.first == b.first;
-                        }),
-            this->coarse_cell_vertices.end());
-        }
-
-        // make cells unique
-        for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
-          {
-            std::sort(this->cell_infos[i].begin(),
-                      this->cell_infos[i].end(),
-                      [](const auto &a, const auto &b) { return a.id < b.id; });
-            this->cell_infos[i].erase(std::unique(this->cell_infos[i].begin(),
-                                                  this->cell_infos[i].end(),
-                                                  [](const auto &a,
-                                                     const auto &b) {
-                                                    return a.id == b.id;
-                                                  }),
-                                      this->cell_infos[i].end());
-          }
-      }
-
-      Description<dim, spacedim>
-      convert(const MPI_Comm comm,
-              const typename Triangulation<dim, spacedim>::MeshSmoothing
-                mesh_smoothing)
-      {
-        Description<dim, spacedim> description;
-
-        // copy communicator
-        description.comm = comm;
-
-        // set settings (no multigrid levels are set for now)
-        description.settings =
-          TriangulationDescription::Settings::default_setting;
-
-        // use mesh smoothing from base triangulation
-        description.smoothing = mesh_smoothing;
-
-        std::map<unsigned int, unsigned int> map;
-
-        for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
-          {
-            description.coarse_cell_vertices.push_back(
-              this->coarse_cell_vertices[i].second);
-            map[this->coarse_cell_vertices[i].first] = i;
-          }
-
-        description.coarse_cells = this->coarse_cells;
-
-        for (auto &cell : description.coarse_cells)
-          for (unsigned int v = 0; v < cell.vertices.size(); ++v)
-            cell.vertices[v] = map[cell.vertices[v]];
-
-        description.coarse_cell_index_to_coarse_cell_id =
-          this->coarse_cell_index_to_coarse_cell_id;
-        description.cell_infos = this->cell_infos;
-
-        return description;
-      }
-
-      std::vector<dealii::CellData<dim>> coarse_cells;
-
-      std::vector<std::pair<unsigned int, Point<spacedim>>>
-        coarse_cell_vertices;
-
-      std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
-
-      std::vector<std::vector<CellData<dim>>> cell_infos;
-    };
-
-
-
-    template <int dim, int spacedim>
-    void
-    fill_cell_infos(const TriaIterator<CellAccessor<dim, spacedim>> &cell,
-                    std::vector<std::vector<CellData<dim>>> &        cell_infos,
-                    const LinearAlgebra::distributed::Vector<double> &partition)
-    {
-      CellData<dim> cell_info;
-
-      // save coarse-cell id
-      cell_info.id = cell->id().template to_binary<dim>();
-
-      // save boundary_ids of each face of this cell
-      for (const auto f : cell->face_indices())
-        {
-          types::boundary_id boundary_ind = cell->face(f)->boundary_id();
-          if (boundary_ind != numbers::internal_face_boundary_id)
-            cell_info.boundary_ids.emplace_back(f, boundary_ind);
-        }
-
-      // save manifold id
-      {
-        // ... of cell
-        cell_info.manifold_id = cell->manifold_id();
-
-        // ... of lines
-        if (dim >= 2)
-          for (const auto line : cell->line_indices())
-            cell_info.manifold_line_ids[line] = cell->line(line)->manifold_id();
-
-        // ... of quads
-        if (dim == 3)
-          for (const auto f : cell->face_indices())
-            cell_info.manifold_quad_ids[f] = cell->quad(f)->manifold_id();
-      }
-
-      // subdomain and level subdomain id
-      if (cell->is_active() && (cell->is_artificial() == false))
-        cell_info.subdomain_id = static_cast<unsigned int>(
-          partition[cell->global_active_cell_index()]);
-      else
-        cell_info.subdomain_id = numbers::artificial_subdomain_id;
-
-      cell_info.level_subdomain_id = numbers::artificial_subdomain_id;
-
-      cell_infos[cell->level()].emplace_back(cell_info);
-    }
-
 
 
     template <int dim, int spacedim>
@@ -741,8 +880,10 @@ namespace TriangulationDescription
 #endif
 
       if (partition.size() == 0)
-        return create_description_from_triangulation(tria,
-                                                     tria.get_communicator());
+        {
+          return create_description_from_triangulation(tria,
+                                                       tria.get_communicator());
+        }
 
       partition.update_ghost_values();
 
@@ -758,164 +899,52 @@ namespace TriangulationDescription
                                          relevant_processes.end());
       }();
 
-      // 2) determine coinciding vertices (important for periodicity)
-      std::map<unsigned int, std::vector<unsigned int>>
-                                           coinciding_vertex_groups;
-      std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
-
-      GridTools::collect_coinciding_vertices(tria,
-                                             coinciding_vertex_groups,
-                                             vertex_to_coinciding_vertex_group);
-
-      const auto add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
-        [&coinciding_vertex_groups, &vertex_to_coinciding_vertex_group](
-          const auto &       cell,
-          std::vector<bool> &vertices_owned_by_locally_owned_cells) {
-          for (const auto v : cell->vertex_indices())
-            {
-              vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] =
-                true;
-              const auto coinciding_vertex_group =
-                vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
-              if (coinciding_vertex_group !=
-                  vertex_to_coinciding_vertex_group.end())
-                for (const auto &co_vertex : coinciding_vertex_groups.at(
-                       coinciding_vertex_group->second))
-                  vertices_owned_by_locally_owned_cells[co_vertex] = true;
-            }
-        };
-
-      // 3) create a description (locally owned cell and a layer of ghost cells
-      // and all their parents)
-      std::vector<DescriptionTemp<dim, spacedim>> description_temp(
-        relevant_processes.size());
+      TriangulationDescription::Settings settings =
+        TriangulationDescription::Settings::default_setting;
 
-      std::vector<bool> old_user_flags;
-      if (description_temp.size() > 0)
-        tria.save_user_flags(old_user_flags);
-
-      for (unsigned int i = 0; i < description_temp.size(); ++i)
-        {
-          const unsigned int proc               = relevant_processes[i];
-          auto &             description_temp_i = description_temp[i];
-          description_temp_i.cell_infos.resize(
-            tria.get_triangulation().n_global_levels());
-
-          // clear user_flags
-          const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
-            .clear_user_flags();
-
-          for (int level = tria.get_triangulation().n_global_levels() - 1;
-               level >= 0;
-               --level)
-            {
-              // mark all vertices attached to locally owned cells
-              std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
-                tria.n_vertices());
-              for (const auto &cell : tria.active_cell_iterators())
-                if (cell->is_active() && cell->is_locally_owned() &&
-                    static_cast<unsigned int>(
-                      partition[cell->global_active_cell_index()]) == proc)
-                  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
-                    cell, vertices_owned_by_locally_owned_cells_on_level);
-
-              // helper function if a cell is locally relevant (active and
-              // connected to cell via a vertex)
-              const auto is_locally_relevant_on_level = [&](const auto &cell) {
-                for (const auto v : cell->vertex_indices())
-                  if (vertices_owned_by_locally_owned_cells_on_level
-                        [cell->vertex_index(v)])
-                    return true;
-                return false;
-              };
-
-              // collect locally relevant cells (including their parents)
-              for (const auto &cell : tria.cell_iterators_on_level(level))
-                if (is_locally_relevant_on_level(cell))
-                  set_user_flag_and_of_its_parents(cell);
-            }
-
-          for (unsigned int level = 0; level < tria.n_global_levels(); ++level)
-            for (const auto &cell : tria.cell_iterators_on_level(level))
-              if (cell->user_flag_set())
-                fill_cell_infos(cell, description_temp_i.cell_infos, partition);
-
-          // collect coarse-grid cells
-          std::vector<bool> vertices_locally_relevant(tria.n_vertices(), false);
+      const auto subdomain_id_function = [&partition](const auto &cell) {
+        if ((cell->is_active() && (cell->is_artificial() == false)))
+          return static_cast<unsigned int>(
+            partition[cell->global_active_cell_index()]);
+        else
+          return numbers::artificial_subdomain_id;
+      };
 
-          for (const auto &cell : tria.cell_iterators_on_level(0))
-            {
-              if (!cell->user_flag_set())
-                continue;
+      const auto level_subdomain_id_function = [](const auto & /*cell*/) {
+        return numbers::artificial_subdomain_id;
+      };
 
-              // extract cell definition
-              dealii::CellData<dim> cell_data(cell->n_vertices());
-              cell_data.material_id = cell->material_id();
-              cell_data.manifold_id = cell->manifold_id();
-              for (const auto v : cell->vertex_indices())
-                cell_data.vertices[v] = cell->vertex_index(v);
-              description_temp_i.coarse_cells.push_back(cell_data);
-
-              // mark cell vertices as relevant
-              for (const auto v : cell->vertex_indices())
-                vertices_locally_relevant[cell->vertex_index(v)] = true;
-
-              // save translation for corase grid: lid -> gid
-              description_temp_i.coarse_cell_index_to_coarse_cell_id.push_back(
-                cell->id().get_coarse_cell_id());
-            }
+      CreateDescriptionFromTriangulationHelper<dim, spacedim> helper(
+        tria,
+        subdomain_id_function,
+        level_subdomain_id_function,
+        tria.get_communicator(),
+        settings);
 
-          // collect coarse-grid vertices
-          for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
-            if (vertices_locally_relevant[i])
-              description_temp_i.coarse_cell_vertices.emplace_back(
-                i, tria.get_vertices()[i]);
-        }
+      // create a description (locally owned cell and a layer of ghost cells
+      // and all their parents)
+      std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
+      descriptions_per_rank.reserve(relevant_processes.size());
 
-      // restore flags
-      if (description_temp.size() > 0)
-        const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
-          .load_user_flags(old_user_flags);
+      for (const auto rank : relevant_processes)
+        descriptions_per_rank.emplace_back(
+          helper.template create_description_for_rank<
+            DescriptionTemp<dim, spacedim>>(rank));
 
       // collect description from all processes that used to own locally-owned
       // active cells of this process in a single description
       DescriptionTemp<dim, spacedim> description_merged;
-
-      dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess<char, char>
-        process([&]() { return relevant_processes; },
-                [&](const unsigned int other_rank,
-                    std::vector<char> &send_buffer) {
-                  const auto ptr = std::find(relevant_processes.begin(),
-                                             relevant_processes.end(),
-                                             other_rank);
-
-                  Assert(ptr != relevant_processes.end(), ExcInternalError());
-
-                  const auto other_rank_index =
-                    std::distance(relevant_processes.begin(), ptr);
-
-                  send_buffer =
-                    dealii::Utilities::pack(description_temp[other_rank_index],
-                                            false);
-                },
-                [&](const unsigned int &,
-                    const std::vector<char> &recv_buffer,
-                    std::vector<char> &) {
-                  description_merged.merge(
-                    dealii::Utilities::unpack<DescriptionTemp<dim, spacedim>>(
-                      recv_buffer, false));
-                });
-
-      dealii::Utilities::MPI::ConsensusAlgorithms::Selector<char, char>(
-        process, partition.get_mpi_communicator())
-        .run();
+      description_merged.collect(relevant_processes,
+                                 descriptions_per_rank,
+                                 partition.get_mpi_communicator());
 
       // remove redundant entries
       description_merged.reduce();
 
       // convert to actual description
       return description_merged.convert(partition.get_mpi_communicator(),
-                                        tria.get_mesh_smoothing());
+                                        tria.get_mesh_smoothing(),
+                                        settings);
     }
 
   } // namespace Utilities