From 57d5067cba3a44f1d0ca904c64d1839306103338 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Wed, 5 Dec 2018 13:52:50 +0100
Subject: [PATCH] Find active cell around point now takes also an rtree of
 vertices.

---
 doc/news/changes/minor/20181205LucaHeltai |  9 ++++++
 include/deal.II/grid/grid_tools.h         | 12 ++++++--
 source/grid/grid_tools.cc                 | 36 ++++++++++++++++++++---
 source/grid/grid_tools.inst.in            |  3 +-
 4 files changed, 52 insertions(+), 8 deletions(-)
 create mode 100644 doc/news/changes/minor/20181205LucaHeltai

diff --git a/doc/news/changes/minor/20181205LucaHeltai b/doc/news/changes/minor/20181205LucaHeltai
new file mode 100644
index 0000000000..126f5cb3a2
--- /dev/null
+++ b/doc/news/changes/minor/20181205LucaHeltai
@@ -0,0 +1,9 @@
+New: GridTools::find_active_cell_around_point now allows you to specify an (optional) rtree, constructed
+from the used vertices of the triangulation. Once you have built a tree, querying for a nearest
+vertex is an O(log(N)) operation, where N is the number of used vertices. You can ask a GridTools::Cache
+object to return a tree that is compatible with the new function signature. 
+The previous version of this function had a cost that was O(N^2) when the point was not in the cell_hint
+object. 
+<br>
+(Luca Heltai, 2018/12/05)
+
diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 9cdb1543d1..09f39d5293 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -36,6 +36,8 @@
 
 #  include <deal.II/lac/sparsity_tools.h>
 
+#  include <deal.II/numerics/rtree.h>
+
 #  include <boost/archive/binary_iarchive.hpp>
 #  include <boost/archive/binary_oarchive.hpp>
 #  include <boost/optional.hpp>
@@ -1125,8 +1127,10 @@ namespace GridTools
    * A version of the previous function that exploits an already existing
    * map between vertices and cells, constructed using the function
    * GridTools::vertex_to_cell_map, a map of vertex_to_cell_centers, obtained
-   * through GridTools::vertex_to_cell_centers_directions, and a guess
-   * `cell_hint`.
+   * through GridTools::vertex_to_cell_centers_directions, a guess
+   * `cell_hint`, and optionally an RTree constructed from the used
+   * vertices of the Triangulation. All of these structures can be queried
+   * from a GridTools::Cache object.
    *
    * @author Luca Heltai, Rene Gassmoeller, 2017
    */
@@ -1148,7 +1152,9 @@ namespace GridTools
     const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
     const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
       typename MeshType<dim, spacedim>::active_cell_iterator(),
-    const std::vector<bool> &marked_vertices = {});
+    const std::vector<bool> &                              marked_vertices = {},
+    const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
+      {});
 
   /**
    * A version of the previous function where we use that mapping on a given
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index b9b1419c95..ac17460dac 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -1683,7 +1683,8 @@ namespace GridTools
       &                                                  vertex_to_cells,
     const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
     const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
-    const std::vector<bool> &marked_vertices)
+    const std::vector<bool> &                              marked_vertices,
+    const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree)
   {
     std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
               Point<dim>>
@@ -1714,8 +1715,33 @@ namespace GridTools
           }
         else
           {
-            closest_vertex_index =
-              GridTools::find_closest_vertex(mesh, p, marked_vertices);
+            if (!used_vertices_rtree.empty())
+              {
+                // If we have an rtree at our disposal, use it.
+                using ValueType = std::pair<Point<spacedim>, unsigned int>;
+                std::function<bool(const ValueType &)> marked;
+                if (marked_vertices.size() == mesh.n_vertices())
+                  marked = [&marked_vertices](const ValueType &value) -> bool {
+                    return marked_vertices[value.second];
+                  };
+                else
+                  marked = [](const ValueType &) -> bool { return true; };
+
+                std::vector<std::pair<Point<spacedim>, unsigned int>> res;
+                used_vertices_rtree.query(
+                  boost::geometry::index::nearest(p, 1) &&
+                    boost::geometry::index::satisfies(marked),
+                  std::back_inserter(res));
+
+                // We should have one and only one result
+                AssertDimension(res.size(), 1);
+                closest_vertex_index = res[0].second;
+              }
+            else
+              {
+                closest_vertex_index =
+                  GridTools::find_closest_vertex(mesh, p, marked_vertices);
+              }
             vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
           }
 
@@ -4963,6 +4989,7 @@ namespace GridTools
     const auto &vertex_to_cells = cache.get_vertex_to_cell_map();
     const auto &vertex_to_cell_centers =
       cache.get_vertex_to_cell_centers_directions();
+    const auto &used_vertices_rtree = cache.get_used_vertices_rtree();
 
     return find_active_cell_around_point(mapping,
                                          mesh,
@@ -4970,7 +4997,8 @@ namespace GridTools
                                          vertex_to_cells,
                                          vertex_to_cell_centers,
                                          cell_hint,
-                                         marked_vertices);
+                                         marked_vertices,
+                                         used_vertices_rtree);
   }
 
   template <int spacedim>
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 03bb9c6d38..d8c1a41e5f 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -36,7 +36,8 @@ for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS;
         const dealii::internal::ActiveCellIterator<deal_II_dimension,
                                                    deal_II_space_dimension,
                                                    X>::type &,
-        const std::vector<bool> &);
+        const std::vector<bool> &,
+        const RTree<std::pair<Point<deal_II_space_dimension>, unsigned int>> &);
 
       template std::vector<BoundingBox<deal_II_space_dimension>>
       compute_mesh_predicate_bounding_box<X>(
-- 
2.39.5