Reference documentation for deal.II version 9.5.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// Copyright (C) 2021 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#include <deal.II/base/config.h>
17
21
23
24#include <deal.II/fe/mapping.h>
25
29#include <deal.II/grid/tria.h>
30
32
33
34namespace Utilities
35{
36 namespace MPI
37 {
38 template <int dim, int spacedim>
40 const double tolerance,
41 const bool enforce_unique_mapping,
42 const unsigned int rtree_level,
43 const std::function<std::vector<bool>()> &marked_vertices)
44 : tolerance(tolerance)
45 , enforce_unique_mapping(enforce_unique_mapping)
46 , rtree_level(rtree_level)
47 , marked_vertices(marked_vertices)
48 , ready_flag(false)
49 {}
50
51
52
53 template <int dim, int spacedim>
55 {
56 if (tria_signal.connected())
57 tria_signal.disconnect();
58 }
59
60
61
62 template <int dim, int spacedim>
63 void
65 const std::vector<Point<spacedim>> &points,
67 const Mapping<dim, spacedim> & mapping)
68 {
69#ifndef DEAL_II_WITH_MPI
70 Assert(false, ExcNeedsMPI());
71 (void)points;
72 (void)tria;
73 (void)mapping;
74#else
75 if (tria_signal.connected())
76 tria_signal.disconnect();
77
78 tria_signal =
79 tria.signals.any_change.connect([&]() { this->ready_flag = false; });
81 std::vector<BoundingBox<spacedim>> local_boxes;
82 for (const auto &cell :
83 tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
84 local_boxes.push_back(mapping.get_bounding_box(cell));
85
86 // create r-tree of bounding boxes
87 const auto local_tree = pack_rtree(local_boxes);
88
89 // compress r-tree to a minimal set of bounding boxes
90 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(1);
91 global_bboxes[0] = extract_rtree_level(local_tree, rtree_level);
92
93 const GridTools::Cache<dim, spacedim> cache(tria, mapping);
95 const auto data =
97 cache,
98 points,
99 global_bboxes,
100 marked_vertices ? marked_vertices() : std::vector<bool>(),
101 tolerance,
102 true,
103 enforce_unique_mapping);
104
105 this->reinit(data, tria, mapping);
106#endif
107 }
109
110
111 template <int dim, int spacedim>
112 void
114 const GridTools::internal::
115 DistributedComputePointLocationsInternal<dim, spacedim> &data,
116 const Triangulation<dim, spacedim> & tria,
117 const Mapping<dim, spacedim> & mapping)
118 {
119 this->tria = &tria;
120 this->mapping = &mapping;
121
122 this->recv_ranks = data.recv_ranks;
123 this->recv_ptrs = data.recv_ptrs;
124
125 this->send_ranks = data.send_ranks;
126 this->send_ptrs = data.send_ptrs;
127
128 this->recv_permutation = {};
129 this->recv_permutation.resize(data.recv_components.size());
130 this->point_ptrs.assign(data.n_searched_points + 1, 0);
131 for (unsigned int i = 0; i < data.recv_components.size(); ++i)
132 {
133 AssertIndexRange(std::get<2>(data.recv_components[i]),
134 this->recv_permutation.size());
135 this->recv_permutation[std::get<2>(data.recv_components[i])] = i;
136
137 AssertIndexRange(std::get<1>(data.recv_components[i]) + 1,
138 this->point_ptrs.size());
139 this->point_ptrs[std::get<1>(data.recv_components[i]) + 1]++;
141
142 std::tuple<unsigned int, unsigned int> n_owning_processes_default{
144 std::tuple<unsigned int, unsigned int> n_owning_processes_local =
145 n_owning_processes_default;
146
147 for (unsigned int i = 0; i < data.n_searched_points; ++i)
148 {
149 std::get<0>(n_owning_processes_local) =
150 std::min(std::get<0>(n_owning_processes_local),
151 this->point_ptrs[i + 1]);
152 std::get<1>(n_owning_processes_local) =
153 std::max(std::get<1>(n_owning_processes_local),
154 this->point_ptrs[i + 1]);
155
156 this->point_ptrs[i + 1] += this->point_ptrs[i];
157 }
158
159 const auto n_owning_processes_global =
160 Utilities::MPI::all_reduce<std::tuple<unsigned int, unsigned int>>(
161 n_owning_processes_local,
162 tria.get_communicator(),
163 [&](const auto &a,
164 const auto &b) -> std::tuple<unsigned int, unsigned int> {
165 if (a == n_owning_processes_default)
166 return b;
167
168 if (b == n_owning_processes_default)
169 return a;
170
171 return std::tuple<unsigned int, unsigned int>{
172 std::min(std::get<0>(a), std::get<0>(b)),
173 std::max(std::get<1>(a), std::get<1>(b))};
174 });
175
176 if (n_owning_processes_global == n_owning_processes_default)
177 {
178 unique_mapping = true;
180 }
181 else
182 {
183 unique_mapping = (std::get<0>(n_owning_processes_global) == 1) &&
184 (std::get<1>(n_owning_processes_global) == 1);
185 all_points_found_flag = std::get<0>(n_owning_processes_global) > 0;
186 }
187
190
191 cell_data = {};
192 send_permutation = {};
193
194 std::pair<int, int> dummy{-1, -1};
195 for (const auto &i : data.send_components)
196 {
197 if (dummy != std::get<0>(i))
198 {
199 dummy = std::get<0>(i);
200 cell_data.cells.emplace_back(dummy);
201 cell_data.reference_point_ptrs.emplace_back(
203 }
204
205 cell_data.reference_point_values.emplace_back(std::get<3>(i));
206 send_permutation.emplace_back(std::get<5>(i));
207 }
208
209 cell_data.reference_point_ptrs.emplace_back(
211
212 this->ready_flag = true;
213 }
214
215
216
217 template <int dim, int spacedim>
220 {
221 return cell_data;
222 }
223
224
225
226 template <int dim, int spacedim>
227 const std::vector<unsigned int> &
229 {
230 return point_ptrs;
231 }
232
233
234
235 template <int dim, int spacedim>
236 bool
238 {
239 return unique_mapping;
240 }
241
242
243
244 template <int dim, int spacedim>
245 bool
247 {
248 return all_points_found_flag;
249 }
250
251
252
253 template <int dim, int spacedim>
254 bool
256 const unsigned int i) const
257 {
258 AssertIndexRange(i, point_ptrs.size() - 1);
259
260 if (all_points_found_flag)
261 return true;
262 else
263 return (point_ptrs[i + 1] - point_ptrs[i]) > 0;
264 }
265
266
267
268 template <int dim, int spacedim>
271 {
272 return *tria;
273 }
274
275
276
277 template <int dim, int spacedim>
280 {
281 return *mapping;
282 }
283
284
285
286 template <int dim, int spacedim>
287 bool
289 {
290 return ready_flag;
291 }
292
293 } // end of namespace MPI
294} // end of namespace Utilities
295
296#include "mpi_remote_point_evaluation.inst"
297
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
Signals signals
Definition tria.h:2479
RemotePointEvaluation(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={})
const Triangulation< dim, spacedim > & get_triangulation() const
SmartPointer< const Mapping< dim, spacedim > > mapping
const std::vector< unsigned int > & get_point_ptrs() 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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)
static const unsigned int invalid_unsigned_int
Definition types.h:213
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 > &)
std::vector< BoundingBox< boost::geometry::dimension< typename Rtree::indexable_type >::value > > extract_rtree_level(const Rtree &tree, const unsigned int level)
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
const ::Triangulation< dim, spacedim > & tria