Make mapping in GridTools::distributed_compute_point_locations() unique

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 7b973c2f3688c1dd27035e6179483288de8e5549..240c194d18757f90247c10e2d60f5822a74632ce 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -5640,7 +5640,7 @@ namespace GridTools
   {
     // run internal function ...
     const auto all = internal::distributed_compute_point_locations(
-                       cache, points, global_bboxes, tolerance, false)
+                       cache, points, global_bboxes, tolerance, false, true)
                        .send_components;
 
     // ... and reshuffle the data
@@ -5749,7 +5749,8 @@ namespace GridTools
       const Point<spacedim> &                                      point,
       typename Triangulation<dim, spacedim>::active_cell_iterator &cell_hint,
       const std::vector<bool> &marked_vertices,
-      const double             tolerance)
+      const double             tolerance,
+      const bool               enforce_unique_mapping)
     {
       std::vector<
         std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
@@ -5768,6 +5769,9 @@ namespace GridTools
         tolerance,
         &cache.get_locally_owned_cell_bounding_boxes_rtree());
 
+      const unsigned int my_rank = Utilities::MPI::this_mpi_process(
+        cache.get_triangulation().get_communicator());
+
       cell_hint = first_cell.first;
       if (cell_hint.state() == IteratorState::valid)
         {
@@ -5779,6 +5783,20 @@ namespace GridTools
               tolerance,
               first_cell);
 
+          if (enforce_unique_mapping)
+            {
+              // check if the rank of this process is the lowest of all cells
+              // if not, the other process will handle this cell and we don't
+              // have to do here anything in the case of unique mapping
+              unsigned int lowes_rank = numbers::invalid_unsigned_int;
+
+              for (const auto &cell : active_cells_around_point)
+                lowes_rank = std::min(lowes_rank, cell.first->subdomain_id());
+
+              if (lowes_rank != my_rank)
+                return {};
+            }
+
           locally_owned_active_cells_around_point.reserve(
             active_cells_around_point.size());
 
@@ -5791,7 +5809,12 @@ namespace GridTools
                 locally_owned_active_cells_around_point.end(),
                 [](const auto &a, const auto &b) { return a.first < b.first; });
 
-      return locally_owned_active_cells_around_point;
+      if (enforce_unique_mapping &&
+          locally_owned_active_cells_around_point.size() > 1)
+        // in the case of unique mapping, we only need a single cell
+        return {locally_owned_active_cells_around_point.front()};
+      else
+        return locally_owned_active_cells_around_point;
     }
 
 
@@ -5806,8 +5829,6 @@ namespace GridTools
       const bool                                             perform_handshake,
       const bool enforce_unique_mapping)
     {
-      Assert(!enforce_unique_mapping || perform_handshake, ExcInternalError());
-
       DistributedComputePointLocationsInternal<dim, spacedim> result;
 
       auto &send_components = result.send_components;
@@ -5879,7 +5900,8 @@ namespace GridTools
                   index_and_point.second,
                   cell_hint,
                   marked_vertices,
-                  tolerance);
+                  tolerance,
+                  enforce_unique_mapping);
 
               for (const auto &cell_and_reference_position :
                    cells_and_reference_positions)
@@ -5893,18 +5915,9 @@ namespace GridTools
                     cell_and_reference_position.second,
                     index_and_point.second,
                     numbers::invalid_unsigned_int);
-
-                  if (enforce_unique_mapping)
-                    break; // in the case of unique mapping, we only need a
-                           // single cell (we take the first)
                 }
 
-              if (perform_handshake)
-                request_buffer_temp[i] =
-                  enforce_unique_mapping ?
-                    std::min<unsigned int>(
-                      1, cells_and_reference_positions.size()) :
-                    cells_and_reference_positions.size();
+              request_buffer_temp[i] = cells_and_reference_positions.size();
             }
 
           if (perform_handshake)
@@ -5946,143 +5959,6 @@ namespace GridTools
         process, cache.get_triangulation().get_communicator())
         .run();
 
