From: Johannes Heinz <johannes.heinz@tuwien.ac.at>
Date: Fri, 28 Apr 2023 09:02:36 +0000 (+0200)
Subject: make it possible to fill RPE with GT::DistributedComputePointLocationsInternal
X-Git-Tag: v9.5.0-rc1~254^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=042dce82a5f9de0e30cc722de84744ddcabe692b;p=dealii.git

make it possible to fill RPE with GT::DistributedComputePointLocationsInternal
---

diff --git a/include/deal.II/base/mpi_remote_point_evaluation.h b/include/deal.II/base/mpi_remote_point_evaluation.h
index 7dfbf220f9..b357d28f6d 100644
--- a/include/deal.II/base/mpi_remote_point_evaluation.h
+++ b/include/deal.II/base/mpi_remote_point_evaluation.h
@@ -25,6 +25,15 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+// Forward declarations
+namespace GridTools
+{
+  namespace internal
+  {
+    template <int dim, int spacedim>
+    struct DistributedComputePointLocationsInternal;
+  }
+} // namespace GridTools
 
 namespace Utilities
 {
@@ -86,6 +95,21 @@ namespace Utilities
              const Triangulation<dim, spacedim> &tria,
              const Mapping<dim, spacedim> &      mapping);
 
+      /**
+       * Set up internal data structures and communication pattern based on
+       * GridTools::internal::DistributedComputePointLocationsInternal.
+       *
+       * This function is called internally by the reinit() function above.
+       * Having it as a seperate function makes it possible to setup the class
+       * if it is known in which cells corresponding reference points are
+       * located (e.g. if intersections of cells are known).
+       */
+      void
+      reinit(const GridTools::internal::
+               DistributedComputePointLocationsInternal<dim, spacedim> &data,
+             const Triangulation<dim, spacedim> &                       tria,
+             const Mapping<dim, spacedim> &mapping);
+
       /**
        * Data of points positioned in a cell.
        */
diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 7bf25de6f9..d48e80213d 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -1213,6 +1213,22 @@ namespace GridTools
     template <int dim, int spacedim>
     struct DistributedComputePointLocationsInternal
     {
+      DistributedComputePointLocationsInternal();
+
+      /**
+       * Function which sets up @p send_ranks, @p send_ptrs, @p recv_ranks,
+       * and @p recv_ptrs from @p send_components, @p recv_components,
+       * and @p n_searched_points. Internally @p send_components and @p recv_components
+       * are sorted and enumerated.
+       */
+      void
+      finalize_setup();
+
+      /**
+       * Number of searched point locations.
+       */
+      unsigned int n_searched_points;
+
       /**
        * Information of each point on sending/evaluation side. The elements of
        * the tuple are as follows: 0) cell level and index, 1) rank of the
diff --git a/source/base/mpi_remote_point_evaluation.cc b/source/base/mpi_remote_point_evaluation.cc
index 656be3e2a6..43565aa536 100644
--- a/source/base/mpi_remote_point_evaluation.cc
+++ b/source/base/mpi_remote_point_evaluation.cc
@@ -78,9 +78,6 @@ namespace Utilities
       tria_signal =
         tria.signals.any_change.connect([&]() { this->ready_flag = false; });
 
-      this->tria    = &tria;
-      this->mapping = &mapping;
-
       std::vector<BoundingBox<spacedim>> local_boxes;
       for (const auto &cell :
            tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
@@ -105,6 +102,23 @@ namespace Utilities
           true,
           enforce_unique_mapping);
 
+      this->reinit(data, tria, mapping);
+#endif
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    RemotePointEvaluation<dim, spacedim>::reinit(
+      const GridTools::internal::
+        DistributedComputePointLocationsInternal<dim, spacedim> &data,
+      const Triangulation<dim, spacedim> &                       tria,
+      const Mapping<dim, spacedim> &                             mapping)
+    {
+      this->tria    = &tria;
+      this->mapping = &mapping;
+
       this->recv_ranks = data.recv_ranks;
       this->recv_ptrs  = data.recv_ptrs;
 
@@ -113,7 +127,7 @@ namespace Utilities
 
       this->recv_permutation = {};
       this->recv_permutation.resize(data.recv_components.size());
-      this->point_ptrs.assign(points.size() + 1, 0);
+      this->point_ptrs.assign(data.n_searched_points + 1, 0);
       for (unsigned int i = 0; i < data.recv_components.size(); ++i)
         {
           AssertIndexRange(std::get<2>(data.recv_components[i]),
@@ -130,7 +144,7 @@ namespace Utilities
       std::tuple<unsigned int, unsigned int> n_owning_processes_local =
         n_owning_processes_default;
 
-      for (unsigned int i = 0; i < points.size(); ++i)
+      for (unsigned int i = 0; i < data.n_searched_points; ++i)
         {
           std::get<0>(n_owning_processes_local) =
             std::min(std::get<0>(n_owning_processes_local),
@@ -196,7 +210,6 @@ namespace Utilities
         cell_data.reference_point_values.size());
 
       this->ready_flag = true;
-#endif
     }
 
 
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index dfce19ff5d..8afe92502e 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -6168,6 +6168,127 @@ namespace GridTools
         return locally_owned_active_cells_around_point;
     }
 
+    template <int dim, int spacedim>
+    DistributedComputePointLocationsInternal<dim, spacedim>::
+      DistributedComputePointLocationsInternal()
+      : n_searched_points(numbers::invalid_unsigned_int)
+    {}
+
+    template <int dim, int spacedim>
+    void
+    DistributedComputePointLocationsInternal<dim, spacedim>::finalize_setup()
+    {
+      // before reshuffeling the data check if data.recv_components and
+      // n_searched_points are in a valid state.
+      Assert(n_searched_points != numbers::invalid_unsigned_int,
+             ExcInternalError());
+      Assert(recv_components.empty() ||
+               std::get<1>(*std::max_element(recv_components.begin(),
+                                             recv_components.end(),
+                                             [](const auto &a, const auto &b) {
+                                               return std::get<1>(a) <
+                                                      std::get<1>(b);
+                                             })) < n_searched_points,
+             ExcInternalError());
+
+      send_ranks.clear();
+      recv_ranks.clear();
+      send_ptrs.clear();
+      recv_ptrs.clear();
+
+      if (true)
+        {
+          // sort according to rank (and point index and cell) -> make
+          // deterministic
+          std::sort(send_components.begin(),
+                    send_components.end(),
+                    [&](const auto &a, const auto &b) {
+                      if (std::get<1>(a) != std::get<1>(b)) // rank
+                        return std::get<1>(a) < std::get<1>(b);
+
+                      if (std::get<2>(a) != std::get<2>(b)) // point index
+                        return std::get<2>(a) < std::get<2>(b);
+
+                      return std::get<0>(a) < std::get<0>(b); // cell
+                    });
+
+          // perform enumeration and extract rank information
+          for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
+               i < send_components.size();
+               ++i)
+            {
+              std::get<5>(send_components[i]) = i;
+
+              if (dummy != std::get<1>(send_components[i]))
+                {
+                  dummy = std::get<1>(send_components[i]);
+                  send_ranks.push_back(dummy);
+                  send_ptrs.push_back(i);
+                }
+            }
+          send_ptrs.push_back(send_components.size());
+
+          // sort according to cell, rank, point index (while keeping
+          // partial ordering)
+          std::sort(send_components.begin(),
+                    send_components.end(),
+                    [&](const auto &a, const auto &b) {
+                      if (std::get<0>(a) != std::get<0>(b))
+                        return std::get<0>(a) < std::get<0>(b); // cell
+
+                      if (std::get<1>(a) != std::get<1>(b))
+                        return std::get<1>(a) < std::get<1>(b); // rank
+
+                      if (std::get<2>(a) != std::get<2>(b))
+                        return std::get<2>(a) < std::get<2>(b); // point index
+
+                      return std::get<5>(a) < std::get<5>(b); // enumeration
+                    });
+        }
+
+      if (recv_components.size() > 0)
+        {
+          // sort according to rank (and point index) -> make deterministic
+          std::sort(recv_components.begin(),
+                    recv_components.end(),
+                    [&](const auto &a, const auto &b) {
+                      if (std::get<0>(a) != std::get<0>(b))
+                        return std::get<0>(a) < std::get<0>(b); // rank
+
+                      return std::get<1>(a) < std::get<1>(b); // point index
+                    });
+
+          // perform enumeration and extract rank information
+          for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
+               i < recv_components.size();
+               ++i)
+            {
+              std::get<2>(recv_components[i]) = i;
+
+              if (dummy != std::get<0>(recv_components[i]))
+                {
+                  dummy = std::get<0>(recv_components[i]);
+                  recv_ranks.push_back(dummy);
+                  recv_ptrs.push_back(i);
+                }
+            }
+          recv_ptrs.push_back(recv_components.size());
+
+          // sort according to point index and rank (while keeping partial
+          // ordering)
+          std::sort(recv_components.begin(),
+                    recv_components.end(),
+                    [&](const auto &a, const auto &b) {
+                      if (std::get<1>(a) != std::get<1>(b))
+                        return std::get<1>(a) < std::get<1>(b); // point index
+
+                      if (std::get<0>(a) != std::get<0>(b))
+                        return std::get<0>(a) < std::get<0>(b); // rank
+
+                      return std::get<2>(a) < std::get<2>(b); // enumeration
+                    });
+        }
+    }
 
 
     template <int dim, int spacedim>
@@ -6182,13 +6303,10 @@ namespace GridTools
       const bool enforce_unique_mapping)
     {
       DistributedComputePointLocationsInternal<dim, spacedim> result;
+      result.n_searched_points = points.size();
 
       auto &send_components = result.send_components;
-      auto &send_ranks      = result.send_ranks;
-      auto &send_ptrs       = result.send_ptrs;
       auto &recv_components = result.recv_components;
-      auto &recv_ranks      = result.recv_ranks;
-      auto &recv_ptrs       = result.recv_ptrs;
 
       const auto comm = cache.get_triangulation().get_communicator();
 
@@ -6330,98 +6448,7 @@ namespace GridTools
         process_answer,
         comm);
 
-      if (true)
-        {
-          // sort according to rank (and point index and cell) -> make
-          // deterministic
-          std::sort(send_components.begin(),
-                    send_components.end(),
-                    [&](const auto &a, const auto &b) {
-                      if (std::get<1>(a) != std::get<1>(b)) // rank
-                        return std::get<1>(a) < std::get<1>(b);
-
-                      if (std::get<2>(a) != std::get<2>(b)) // point index
-                        return std::get<2>(a) < std::get<2>(b);
-
-                      return std::get<0>(a) < std::get<0>(b); // cell
-                    });
-
-          // perform enumeration and extract rank information
-          for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
-               i < send_components.size();
-               ++i)
-            {
-              std::get<5>(send_components[i]) = i;
-
-              if (dummy != std::get<1>(send_components[i]))
-                {
-                  dummy = std::get<1>(send_components[i]);
-                  send_ranks.push_back(dummy);
-                  send_ptrs.push_back(i);
-                }
-            }
-          send_ptrs.push_back(send_components.size());
-
-          // sort according to cell, rank, point index (while keeping
-          // partial ordering)
-          std::sort(send_components.begin(),
-                    send_components.end(),
-                    [&](const auto &a, const auto &b) {
-                      if (std::get<0>(a) != std::get<0>(b))
-                        return std::get<0>(a) < std::get<0>(b); // cell
-
-                      if (std::get<1>(a) != std::get<1>(b))
-                        return std::get<1>(a) < std::get<1>(b); // rank
-
-                      if (std::get<2>(a) != std::get<2>(b))
-                        return std::get<2>(a) < std::get<2>(b); // point index
-
-                      return std::get<5>(a) < std::get<5>(b); // enumeration
-                    });
-        }
-
-      if (perform_handshake)
-        {
-          // sort according to rank (and point index) -> make deterministic
-          std::sort(recv_components.begin(),
-                    recv_components.end(),
-                    [&](const auto &a, const auto &b) {
-                      if (std::get<0>(a) != std::get<0>(b))
-                        return std::get<0>(a) < std::get<0>(b); // rank
-
-                      return std::get<1>(a) < std::get<1>(b); // point index
-                    });
-
-          // perform enumeration and extract rank information
-          for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
-               i < recv_components.size();
-               ++i)
-            {
-              std::get<2>(recv_components[i]) = i;
-
-              if (dummy != std::get<0>(recv_components[i]))
-                {
-                  dummy = std::get<0>(recv_components[i]);
-                  recv_ranks.push_back(dummy);
-                  recv_ptrs.push_back(i);
-                }
-            }
-          recv_ptrs.push_back(recv_components.size());
-
-          // sort according to point index and rank (while keeping partial
-          // ordering)
-          std::sort(recv_components.begin(),
-                    recv_components.end(),
-                    [&](const auto &a, const auto &b) {
-                      if (std::get<1>(a) != std::get<1>(b))
-                        return std::get<1>(a) < std::get<1>(b); // point index
-
-                      if (std::get<0>(a) != std::get<0>(b))
-                        return std::get<0>(a) < std::get<0>(b); // rank
-
-                      return std::get<2>(a) < std::get<2>(b); // enumeration
-                    });
-        }
+      result.finalize_setup();
 
       return result;
     }
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index ec06814639..9e5d8fe338 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -141,6 +141,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         const std::vector<bool> &marked_vertices,
         const bool               enforce_unique_mapping);
 
+      template struct internal::DistributedComputePointLocationsInternal<
+        deal_II_dimension,
+        deal_II_space_dimension>;
+
       template internal::DistributedComputePointLocationsInternal<
         deal_II_dimension,
         deal_II_space_dimension>