Reference documentation for deal.II version 9.3.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\}}\)
mpi_remote_point_evaluation.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 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 
20 #include <deal.II/base/mpi_consensus_algorithms.templates.h>
22 
24 
25 #include <deal.II/fe/mapping.h>
26 
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  : tolerance(tolerance)
43  , enforce_unique_mapping(enforce_unique_mapping)
44  , ready_flag(false)
45  {}
46 
47 
48 
49  template <int dim, int spacedim>
51  {
52  if (tria_signal.connected())
53  tria_signal.disconnect();
54  }
55 
56 
57 
58  template <int dim, int spacedim>
59  void
61  const std::vector<Point<spacedim>> &points,
64  {
65 #ifndef DEAL_II_WITH_MPI
66  Assert(false, ExcNeedsMPI());
67  (void)points;
68  (void)tria;
69  (void)mapping;
70 #else
71  if (tria_signal.connected())
72  tria_signal.disconnect();
73 
74  tria_signal =
75  tria.signals.any_change.connect([&]() { this->ready_flag = false; });
76 
77  this->tria = &tria;
78  this->mapping = &mapping;
79 
80  std::vector<BoundingBox<spacedim>> local_boxes;
81  for (const auto &cell : tria.active_cell_iterators())
82  if (cell->is_locally_owned())
83  local_boxes.push_back(mapping.get_bounding_box(cell));
84 
85  // create r-tree of bounding boxes
86  const auto local_tree = pack_rtree(local_boxes);
87 
88  // compress r-tree to a minimal set of bounding boxes
89  const auto local_reduced_box = extract_rtree_level(local_tree, 0);
90 
91  // gather bounding boxes of other processes
92  const auto global_bboxes =
93  Utilities::MPI::all_gather(tria.get_communicator(), local_reduced_box);
94 
95  const GridTools::Cache<dim, spacedim> cache(tria, mapping);
96 
97  const auto data =
99  cache,
100  points,
101  global_bboxes,
102  tolerance,
103  true,
105 
106  this->recv_ranks = data.recv_ranks;
107  this->recv_ptrs = data.recv_ptrs;
108 
109  this->send_ranks = data.send_ranks;
110  this->send_ptrs = data.send_ptrs;
111 
112  this->recv_permutation = {};
113  this->recv_permutation.resize(data.recv_components.size());
114  this->point_ptrs.assign(points.size() + 1, 0);
115  for (unsigned int i = 0; i < data.recv_components.size(); ++i)
116  {
117  AssertIndexRange(std::get<2>(data.recv_components[i]),
118  this->recv_permutation.size());
119  this->recv_permutation[std::get<2>(data.recv_components[i])] = i;
120 
121  AssertIndexRange(std::get<1>(data.recv_components[i]) + 1,
122  this->point_ptrs.size());
123  this->point_ptrs[std::get<1>(data.recv_components[i]) + 1]++;
124  }
125 
126  unique_mapping = true;
127  for (unsigned int i = 0; i < points.size(); ++i)
128  {
129  if (unique_mapping && this->point_ptrs[i + 1] != 1)
130  unique_mapping = false;
131  this->point_ptrs[i + 1] += this->point_ptrs[i];
132  }
133 
134  cell_data = {};
135  send_permutation = {};
136 
137  std::pair<int, int> dummy{-1, -1};
138  for (const auto &i : data.send_components)
139  {
140  if (dummy != std::get<0>(i))
141  {
142  dummy = std::get<0>(i);
143  cell_data.cells.emplace_back(dummy);
144  cell_data.reference_point_ptrs.emplace_back(
146  }
147 
148  cell_data.reference_point_values.emplace_back(std::get<3>(i));
149  send_permutation.emplace_back(std::get<5>(i));
150  }
151 
152  cell_data.reference_point_ptrs.emplace_back(
154 
155  this->ready_flag = true;
156 #endif
157  }
158 
159 
160  template <int dim, int spacedim>
161  const std::vector<unsigned int> &
163  {
164  return point_ptrs;
165  }
166 
167 
168 
169  template <int dim, int spacedim>
170  bool
172  {
173  return unique_mapping;
174  }
175 
176 
177 
178  template <int dim, int spacedim>
181  {
182  return *tria;
183  }
184 
185 
186 
187  template <int dim, int spacedim>
188  const Mapping<dim, spacedim> &
190  {
191  return *mapping;
192  }
193 
194 
195 
196  template <int dim, int spacedim>
197  bool
199  {
200  return ready_flag;
201  }
202 
203  } // end of namespace MPI
204 } // end of namespace Utilities
205 
206 #include "mpi_remote_point_evaluation.inst"
207 
static ::ExceptionBase & ExcNeedsMPI()
virtual MPI_Comm get_communicator() const
Definition: tria.cc:10124
void reinit(const std::vector< Point< spacedim >> &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
RemotePointEvaluation(const double tolerance=1e-6, const bool enforce_unique_mapping=false)
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12154
std::vector< BoundingBox< boost::geometry::dimension< typename Rtree::indexable_type >::value > > extract_rtree_level(const Rtree &tree, const unsigned int level)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
const std::vector< unsigned int > & get_point_ptrs() const
SmartPointer< const Triangulation< dim, spacedim > > tria
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
SmartPointer< const Mapping< dim, spacedim > > mapping
#define Assert(cond, exc)
Definition: exceptions.h:1465
Signals signals
Definition: tria.h:2295
Abstract base class for mapping classes.
Definition: mapping.h:303
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
const Mapping< dim, spacedim > & get_mapping() const
const Triangulation< dim, spacedim > & get_triangulation() const
virtual BoundingBox< spacedim > get_bounding_box(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: mapping.cc:74
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
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 double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition: grid_tools.cc:5779
std::vector< T > all_gather(const MPI_Comm &comm, const T &object_to_send)