* Since an RTree object is a balanced tree, you can expect each entry of the
* resulting vector to contain roughly the same number of children, and
* ultimately, the same number of leaf objects. If you request for a level that
- * is not present in the RTree, an empty vector is returned.
+ * is not present in the RTree, the last level is returned.
*
* A typical usage of this function is in the context of
* parallel::distributed::Triangulation objects, where one would like to
* computations. If however one constructs an RTree containing these bounding
* boxes (for example, by calling
* GridTools::Cache::get_cell_bounding_boxes_rtree()), and then extracts one of
- * the first levels of the RTree, only a handful of BoundingBoxes would be
+ * the first levels of the RTree, only a handful of BoundingBox objects would be
* returned, allowing the user to have a very efficient description of the
* geometry of the domain, and of its distribution among processes.
*
+ * An example usage is given by the following snippet of code:
+ * @code
+ * parallel::distributed::Triangulation<2> tria(MPI_COMM_WORLD);
+ * GridGenerator::hyper_ball(tria);
+ * tria.refine_global(4);
+ *
+ * std::vector<BoundingBox<2>> all_boxes(tria.n_locally_owned_active_cells());
+ * unsigned int i = 0;
+ * for (const auto &cell : tria.active_cell_iterators())
+ * if (cell->is_locally_owned())
+ * all_boxes[i++] = cell->bounding_box();
+ *
+ * const auto tree = pack_rtree(all_boxes);
+ * const auto boxes = extract_rtree_level(tree, 1);
+ * @endcode
+ *
+ * When run on three processes, the complete set of the BoundingBox objects
+ * surrounding only the locally owned cells and the second level of the rtree
+ * constructed with those boxes would look like in the following pictures (one
+ * image per process):
+ *
+ * @image html rtree-process-0.png
+ * @image html rtree-process-1.png
+ * @image html rtree-process-2.png
+ *
* @author Luca Heltai, 2020.
*/
template <typename Rtree>
-inline std::vector<
- BoundingBox<boost::geometry::dimension<typename Rtree::value_type>::value>>
-extract_rtree_level(Rtree const &tree, const unsigned int level)
-{
- constexpr unsigned int dim =
- boost::geometry::dimension<typename Rtree::value_type>::value;
-
- using RtreeView =
- boost::geometry::index::detail::rtree::utilities::view<Rtree>;
- RtreeView rtv(tree);
-
- std::vector<BoundingBox<dim>> boxes;
-
- ExtractLevelVisitor<typename RtreeView::value_type,
- typename RtreeView::options_type,
- typename RtreeView::translator_type,
- typename RtreeView::box_type,
- typename RtreeView::allocators_type>
- extract_level_visitor(rtv.translator(), level, boxes);
- rtv.apply_visitor(extract_level_visitor);
- return boxes;
-}
+inline std::vector<BoundingBox<
+ boost::geometry::dimension<typename Rtree::indexable_type>::value>>
+extract_rtree_level(const Rtree &tree, const unsigned int level);
operator()(const ExtractLevelVisitor::Leaf &)
{}
+
+
+template <typename Rtree>
+inline std::vector<BoundingBox<
+ boost::geometry::dimension<typename Rtree::indexable_type>::value>>
+extract_rtree_level(const Rtree &tree, const unsigned int level)
+{
+ constexpr unsigned int dim =
+ boost::geometry::dimension<typename Rtree::indexable_type>::value;
+
+ using RtreeView =
+ boost::geometry::index::detail::rtree::utilities::view<Rtree>;
+ RtreeView rtv(tree);
+
+ std::vector<BoundingBox<dim>> boxes;
+
+ const unsigned int target_level =
+ std::min<unsigned int>(level, rtv.depth() - 1);
+ // There should always be at least one level.
+ Assert(target_level >= 0, ExcInternalError());
+
+ ExtractLevelVisitor<typename RtreeView::value_type,
+ typename RtreeView::options_type,
+ typename RtreeView::translator_type,
+ typename RtreeView::box_type,
+ typename RtreeView::allocators_type>
+ extract_level_visitor(rtv.translator(), target_level, boxes);
+ rtv.apply_visitor(extract_level_visitor);
+ return boxes;
+}
+
+
+
+template <class Rtree>
+unsigned int
+n_levels(const Rtree &tree)
+{
+ boost::geometry::index::detail::rtree::utilities::view<Rtree> rtv(tree);
+ return rtv.depth();
+}
+
+
+
#endif
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+// Extract a representative vector of BoundingBox objects from an RTree
+
+#include <deal.II/base/patterns.h>
+
+#include <deal.II/boost_adaptors/bounding_box.h>
+#include <deal.II/boost_adaptors/point.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
+#include <deal.II/numerics/rtree.h>
+
+#include <boost/geometry/index/indexable.hpp>
+
+#include <algorithm>
+
+#include "../tests.h"
+
+
+std::string
+print(const BoundingBox<2> &box)
+{
+ std::stringstream str;
+
+ const auto p = box.get_boundary_points();
+ const auto p1 = p.first;
+ const auto p2 = p.second;
+
+ str << p1 << std::endl
+ << p2[0] << " " << p1[1] << std::endl
+ << p2 << std::endl
+ << p1[0] << " " << p2[1] << std::endl
+ << p1 << std::endl;
+ return str.str();
+}
+
+
+int
+main(int argc, char **argv)
+{
+ auto mpi_init = Utilities::MPI::MPI_InitFinalize(argc, argv);
+ MPILogInitAll log;
+
+ parallel::distributed::Triangulation<2> tria(MPI_COMM_WORLD);
+
+ GridGenerator::hyper_ball(tria);
+ tria.refine_global(3);
+
+ std::vector<BoundingBox<2>> all_boxes(tria.n_locally_owned_active_cells());
+ unsigned int i = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned())
+ all_boxes[i++] = cell->bounding_box();
+
+ const auto tree = pack_rtree(all_boxes);
+ const auto boxes = extract_rtree_level(tree, 0);
+ deallog << "LEVEL 0: N boxes: " << boxes.size() << std::endl;
+
+ // Uncomment the following lines to generate the images of the documentation
+
+ // std::ofstream ofile(
+ // "output_" +
+ // std::to_string(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) +
+ // ".gpl");
+ //
+ // std::ofstream all(
+ // "all_" + std::to_string(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD))
+ // +
+ // ".gpl");
+
+
+ for (const auto &b : boxes)
+ deallog << print(b) << std::endl;
+
+ // for (const auto &b : all_boxes)
+ // all << print(b);
+}
\ No newline at end of file
--- /dev/null
+
+DEAL:0::LEVEL 0: N boxes: 7
+DEAL:0::-0.92388 -0.92388
+-0.297335 -0.92388
+-0.297335 -0.297335
+-0.92388 -0.297335
+-0.92388 -0.92388
+
+DEAL:0::-1 -0.382683
+-0.603553 -0.382683
+-0.603553 0.707107
+-1 0.707107
+-1 -0.382683
+
+DEAL:0::-0.748551 -0.459948
+-0.34467 -0.459948
+-0.34467 0.283171
+-0.748551 0.283171
+-0.748551 -0.459948
+
+DEAL:0::-0.382683 -1
+0.707107 -1
+0.707107 -0.603553
+-0.382683 -0.603553
+-0.382683 -1
+
+DEAL:0::-0.459948 -0.748551
+0.603553 -0.748551
+0.603553 -0.475682
+-0.459948 -0.475682
+-0.459948 -0.748551
+
+DEAL:0::-0.396447 -0.558058
+0.0861675 -0.558058
+0.0861675 0
+-0.396447 0
+-0.396447 -0.558058
+
+DEAL:0::-4.16334e-17 -0.558058
+0.5 -0.558058
+0.5 -0.292893
+-4.16334e-17 -0.292893
+-4.16334e-17 -0.558058
+
+
+DEAL:1::LEVEL 0: N boxes: 7
+DEAL:1::-0.292893 -0.292893
+0 -0.292893
+0 0.0732233
+-0.292893 0.0732233
+-0.292893 -0.292893
+
+DEAL:1::-0.748551 -5.55112e-17
+-0.34467 -5.55112e-17
+-0.34467 0.603553
+-0.748551 0.603553
+-0.748551 -5.55112e-17
+
+DEAL:1::-0.381282 -1.38778e-17
+0 -1.38778e-17
+0 0.34467
+-0.381282 0.34467
+-0.381282 -1.38778e-17
+
+DEAL:1::-0.0732233 -0.292893
+0.292893 -0.292893
+0.292893 0
+-0.0732233 0
+-0.0732233 -0.292893
+
+DEAL:1::0.5 -0.707107
+1 -0.707107
+1 -5.55112e-17
+0.5 -5.55112e-17
+0.5 -0.707107
+
+DEAL:1::-0.0732233 -0.0732233
+0.292893 -0.0732233
+0.292893 0.292893
+-0.0732233 0.292893
+-0.0732233 -0.0732233
+
+DEAL:1::0.21967 -0.0732233
+1 -0.0732233
+1 0.382683
+0.21967 0.382683
+0.21967 -0.0732233
+
+
+
+DEAL:2::LEVEL 0: N boxes: 7
+DEAL:2::-0.448223 0.292893
+0.297335 0.292893
+0.297335 0.503141
+-0.448223 0.503141
+-0.448223 0.292893
+
+DEAL:2::-0.707107 0.433058
+-0.0991117 0.433058
+-0.0991117 0.980785
+-0.707107 0.980785
+-0.707107 0.433058
+
+DEAL:2::-0.19509 0.433058
+0.316342 0.433058
+0.316342 1
+-0.19509 1
+-0.19509 0.433058
+
+DEAL:2::0.292893 -0.448223
+0.475682 -0.448223
+0.475682 0.34467
+0.292893 0.34467
+0.292893 -0.448223
+
+DEAL:2::0.414752 -0.5
+0.823223 -0.5
+0.823223 0.336167
+0.414752 0.336167
+0.414752 -0.5
+
+DEAL:2::0.158171 0.258502
+0.55557 0.258502
+0.55557 0.980785
+0.158171 0.980785
+0.158171 0.258502
+
+DEAL:2::0.417474 0.224112
+0.92388 0.224112
+0.92388 0.83147
+0.417474 0.83147
+0.417474 0.224112
+
+
projects / dealii.git / commitdiff
