]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New function GridTools::guess_point_owner which takes as input an rtree and a vector... 7752/head
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Mon, 25 Feb 2019 15:52:37 +0000 (16:52 +0100)
committerGiovanni Alzetta <giovannialzetta@hotmail.it>
Fri, 10 May 2019 07:42:37 +0000 (09:42 +0200)
doc/news/changes/minor/20190225GiovanniAlzetta [new file with mode: 0644]
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
tests/grid/grid_tools_guess_pt_owner_01.cc [moved from tests/grid/grid_tools_guess_pt_owner_1.cc with 100% similarity]
tests/grid/grid_tools_guess_pt_owner_01.output [moved from tests/grid/grid_tools_guess_pt_owner_1.output with 100% similarity]
tests/grid/grid_tools_guess_pt_owner_02.cc [new file with mode: 0644]
tests/grid/grid_tools_guess_pt_owner_02.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190225GiovanniAlzetta b/doc/news/changes/minor/20190225GiovanniAlzetta
new file mode 100644 (file)
index 0000000..e5c6533
--- /dev/null
@@ -0,0 +1,4 @@
+New: Function GridTools::guess_point_owner(), which uses a covering rtree
+to guess which processes own the points of the given vector of points.
+<br>
+(Giovanni Alzetta, 2019/02/25)
index 00e434cc3af2996e76342b18c2319fefc2985ebc..85e422ea67e93af32d26fd5973ebcc8024c64d89 100644 (file)
@@ -43,6 +43,7 @@
 #  include <boost/archive/binary_iarchive.hpp>
 #  include <boost/archive/binary_oarchive.hpp>
 #  include <boost/geometry/index/detail/serialization.hpp>
+#  include <boost/geometry/index/rtree.hpp>
 #  include <boost/optional.hpp>
 #  include <boost/serialization/array.hpp>
 #  include <boost/serialization/vector.hpp>
@@ -1641,6 +1642,51 @@ namespace GridTools
     const std::vector<Point<spacedim>> &                   points);
 
 
