From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 29 May 2021 20:58:31 +0000 (+0200)
Subject: GridTools::find_active_cell_around_point(): improve detection if point is outside... 
X-Git-Tag: v9.4.0-rc1~1308^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a8c270f99d0e0aedfeac07ecfb5bb6d75de40fcd;p=dealii.git

GridTools::find_active_cell_around_point(): improve detection if point is outside of domain
---

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index 5b3e934f4a..2fa4e1413d 100644
--- 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 1370346917..fcbd2772d1 100644
--- 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 e8fefe6909..57d21677f4 100644
--- 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>(