p::f::T: allow repartitioning

include/deal.II/distributed/fully_distributed_tria.h

@@ -19,6 +19,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/distributed/repartitioning_policy_tools.h>
 #include <deal.II/distributed/tria_base.h>
 
 #include <deal.II/grid/grid_tools.h>
@@ -196,6 +197,22 @@ namespace parallel
                                  const unsigned int)> &partitioner,
         const TriangulationDescription::Settings &     settings);
 
+      /**
+       * Register a partitioner, which is used within the method
+       * repartition().
+       */
+      void
+      set_partitioner(
+        const RepartitioningPolicyTools::Base<dim, spacedim> &partitioner,
+        const TriangulationDescription::Settings &            settings);
+
+      /**
+       * Execute repartitioning and use the partitioner attached by the
+       * method set_partitioner();
+       */
+      void
+      repartition();
+
       /**
        * Coarsen and refine the mesh according to refinement and coarsening
        * flags set.
@@ -301,6 +318,12 @@ namespace parallel
                          const unsigned int)>
         partitioner;
 
+      /**
+       * Partitioner used during repartition().
+       */
+      SmartPointer<const RepartitioningPolicyTools::Base<dim, spacedim>>
+        partitioner_distributed;
+
       /**
        * Sorted list of pairs of coarse-cell ids and their indices.
        */

include/deal.II/distributed/repartitioning_policy_tools.h

@@ -46,7 +46,7 @@ namespace RepartitioningPolicyTools
    * See the description of RepartitioningPolicyTools for more information.
    */
   template <int dim, int spacedim = dim>
