Merge pull request #12602 from peterrum/rpe_unique_fix

Fix tolerances in GridTools::internal::distributed_compute_point_locations()

source/base/mpi_remote_point_evaluation.cc

@@ -134,6 +134,9 @@ namespace Utilities
           this->point_ptrs[i + 1] += this->point_ptrs[i];
         }
 
+      Assert(enforce_unique_mapping == false || unique_mapping,
+             ExcInternalError());
+
       cell_data        = {};
       send_permutation = {};
 

source/grid/grid_tools.cc

@@ -2786,8 +2786,20 @@ namespace GridTools
     if (relevant_cell_bounding_boxes_rtree != nullptr &&
         !relevant_cell_bounding_boxes_rtree->empty())
       {
+        // create a bounding box around point p with 2*tolerance as side length.
+        auto p1 = p;
+        auto p2 = p;
+
+        for (int d = 0; d < spacedim; ++d)
+          {
+            p1[d] = p1[d] - tolerance;
+            p2[d] = p2[d] + tolerance;
+          }
+
+        BoundingBox<spacedim> bb({p1, p2});
+
         if (relevant_cell_bounding_boxes_rtree->qbegin(
-              boost::geometry::index::intersects(p)) ==
+              boost::geometry::index::intersects(bb)) ==
             relevant_cell_bounding_boxes_rtree->qend())
           return cell_and_position;
       }
@@ -5731,7 +5743,8 @@ namespace GridTools
                std::vector<unsigned int>>
     guess_point_owner(
       const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
-      const std::vector<Point<spacedim>> &                   points)
+      const std::vector<Point<spacedim>> &                   points,
+      const double                                           tolerance)
     {
       std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
       ranks_and_indices.reserve(points.size());
@@ -5741,7 +5754,7 @@ namespace GridTools
           const auto &point = points[i];
           for (unsigned rank = 0; rank < global_bboxes.size(); ++rank)
             for (const auto &box : global_bboxes[rank])
-              if (box.point_inside(point))
+              if (box.point_inside(point, tolerance))
                 {
                   ranks_and_indices.emplace_back(rank, i);
                   break;
@@ -5876,7 +5889,7 @@ namespace GridTools
       auto &recv_ptrs       = result.recv_ptrs;
 
       const auto potential_owners =
-        internal::guess_point_owner(global_bboxes, points);
+        internal::guess_point_owner(global_bboxes, points, tolerance);
 
       const auto &potential_owners_ranks   = std::get<0>(potential_owners);
       const auto &potential_owners_ptrs    = std::get<1>(potential_owners);

tests/remote_point_evaluation/mapping_03.cc

@@ -0,0 +1,93 @@
+// ---------------------------------------------------------------------
+//
+// 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 Utilities::MPI::RemotePointEvaluation with and without unique mapping
+// for a within a tolerance.
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_q.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(const bool enforce_unique_map)
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_cube(tria, 2);
+  tria.refine_global(1);
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> vec(dof_handler.n_dofs());
+
+  vec = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) + 1;
+
+  MappingQ1<dim> mapping;
+
+  std::vector<Point<dim>> evaluation_points;
+
+  evaluation_points.emplace_back(0.5, 0.5 + 1e-7); // should be assigned to rank
+                                                   // 0 for unique mapping
+
+  Utilities::MPI::RemotePointEvaluation<dim> eval(1e-6, enforce_unique_map);
+
+  const auto result_avg =
+    VectorTools::point_values<1>(mapping,
+                                 dof_handler,
+                                 vec,
+                                 evaluation_points,
+                                 eval,
+                                 VectorTools::EvaluationFlags::avg);
+  const auto result_min = VectorTools::point_values<1>(
+    eval, dof_handler, vec, VectorTools::EvaluationFlags::min);
+  const auto result_max = VectorTools::point_values<1>(
+    eval, dof_handler, vec, VectorTools::EvaluationFlags::max);
+  const auto result_insert = VectorTools::point_values<1>(
+    eval, dof_handler, vec, VectorTools::EvaluationFlags::insert);
+
+  for (unsigned int i = 0; i < evaluation_points.size(); ++i)
+    deallog << i << " " << result_avg[i] << " " << result_min[i] << " "
+            << result_max[i] << " " << result_insert[i] << " "
+            << eval.get_point_ptrs()[i + 1] - eval.get_point_ptrs()[i]
+            << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  test<2>(false);
+  test<2>(true);
+}

tests/remote_point_evaluation/mapping_03.with_p4est=true.mpirun=2.output

@@ -0,0 +1,7 @@
+
+DEAL:0::0 1.50000 1.00000 2.00000 1.00000 4
+DEAL:0::0 1.00000 1.00000 1.00000 1.00000 1
+
+DEAL:1::0 1.50000 1.00000 2.00000 1.00000 4
+DEAL:1::0 1.00000 1.00000 1.00000 1.00000 1
+