+  /**
+   * Given a covering rtree (see GridTools::Cache::get_covering_rtree()), and an
+   * array of points, find a superset of processes which, individually,
+   * may own the cell containing the point.
+   *
+   * For further details see GridTools::guess_point_owner; here only
+   * different input/output types are reported:
+   *
+   * @param[in] covering_rtree RTRee which enables us to identify which
+   * process(es) in a parallel computation may own the cell that
+   * surrounds a given point.
+   *
+   * @return A tuple containing the following information:
+   *  - A map indexed by processor ranks. For each rank it contains
+   *   a vector of the indices of points it might own.
+   *  - A map from the index <code>unsigned int</code> of the point in @p points
+   *   to the rank of the owner; these are points for which a single possible
+   *   owner was found.
+   *  - A map from the index <code>unsigned int</code> of the point in @p points
+   *   to the ranks of the guessed owners; these are points for which multiple
+   *   possible owners were found.
+   *
+   * @note The actual return type of this function, i.e., the type referenced
+   * above as @p return_type, is
+   * @code
+   * std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
+   *            std::map<unsigned int, unsigned int>,
+   *            std::map<unsigned int, std::vector<unsigned int>>>
+   * @endcode
+   * The type is abbreviated in the online documentation to improve readability
+   * of this page.
+   */
+  template <int spacedim>
+#  ifndef DOXYGEN
+  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
+             std::map<unsigned int, unsigned int>,
+             std::map<unsigned int, std::vector<unsigned int>>>
+#  else
+  return_type
+#  endif
+  guess_point_owner(
+    const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
+    const std::vector<Point<spacedim>> &                         points);
+
+
   /**
    * Return the adjacent cells of all the vertices. If a vertex is also a
    * hanging node, the associated coarse cell is also returned. The vertices
index 58d12d3a94fca95bb4d214afa9f941cc1a79435a..fd42a00dd75456bd2f2858de3381ae6dcbb73b55 100644 (file)
@@ -2103,6 +2103,59 @@ namespace GridTools
     return output_tuple;
   }
 
+  template <int spacedim>
+#ifndef DOXYGEN
+  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
+             std::map<unsigned int, unsigned int>,
+             std::map<unsigned int, std::vector<unsigned int>>>
+#else
+  return_type
+#endif
+  guess_point_owner(
+    const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
+    const std::vector<Point<spacedim>> &                         points)
+  {
+    std::map<unsigned int, std::vector<unsigned int>> point_owners;
+    std::map<unsigned int, unsigned int>              map_owners_found;
+    std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
+    std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> search_result;
+
+    unsigned int n_points = points.size();
+    for (unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
+      {
+        search_result.clear(); // clearing last output
+
+        // Running tree search
+        covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
+                             std::back_inserter(search_result));
+
+        // Keep track of how many processes we guess to own the point
+        std::set<unsigned int> owners_found;
+        // Check in which other processes the point might be
+        for (const auto &rank_bbox : search_result)
+          {
+            // Try to add the owner to the owners found,
+            // and check if it was already present
+            const bool pt_inserted = owners_found.insert(pt_n).second;
+            if (pt_inserted)
+              point_owners[rank_bbox.second].emplace_back(pt_n);
+          }
+        Assert(owners_found.size() > 0,
+               ExcMessage("No owners found for the point " +
+                          std::to_string(pt_n)));
+        if (owners_found.size() == 1)
+          map_owners_found[pt_n] = *owners_found.begin();
+        else
+          // Multiple owners
+          std::copy(owners_found.begin(),
+                    owners_found.end(),
+                    std::back_inserter(map_owners_guessed[pt_n]));
+      }
+
+    return std::make_tuple(std::move(point_owners),
+                           std::move(map_owners_found),
+                           std::move(map_owners_guessed));
+  }
 
 
   template <int dim, int spacedim>
index 23dd04b907b6debaae230217d13a17a10d5c1b03..5d410d900affe5e8746b15de280ac7b41a7ad36f 100644 (file)
@@ -163,6 +163,14 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS)
       const std::vector<std::vector<BoundingBox<deal_II_space_dimension>>> &,
       const std::vector<Point<deal_II_space_dimension>> &);
 
+    template std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
+                        std::map<unsigned int, unsigned int>,
+                        std::map<unsigned int, std::vector<unsigned int>>>
+    GridTools::guess_point_owner(
+      const RTree<std::pair<BoundingBox<deal_II_space_dimension>, unsigned int>>
+        &,
+      const std::vector<Point<deal_II_space_dimension>> &);
+
     template RTree<
       std::pair<BoundingBox<deal_II_space_dimension>, unsigned int>>
     GridTools::build_global_description_tree(
diff --git a/tests/grid/grid_tools_guess_pt_owner_02.cc b/tests/grid/grid_tools_guess_pt_owner_02.cc
new file mode 100644 (file)
index 0000000..4831361
--- /dev/null
@@ -0,0 +1,201 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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 for GridTools::guess_point_owner, tree version:
+// simulating different number of processes with different
+// number of bounding boxes
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/point.h>
+
+#include <deal.II/grid/grid_tools.h>
+
+#include <tuple>
+
+#include "../tests.h"
+
+template <int spacedim>
+void
+test_point_owner(unsigned int n_procs)
+{
+  // Step 1: creating the vector of bounding boxes
+  // We want to simulate n_procs bounding boxes:
+  // each process of rank rk has rk + 1 bounding boxes
+  // the bounding boxes describe a mesh which is a
+  // segment (in dim 1), a rectangle (in dim 2),
+  // a parallelogram (in dim 3)
+  // and using the center point of each bounding box
+
+  std::vector<std::pair<BoundingBox<spacedim>, unsigned int>> global_bboxes;
+
+  unsigned int tot_bbox = 0;
+
+  for (unsigned int rk = 0; rk < n_procs; ++rk)
+    for (unsigned int box = 0; box < rk; ++box)
+      {
+        std::pair<Point<spacedim>, Point<spacedim>> boundaries;
+        boundaries.first[0]  = tot_bbox;
+        boundaries.second[0] = tot_bbox + 1;
+        for (int i = 1; i < spacedim; i++)
+          boundaries.second[i] = 1;
+
+        BoundingBox<spacedim> new_box(boundaries);
+        global_bboxes.emplace_back(std::make_pair(new_box, rk));
+        ++tot_bbox;
+      }
+
+  // Creating the vector of center points
+  std::vector<Point<spacedim>> points(tot_bbox);
+  for (unsigned int i = 0; i < tot_bbox; ++i)
+    {
+      points[i][0] = i + 0.5;
+      for (unsigned int j = 1; j < spacedim; j++)
+        points[i][j] = 0.5;
+    }
+
+  // Step 2: building the rtree
+  RTree<std::pair<BoundingBox<spacedim>, unsigned int>> covering_rtree(
+    global_bboxes.begin(), global_bboxes.end());
+
+  // Step 2: testing the function
+  auto output_tp = GridTools::guess_point_owner(covering_rtree, points);
+
+  bool test_passed = true;
+
+  if (std::get<2>(output_tp).size() != 0)
+    {
+      test_passed = false;
+      deallog << " Test 1 failed: multiple owners" << std::endl;
+    }
+  else
+    {
+      // Check the points are all in the correct rank
+      unsigned int tot_pt = 0;
+      for (unsigned int rk = 0; rk < n_procs; ++rk)
+        {
+          const auto & rk_points = std::get<0>(output_tp)[rk];
+          unsigned int first_el  = rk * (rk + 1) / 2 - 1;
+          unsigned int last_el   = (rk + 1) * (rk + 2) / 2;
+          for (unsigned int box = 0; box < rk; ++box)
+            {
+              if (std::find(rk_points.begin() + first_el,
+                            rk_points.begin() + last_el,
+                            tot_pt) == rk_points.end())
+                {
+                  deallog << "Point " << tot_pt << " not found in rank " << rk
+                          << std::endl;
+                  test_passed = false;
+                }
+              else
+                ++tot_pt;
+            }
+        }
+    }
+
+  // Step 3: creating random points and adding bounding boxes
+
+  // Creating the random points
+  points.clear();
+
+  unsigned int n_points = 2 * tot_bbox;
+  for (std::size_t i = 0; i < n_points; ++i)
+    {
+      Point<spacedim> p;
+      p[0] = tot_bbox * double(Testing::rand()) / RAND_MAX;
+      for (unsigned int d = 1; d < spacedim; ++d)
+        p[d] = double(Testing::rand()) / RAND_MAX;
+      points.push_back(p);
+    }
+
+  // Adding a bounding box to some processes which covers
+  // part of the domain already covered by other processes
+
+  for (unsigned int rk = 1; rk < n_procs; rk += 2)
+    {
+      unsigned int first_el = rk * (rk + 1) / 2 - 1;
+      unsigned int last_el  = (rk + 1) * (rk + 2) / 2;
+      std::pair<Point<spacedim>, Point<spacedim>> boundaries =
+        global_bboxes[first_el].first.get_boundary_points();
+
+      boundaries.first[0] -= 0.5 * double(Testing::rand()) / RAND_MAX;
+      boundaries.second[0] -= 0.5 * double(Testing::rand()) / RAND_MAX;
+
+      BoundingBox<spacedim> new_box(boundaries);
+      global_bboxes.emplace_back(std::make_pair(new_box, rk));
+    }
+
+  // Building a new tree
+  RTree<std::pair<BoundingBox<spacedim>, unsigned int>> covering_rtree_2(
+    global_bboxes.begin(), global_bboxes.end());
+
+  // Running the function again and checking the output
+  output_tp = GridTools::guess_point_owner(covering_rtree_2, points);
+
+  for (unsigned int rk = 0; rk < n_procs; ++rk)
+    {
+      for (const auto &pt : std::get<0>(output_tp)[rk])
+        {
+          bool found = false;
+          for (const auto &bbox : global_bboxes)
+            if (bbox.second == rk && bbox.first.point_inside(points[pt]))
+              {
+                found = true;
+                break;
+              }
+
+          if (!found)
+            {
+              test_passed = false;
+              deallog << "Couldn't find point " << pt << " in process " << rk
+                      << std::endl;
+            }
+        }
+    }
+
+  if (test_passed)
+    deallog << "TEST PASSED" << std::endl;
+  else
+    deallog << "TEST FAILED" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog << "Test for GridTools::guess_point_owner " << std::endl;
+  deallog << std::endl << "Test for dimension 1" << std::endl;
+  for (unsigned int d = 1; d < 4; ++d)
+    {
+      deallog << "Simulating " << d << " processes" << std::endl;
+      test_point_owner<1>(d);
+    }
+
+  deallog << std::endl << "Test for dimension 2" << std::endl;
+  for (unsigned int d = 2; d < 10; ++d)
+    {
+      deallog << "Simulating " << d << " processes" << std::endl;
+      test_point_owner<2>(d);
+    }
+
+
+  deallog << std::endl << "Test for dimension 3" << std::endl;
+  for (unsigned int d = 1; d < 17; d += 2)
+    {
+      deallog << "Simulating " << d << " processes" << std::endl;
+      test_point_owner<3>(d);
+    }
+}
diff --git a/tests/grid/grid_tools_guess_pt_owner_02.output b/tests/grid/grid_tools_guess_pt_owner_02.output
new file mode 100644 (file)
index 0000000..816e811
--- /dev/null
@@ -0,0 +1,47 @@
+
+DEAL::Test for GridTools::guess_point_owner 
+DEAL::
+DEAL::Test for dimension 1
+DEAL::Simulating 1 processes
+DEAL::TEST PASSED
+DEAL::Simulating 2 processes
+DEAL::TEST PASSED
+DEAL::Simulating 3 processes
+DEAL::TEST PASSED
+DEAL::
+DEAL::Test for dimension 2
+DEAL::Simulating 2 processes
+DEAL::TEST PASSED
+DEAL::Simulating 3 processes
+DEAL::TEST PASSED
+DEAL::Simulating 4 processes
+DEAL::TEST PASSED
+DEAL::Simulating 5 processes
+DEAL::TEST PASSED
+DEAL::Simulating 6 processes
+DEAL::TEST PASSED
+DEAL::Simulating 7 processes
+DEAL::TEST PASSED
+DEAL::Simulating 8 processes
+DEAL::TEST PASSED
+DEAL::Simulating 9 processes
+DEAL::TEST PASSED
+DEAL::
+DEAL::Test for dimension 3
+DEAL::Simulating 1 processes
+DEAL::TEST PASSED
+DEAL::Simulating 3 processes
+DEAL::TEST PASSED
+DEAL::Simulating 5 processes
+DEAL::Point 2 not found in rank 2
+DEAL::TEST FAILED
+DEAL::Simulating 7 processes
+DEAL::TEST PASSED
+DEAL::Simulating 9 processes
+DEAL::TEST PASSED
+DEAL::Simulating 11 processes
+DEAL::TEST PASSED
+DEAL::Simulating 13 processes
+DEAL::TEST PASSED
+DEAL::Simulating 15 processes
+DEAL::TEST PASSED

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.