From 57d5067cba3a44f1d0ca904c64d1839306103338 Mon Sep 17 00:00:00 2001 From: Luca Heltai 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. +
+(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 +# include + # include # include # include @@ -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>> &vertex_to_cell_centers, const typename MeshType::active_cell_iterator &cell_hint = typename MeshType::active_cell_iterator(), - const std::vector &marked_vertices = {}); + const std::vector & marked_vertices = {}, + const RTree, 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>> &vertex_to_cell_centers, const typename MeshType::active_cell_iterator &cell_hint, - const std::vector &marked_vertices) + const std::vector & marked_vertices, + const RTree, unsigned int>> &used_vertices_rtree) { std::pair::active_cell_iterator, Point> @@ -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, unsigned int>; + std::function 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, 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 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::type &, - const std::vector &); + const std::vector &, + const RTree, unsigned int>> &); template std::vector> compute_mesh_predicate_bounding_box( -- 2.39.5