GridTools::find_active_cell_around_point(): improve detection if point is outside...

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 5b3e934f4a1df722214393664855c97fd743d2a0..2fa4e1413d66c4ac126f895ae250db37076cf6bf 100644 (file)
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -1501,7 +1501,11 @@ namespace GridTools
     const std::vector<bool> &                              marked_vertices = {},
     const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
       RTree<std::pair<Point<spacedim>, unsigned int>>{},
-    const double tolerance = 1.e-10);
+    const double tolerance = 1.e-10,
+    const RTree<
+      std::pair<BoundingBox<spacedim>,
+                typename Triangulation<dim, spacedim>::active_cell_iterator>>
+      *relevant_cell_bounding_boxes_rtree = nullptr);
 
   /**
    * As compared to the functions above, this function identifies all active

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 1370346917e9869b2d66f6f962e1a6217dfcf3af..fcbd2772d1bb05240f3275d7ee60a439efe1ba73 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -2767,7 +2767,11 @@ namespace GridTools
     const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
     const std::vector<bool> &                              marked_vertices,
     const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree,
-    const double                                           tolerance)
+    const double                                           tolerance,
+    const RTree<
+      std::pair<BoundingBox<spacedim>,
+                typename Triangulation<dim, spacedim>::active_cell_iterator>>
+      *relevant_cell_bounding_boxes_rtree)
   {
     std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
               Point<dim>>
@@ -2778,6 +2782,21 @@ namespace GridTools
               Point<dim>>
       cell_and_position_approx;
 
+    if (relevant_cell_bounding_boxes_rtree != nullptr &&
+        !relevant_cell_bounding_boxes_rtree->empty())
+      {
+        using point_t = boost::geometry::model::
+          point<double, spacedim, boost::geometry::cs::cartesian>;
+        const auto pt =
+          spacedim == 1 ?
+            point_t(p[0]) :
+            (spacedim == 2 ? point_t(p[0], p[1]) : point_t(p[0], p[1], p[2]));
+        if (relevant_cell_bounding_boxes_rtree->qbegin(
+              boost::geometry::index::intersects(pt)) ==
+            relevant_cell_bounding_boxes_rtree->qend())
+          return cell_and_position;
+      }
+
     bool found_cell  = false;
     bool approx_cell = false;
 
@@ -5744,7 +5763,16 @@ namespace GridTools
         locally_owned_active_cells_around_point;
 
       const auto first_cell = GridTools::find_active_cell_around_point(
-        cache, point, cell_hint, marked_vertices, tolerance);
+        cache.get_mapping(),
+        cache.get_triangulation(),
+        point,
+        cache.get_vertex_to_cell_map(),
+        cache.get_vertex_to_cell_centers_directions(),
+        cell_hint,
+        marked_vertices,
+        cache.get_used_vertices_rtree(),
+        tolerance,
+        &cache.get_locally_owned_cell_bounding_boxes_rtree());
 
       cell_hint = first_cell.first;
       if (cell_hint.state() == IteratorState::valid)

diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index e8fefe6909e0ed783b99489ddd43af73be500090..57d21677f4e113ed19a38f5d2c78a3e74b0b58da 100644 (file)
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -38,7 +38,11 @@ for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS;
                                                    X>::type &,
         const std::vector<bool> &,
         const RTree<std::pair<Point<deal_II_space_dimension>, unsigned int>> &,
-        const double);
+        const double,
+        const RTree<std::pair<
+          BoundingBox<deal_II_space_dimension>,
+          typename Triangulation<deal_II_dimension, deal_II_space_dimension>::
+            active_cell_iterator>> *);
 
       template std::vector<BoundingBox<deal_II_space_dimension>>
       compute_mesh_predicate_bounding_box<X>(