-  class Base
+  class Base : public Subscriptor
   {
   public:
     /**

source/distributed/fully_distributed_tria.cc

@@ -18,6 +18,7 @@
 #include <deal.II/base/mpi.h>
 
 #include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/repartitioning_policy_tools.h>
 
 #include <deal.II/grid/grid_tools.h>
 
@@ -310,6 +311,46 @@ namespace parallel
 
 
 
+    template <int dim, int spacedim>
+    void
+    Triangulation<dim, spacedim>::set_partitioner(
+      const RepartitioningPolicyTools::Base<dim, spacedim> &partitioner,
+      const TriangulationDescription::Settings &            settings)
+    {
+      this->partitioner_distributed = &partitioner;
+      this->settings                = settings;
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    Triangulation<dim, spacedim>::repartition()
+    {
+      // signal that repartitioning has started
+      this->signals.pre_distributed_repartition();
+
+      // create construction_data with the help of the partitioner
+      const auto construction_data = TriangulationDescription::Utilities::
+        create_description_from_triangulation(
+          *this,
+          this->partitioner_distributed->partition(*this),
+          this->settings);
+
+      // clear old content
+      this->clear();
+      this->coarse_cell_id_to_coarse_cell_index_vector.clear();
+      this->coarse_cell_index_to_coarse_cell_id_vector.clear();
+
+      // use construction_data to set up new triangulation
+      this->create_triangulation(construction_data);
+
+      // signal that repartitioning has completed
+      this->signals.post_distributed_repartition();
+    }
+
+
+
     template <int dim, int spacedim>
     void
     Triangulation<dim, spacedim>::execute_coarsening_and_refinement()

source/grid/tria.cc

@@ -11871,11 +11871,8 @@ Triangulation<dim, spacedim>::create_triangulation(
 
           // boundary ids
           for (auto pair : cell_info->boundary_ids)
-            {
-              Assert(cell->at_boundary(pair.first),
-                     ExcMessage("Cell face is not on the boundary!"));
+            if (cell->face(pair.first)->at_boundary())
               cell->face(pair.first)->set_boundary_id(pair.second);
-            }
         }
     }
 }

source/grid/tria_description.cc

@@ -17,6 +17,7 @@
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/mpi_consensus_algorithms.h>
 
+#include <deal.II/distributed/fully_distributed_tria.h>
 #include <deal.II/distributed/tria.h>
 
 #include <deal.II/dofs/dof_accessor.h>
@@ -109,7 +110,8 @@ namespace TriangulationDescription
         collect(
           const std::vector<unsigned int> &                  relevant_processes,
           const std::vector<DescriptionTemp<dim, spacedim>> &description_temp,
-          const MPI_Comm &                                   comm)
+          const MPI_Comm &                                   comm,
+          const bool vertices_have_unique_ids)
         {
           dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess<char,
                                                                         char>
@@ -135,7 +137,8 @@ namespace TriangulationDescription
                   std::vector<char> &) {
                 this->merge(
                   dealii::Utilities::unpack<DescriptionTemp<dim, spacedim>>(
-                    recv_buffer, false));
+                    recv_buffer, false),
+                  vertices_have_unique_ids);
               });
 
           dealii::Utilities::MPI::ConsensusAlgorithms::Selector<char, char>(
@@ -149,17 +152,75 @@ namespace TriangulationDescription
          * This done by the reduce() function.
          */
         void
-        merge(const DescriptionTemp<dim, spacedim> &other)
+        merge(const DescriptionTemp<dim, spacedim> &other,
+              const bool                            vertices_have_unique_ids)
         {
           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());
+          if (vertices_have_unique_ids == false) // need to compare points
+            {
+              // comparator of points
+              const auto comp = [](const auto &a, const auto &b) {
+                const double tolerance = 1e-10;
+
+                for (unsigned int i = 0; i < spacedim; ++i)
+                  if (std::abs(a[i] - b[i]) > tolerance)
+                    return a[i] < b[i];
+
+                return false;
+              };
+
+              // map point to local vertex index
+              std::map<Point<spacedim>, unsigned int, decltype(comp)>
+                map_point_to_local_vertex_index(comp);
+
+              // ... initialize map with existing points
+              for (unsigned int i = 0; i < this->coarse_cell_vertices.size();
+                   ++i)
+                map_point_to_local_vertex_index[coarse_cell_vertices[i]
+                                                  .second] = i;
+
+              // map local vertex indices within other to the new local indices
+              std::map<unsigned int, unsigned int>
+                map_old_to_new_local_vertex_index;
+
+              // 1) renumerate vertices in other and insert into maps
+              unsigned int counter = coarse_cell_vertices.size();
+              for (const auto &p : other.coarse_cell_vertices)
+                if (map_point_to_local_vertex_index.find(p.second) ==
+                    map_point_to_local_vertex_index.end())
+                  {
+                    this->coarse_cell_vertices.emplace_back(counter, p.second);
+                    map_point_to_local_vertex_index[p.second] =
+                      map_old_to_new_local_vertex_index[p.first] = counter++;
+                  }
+                else
+                  map_old_to_new_local_vertex_index[p.first] =
+                    map_point_to_local_vertex_index[p.second];
+
+              // 2) renumerate vertices of cells
+              auto other_coarse_cells_copy = other.coarse_cells;
+
+              for (auto &cell : other_coarse_cells_copy)
+                for (auto &v : cell.vertices)
+                  v = map_old_to_new_local_vertex_index[v];
+
+              this->coarse_cells.insert(this->coarse_cells.end(),
+                                        other_coarse_cells_copy.begin(),
+                                        other_coarse_cells_copy.end());
+            }
+          else
+            {
+              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(),
@@ -228,7 +289,13 @@ namespace TriangulationDescription
               std::unique(this->coarse_cell_vertices.begin(),
                           this->coarse_cell_vertices.end(),
                           [](const auto &a, const auto &b) {
-                            return a.first == b.first;
+                            if (a.first == b.first)
+                              {
+                                Assert(a.second.distance(b.second) < 10e-8,
+                                       ExcInternalError());
+                                return true;
+                              }
+                            return false;
                           }),
               this->coarse_cell_vertices.end());
           }
@@ -1025,9 +1092,13 @@ namespace TriangulationDescription
       // 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;
-      description_merged.collect(relevant_processes,
-                                 descriptions_per_rank,
-                                 partition.get_mpi_communicator());
+      description_merged.collect(
+        relevant_processes,
+        descriptions_per_rank,
+        partition.get_mpi_communicator(),
+        dynamic_cast<
+          const parallel::fullydistributed::Triangulation<dim, spacedim> *>(
+          &tria) == nullptr);
 
       // remove redundant entries
       description_merged.reduce();