From a8c270f99d0e0aedfeac07ecfb5bb6d75de40fcd Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 29 May 2021 22:58:31 +0200 Subject: [PATCH] GridTools::find_active_cell_around_point(): improve detection if point is outside of domain --- include/deal.II/grid/grid_tools.h | 6 +++++- source/grid/grid_tools.cc | 32 +++++++++++++++++++++++++++++-- source/grid/grid_tools.inst.in | 6 +++++- 3 files changed, 40 insertions(+), 4 deletions(-) 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 & marked_vertices = {}, const RTree, unsigned int>> &used_vertices_rtree = RTree, unsigned int>>{}, - const double tolerance = 1.e-10); + const double tolerance = 1.e-10, + const RTree< + std::pair, + typename Triangulation::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::active_cell_iterator &cell_hint, const std::vector & marked_vertices, const RTree, unsigned int>> &used_vertices_rtree, - const double tolerance) + const double tolerance, + const RTree< + std::pair, + typename Triangulation::active_cell_iterator>> + *relevant_cell_bounding_boxes_rtree) { std::pair::active_cell_iterator, Point> @@ -2778,6 +2782,21 @@ namespace GridTools Point> cell_and_position_approx; + if (relevant_cell_bounding_boxes_rtree != nullptr && + !relevant_cell_bounding_boxes_rtree->empty()) + { + using point_t = boost::geometry::model:: + point; + 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 &, const RTree, unsigned int>> &, - const double); + const double, + const RTree, + typename Triangulation:: + active_cell_iterator>> *); template std::vector> compute_mesh_predicate_bounding_box( -- 2.39.5