--- /dev/null
+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)
# 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>
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
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>
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(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+ }
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff