From e5a032b87b499101d40ca79d9299acbbac5adf2a Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 6 Mar 2021 09:26:39 +0100
Subject: [PATCH] Enable set up of active_cell_index_partitioner for p:s:T

---
 include/deal.II/grid/tria.h     |  7 +++++-
 source/distributed/tria_base.cc | 43 ++++++++++++++++++++++++++++++---
 2 files changed, 45 insertions(+), 5 deletions(-)

diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index a840993ab3..edb113f62a 100644
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -4234,8 +4234,13 @@ Triangulation<dim, spacedim>::load(Archive &ar, const unsigned int)
   // this here. don't forget to first resize the fields appropriately
   {
     for (auto &level : levels)
-      level->active_cell_indices.resize(level->refine_flags.size());
+      {
+        level->active_cell_indices.resize(level->refine_flags.size());
+        level->global_active_cell_indices.resize(level->refine_flags.size());
+        level->global_level_cell_indices.resize(level->refine_flags.size());
+      }
     reset_active_cell_indices();
+    reset_global_cell_indices();
   }
 
 
diff --git a/source/distributed/tria_base.cc b/source/distributed/tria_base.cc
index b379bb8ffa..0abf54bfc5 100644
--- a/source/distributed/tria_base.cc
+++ b/source/distributed/tria_base.cc
@@ -19,6 +19,7 @@
 #include <deal.II/base/mpi.templates.h>
 #include <deal.II/base/utilities.h>
 
+#include <deal.II/distributed/shared_tria.h>
 #include <deal.II/distributed/tria_base.h>
 
 #include <deal.II/grid/grid_tools.h>
@@ -380,11 +381,45 @@ namespace parallel
 #ifndef DEAL_II_WITH_MPI
     Assert(false, ExcNeedsMPI());
 #else
+    if (const auto pst =
+          dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
+            this))
+      if (pst->with_artificial_cells() == false)
+        {
+          // Specialization for parallel::shared::Triangulation without
+          // artificial cells. The code below only works if a halo of a single
+          // ghost cells is needed.
+
+          std::vector<unsigned int> cell_counter(n_subdomains + 1);
+
+          // count number of cells of each process
+          for (const auto &cell : this->active_cell_iterators())
+            cell_counter[cell->subdomain_id() + 1]++;
+
+          // take prefix sum to obtain offset of each process
+          for (unsigned int i = 0; i < n_subdomains; ++i)
+            cell_counter[i + 1] += cell_counter[i];
+
+          // create partitioners
+          IndexSet is_local(cell_counter.back());
+          is_local.add_range(cell_counter[my_subdomain],
+                             cell_counter[my_subdomain + 1]);
+          IndexSet is_ghost(cell_counter.back());
+          number_cache.active_cell_index_partitioner =
+            Utilities::MPI::Partitioner(is_local,
+                                        is_ghost,
+                                        this->mpi_communicator);
+
+          // set global active cell indices and increment process-local counters
+          for (const auto &cell : this->active_cell_iterators())
+            cell->set_global_active_cell_index(
+              cell_counter[cell->subdomain_id()]++);
 
-    // currently only implemented for distributed triangulations
-    if (dynamic_cast<const parallel::DistributedTriangulationBase<dim, spacedim>
-                       *>(this) == nullptr)
-      return;
+          Assert(this->is_multilevel_hierarchy_constructed() == false,
+                 ExcNotImplemented());
+
+          return;
+        }
 
     // 1) determine number of active locally-owned cells
     const types::global_cell_index n_locally_owned_cells =
-- 
2.39.5