Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10:02+00:00
\(\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\}}\)
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 
34 namespace 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; });
80 
81  std::vector<BoundingBox<spacedim>> local_boxes;
82  for (const auto &cell :
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);
94 
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  }
108 
109 
110 
111  template <int dim, int spacedim>
112  void
114  const GridTools::internal::
115  DistributedComputePointLocationsInternal<dim, spacedim> &data,
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]++;
140  }
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,
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;
179  all_points_found_flag = 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 
188  Assert(enforce_unique_mapping == false || unique_mapping,
189  ExcInternalError());
190 
191  cell_data = std::make_unique<CellData>(tria);
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(
202  cell_data->reference_point_values.size());
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(
210  cell_data->reference_point_values.size());
211 
212  this->ready_flag = true;
213  }
214 
215 
216 
217  template <int dim, int spacedim>
221  {}
222 
223 
224 
225  template <int dim, int spacedim>
228  {
229  return {0, static_cast<unsigned int>(cells.size())};
230  }
231 
232 
233 
234  template <int dim, int spacedim>
237  const unsigned int cell) const
238  {
239  AssertIndexRange(cell, cells.size());
240  return {&triangulation, cells[cell].first, cells[cell].second};
241  }
242 
243 
244 
245  template <int dim, int spacedim>
248  const unsigned int cell) const
249  {
250  AssertIndexRange(cell, cells.size());
251  return {reference_point_values.data() + reference_point_ptrs[cell],
252  reference_point_ptrs[cell + 1] - reference_point_ptrs[cell]};
253  }
254 
255 
256 
257  template <int dim, int spacedim>
260  {
261  return *cell_data;
262  }
263 
264 
265 
266  template <int dim, int spacedim>
267  const std::vector<unsigned int> &
269  {
270  return point_ptrs;
271  }
272 
273 
274 
275  template <int dim, int spacedim>
276  bool
278  {
280  }
281 
282 
283 
284  template <int dim, int spacedim>
285  bool
287  {
288  return all_points_found_flag;
289  }
290 
291 
292 
293  template <int dim, int spacedim>
294  bool
296  const unsigned int i) const
297  {
298  AssertIndexRange(i, point_ptrs.size() - 1);
299 
301  return true;
302  else
303  return (point_ptrs[i + 1] - point_ptrs[i]) > 0;
304  }
305 
306 
307 
308  template <int dim, int spacedim>
311  {
312  return *tria;
313  }
314 
315 
316 
317  template <int dim, int spacedim>
318  const Mapping<dim, spacedim> &
320  {
321  return *mapping;
322  }
323 
324 
325 
326  template <int dim, int spacedim>
327  bool
329  {
330  return ready_flag;
331  }
332 
333  } // end of namespace MPI
334 } // end of namespace Utilities
335 
336 #include "mpi_remote_point_evaluation.inst"
337 
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual MPI_Comm get_communicator() const
Signals signals
Definition: tria.h:2528
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
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
SmartPointer< const Triangulation< dim, spacedim > > tria
bool point_found(const unsigned int i) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
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)
Definition: grid_tools.cc:3942
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static const unsigned int invalid_unsigned_int
Definition: types.h:221
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
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)
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
const ::Triangulation< dim, spacedim > & tria