+/**
+ * Given an Rtree object @p tree and a target level @p level, this function returns
+ * the bounding boxes associated to the children on level l+1 and stores them in
+ * a vector. The resulting type is hence a vector of vectors of BoundingBox.
+ * If @p v is such a vector, then @p v has the following property:
+ * @p v[i] is a vector with all of the bounding boxes associated to the children
+ * of the i-th node of the tree.
+ */
+template <typename Rtree>
+std::vector<std::vector<BoundingBox<
+ boost::geometry::dimension<typename Rtree::indexable_type>::value>>>
+extract_children_of_level(const Rtree &tree, const unsigned int level);
+
+
+
// Inline and template functions
#ifndef DOXYGEN
}
+template <typename Value,
+ typename Options,
+ typename Translator,
+ typename Box,
+ typename Allocators>
+struct NodeVisitor : public boost::geometry::index::detail::rtree::visitor<
+ Value,
+ typename Options::parameters_type,
+ Box,
+ Allocators,
+ typename Options::node_tag,
+ true>::type
+{
+ inline NodeVisitor(
+ const Translator &translator,
+ unsigned int target_level,
+ std::vector<
+ std::vector<BoundingBox<boost::geometry::dimension<Box>::value>>> &boxes);
+
+
+ /**
+ * An alias that identifies an InternalNode of the tree.
+ */
+ using InternalNode =
+ typename boost::geometry::index::detail::rtree::internal_node<
+ Value,
+ typename Options::parameters_type,
+ Box,
+ Allocators,
+ typename Options::node_tag>::type;
+
+ /**
+ * An alias that identifies a Leaf of the tree.
+ */
+ using Leaf = typename boost::geometry::index::detail::rtree::leaf<
+ Value,
+ typename Options::parameters_type,
+ Box,
+ Allocators,
+ typename Options::node_tag>::type;
+
+ /**
+ * Implements the visitor interface for InternalNode objects. If the node
+ * belongs to the level next to @p target_level, then fill the bounding box vector for that node.
+ */
+ inline void
+ operator()(InternalNode const &node);
+
+ /**
+ * Implements the visitor interface for Leaf objects.
+ */
+ inline void
+ operator()(Leaf const &);
+
+ /**
+ * Translator interface, required by the boost implementation of the rtree.
+ */
+ Translator const &translator;
+
+ /**
+ * Store the level we are currently visiting.
+ */
+ size_t level;
+
+ /**
+ * Index used to keep track of the number of different visited nodes during
+ * recursion/
+ */
+ size_t node_counter;
+
+ /**
+ * The level where children are living.
+ * Before: "we want to extract from the RTree object."
+ */
+ const size_t target_level;
+
+ /**
+ * A reference to the input vector of vector of BoundingBox objects. This
+ * vector v has the following property: v[i] = vector with all
+ * of the BoundingBox bounded by the i-th node of the Rtree.
+ */
+ std::vector<std::vector<BoundingBox<boost::geometry::dimension<Box>::value>>>
+ &boxes_in_boxes;
+};
+
+
+
+template <typename Value,
+ typename Options,
+ typename Translator,
+ typename Box,
+ typename Allocators>
+NodeVisitor<Value, Options, Translator, Box, Allocators>::NodeVisitor(
+ const Translator & translator,
+ const unsigned int target_level,
+ std::vector<std::vector<BoundingBox<boost::geometry::dimension<Box>::value>>>
+ &bb_in_boxes)
+ : translator(translator)
+ , level(0)
+ , node_counter(0)
+ , target_level(target_level)
+ , boxes_in_boxes(bb_in_boxes)
+{}
+
+
+
+template <typename Value,
+ typename Options,
+ typename Translator,
+ typename Box,
+ typename Allocators>
+void
+NodeVisitor<Value, Options, Translator, Box, Allocators>::operator()(
+ const NodeVisitor::InternalNode &node)
+{
+ using elements_type =
+ typename boost::geometry::index::detail::rtree::elements_type<
+ InternalNode>::type; // pairs of bounding box and pointer to child node
+ const elements_type &elements =
+ boost::geometry::index::detail::rtree::elements(node);
+
+ if (level == target_level)
+ {
+ const unsigned int n_children = elements.size();
+ const auto offset = boxes_in_boxes.size();
+ boxes_in_boxes.resize(offset + n_children);
+ }
+
+ if (level == target_level + 1)
+ {
+ // I have now access to children of level target_level
+ boxes_in_boxes[node_counter].resize(elements.size()); // number of bboxes
+ unsigned int i = 0;
+ for (typename elements_type::const_iterator it = elements.begin();
+ it != elements.end();
+ ++it)
+ {
+ boost::geometry::convert(it->first, boxes_in_boxes[node_counter][i]);
+ ++i;
+ }
+ // Children have been stored, go to the next parent.
+ ++node_counter;
+ return;
+ }
+
+ size_t level_backup = level;
+ ++level;
+
+ for (typename elements_type::const_iterator it = elements.begin();
+ it != elements.end();
+ ++it)
+ {
+ boost::geometry::index::detail::rtree::apply_visitor(*this, *it->second);
+ }
+
+ level = level_backup;
+}
+
+
+
+template <typename Value,
+ typename Options,
+ typename Translator,
+ typename Box,
+ typename Allocators>
+void
+NodeVisitor<Value, Options, Translator, Box, Allocators>::operator()(
+ const NodeVisitor::Leaf &)
+{
+ // No children for leaf nodes.
+ boxes_in_boxes.clear();
+}
+
+template <typename Rtree>
+inline std::vector<std::vector<BoundingBox<
+ boost::geometry::dimension<typename Rtree::indexable_type>::value>>>
+extract_children_of_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<std::vector<BoundingBox<dim>>> boxes_in_boxes;
+
+ if (rtv.depth() == 0)
+ {
+ // The below algorithm does not work for `rtv.depth()==0`, which might
+ // happen if the number entries in the tree is too small.
+ // In this case, simply return a single bounding box.
+ boxes_in_boxes.resize(1);
+ boxes_in_boxes[0].resize(1);
+ boost::geometry::convert(tree.bounds(), boxes_in_boxes[0][0]);
+ }
+ else
+ {
+ const unsigned int target_level =
+ std::min<unsigned int>(level, rtv.depth() - 1);
+
+ NodeVisitor<typename RtreeView::value_type,
+ typename RtreeView::options_type,
+ typename RtreeView::translator_type,
+ typename RtreeView::box_type,
+ typename RtreeView::allocators_type>
+ node_visitor(rtv.translator(), target_level, boxes_in_boxes);
+ rtv.apply_visitor(node_visitor);
+ }
+
+ return boxes_in_boxes;
+}
+
+
#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+// Check that the number of total children of a given level is correct.
+
+#include <deal.II/base/patterns.h>
+
+#include <deal.II/boost_adaptors/bounding_box.h>
+#include <deal.II/boost_adaptors/point.h>
+
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/rtree.h>
+
+#include <algorithm>
+
+#include "../tests.h"
+
+int
+main()
+{
+ initlog(true);
+
+ Triangulation<2> tria;
+ MappingQ1<2> mapping;
+ GridGenerator::hyper_ball(tria);
+ tria.refine_global(6);
+
+ namespace bgi = boost::geometry::index;
+ std::vector<
+ std::pair<BoundingBox<2>, typename Triangulation<2>::active_cell_iterator>>
+ boxes(tria.n_active_cells());
+ unsigned int i = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell);
+
+ const auto tree = pack_rtree<bgi::linear<16>>(boxes);
+ for (const unsigned int level : {0, 1, 2})
+ {
+ const auto bboxes = extract_rtree_level(tree, level + 1);
+ deallog << "LEVEL " + std::to_string(level + 1) + " N boxes: "
+ << bboxes.size() << std::endl;
+
+ const auto boxes_in_boxes = extract_children_of_level(tree, level);
+ unsigned int total_bboxes = 0;
+ for (unsigned int i = 0; i < boxes_in_boxes.size(); ++i)
+ total_bboxes += boxes_in_boxes[i].size();
+
+ if (level == 2)
+ {
+ Assert(boxes_in_boxes.size() == 0,
+ ExcMessage("Leafs have no children."));
+ }
+ else
+ {
+ Assert(total_bboxes == bboxes.size(),
+ ExcMessage(
+ "The number of total children of level " +
+ std::to_string(level) +
+ " should be equal to the number of boxes on level " +
+ std::to_string(level + 1)));
+ }
+
+ deallog << "OK" << std::endl;
+ }
+
+ // for (unsigned int i = 0; i < boxes_in_boxes.size(); ++i)
+ // for (const auto &b : boxes_in_boxes[i])
+ // deallog << "Box: " <<
+ // Patterns::Tools::to_string(b.get_boundary_points())
+ // << std::endl;
+}
--- /dev/null
+
+DEAL::LEVEL 1 N boxes: 80
+DEAL::OK
+DEAL::LEVEL 2 N boxes: 1280
+DEAL::OK
+DEAL::LEVEL 3 N boxes: 1280
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+// Use information from NodeVisitor to print show explicitely the boxes
+// associated to each parent node on the previous level
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/bounding_box_data_out.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <algorithm>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace bg = boost::geometry;
+namespace bgi = boost::geometry::index;
+
+
+
+template <int dim, int spacedim, unsigned int max_elem_per_node>
+void
+test(const unsigned int ref = 6, const unsigned int level = 0)
+{
+ Triangulation<dim, spacedim> tria;
+ MappingQ<dim> mapping(1);
+ GridGenerator::hyper_ball(tria);
+ tria.refine_global(ref);
+
+ std::vector<
+ std::pair<BoundingBox<spacedim>,
+ typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ boxes(tria.n_active_cells());
+ unsigned int i = 0;
+ for (const auto &cell : tria.active_cell_iterators())
+ boxes[i++] = std::make_pair(mapping.get_bounding_box(cell), cell);
+
+ const auto tree = pack_rtree<bgi::linear<max_elem_per_node>>(boxes);
+ Assert(n_levels(tree) == 2,
+ ExcMessage("Two levels are needed for this test."));
+ const auto boxes_in_level = extract_rtree_level(tree, level);
+ deallog << "N Boxes in level " + std::to_string(level) + " "
+ << boxes_in_level.size() << std::endl;
+ const auto boxes_in_boxes = extract_children_of_level(tree, level);
+
+ unsigned int total_number_of_boxes = 0;
+ for (unsigned int i = 0; i < boxes_in_boxes.size(); ++i)
+ {
+ deallog << "boxes_in_boxes[" << i
+ << "] has size: " << boxes_in_boxes[i].size() << std::endl;
+ total_number_of_boxes += boxes_in_boxes[i].size();
+
+ for (const auto &b : boxes_in_boxes[i])
+ deallog << "Box: "
+ << Patterns::Tools::to_string(b.get_boundary_points())
+ << std::endl;
+ }
+
+ // Uncomment to display the partition
+ // {
+ // i = 0;
+ // std::vector<
+ // std::pair<BoundingBox<spacedim>,
+ // typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ // cells;
+
+ // for (unsigned int k = 0; k < boxes_in_boxes.size(); ++k)
+ // {
+ // deallog << "boxes_in_boxes[" << k << "]=" << boxes_in_boxes[k].size()
+ // << std::endl;
+ // for (const auto &box : boxes_in_boxes[k])
+ // {
+ // tree.query(bgi::within(box), std::back_inserter(cells));
+ // deallog << "Number of cells inside " << std::to_string(k)
+ // << "-th bounding box: " << cells.size() << std::endl;
+ // for (const auto &my_pair : cells)
+ // {
+ // deallog << my_pair.second->active_cell_index() << std::endl;
+ // my_pair.second->set_subdomain_id(k);
+ // }
+
+ // cells.clear();
+ // }
+ // }
+
+
+ // for (unsigned int j = 0; j < i; ++j)
+ // deallog << GridTools::count_cells_with_subdomain_association(tria, j)
+ // << " cells have subdomain " + std::to_string(j) << std::endl;
+ // GridOut grid_out_svg;
+ // GridOutFlags::Svg svg_flags;
+ // svg_flags.label_subdomain_id = true;
+ // svg_flags.coloring = GridOutFlags::Svg::subdomain_id;
+ // grid_out_svg.set_flags(svg_flags);
+ // std::ofstream out("grid_ball_agglomerates" +
+ // std::to_string(level) + ".svg");
+ // grid_out_svg.write_svg(tria, out);
+ // }
+
+ // Print hierarchy using boost utilities
+ // bgi::detail::rtree::utilities::print(std::cout, tree);
+}
+
+int
+main()
+{
+ initlog(true);
+
+ static constexpr unsigned int global_refinements = 4; // 6
+ static constexpr unsigned int max_per_node = 16;
+ const unsigned int level = 0; // 1
+ test<2, 2, max_per_node>(global_refinements, level);
+}
--- /dev/null
+
+DEAL::N Boxes in level 0 5
+DEAL::boxes_in_boxes[0] has size: 16
+DEAL::Box: -0.881921, -0.83147 : -0.531665, -0.450016
+DEAL::Box: -0.95694, -0.53175 : -0.536612, -0.276678
+DEAL::Box: -0.607969, -0.92388 : -0.366098, -0.518306
+DEAL::Box: -0.579054, -0.551777 : -0.369543, -0.29632
+DEAL::Box: -0.980785, -0.369233 : -0.591529, -0.125
+DEAL::Box: -1, -0.19509 : -0.609835, -5.55112e-17
+DEAL::Box: -0.635452, -0.355584 : -0.384288, -0.138959
+DEAL::Box: -0.646447, -0.1875 : -0.376705, -1.38778e-17
+DEAL::Box: -0.450016, -0.95694 : -0.224817, -0.522164
+DEAL::Box: -0.4375, -0.554917 : -0.224112, -0.324238
+DEAL::Box: -0.299756, -0.995185 : -0.0837004, -0.538182
+DEAL::Box: -0.28014, -0.570217 : -0.0991117, -0.323358
+DEAL::Box: -0.414752, -0.370558 : -0.256282, -0.146447
+DEAL::Box: -0.418611, -0.185279 : -0.256282, 0
+DEAL::Box: -0.301586, -0.353823 : -0.109835, -0.183058
+DEAL::Box: -0.292893, -0.183058 : -0.146447, 0.0366117
+DEAL::boxes_in_boxes[1] has size: 16
+DEAL::Box: -1, -1.11022e-16 : -0.609835, 0.19509
+DEAL::Box: -0.980785, 0.125 : -0.591529, 0.369233
+DEAL::Box: -0.646447, -5.55112e-17 : -0.376705, 0.1875
+DEAL::Box: -0.635452, 0.138959 : -0.384288, 0.355584
+DEAL::Box: -0.95694, 0.276678 : -0.536612, 0.53175
+DEAL::Box: -0.881921, 0.450016 : -0.531665, 0.83147
+DEAL::Box: -0.579054, 0.29632 : -0.369543, 0.551777
+DEAL::Box: -0.607969, 0.518306 : -0.366098, 0.92388
+DEAL::Box: -0.418611, -1.38778e-17 : -0.256282, 0.215419
+DEAL::Box: -0.414752, 0.172335 : -0.256282, 0.370558
+DEAL::Box: -0.292893, 0 : -0.109835, 0.183058
+DEAL::Box: -0.301586, 0.146447 : -0.109835, 0.353823
+DEAL::Box: -0.4375, 0.324238 : -0.224112, 0.554917
+DEAL::Box: -0.450016, 0.522164 : -0.224817, 0.95694
+DEAL::Box: -0.28014, 0.323358 : -0.0991117, 0.570217
+DEAL::Box: -0.299756, 0.538182 : -0.0837004, 0.995185
+DEAL::boxes_in_boxes[2] has size: 16
+DEAL::Box: -0.167401, -1 : 0.0980171, -0.672244
+DEAL::Box: -0.141585, -0.734835 : 0.0707927, -0.460517
+DEAL::Box: 0.0666464, -0.995185 : 0.316342, -0.635452
+DEAL::Box: 0.0495558, -0.679374 : 0.29632, -0.451364
+DEAL::Box: -0.148667, -0.502423 : 0.0926396, -0.292893
+DEAL::Box: -0.146447, -0.337087 : 0.0796954, -0.21967
+DEAL::Box: 0.0463198, -0.490982 : 0.297335, -0.323358
+DEAL::Box: 0.0732233, -0.353823 : 0.301586, -0.21967
+DEAL::Box: 0.237256, -0.95694 : 0.531665, -0.594264
+DEAL::Box: 0.237056, -0.641192 : 0.508188, -0.422335
+DEAL::Box: 0.438711, -0.881921 : 0.83147, -0.531665
+DEAL::Box: 0.4375, -0.629442 : 0.881921, -0.40839
+DEAL::Box: 0.247779, -0.474112 : 0.490129, -0.277919
+DEAL::Box: 0.256282, -0.34467 : 0.503141, -0.185279
+DEAL::Box: 0.461953, -0.484625 : 0.95694, -0.276678
+DEAL::Box: 0.475682, -0.369233 : 0.980785, -0.167401
+DEAL::boxes_in_boxes[3] has size: 16
+DEAL::Box: -0.146447, -0.21967 : -0.0366117, 0
+DEAL::Box: -0.0732233, -0.21967 : 0.0732233, 0
+DEAL::Box: -0.146447, 0 : -0.0366117, 0.21967
+DEAL::Box: -0.0732233, -0.0366117 : 0.0732233, 0.21967
+DEAL::Box: 0.0366117, -0.21967 : 0.183058, 0
+DEAL::Box: 0.146447, -0.21967 : 0.256282, 0
+DEAL::Box: 0.0366117, -0.0366117 : 0.183058, 0.21967
+DEAL::Box: 0.146447, 0 : 0.256282, 0.21967
+DEAL::Box: -0.148667, 0.21967 : -1.38778e-17, 0.502423
+DEAL::Box: -0.0430837, 0.183058 : 0.0861675, 0.46967
+DEAL::Box: -0.141585, 0.460517 : 0.0707927, 0.734835
+DEAL::Box: -0.167401, 0.672244 : 0.0980171, 1
+DEAL::Box: 0.0430837, 0.183058 : 0.185279, 0.460517
+DEAL::Box: 0.138959, 0.21967 : 0.297335, 0.479541
+DEAL::Box: 0.0495558, 0.442211 : 0.29632, 0.672244
+DEAL::Box: 0.0666464, 0.617055 : 0.299756, 0.995185
+DEAL::boxes_in_boxes[4] has size: 16
+DEAL::Box: 0.256282, -0.239086 : 0.51687, -0.0861675
+DEAL::Box: 0.256282, -0.129251 : 0.513864, 0.0366117
+DEAL::Box: 0.256282, -3.46945e-17 : 0.513864, 0.148667
+DEAL::Box: 0.256282, 0.0991117 : 0.5306, 0.247779
+DEAL::Box: 0.479541, -0.25 : 0.716348, -3.46945e-17
+DEAL::Box: 0.672244, -0.250164 : 1, -6.245e-17
+DEAL::Box: 0.490982, -6.93889e-17 : 0.734835, 0.25
+DEAL::Box: 0.679374, -1.11022e-16 : 1, 0.263072
+DEAL::Box: 0.239086, 0.199238 : 0.396447, 0.456658
+DEAL::Box: 0.34467, 0.198223 : 0.503141, 0.474112
+DEAL::Box: 0.237056, 0.422335 : 0.508188, 0.641192
+DEAL::Box: 0.224817, 0.594264 : 0.507759, 0.95694
+DEAL::Box: 0.461953, 0.1875 : 0.679374, 0.5
+DEAL::Box: 0.63361, 0.176631 : 0.980785, 0.483853
+DEAL::Box: 0.428636, 0.4375 : 0.65533, 0.881921
+DEAL::Box: 0.577665, 0.40839 : 0.881921, 0.77301
projects / dealii.git / commitdiff