From d7004e60dce2b93a9881ec031324f4212be42893 Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Mon, 25 Feb 2019 16:52:37 +0100 Subject: [PATCH] New function GridTools::guess_point_owner which takes as input an rtree and a vector of points --- .../changes/minor/20190225GiovanniAlzetta | 4 + include/deal.II/grid/grid_tools.h | 46 ++++ source/grid/grid_tools.cc | 53 +++++ source/grid/grid_tools.inst.in | 8 + ...r_1.cc => grid_tools_guess_pt_owner_01.cc} | 0 ...ut => grid_tools_guess_pt_owner_01.output} | 0 tests/grid/grid_tools_guess_pt_owner_02.cc | 201 ++++++++++++++++++ .../grid/grid_tools_guess_pt_owner_02.output | 47 ++++ 8 files changed, 359 insertions(+) create mode 100644 doc/news/changes/minor/20190225GiovanniAlzetta rename tests/grid/{grid_tools_guess_pt_owner_1.cc => grid_tools_guess_pt_owner_01.cc} (100%) rename tests/grid/{grid_tools_guess_pt_owner_1.output => grid_tools_guess_pt_owner_01.output} (100%) create mode 100644 tests/grid/grid_tools_guess_pt_owner_02.cc create mode 100644 tests/grid/grid_tools_guess_pt_owner_02.output diff --git a/doc/news/changes/minor/20190225GiovanniAlzetta b/doc/news/changes/minor/20190225GiovanniAlzetta new file mode 100644 index 0000000000..e5c653343e --- /dev/null +++ b/doc/news/changes/minor/20190225GiovanniAlzetta @@ -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. +
+(Giovanni Alzetta, 2019/02/25) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 00e434cc3a..85e422ea67 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -43,6 +43,7 @@ # include # include # include +# include # include # include # include @@ -1641,6 +1642,51 @@ namespace GridTools const std::vector> & 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 unsigned int 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 unsigned int 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, + * std::map>> + * @endcode + * The type is abbreviated in the online documentation to improve readability + * of this page. + */ + template +# ifndef DOXYGEN + std::tuple>, + std::map, + std::map>> +# else + return_type +# endif + guess_point_owner( + const RTree, unsigned int>> &covering_rtree, + const std::vector> & 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 diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 58d12d3a94..fd42a00dd7 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2103,6 +2103,59 @@ namespace GridTools return output_tuple; } + template +#ifndef DOXYGEN + std::tuple>, + std::map, + std::map>> +#else + return_type +#endif + guess_point_owner( + const RTree, unsigned int>> &covering_rtree, + const std::vector> & points) + { + std::map> point_owners; + std::map map_owners_found; + std::map> map_owners_guessed; + std::vector, 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 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 diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 23dd04b907..5d410d900a 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -163,6 +163,14 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS) const std::vector>> &, const std::vector> &); + template std::tuple>, + std::map, + std::map>> + GridTools::guess_point_owner( + const RTree, unsigned int>> + &, + const std::vector> &); + template RTree< std::pair, unsigned int>> GridTools::build_global_description_tree( diff --git a/tests/grid/grid_tools_guess_pt_owner_1.cc b/tests/grid/grid_tools_guess_pt_owner_01.cc similarity index 100% rename from tests/grid/grid_tools_guess_pt_owner_1.cc rename to tests/grid/grid_tools_guess_pt_owner_01.cc diff --git a/tests/grid/grid_tools_guess_pt_owner_1.output b/tests/grid/grid_tools_guess_pt_owner_01.output similarity index 100% rename from tests/grid/grid_tools_guess_pt_owner_1.output rename to tests/grid/grid_tools_guess_pt_owner_01.output 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 index 0000000000..483136152f --- /dev/null +++ b/tests/grid/grid_tools_guess_pt_owner_02.cc @@ -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 +#include + +#include + +#include + +#include "../tests.h" + +template +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, 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> boundaries; + boundaries.first[0] = tot_bbox; + boundaries.second[0] = tot_bbox + 1; + for (int i = 1; i < spacedim; i++) + boundaries.second[i] = 1; + + BoundingBox new_box(boundaries); + global_bboxes.emplace_back(std::make_pair(new_box, rk)); + ++tot_bbox; + } + + // Creating the vector of center points + std::vector> 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, 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 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> 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 new_box(boundaries); + global_bboxes.emplace_back(std::make_pair(new_box, rk)); + } + + // Building a new tree + RTree, 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 index 0000000000..816e8115e0 --- /dev/null +++ b/tests/grid/grid_tools_guess_pt_owner_02.output @@ -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 -- 2.39.5