Reference documentation for deal.II version Git d780fd0cc4 2019-11-15 10:11:22 -0500
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tria_base.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2019 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 
17 #include <deal.II/base/logstream.h>
18 #include <deal.II/base/memory_consumption.h>
19 #include <deal.II/base/utilities.h>
20 
21 #include <deal.II/distributed/tria_base.h>
22 
23 #include <deal.II/grid/grid_tools.h>
24 #include <deal.II/grid/tria.h>
25 #include <deal.II/grid/tria_accessor.h>
26 #include <deal.II/grid/tria_iterator.h>
27 
28 #include <deal.II/lac/sparsity_pattern.h>
29 #include <deal.II/lac/sparsity_tools.h>
30 #include <deal.II/lac/vector_memory.h>
31 
32 #include <algorithm>
33 #include <fstream>
34 #include <iostream>
35 #include <numeric>
36 
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 namespace parallel
41 {
42  template <int dim, int spacedim>
44  MPI_Comm mpi_communicator,
45  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
46  smooth_grid,
47  const bool check_for_distorted_cells)
48  : ::Triangulation<dim, spacedim>(smooth_grid,
49  check_for_distorted_cells)
50  , mpi_communicator(mpi_communicator)
51  , my_subdomain(Utilities::MPI::this_mpi_process(this->mpi_communicator))
52  , n_subdomains(Utilities::MPI::n_mpi_processes(this->mpi_communicator))
53  {
54 #ifndef DEAL_II_WITH_MPI
55  Assert(false,
56  ExcMessage("You compiled deal.II without MPI support, for "
57  "which parallel::TriangulationBase is not available."));
58 #endif
59  }
60 
61 
62 
63  template <int dim, int spacedim>
64  void
66  const ::Triangulation<dim, spacedim> &other_tria)
67  {
68 #ifndef DEAL_II_WITH_MPI
69  (void)other_tria;
70  Assert(false,
71  ExcMessage("You compiled deal.II without MPI support, for "
72  "which parallel::TriangulationBase is not available."));
73 #else
75 
76  if (const ::parallel::TriangulationBase<dim, spacedim> *other_tria_x =
77  dynamic_cast<const ::parallel::TriangulationBase<dim, spacedim>
78  *>(&other_tria))
79  {
80  // release unused vector memory because we will have very different
81  // vectors now
82  ::internal::GrowingVectorMemoryImplementation::
83  release_all_unused_memory();
84  }
85 #endif
86  }
87 
88 
89 
90  template <int dim, int spacedim>
91  std::size_t
93  {
94  std::size_t mem =
99  number_cache.n_global_active_cells) +
101  return mem;
102  }
103 
104  template <int dim, int spacedim>
106  {
107  // release unused vector memory because the vector layout is going to look
108  // very different now
109  ::internal::GrowingVectorMemoryImplementation::
110  release_all_unused_memory();
111  }
112 
113  template <int dim, int spacedim>
116  , n_global_levels(0)
117  {}
118 
119  template <int dim, int spacedim>
120  unsigned int
122  {
123  return number_cache.n_locally_owned_active_cells;
124  }
125 
126  template <int dim, int spacedim>
127  unsigned int
129  {
130  return number_cache.n_global_levels;
131  }
132 
133  template <int dim, int spacedim>
136  {
137  return number_cache.n_global_active_cells;
138  }
139 
140  template <int dim, int spacedim>
141  std::vector<unsigned int>
144  {
145  ;
146 #ifdef DEAL_II_WITH_MPI
147  std::vector<unsigned int> n_locally_owned_active_cells_per_processor(
149 
150  if (this->n_levels() > 0)
151  {
152  const int ierr =
153  MPI_Allgather(&number_cache.n_locally_owned_active_cells,
154  1,
155  MPI_UNSIGNED,
156  n_locally_owned_active_cells_per_processor.data(),
157  1,
158  MPI_UNSIGNED,
159  this->mpi_communicator);
160  AssertThrowMPI(ierr);
161  }
162 
163  return n_locally_owned_active_cells_per_processor;
164 #else
165  return {number_cache.n_locally_owned_active_cells};
166 #endif
167  }
168 
169  template <int dim, int spacedim>
170  const MPI_Comm &
172  {
173  return mpi_communicator;
174  }
175 
176 #ifdef DEAL_II_WITH_MPI
177  template <int dim, int spacedim>
178  void
180  {
181  number_cache.ghost_owners.clear();
182  number_cache.level_ghost_owners.clear();
183  number_cache.n_locally_owned_active_cells = 0;
184 
185  if (this->n_levels() == 0)
186  {
187  // Skip communication done below if we do not have any cells
188  // (meaning the Triangulation is empty on all processors). This will
189  // happen when called from the destructor of Triangulation, which
190  // can get called during exception handling causing a hang in this
191  // function.
192  number_cache.n_global_active_cells = 0;
193  number_cache.n_global_levels = 0;
194  return;
195  }
196 
197 
198  {
199  // find ghost owners
201  this->begin_active();
202  cell != this->end();
203  ++cell)
204  if (cell->is_ghost())
205  number_cache.ghost_owners.insert(cell->subdomain_id());
206 
207  Assert(number_cache.ghost_owners.size() <
209  ExcInternalError());
210  }
211 
212  if (this->n_levels() > 0)
214  this->begin_active();
215  cell != this->end();
216  ++cell)
217  if (cell->subdomain_id() == my_subdomain)
218  ++number_cache.n_locally_owned_active_cells;
219 
220  number_cache.n_global_active_cells =
222  this->mpi_communicator);
223  number_cache.n_global_levels =
225  }
226 
227 
228 
229  template <int dim, int spacedim>
230  void
232  {
233  number_cache.level_ghost_owners.clear();
234 
235  // if there is nothing to do, then do nothing
236  if (this->n_levels() == 0)
237  return;
238 
239  // find level ghost owners
241  this->begin();
242  cell != this->end();
243  ++cell)
244  if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
245  cell->level_subdomain_id() != this->locally_owned_subdomain())
246  this->number_cache.level_ghost_owners.insert(
247  cell->level_subdomain_id());
248 
249 # ifdef DEBUG
250  // Check that level_ghost_owners is symmetric by sending a message to
251  // everyone
252  {
253  int ierr = MPI_Barrier(this->mpi_communicator);
254  AssertThrowMPI(ierr);
255 
256  const int mpi_tag = Utilities::MPI::internal::Tags::
258 
259  // important: preallocate to avoid (re)allocation:
260  std::vector<MPI_Request> requests(
261  this->number_cache.level_ghost_owners.size());
262  unsigned int dummy = 0;
263  unsigned int req_counter = 0;
264 
265  for (std::set<types::subdomain_id>::iterator it =
266  this->number_cache.level_ghost_owners.begin();
267  it != this->number_cache.level_ghost_owners.end();
268  ++it, ++req_counter)
269  {
270  ierr = MPI_Isend(&dummy,
271  1,
272  MPI_UNSIGNED,
273  *it,
274  mpi_tag,
275  this->mpi_communicator,
276  &requests[req_counter]);
277  AssertThrowMPI(ierr);
278  }
279 
280  for (std::set<types::subdomain_id>::iterator it =
281  this->number_cache.level_ghost_owners.begin();
282  it != this->number_cache.level_ghost_owners.end();
283  ++it)
284  {
285  unsigned int dummy;
286  ierr = MPI_Recv(&dummy,
287  1,
288  MPI_UNSIGNED,
289  *it,
290  mpi_tag,
291  this->mpi_communicator,
292  MPI_STATUS_IGNORE);
293  AssertThrowMPI(ierr);
294  }
295 
296  if (requests.size() > 0)
297  {
298  ierr =
299  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
300  AssertThrowMPI(ierr);
301  }
302 
303  ierr = MPI_Barrier(this->mpi_communicator);
304  AssertThrowMPI(ierr);
305  }
306 # endif
307 
308  Assert(this->number_cache.level_ghost_owners.size() <
310  ExcInternalError());
311  }
312 
313 #else
314 
315  template <int dim, int spacedim>
316  void
318  {
319  Assert(false, ExcNotImplemented());
320  }
321 
322  template <int dim, int spacedim>
323  void
325  {
326  Assert(false, ExcNotImplemented());
327  }
328 
329 #endif
330 
331  template <int dim, int spacedim>
334  {
335  return my_subdomain;
336  }
337 
338 
339 
340  template <int dim, int spacedim>
341  const std::set<types::subdomain_id> &
343  {
344  return number_cache.ghost_owners;
345  }
346 
347 
348 
349  template <int dim, int spacedim>
350  const std::set<types::subdomain_id> &
352  {
353  return number_cache.level_ghost_owners;
354  }
355 
356 
357 
358  template <int dim, int spacedim>
359  std::map<unsigned int, std::set<::types::subdomain_id>>
361  const
362  {
363  // 1) collect for each vertex on periodic faces all vertices it coincides
364  // with
365  std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
366  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
367 
369  coinciding_vertex_groups,
370  vertex_to_coinciding_vertex_group);
371 
372  // 2) collect vertices belonging to local cells
373  std::vector<bool> vertex_of_own_cell(this->n_vertices(), false);
374  for (const auto &cell : this->active_cell_iterators())
375  if (cell->is_locally_owned())
376  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
377  vertex_of_own_cell[cell->vertex_index(v)] = true;
378 
379  // 3) for each vertex belonging to a locally owned cell all ghost
380  // neighbors (including the periodic own)
381  std::map<unsigned int, std::set<types::subdomain_id>> result;
382 
383  // loop over all active ghost cells
384  for (const auto &cell : this->active_cell_iterators())
385  if (cell->is_ghost())
386  {
387  const types::subdomain_id owner = cell->subdomain_id();
388 
389  // loop over all its vertices
390  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
391  ++v)
392  {
393  // set owner if vertex belongs to a local cell
394  if (vertex_of_own_cell[cell->vertex_index(v)])
395  result[cell->vertex_index(v)].insert(owner);
396 
397  // mark also nodes coinciding due to periodicity
398  auto coinciding_vertex_group =
399  vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
400  if (coinciding_vertex_group !=
401  vertex_to_coinciding_vertex_group.end())
402  for (auto coinciding_vertex :
403  coinciding_vertex_groups[coinciding_vertex_group->second])
404  if (vertex_of_own_cell[coinciding_vertex])
405  result[coinciding_vertex].insert(owner);
406  }
407  }
408 
409  return result;
410  }
411 
412 
413 
414  template <int dim, int spacedim>
416  MPI_Comm mpi_communicator,
417  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
418  smooth_grid,
419  const bool check_for_distorted_cells)
420  : ::parallel::TriangulationBase<dim, spacedim>(
421  mpi_communicator,
422  smooth_grid,
423  check_for_distorted_cells)
424  {}
425 
426 } // end namespace parallel
427 
428 
429 
430 /*-------------- Explicit Instantiations -------------------------------*/
431 #include "tria_base.inst"
432 
433 DEAL_II_NAMESPACE_CLOSE
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
Definition: tria.cc:10481
unsigned int n_vertices() const
const std::set< types::subdomain_id > & level_ghost_owners() const
Definition: tria_base.cc:351
const std::set< types::subdomain_id > & ghost_owners() const
Definition: tria_base.cc:342
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition: tria_base.cc:65
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12115
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11943
TriangulationBase<dim, spacedim>::fill_level_ghost_owners()
Definition: mpi_tags.h:103
virtual ~TriangulationBase() override
Definition: tria_base.cc:105
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11923
virtual types::global_dof_index n_global_active_cells() const override
Definition: tria_base.cc:135
types::global_dof_index n_global_active_cells
Definition: tria_base.h:240
unsigned int n_levels() const
cell_iterator end() const
Definition: tria.cc:12009
std::set< types::subdomain_id > ghost_owners
Definition: tria_base.h:251
virtual std::size_t memory_consumption() const override
Definition: tria_base.cc:92
static ::ExceptionBase & ExcMessage(std::string arg1)
const bool check_for_distorted_cells
Definition: tria.h:4048
unsigned int subdomain_id
Definition: types.h:43
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1411
virtual std::size_t memory_consumption() const
Definition: tria.cc:15109
types::subdomain_id locally_owned_subdomain() const override
Definition: tria_base.cc:333
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:5663
DistributedTriangulationBase(MPI_Comm mpi_communicator, const typename ::Triangulation< dim, spacedim >::MeshSmoothing smooth_grid=(::Triangulation< dim, spacedim >::none), const bool check_for_distorted_cells=false)
Definition: tria_base.cc:415
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:75
virtual const MPI_Comm & get_communicator() const
Definition: tria_base.cc:171
TriangulationBase(MPI_Comm mpi_communicator, const typename ::Triangulation< dim, spacedim >::MeshSmoothing smooth_grid=(::Triangulation< dim, spacedim >::none), const bool check_for_distorted_cells=false)
Definition: tria_base.cc:43
std::set< types::subdomain_id > level_ghost_owners
Definition: tria_base.h:256
virtual std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors() const
Definition: tria_base.cc:360
unsigned int global_dof_index
Definition: types.h:89
const types::subdomain_id artificial_subdomain_id
Definition: types.h:296
virtual void update_number_cache()
Definition: tria_base.cc:179
Definition: cuda.h:31
unsigned int n_locally_owned_active_cells() const
Definition: tria_base.cc:121
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1699
std::vector< unsigned int > compute_n_locally_owned_active_cells_per_processor() const
Definition: tria_base.cc:143
virtual unsigned int n_global_levels() const override
Definition: tria_base.cc:128
MeshSmoothing smooth_grid
Definition: tria.h:3629
static ::ExceptionBase & ExcNotImplemented()
types::subdomain_id my_subdomain
Definition: tria_base.h:219
T max(const T &t, const MPI_Comm &mpi_communicator)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()