]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce NodeVisitor
authorMarco Feder <marco.feder@sissa.it>
Wed, 26 Jul 2023 22:03:15 +0000 (22:03 +0000)
committerMarco Feder <marco.feder@sissa.it>
Fri, 28 Jul 2023 13:54:11 +0000 (13:54 +0000)
include/deal.II/numerics/rtree.h
tests/boost/node_visitor_01.cc [new file with mode: 0644]
tests/boost/node_visitor_01.output [new file with mode: 0644]
tests/boost/node_visitor_02.cc [new file with mode: 0644]
tests/boost/node_visitor_02.output [new file with mode: 0644]

index cbdaefe186d6ed390e104c1d8bb2a6d0b2840839..e25bcebddd8498ede2909d57944f70a7e0312beb 100644 (file)
@@ -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 <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
 
@@ -597,6 +612,220 @@ n_levels(const Rtree &tree)
 }
 
 
+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
 
diff --git a/tests/boost/node_visitor_01.cc b/tests/boost/node_visitor_01.cc
new file mode 100644 (file)
index 0000000..cfcd715
--- /dev/null
@@ -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 <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;
+}
diff --git a/tests/boost/node_visitor_01.output b/tests/boost/node_visitor_01.output
new file mode 100644 (file)
index 0000000..3bbd1a7
--- /dev/null
@@ -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 (file)
index 0000000..006a967
--- /dev/null
@@ -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 <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);
+}
diff --git a/tests/boost/node_visitor_02.output b/tests/boost/node_visitor_02.output
new file mode 100644 (file)
index 0000000..41af87d
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.