Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mpi_remote_point_evaluation.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
20
22
23#include <deal.II/fe/mapping.h>
24
28#include <deal.II/grid/tria.h>
29
31
32
33namespace Utilities
34{
35 namespace MPI
36 {
37 template <int dim, int spacedim>
39 const double tolerance,
40 const bool enforce_unique_mapping,
41 const unsigned int rtree_level,
42 const std::function<std::vector<bool>()> &marked_vertices)
43 : tolerance(tolerance)
44 , enforce_unique_mapping(enforce_unique_mapping)
45 , rtree_level(rtree_level)
46 , marked_vertices(marked_vertices)
47 {}
48
49
50
51 template <int dim, int spacedim>
57
58
59
60 template <int dim, int spacedim>
62 const double tolerance,
63 const bool enforce_unique_mapping,
64 const unsigned int rtree_level,
65 const std::function<std::vector<bool>()> &marked_vertices)
66 : additional_data(tolerance,
67 enforce_unique_mapping,
68 rtree_level,
69 marked_vertices)
70 , ready_flag(false)
71 {}
72
73
74
75 template <int dim, int spacedim>
81
82
83
84 template <int dim, int spacedim>
85 void
87 const std::vector<Point<spacedim>> &points,
90 {
92
93 this->reinit(cache, points);
94 }
95
96
97
98 template <int dim, int spacedim>
99 void
102 const std::vector<Point<spacedim>> &points)
103 {
104#ifndef DEAL_II_WITH_MPI
105 Assert(false, ExcNeedsMPI());
106 (void)cache;
107 (void)points;
108#else
109 if (tria_signal.connected())
110 tria_signal.disconnect();
111
112 tria_signal = cache.get_triangulation().signals.any_change.connect(
113 [&]() { this->ready_flag = false; });
114
115 // compress r-tree to a minimal set of bounding boxes
116 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes;
117 global_bboxes.emplace_back(
120
121 const auto data =
123 cache,
124 points,
125 global_bboxes,
127 std::vector<bool>(),
128 additional_data.tolerance,
129 true,
130 additional_data.enforce_unique_mapping);
131
132 this->reinit(data, cache.get_triangulation(), cache.get_mapping());
133#endif
134 }
135
136
137
138 template <int dim, int spacedim>
139 void
141 const GridTools::internal::
142 DistributedComputePointLocationsInternal<dim, spacedim> &data,
143 const Triangulation<dim, spacedim> &tria,
144 const Mapping<dim, spacedim> &mapping)
145 {
146 this->tria = &tria;
147 this->mapping = &mapping;
148
149 this->recv_ranks = data.recv_ranks;
150 this->recv_ptrs = data.recv_ptrs;
151
152 this->send_ranks = data.send_ranks;
153 this->send_ptrs = data.send_ptrs;
154
155 this->recv_permutation = {};
156 this->recv_permutation.resize(data.recv_components.size());
157 this->point_ptrs.assign(data.n_searched_points + 1, 0);
158 for (unsigned int i = 0; i < data.recv_components.size(); ++i)
159 {
160 AssertIndexRange(std::get<2>(data.recv_components[i]),
161 this->recv_permutation.size());
162 this->recv_permutation[std::get<2>(data.recv_components[i])] = i;
163
164 AssertIndexRange(std::get<1>(data.recv_components[i]) + 1,
165 this->point_ptrs.size());
166 this->point_ptrs[std::get<1>(data.recv_components[i]) + 1]++;
167 }
169 std::tuple<unsigned int, unsigned int> n_owning_processes_default{
171 std::tuple<unsigned int, unsigned int> n_owning_processes_local =
172 n_owning_processes_default;
173
174 for (unsigned int i = 0; i < data.n_searched_points; ++i)
175 {
176 std::get<0>(n_owning_processes_local) =
177 std::min(std::get<0>(n_owning_processes_local),
178 this->point_ptrs[i + 1]);
179 std::get<1>(n_owning_processes_local) =
180 std::max(std::get<1>(n_owning_processes_local),
181 this->point_ptrs[i + 1]);
182
183 this->point_ptrs[i + 1] += this->point_ptrs[i];
184 }
185
186 const auto n_owning_processes_global =
188 n_owning_processes_local,
189 tria.get_communicator(),
190 [&](const auto &a,
191 const auto &b) -> std::tuple<unsigned int, unsigned int> {
192 if (a == n_owning_processes_default)
193 return b;
194
195 if (b == n_owning_processes_default)
196 return a;
197
198 return std::tuple<unsigned int, unsigned int>{
199 std::min(std::get<0>(a), std::get<0>(b)),
200 std::max(std::get<1>(a), std::get<1>(b))};
201 });
202
203 if (n_owning_processes_global == n_owning_processes_default)
204 {
205 unique_mapping = true;
207 }
208 else
209 {
210 unique_mapping = (std::get<0>(n_owning_processes_global) == 1) &&
211 (std::get<1>(n_owning_processes_global) == 1);
212 all_points_found_flag = std::get<0>(n_owning_processes_global) > 0;
213 }
214
217
218 cell_data = std::make_unique<CellData>(tria);
219 send_permutation = {};
220
221 std::pair<int, int> dummy{-1, -1};
222 for (const auto &i : data.send_components)
223 {
224 if (dummy != std::get<0>(i))
225 {
226 dummy = std::get<0>(i);
227 cell_data->cells.emplace_back(dummy);
228 cell_data->reference_point_ptrs.emplace_back(
229 cell_data->reference_point_values.size());
230 }
231
232 cell_data->reference_point_values.emplace_back(std::get<3>(i));
233 send_permutation.emplace_back(std::get<5>(i));
234 }
235
236 cell_data->reference_point_ptrs.emplace_back(
237 cell_data->reference_point_values.size());
239 unsigned int max_size_recv = 0;
240 for (unsigned int i = 0; i < recv_ranks.size(); ++i)
241 max_size_recv =
242 std::max(max_size_recv, recv_ptrs[i + 1] - recv_ptrs[i]);
243
244 unsigned int max_size_send = 0;
245 for (unsigned int i = 0; i < send_ranks.size(); ++i)
246 max_size_send =
247 std::max(max_size_send, send_ptrs[i + 1] - send_ptrs[i]);
248
250 std::max(send_permutation.size() * 2 + max_size_recv,
251 point_ptrs.back() + send_permutation.size() + max_size_send);
252
254
255 // invert permutation matrices
257 for (unsigned int c = 0; c < recv_permutation.size(); ++c)
259
261 for (unsigned int c = 0; c < send_permutation.size(); ++c)
263
264 this->ready_flag = true;
265 }
266
267
268
269 template <int dim, int spacedim>
274
275
276
277 template <int dim, int spacedim>
280 {
281 return {0, static_cast<unsigned int>(cells.size())};
282 }
283
284
285
286 template <int dim, int spacedim>
289 const unsigned int cell) const
290 {
291 AssertIndexRange(cell, cells.size());
292 return {&triangulation, cells[cell].first, cells[cell].second};
293 }
294
295
296
297 template <int dim, int spacedim>
300 const unsigned int cell) const
301 {
302 AssertIndexRange(cell, cells.size());
303 return {reference_point_values.data() + reference_point_ptrs[cell],
304 reference_point_ptrs[cell + 1] - reference_point_ptrs[cell]};
305 }
306
307
308
309 template <int dim, int spacedim>
315
316
317
318 template <int dim, int spacedim>
319 const std::vector<unsigned int> &
324
325
326
327 template <int dim, int spacedim>
328 bool
333
334
335
336 template <int dim, int spacedim>
337 bool
342
343
344
345 template <int dim, int spacedim>
346 bool
348 const unsigned int i) const
349 {
350 AssertIndexRange(i, point_ptrs.size() - 1);
351
353 return true;
354 else
355 return (point_ptrs[i + 1] - point_ptrs[i]) > 0;
356 }
357
358
360 template <int dim, int spacedim>
366
367
369 template <int dim, int spacedim>
375
376
377
378 template <int dim, int spacedim>
379 bool
384
385
387 template <int dim, int spacedim>
388 const std::vector<unsigned int> &
393
394
395
396 template <int dim, int spacedim>
397 const std::vector<unsigned int> &
399 {
401 }
402
403 } // end of namespace MPI
404} // end of namespace Utilities
405
406#include "mpi_remote_point_evaluation.inst"
407
const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator > > & get_locally_owned_cell_bounding_boxes_rtree() const
const Mapping< dim, spacedim > & get_mapping() const
const Triangulation< dim, spacedim > & get_triangulation() const
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
Signals signals
Definition tria.h:2513
std_cxx20::ranges::iota_view< unsigned int, unsigned int > cell_indices() const
ArrayView< const Point< dim > > get_unit_points(const unsigned int cell) const
Triangulation< dim, spacedim >::active_cell_iterator get_active_cell_iterator(const unsigned int cell) const
CellData(const Triangulation< dim, spacedim > &triangulation)
const std::vector< unsigned int > & get_send_permutation() const
RemotePointEvaluation(const AdditionalData &additional_data=AdditionalData())
const Triangulation< dim, spacedim > & get_triangulation() const
SmartPointer< const Mapping< dim, spacedim > > mapping
const std::vector< unsigned int > & get_point_ptrs() const
const std::vector< unsigned int > & get_inverse_recv_permutation() const
const Mapping< dim, spacedim > & get_mapping() const
void reinit(const std::vector< Point< spacedim > > &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
SmartPointer< const Triangulation< dim, spacedim > > tria
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim > > &points, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
T all_reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner)
static const unsigned int invalid_unsigned_int
Definition types.h:220
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< BoundingBox< boost::geometry::dimension< typename Rtree::indexable_type >::value > > extract_rtree_level(const Rtree &tree, const unsigned int level)
AdditionalData(const double tolerance=1e-6, const bool enforce_unique_mapping=false, const unsigned int rtree_level=0, const std::function< std::vector< bool >()> &marked_vertices={})