-      // for unique mapping, we need to modify recv_components and
-      // send_components consistently
-      if (enforce_unique_mapping)
-        {
-          std::vector<unsigned int> mask_recv(recv_components.size());
-          std::vector<unsigned int> mask_send(send_components.size());
-
-          // set up new recv_components and monitor which entries (we keep
-          // the entry of the lowest rank) have been eliminated so that we can
-          // communicate it
-          auto recv_components_copy = recv_components;
-          recv_components.clear();
-
-          for (unsigned int i = 0; i < recv_components_copy.size(); ++i)
-            std::get<2>(recv_components_copy[i]) = i;
-
-          std::sort(recv_components_copy.begin(),
-                    recv_components_copy.end(),
-                    [&](const auto &a, const auto &b) {
-                      if (std::get<0>(a) != std::get<0>(b)) // rank
-                        return std::get<0>(a) < std::get<0>(b);
-
-                      return std::get<2>(a) < std::get<2>(b); // enumeration
-                    });
-
-          std::vector<bool> unique(points.size(), false);
-
-          std::vector<unsigned int> recv_ranks;
-          std::vector<unsigned int> recv_ptrs;
-
-          for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
-               i < recv_components_copy.size();
-               ++i)
-            {
-              if (dummy != std::get<0>(recv_components_copy[i]))
-                {
-                  dummy = std::get<0>(recv_components_copy[i]);
-                  recv_ranks.push_back(dummy);
-                  recv_ptrs.push_back(i);
-                }
-
-              if (unique[std::get<1>(recv_components_copy[i])] == false)
-                {
-                  recv_components.emplace_back(recv_components_copy[i]);
-                  mask_recv[i]                                 = 1;
-                  unique[std::get<1>(recv_components_copy[i])] = true;
-                }
-              else
-                {
-                  mask_recv[i] = 0;
-                }
-            }
-          recv_ptrs.push_back(recv_components_copy.size());
-
-          Assert(std::all_of(unique.begin(),
-                             unique.end(),
-                             [](const auto &v) { return v; }),
-                 ExcInternalError());
-
-
-          // prepare send_components so that not needed entries can be
-          // eliminated later on
-          auto send_components_copy = send_components;
-          send_components.clear();
-
-          for (unsigned int i = 0; i < send_components_copy.size(); ++i)
-            std::get<5>(send_components_copy[i]) = i;
-
-          std::sort(send_components_copy.begin(),
-                    send_components_copy.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);
-
-                      return std::get<5>(a) < std::get<5>(b); // enumeration
-                    });
-
-          std::vector<unsigned int> send_ranks;
-          std::vector<unsigned int> send_ptrs;
-
-          for (unsigned int i = 0, dummy = numbers::invalid_unsigned_int;
-               i < send_components_copy.size();
-               ++i)
-            {
-              if (dummy != std::get<1>(send_components_copy[i]))
-                {
-                  dummy = std::get<1>(send_components_copy[i]);
-                  send_ranks.push_back(dummy);
-                  send_ptrs.push_back(i);
-                }
-            }
-          send_ptrs.push_back(send_components_copy.size());
-
-          // perform communication
-#ifdef DEAL_II_WITH_MPI
-          std::vector<MPI_Request> req(send_ranks.size() + recv_ranks.size());
-
-          for (unsigned int i = 0; i < send_ranks.size(); ++i)
-            {
-              const auto ierr =
-                MPI_Irecv(mask_send.data() + send_ptrs[i],
-                          send_ptrs[i + 1] - send_ptrs[i],
-                          MPI_UNSIGNED,
-                          send_ranks[i],
-                          Utilities::MPI::internal::Tags::
-                            distributed_compute_point_locations,
-                          cache.get_triangulation().get_communicator(),
-                          &req[i]);
-              AssertThrowMPI(ierr);
-            }
-
-          for (unsigned int i = 0; i < recv_ranks.size(); ++i)
-            {
-              const auto ierr =
-                MPI_Isend(mask_recv.data() + recv_ptrs[i],
-                          recv_ptrs[i + 1] - recv_ptrs[i],
-                          MPI_UNSIGNED,
-                          recv_ranks[i],
-                          Utilities::MPI::internal::Tags::
-                            distributed_compute_point_locations,
-                          cache.get_triangulation().get_communicator(),
-                          &req[i] + send_ranks.size());
-              AssertThrowMPI(ierr);
-            }
-
-          auto ierr = MPI_Waitall(req.size(), req.data(), MPI_STATUSES_IGNORE);
-          AssertThrowMPI(ierr);
-#else
-          mask_send = mask_recv;
-#endif
-
-          // eliminate not needed entries
-          for (unsigned int i = 0; i < send_components_copy.size(); ++i)
-            if (mask_send[i] == 1)
-              send_components.emplace_back(send_components_copy[i]);
-        }
-
       if (true)
         {
           // sort according to rank (and point index and cell) -> make

diff --git a/tests/remote_point_evaluation/mapping_02.cc b/tests/remote_point_evaluation/mapping_02.cc
new file mode 100644 (file)
index 0000000..29f4e9d
--- /dev/null
+++ b/tests/remote_point_evaluation/mapping_02.cc
@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Test GridTools::distributed_compute_point_locations() that it works with
+// unique mapping.
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_vector.h>
+
+#include <deal.II/numerics/vector_tools_evaluate.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+
+template <int dim>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_cube(tria, 2);
+  tria.refine_global(1);
+
+  FE_DGQ<dim>     fe(0);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> vec(dof_handler.n_dofs());
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    if (cell->is_locally_owned())
+      vec[cell->global_active_cell_index()] = cell->global_active_cell_index();
+
+  MappingQ1<dim> mapping;
+
+  std::vector<BoundingBox<dim>> local_boxes;
+  for (const auto &cell : tria.cell_iterators_on_level(0))
+    local_boxes.push_back(mapping.get_bounding_box(cell));
+
+  // create r-tree of bounding boxes
+  const auto local_tree = pack_rtree(local_boxes);
+
+  // compress r-tree to a minimal set of bounding boxes
+  const auto local_reduced_box = extract_rtree_level(local_tree, 0);
+
+  // gather bounding boxes of other processes
+  const auto global_bboxes =
+    Utilities::MPI::all_gather(tria.get_communicator(), local_reduced_box);
+
+  const GridTools::Cache<dim> cache(tria, mapping);
+
+  std::vector<Point<dim>> local_points;
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    for (unsigned int j = 0; j <= 4; ++j)
+      for (unsigned int i = 0; i <= 4; ++i)
+        local_points.emplace_back(i * 0.25, j * 0.25);
+
+  const auto temp =
+    GridTools::distributed_compute_point_locations(cache,
+                                                   local_points,
+                                                   global_bboxes);
+
+  unsigned int count = 0;
+
+  for (const auto &i : std::get<1>(temp))
+    count += i.size();
+
+  const auto global_count = Utilities::MPI::sum(count, MPI_COMM_WORLD);
+  const auto expected_global_count =
+    Utilities::MPI::sum(local_points.size(), MPI_COMM_WORLD);
+
+  deallog << global_count << " " << expected_global_count << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  initlog();
+
+  test<2>();
+}

diff --git a/tests/remote_point_evaluation/mapping_02.mpirun=1.output b/tests/remote_point_evaluation/mapping_02.mpirun=1.output
new file mode 100644 (file)
index 0000000..62654b9
--- /dev/null
+++ b/tests/remote_point_evaluation/mapping_02.mpirun=1.output
@@ -0,0 +1,2 @@
+
+DEAL::25 25

diff --git a/tests/remote_point_evaluation/mapping_02.mpirun=4.output b/tests/remote_point_evaluation/mapping_02.mpirun=4.output
new file mode 100644 (file)
index 0000000..62654b9
--- /dev/null
+++ b/tests/remote_point_evaluation/mapping_02.mpirun=4.output
@@ -0,0 +1,2 @@
+
+DEAL::25 25