From 95d1987cd96453a43c60700f5b75f77e7e29c9e9 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Wed, 26 Jul 2023 22:03:15 +0000 Subject: [PATCH] Introduce NodeVisitor --- include/deal.II/numerics/rtree.h | 229 +++++++++++++++++++++++++++++ tests/boost/node_visitor_01.cc | 87 +++++++++++ tests/boost/node_visitor_01.output | 7 + tests/boost/node_visitor_02.cc | 131 +++++++++++++++++ tests/boost/node_visitor_02.output | 87 +++++++++++ 5 files changed, 541 insertions(+) create mode 100644 tests/boost/node_visitor_01.cc create mode 100644 tests/boost/node_visitor_01.output create mode 100644 tests/boost/node_visitor_02.cc create mode 100644 tests/boost/node_visitor_02.output diff --git a/include/deal.II/numerics/rtree.h b/include/deal.II/numerics/rtree.h index cbdaefe186..e25bcebddd 100644 --- a/include/deal.II/numerics/rtree.h +++ b/include/deal.II/numerics/rtree.h @@ -427,6 +427,21 @@ extract_rtree_level(const Rtree &tree, const unsigned int level); +/** + * 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 +std::vector::value>>> +extract_children_of_level(const Rtree &tree, const unsigned int level); + + + // Inline and template functions #ifndef DOXYGEN @@ -597,6 +612,220 @@ n_levels(const Rtree &tree) } +template +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::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::value>>> + &boxes_in_boxes; +}; + + + +template +NodeVisitor::NodeVisitor( + const Translator & translator, + const unsigned int target_level, + std::vector::value>>> + &bb_in_boxes) + : translator(translator) + , level(0) + , node_counter(0) + , target_level(target_level) + , boxes_in_boxes(bb_in_boxes) +{} + + + +template +void +NodeVisitor::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 +void +NodeVisitor::operator()( + const NodeVisitor::Leaf &) +{ + // No children for leaf nodes. + boxes_in_boxes.clear(); +} + +template +inline std::vector::value>>> +extract_children_of_level(const Rtree &tree, const unsigned int level) +{ + constexpr unsigned int dim = + boost::geometry::dimension::value; + + using RtreeView = + boost::geometry::index::detail::rtree::utilities::view; + RtreeView rtv(tree); + + std::vector>> 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(level, rtv.depth() - 1); + + NodeVisitor + node_visitor(rtv.translator(), target_level, boxes_in_boxes); + rtv.apply_visitor(node_visitor); + } + + return boxes_in_boxes; +} + + #endif diff --git a/tests/boost/node_visitor_01.cc b/tests/boost/node_visitor_01.cc new file mode 100644 index 0000000000..cfcd7150eb --- /dev/null +++ b/tests/boost/node_visitor_01.cc @@ -0,0 +1,87 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include +#include + +#include + +#include +#include + +#include + +#include + +#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, 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>(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; +} diff --git a/tests/boost/node_visitor_01.output b/tests/boost/node_visitor_01.output new file mode 100644 index 0000000000..3bbd1a7720 --- /dev/null +++ b/tests/boost/node_visitor_01.output @@ -0,0 +1,7 @@ + +DEAL::LEVEL 1 N boxes: 80 +DEAL::OK +DEAL::LEVEL 2 N boxes: 1280 +DEAL::OK +DEAL::LEVEL 3 N boxes: 1280 +DEAL::OK diff --git a/tests/boost/node_visitor_02.cc b/tests/boost/node_visitor_02.cc new file mode 100644 index 0000000000..006a967081 --- /dev/null +++ b/tests/boost/node_visitor_02.cc @@ -0,0 +1,131 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; +namespace bg = boost::geometry; +namespace bgi = boost::geometry::index; + + + +template +void +test(const unsigned int ref = 6, const unsigned int level = 0) +{ + Triangulation tria; + MappingQ mapping(1); + GridGenerator::hyper_ball(tria); + tria.refine_global(ref); + + std::vector< + std::pair, + typename Triangulation::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>(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, + // typename Triangulation::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); +} diff --git a/tests/boost/node_visitor_02.output b/tests/boost/node_visitor_02.output new file mode 100644 index 0000000000..41af87d067 --- /dev/null +++ b/tests/boost/node_visitor_02.output @@ -0,0 +1,87 @@ + +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 -- 2.39.5