From: Peter Munch Date: Wed, 28 Jul 2021 13:07:22 +0000 (+0200) Subject: Fix tolerances in GridTools::internal::distributed_compute_point_locations() X-Git-Tag: v9.4.0-rc1~1118^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12602%2Fhead;p=dealii.git Fix tolerances in GridTools::internal::distributed_compute_point_locations() --- diff --git a/source/base/mpi_remote_point_evaluation.cc b/source/base/mpi_remote_point_evaluation.cc index d6258e1e76..e146270219 100644 --- a/source/base/mpi_remote_point_evaluation.cc +++ b/source/base/mpi_remote_point_evaluation.cc @@ -134,6 +134,9 @@ namespace Utilities this->point_ptrs[i + 1] += this->point_ptrs[i]; } + Assert(enforce_unique_mapping == false || unique_mapping, + ExcInternalError()); + cell_data = {}; send_permutation = {}; diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index f2266b99b5..08c648658f 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2786,8 +2786,20 @@ namespace GridTools if (relevant_cell_bounding_boxes_rtree != nullptr && !relevant_cell_bounding_boxes_rtree->empty()) { + // create a bounding box around point p with 2*tolerance as side length. + auto p1 = p; + auto p2 = p; + + for (int d = 0; d < spacedim; ++d) + { + p1[d] = p1[d] - tolerance; + p2[d] = p2[d] + tolerance; + } + + BoundingBox bb({p1, p2}); + if (relevant_cell_bounding_boxes_rtree->qbegin( - boost::geometry::index::intersects(p)) == + boost::geometry::index::intersects(bb)) == relevant_cell_bounding_boxes_rtree->qend()) return cell_and_position; } @@ -5731,7 +5743,8 @@ namespace GridTools std::vector> guess_point_owner( const std::vector>> &global_bboxes, - const std::vector> & points) + const std::vector> & points, + const double tolerance) { std::vector> ranks_and_indices; ranks_and_indices.reserve(points.size()); @@ -5741,7 +5754,7 @@ namespace GridTools const auto &point = points[i]; for (unsigned rank = 0; rank < global_bboxes.size(); ++rank) for (const auto &box : global_bboxes[rank]) - if (box.point_inside(point)) + if (box.point_inside(point, tolerance)) { ranks_and_indices.emplace_back(rank, i); break; @@ -5876,7 +5889,7 @@ namespace GridTools auto &recv_ptrs = result.recv_ptrs; const auto potential_owners = - internal::guess_point_owner(global_bboxes, points); + internal::guess_point_owner(global_bboxes, points, tolerance); const auto &potential_owners_ranks = std::get<0>(potential_owners); const auto &potential_owners_ptrs = std::get<1>(potential_owners); diff --git a/tests/remote_point_evaluation/mapping_03.cc b/tests/remote_point_evaluation/mapping_03.cc new file mode 100644 index 0000000000..54bcfc811b --- /dev/null +++ b/tests/remote_point_evaluation/mapping_03.cc @@ -0,0 +1,93 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +// Test Utilities::MPI::RemotePointEvaluation with and without unique mapping +// for a within a tolerance. + +#include + +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + + +template +void +test(const bool enforce_unique_map) +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::subdivided_hyper_cube(tria, 2); + tria.refine_global(1); + + FE_Q fe(1); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vec(dof_handler.n_dofs()); + + vec = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) + 1; + + MappingQ1 mapping; + + std::vector> evaluation_points; + + evaluation_points.emplace_back(0.5, 0.5 + 1e-7); // should be assigned to rank + // 0 for unique mapping + + Utilities::MPI::RemotePointEvaluation eval(1e-6, enforce_unique_map); + + const auto result_avg = + VectorTools::point_values<1>(mapping, + dof_handler, + vec, + evaluation_points, + eval, + VectorTools::EvaluationFlags::avg); + const auto result_min = VectorTools::point_values<1>( + eval, dof_handler, vec, VectorTools::EvaluationFlags::min); + const auto result_max = VectorTools::point_values<1>( + eval, dof_handler, vec, VectorTools::EvaluationFlags::max); + const auto result_insert = VectorTools::point_values<1>( + eval, dof_handler, vec, VectorTools::EvaluationFlags::insert); + + for (unsigned int i = 0; i < evaluation_points.size(); ++i) + deallog << i << " " << result_avg[i] << " " << result_min[i] << " " + << result_max[i] << " " << result_insert[i] << " " + << eval.get_point_ptrs()[i + 1] - eval.get_point_ptrs()[i] + << std::endl; +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + MPILogInitAll all; + + test<2>(false); + test<2>(true); +} diff --git a/tests/remote_point_evaluation/mapping_03.with_p4est=true.mpirun=2.output b/tests/remote_point_evaluation/mapping_03.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..9529db7f95 --- /dev/null +++ b/tests/remote_point_evaluation/mapping_03.with_p4est=true.mpirun=2.output @@ -0,0 +1,7 @@ + +DEAL:0::0 1.50000 1.00000 2.00000 1.00000 4 +DEAL:0::0 1.00000 1.00000 1.00000 1.00000 1 + +DEAL:1::0 1.50000 1.00000 2.00000 1.00000 4 +DEAL:1::0 1.00000 1.00000 1.00000 1.00000 1 +