From: Luca Heltai Date: Tue, 26 May 2020 11:10:09 +0000 (+0200) Subject: get_locally_owned_cell_bounding_boxes_rtree and covering. X-Git-Tag: v9.3.0-rc1~1419^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=23a2cbb68bab34ecb7caf9287be3cb3588d5e930;p=dealii.git get_locally_owned_cell_bounding_boxes_rtree and covering. --- diff --git a/doc/news/changes/minor/20200526LucaHeltai b/doc/news/changes/minor/20200526LucaHeltai new file mode 100644 index 0000000000..52eb8d13f7 --- /dev/null +++ b/doc/news/changes/minor/20200526LucaHeltai @@ -0,0 +1,6 @@ +New: GridTools::Cache::get_locally_owned_cell_bounding_boxes_rtree() extracts +a tree of bounding boxes covering the locally owned cells of a triangulation. +This can be used in conjunction with GridTools::Cache::get_covering_rtree() to +make approximate geometrical queries on who owns what spacial region. +
+(Luca Heltai, 2020/05/26) diff --git a/include/deal.II/grid/grid_tools_cache.h b/include/deal.II/grid/grid_tools_cache.h index 03f639f47d..26140d8097 100644 --- a/include/deal.II/grid/grid_tools_cache.h +++ b/include/deal.II/grid/grid_tools_cache.h @@ -127,14 +127,31 @@ namespace GridTools get_used_vertices_rtree() const; /** - * Return the cached RTree object of the cell bounging boxes, constructed - * using the active cell iterators of the stored triangulation. + * Return the cached RTree object of the cell bounding boxes, constructed + * using the active cell iterators of the stored triangulation. For + * parallel::distributed::Triangulation objects, this function will return + * also the bounding boxes of ghost and artificial cells. */ const RTree< std::pair, typename Triangulation::active_cell_iterator>> & get_cell_bounding_boxes_rtree() const; + /** + * Return the cached RTree object of bounding boxes containing locally owned + * active cells, constructed using the active cell iterators of the stored + * triangulation. + * + * In contrast to the previous function, this function builds the RTree + * using only the locally owned cells, i.e., not including ghost or + * artificial cells. The two functions return the same result in serial + * computations. + */ + const RTree< + std::pair, + typename Triangulation::active_cell_iterator>> & + get_locally_owned_cell_bounding_boxes_rtree() const; + /** * Return a reference to the stored triangulation. */ @@ -155,7 +172,9 @@ namespace GridTools * object -- an Rtree -- are pairs of bounding boxes denoting * areas that cover all or parts of the local portion of a * parallel triangulation, and an unsigned int representing - * the process or subdomain that owns these cells. + * the process or subdomain that owns the cells falling within the given + * bounding box. + * * Given a point on a parallel::TriangulationBase, this tree * allows to identify one, or few candidate processes, for * which the point lies on a locally owned cell. @@ -167,12 +186,15 @@ namespace GridTools * Therefore this function must be called by all processes * at the same time. * - * While each box may only cover part of a process's locally - * owned part of the triangulation, the boxes associated with - * each process jointly cover the entire local portion. + * The local bounding boxes are constructed by extracting the + * specified @p level from the rtree object returned by the + * get_locally_owned_cell_bounding_boxes_rtree(). Notice that @p level + * in this context does not refer to the level of the triangulation, but + * refer to the level of the RTree object (see, e.g., + * https://en.wikipedia.org/wiki/R-tree). */ const RTree, unsigned int>> & - get_covering_rtree() const; + get_covering_rtree(const unsigned int level = 0) const; private: /** @@ -207,9 +229,15 @@ namespace GridTools vertex_to_cell_centers; /** - * An rtree object covering the whole mesh. + * A collection of rtree objects covering the whole mesh. + * + * Each entry of the map is constructed from the function + * extract_rtree_level() applied to the rtree returned by the function + * get_locally_owned_cell_bounding_boxes_rtree(), with the specified level + * in input. */ - mutable RTree, unsigned int>> + mutable std::map, unsigned int>>> covering_rtree; /** @@ -232,6 +260,15 @@ namespace GridTools typename Triangulation::active_cell_iterator>> cell_bounding_boxes_rtree; + /** + * Store an RTree object, containing the bounding boxes of the locally owned + * cells of the triangulation. + */ + mutable RTree< + std::pair, + typename Triangulation::active_cell_iterator>> + locally_owned_cell_bounding_boxes_rtree; + /** * Storage for the status of the triangulation signal. */ diff --git a/include/deal.II/grid/grid_tools_cache_update_flags.h b/include/deal.II/grid/grid_tools_cache_update_flags.h index e46315622e..2a89e3236c 100644 --- a/include/deal.II/grid/grid_tools_cache_update_flags.h +++ b/include/deal.II/grid/grid_tools_cache_update_flags.h @@ -74,6 +74,11 @@ namespace GridTools */ update_covering_rtree = 0x040, + /** + * Update an RTree of locally owned cell bounding boxes. + */ + update_locally_owned_cell_bounding_boxes_rtree = 0x080, + /** * Update all objects. */ diff --git a/source/grid/grid_tools_cache.cc b/source/grid/grid_tools_cache.cc index 22e8796b5d..9cd4a8704b 100644 --- a/source/grid/grid_tools_cache.cc +++ b/source/grid/grid_tools_cache.cc @@ -144,43 +144,59 @@ namespace GridTools template - const RTree, unsigned int>> & - Cache::get_covering_rtree() const + const RTree< + std::pair, + typename Triangulation::active_cell_iterator>> & + Cache::get_locally_owned_cell_bounding_boxes_rtree() const { - if (update_flags & update_covering_rtree) + if (update_flags & update_locally_owned_cell_bounding_boxes_rtree) { - // TO DO: the locally owned portion of the domain - // might consists of more separate pieces: a single - // bounding box might not always be the best choice. + std::vector, + typename Triangulation::active_cell_iterator>> + boxes; + boxes.reserve(tria->n_active_cells()); + for (const auto &cell : tria->active_cell_iterators()) + if (cell->is_locally_owned()) + boxes.emplace_back( + std::make_pair(mapping->get_bounding_box(cell), cell)); - // Note: we create the local variable ptr here, because gcc 4.8.5 - // fails to compile if we pass the variable directly. - IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate; - std::function::active_cell_iterator &)> - ptr(locally_owned_cell_predicate); + locally_owned_cell_bounding_boxes_rtree = pack_rtree(boxes); + update_flags = + update_flags & ~update_locally_owned_cell_bounding_boxes_rtree; + } + return locally_owned_cell_bounding_boxes_rtree; + } - const BoundingBox bbox = - GridTools::compute_bounding_box(get_triangulation(), ptr); - std::vector> bbox_v(1, bbox); + template + const RTree, unsigned int>> & + Cache::get_covering_rtree(const unsigned int level) const + { + if (update_flags & update_covering_rtree || + covering_rtree.find(level) == covering_rtree.end()) + { + const auto boxes = + extract_rtree_level(get_locally_owned_cell_bounding_boxes_rtree(), + level); + if (const auto tria_mpi = dynamic_cast *>( &(*tria))) { - covering_rtree = GridTools::build_global_description_tree( - bbox_v, tria_mpi->get_communicator()); + covering_rtree[level] = GridTools::build_global_description_tree( + boxes, tria_mpi->get_communicator()); } else { - covering_rtree = - GridTools::build_global_description_tree(bbox_v, MPI_COMM_SELF); + covering_rtree[level] = + GridTools::build_global_description_tree(boxes, MPI_COMM_SELF); } update_flags = update_flags & ~update_covering_rtree; } - return covering_rtree; + return covering_rtree[level]; } #include "grid_tools_cache.inst" diff --git a/tests/grid/grid_tools_cache_06.cc b/tests/grid/grid_tools_cache_06.cc new file mode 100644 index 0000000000..73772689d9 --- /dev/null +++ b/tests/grid/grid_tools_cache_06.cc @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test that Cache::get_locally_owned_cell_bounding_boxes_rtree() works + +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +namespace bgi = boost::geometry::index; + +template +void +test(unsigned int ref) +{ + const MPI_Comm &mpi_communicator = MPI_COMM_WORLD; + + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_ball(tria); + tria.refine_global(ref); + + GridTools::Cache cache(tria); + + const auto &tree = cache.get_locally_owned_cell_bounding_boxes_rtree(); + const auto box = random_box(); + + deallog << "Query box: " << box.get_boundary_points().first << "; " + << box.get_boundary_points().second << std::endl; + + for (const auto p : tree | bgi::adaptors::queried(bgi::intersects(box))) + deallog << "Cell " << p.second << " intersects box with " + << p.first.get_boundary_points().first << "; " + << p.first.get_boundary_points().second << std::endl; +} + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + test<2>(3); + + return 0; +} diff --git a/tests/grid/grid_tools_cache_06.mpirun=1.with_p4est=true.output b/tests/grid/grid_tools_cache_06.mpirun=1.with_p4est=true.output new file mode 100644 index 0000000000..6a9cf8d81e --- /dev/null +++ b/tests/grid/grid_tools_cache_06.mpirun=1.with_p4est=true.output @@ -0,0 +1,5 @@ + +DEAL:0::Query box: 0.783099 0.394383; 0.840188 0.798440 +DEAL:0::Cell 3.212 intersects box with 0.753761 0.349513; 0.923880 0.555570 +DEAL:0::Cell 3.213 intersects box with 0.655330 0.507759; 0.831470 0.707107 +DEAL:0::Cell 3.214 intersects box with 0.676052 0.316342; 0.836215 0.507759 diff --git a/tests/grid/grid_tools_cache_06.mpirun=2.with_p4est=true.output b/tests/grid/grid_tools_cache_06.mpirun=2.with_p4est=true.output new file mode 100644 index 0000000000..827308ab05 --- /dev/null +++ b/tests/grid/grid_tools_cache_06.mpirun=2.with_p4est=true.output @@ -0,0 +1,8 @@ + +DEAL:0::Query box: 0.783099 0.394383; 0.840188 0.798440 + +DEAL:1::Query box: 0.783099 0.394383; 0.840188 0.798440 +DEAL:1::Cell 3.114 intersects box with 0.676052 0.316342; 0.836215 0.507759 +DEAL:1::Cell 3.112 intersects box with 0.753761 0.349513; 0.923880 0.555570 +DEAL:1::Cell 3.113 intersects box with 0.655330 0.507759; 0.831470 0.707107 + diff --git a/tests/grid/grid_tools_cache_07.cc b/tests/grid/grid_tools_cache_07.cc new file mode 100644 index 0000000000..b8544974c8 --- /dev/null +++ b/tests/grid/grid_tools_cache_07.cc @@ -0,0 +1,71 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// test that Cache::get_covering_rtree() works + +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +namespace bgi = boost::geometry::index; + +template +void +test(unsigned int ref) +{ + const MPI_Comm &mpi_communicator = MPI_COMM_WORLD; + + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_ball(tria); + tria.refine_global(ref); + + GridTools::Cache cache(tria); + + const auto &tree = cache.get_covering_rtree(); + const auto point = random_point(); + + deallog << "Query point: " << point << std::endl; + + for (const auto p : tree | bgi::adaptors::queried(bgi::intersects(point))) + deallog << "Processor " << p.second << " may own the point " << point + << ": it is within " << p.first.get_boundary_points().first << "; " + << p.first.get_boundary_points().second << std::endl; +} + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + test<2>(3); + + return 0; +} diff --git a/tests/grid/grid_tools_cache_07.mpirun=1.with_p4est=true.output b/tests/grid/grid_tools_cache_07.mpirun=1.with_p4est=true.output new file mode 100644 index 0000000000..cc80e90553 --- /dev/null +++ b/tests/grid/grid_tools_cache_07.mpirun=1.with_p4est=true.output @@ -0,0 +1,3 @@ + +DEAL:0::Query point: 0.840188 0.394383 +DEAL:0::Processor 0 may own the point 0.840188 0.394383: it is within 0.316342 -0.923880; 1.00000 0.923880 diff --git a/tests/grid/grid_tools_cache_07.mpirun=2.with_p4est=true.output b/tests/grid/grid_tools_cache_07.mpirun=2.with_p4est=true.output new file mode 100644 index 0000000000..86e8503339 --- /dev/null +++ b/tests/grid/grid_tools_cache_07.mpirun=2.with_p4est=true.output @@ -0,0 +1,9 @@ + +DEAL:0::Query point: 0.840188 0.394383 +DEAL:0::Processor 1 may own the point 0.840188 0.394383: it is within 0.316342 0.283171; 0.923880 0.923880 +DEAL:0::Processor 1 may own the point 0.840188 0.394383: it is within 0.372129 -1.11022e-16; 1.00000 0.417474 + +DEAL:1::Query point: 0.840188 0.394383 +DEAL:1::Processor 1 may own the point 0.840188 0.394383: it is within 0.316342 0.283171; 0.923880 0.923880 +DEAL:1::Processor 1 may own the point 0.840188 0.394383: it is within 0.372129 -1.11022e-16; 1.00000 0.417474 + diff --git a/tests/tests.h b/tests/tests.h index 2fc3e79ab3..fef16c13d4 100644 --- a/tests/tests.h +++ b/tests/tests.h @@ -20,6 +20,7 @@ #include +#include #include #include #include @@ -268,6 +269,20 @@ random_point(const double &min = 0.0, const double &max = 1.0) +// Construct a uniformly distributed random box, with each coordinate +// between min and max +template +inline BoundingBox +random_box(const double &min = 0.0, const double &max = 1.0) +{ + Assert(max >= min, ExcMessage("Make sure max>=min")); + std::vector> p = {random_point(min, max), + random_point(min, max)}; + return BoundingBox(p); +} + + + // given the name of a file, copy it to deallog // and then delete it void