Reference documentation for deal.II version GIT d7aca55de5 2022-08-10 12:50: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\}}\)
shared_tria.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 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/mpi.h>
17 #include <deal.II/base/mpi.templates.h>
18 #include <deal.II/base/utilities.h>
19 
22 
25 #include <deal.II/grid/tria.h>
28 
30 
31 #include <type_traits>
32 
33 
35 
36 #ifdef DEAL_II_WITH_MPI
37 namespace parallel
38 {
39  namespace shared
40  {
41  template <int dim, int spacedim>
43  const MPI_Comm &mpi_communicator,
44  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
45  smooth_grid,
46  const bool allow_artificial_cells,
47  const Settings settings)
48  : ::parallel::TriangulationBase<dim, spacedim>(mpi_communicator,
49  smooth_grid,
50  false)
52  , allow_artificial_cells(allow_artificial_cells)
53  {
54  const auto partition_settings =
57  settings;
58  (void)partition_settings;
59  Assert(partition_settings == partition_auto ||
60  partition_settings == partition_metis ||
61  partition_settings == partition_zoltan ||
62  partition_settings == partition_zorder ||
63  partition_settings == partition_custom_signal,
64  ExcMessage("Settings must contain exactly one type of the active "
65  "cell partitioning scheme."));
66 
69  ExcMessage("construct_multigrid_hierarchy requires "
70  "allow_artificial_cells to be set to true."));
71  }
72 
73 
74 
75  template <int dim, int spacedim>
76  bool
78  {
80  }
81 
82 
83 
84  template <int dim, int spacedim>
85  void
87  {
88 # ifdef DEBUG
89  // Check that all meshes are the same (or at least have the same
90  // total number of active cells):
91  const unsigned int max_active_cells =
92  Utilities::MPI::max(this->n_active_cells(), this->get_communicator());
93  Assert(
94  max_active_cells == this->n_active_cells(),
95  ExcMessage(
96  "A parallel::shared::Triangulation needs to be refined in the same "
97  "way on all processors, but the participating processors don't "
98  "agree on the number of active cells."));
99 # endif
100 
101  auto partition_settings = (partition_zoltan | partition_metis |
102  partition_zorder | partition_custom_signal) &
103  settings;
104  if (partition_settings == partition_auto)
105 # ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
106  partition_settings = partition_zoltan;
107 # elif defined DEAL_II_WITH_METIS
108  partition_settings = partition_metis;
109 # else
110  partition_settings = partition_zorder;
111 # endif
112 
113  if (partition_settings == partition_zoltan)
114  {
115 # ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
116  AssertThrow(false,
117  ExcMessage(
118  "Choosing 'partition_zoltan' requires the library "
119  "to be compiled with support for Zoltan! "
120  "Instead, you might use 'partition_auto' to select "
121  "a partitioning algorithm that is supported "
122  "by your current configuration."));
123 # else
125  this->n_subdomains, *this, SparsityTools::Partitioner::zoltan);
126 # endif
127  }
128  else if (partition_settings == partition_metis)
129  {
130 # ifndef DEAL_II_WITH_METIS
131  AssertThrow(false,
132  ExcMessage(
133  "Choosing 'partition_metis' requires the library "
134  "to be compiled with support for METIS! "
135  "Instead, you might use 'partition_auto' to select "
136  "a partitioning algorithm that is supported "
137  "by your current configuration."));
138 # else
139  GridTools::partition_triangulation(this->n_subdomains,
140  *this,
142 # endif
143  }
144  else if (partition_settings == partition_zorder)
145  {
146  GridTools::partition_triangulation_zorder(this->n_subdomains, *this);
147  }
148  else if (partition_settings == partition_custom_signal)
149  {
150  // User partitions mesh manually
151  }
152  else
153  {
154  AssertThrow(false, ExcInternalError())
155  }
156 
157  // do not partition multigrid levels if user is
158  // defining a custom partition
160  !(settings & partition_custom_signal))
162 
163  true_subdomain_ids_of_cells.resize(this->n_active_cells());
164 
165  // loop over all cells and mark artificial:
166  typename parallel::shared::Triangulation<dim,
167  spacedim>::active_cell_iterator
168  cell = this->begin_active(),
169  endc = this->end();
170 
171  if (allow_artificial_cells)
172  {
173  // get active halo layer of (ghost) cells
174  // parallel::shared::Triangulation<dim>::
175  std::function<bool(
178  predicate = IteratorFilters::SubdomainEqualTo(this->my_subdomain);
179 
180  const std::vector<typename parallel::shared::Triangulation<
181  dim,
182  spacedim>::active_cell_iterator>
183  active_halo_layer_vector =
185  predicate);
186  std::set<typename parallel::shared::Triangulation<dim, spacedim>::
187  active_cell_iterator>
188  active_halo_layer(active_halo_layer_vector.begin(),
189  active_halo_layer_vector.end());
190 
191  for (unsigned int index = 0; cell != endc; cell++, index++)
192  {
193  // store original/true subdomain ids:
194  true_subdomain_ids_of_cells[index] = cell->subdomain_id();
195 
196  if (cell->is_locally_owned() == false &&
197  active_halo_layer.find(cell) == active_halo_layer.end())
198  cell->set_subdomain_id(numbers::artificial_subdomain_id);
199  }
200 
201  // loop over all cells in multigrid hierarchy and mark artificial:
203  {
204  true_level_subdomain_ids_of_cells.resize(this->n_levels());
205 
206  std::function<bool(
208  cell_iterator &)>
210  for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
211  {
212  true_level_subdomain_ids_of_cells[lvl].resize(
213  this->n_cells(lvl));
214 
215  const std::vector<typename parallel::shared::Triangulation<
216  dim,
217  spacedim>::cell_iterator>
218  level_halo_layer_vector =
220  *this, predicate, lvl);
223  level_halo_layer(level_halo_layer_vector.begin(),
224  level_halo_layer_vector.end());
225 
227  cell_iterator cell = this->begin(lvl),
228  endc = this->end(lvl);
229  for (unsigned int index = 0; cell != endc; cell++, index++)
230  {
231  // Store true level subdomain IDs before setting
232  // artificial
233  true_level_subdomain_ids_of_cells[lvl][index] =
234  cell->level_subdomain_id();
235 
236  // for active cells, we must have knowledge of level
237  // subdomain ids of all neighbors to our subdomain, not
238  // just neighbors on the same level. if the cells
239  // subdomain id was not set to artitficial above, we will
240  // also keep its level subdomain id since it is either
241  // owned by this processor or in the ghost layer of the
242  // active mesh.
243  if (cell->is_active() &&
244  cell->subdomain_id() !=
246  continue;
247 
248  // we must have knowledge of our parent in the hierarchy
249  if (cell->has_children())
250  {
251  bool keep_cell = false;
252  for (unsigned int c = 0;
253  c < GeometryInfo<dim>::max_children_per_cell;
254  ++c)
255  if (cell->child(c)->level_subdomain_id() ==
256  this->my_subdomain)
257  {
258  keep_cell = true;
259  break;
260  }
261  if (keep_cell)
262  continue;
263  }
264 
265  // we must have knowledge of our neighbors on the same
266  // level
267  if (!cell->is_locally_owned_on_level() &&
268  level_halo_layer.find(cell) != level_halo_layer.end())
269  continue;
270 
271  // mark all other cells to artificial
272  cell->set_level_subdomain_id(
274  }
275  }
276  }
277  }
278  else
279  {
280  // just store true subdomain ids
281  for (unsigned int index = 0; cell != endc; cell++, index++)
282  true_subdomain_ids_of_cells[index] = cell->subdomain_id();
283  }
284 
285 # ifdef DEBUG
286  {
287  // Assert that each cell is owned by a processor
288  const unsigned int n_my_cells = std::count_if(
289  this->begin_active(),
291  this->end()),
292  [](const auto &i) { return (i.is_locally_owned()); });
293 
294  const unsigned int total_cells =
295  Utilities::MPI::sum(n_my_cells, this->get_communicator());
296  Assert(total_cells == this->n_active_cells(),
297  ExcMessage("Not all cells are assigned to a processor."));
298  }
299 
300  // If running with multigrid, assert that each level
301  // cell is owned by a processor
302  if (settings & construct_multigrid_hierarchy)
303  {
304  const unsigned int n_my_cells =
305  std::count_if(this->begin(), this->end(), [](const auto &i) {
306  return (i.is_locally_owned_on_level());
307  });
308 
309 
310  const unsigned int total_cells =
311  Utilities::MPI::sum(n_my_cells, this->get_communicator());
312  Assert(total_cells == this->n_cells(),
313  ExcMessage("Not all cells are assigned to a processor."));
314  }
315 # endif
316  }
317 
318 
319 
320  template <int dim, int spacedim>
321  bool
323  {
324  return allow_artificial_cells;
325  }
326 
327 
328 
329  template <int dim, int spacedim>
330  const std::vector<types::subdomain_id> &
332  {
333  return true_subdomain_ids_of_cells;
334  }
335 
336 
337 
338  template <int dim, int spacedim>
339  const std::vector<types::subdomain_id> &
341  const unsigned int level) const
342  {
343  Assert(level < true_level_subdomain_ids_of_cells.size(),
344  ExcInternalError());
345  Assert(true_level_subdomain_ids_of_cells[level].size() ==
346  this->n_cells(level),
347  ExcInternalError());
348  return true_level_subdomain_ids_of_cells[level];
349  }
350 
351 
352 
353  template <int dim, int spacedim>
354  void
356  {
357  // make sure that all refinement/coarsening flags are the same on all
358  // processes
359  {
360  // Obtain the type used to store the different possibilities
361  // a cell can be refined. This is a bit awkward because
362  // what `cell->refine_flag_set()` returns is a struct
363  // type, RefinementCase, which internally stores a
364  // std::uint8_t, which actually holds integers of
365  // enum type RefinementPossibilities<dim>::Possibilities.
366  // In the following, use the actual name of the enum, but
367  // make sure that it is in fact a `std::uint8_t` or
368  // equally sized type.
369  using int_type = std::underlying_type_t<
371  static_assert(sizeof(int_type) == sizeof(std::uint8_t),
372  "Internal type mismatch.");
373 
374  std::vector<int_type> refinement_configurations(this->n_active_cells() *
375  2,
376  int_type(0));
377  for (const auto &cell : this->active_cell_iterators())
378  if (cell->is_locally_owned())
379  {
380  refinement_configurations[cell->active_cell_index() * 2 + 0] =
381  static_cast<int_type>(cell->refine_flag_set());
382  refinement_configurations[cell->active_cell_index() * 2 + 1] =
383  static_cast<int_type>(cell->coarsen_flag_set() ? 1 : 0);
384  }
385 
386  Utilities::MPI::max(refinement_configurations,
387  this->get_communicator(),
388  refinement_configurations);
389 
390  for (const auto &cell : this->active_cell_iterators())
391  {
392  cell->clear_refine_flag();
393  cell->clear_coarsen_flag();
394 
395  Assert(
396  (refinement_configurations[cell->active_cell_index() * 2 + 0] >
397  0 ?
398  1 :
399  0) +
400  refinement_configurations[cell->active_cell_index() * 2 +
401  1] <=
402  1,
403  ExcMessage(
404  "Refinement/coarsening flags of cells are not consistent in parallel!"));
405 
406  if (refinement_configurations[cell->active_cell_index() * 2 + 0] !=
407  0)
408  cell->set_refine_flag(RefinementCase<dim>(
409  refinement_configurations[cell->active_cell_index() * 2 + 0]));
410 
411  if (refinement_configurations[cell->active_cell_index() * 2 + 1] >
412  0)
413  cell->set_coarsen_flag();
414  }
415  }
416 
418  partition();
419  this->update_number_cache();
420  }
421 
422 
423 
424  template <int dim, int spacedim>
425  void
427  const std::vector<Point<spacedim>> &vertices,
428  const std::vector<CellData<dim>> & cells,
429  const SubCellData & subcelldata)
430  {
431  try
432  {
434  vertices, cells, subcelldata);
435  }
436  catch (
437  const typename ::Triangulation<dim, spacedim>::DistortedCellList
438  &)
439  {
440  // the underlying triangulation should not be checking for distorted
441  // cells
442  Assert(false, ExcInternalError());
443  }
444  partition();
445  this->update_number_cache();
446  }
447 
448 
449 
450  template <int dim, int spacedim>
451  void
454  &construction_data)
455  {
456  (void)construction_data;
457 
458  Assert(false, ExcInternalError());
459  }
460 
461 
462 
463  template <int dim, int spacedim>
464  void
466  const ::Triangulation<dim, spacedim> &other_tria)
467  {
468  Assert(
469  (dynamic_cast<
470  const ::parallel::DistributedTriangulationBase<dim, spacedim>
471  *>(&other_tria) == nullptr),
472  ExcMessage(
473  "Cannot use this function on parallel::distributed::Triangulation."));
474 
476  other_tria);
477  partition();
478  this->update_number_cache();
479  }
480  } // namespace shared
481 } // namespace parallel
482 
483 #else
484 
485 namespace parallel
486 {
487  namespace shared
488  {
489  template <int dim, int spacedim>
490  bool
492  {
493  Assert(false, ExcNotImplemented());
494  return true;
495  }
496 
497  template <int dim, int spacedim>
498  bool
500  {
501  return false;
502  }
503 
504  template <int dim, int spacedim>
505  const std::vector<unsigned int> &
507  {
508  Assert(false, ExcNotImplemented());
509  return true_subdomain_ids_of_cells;
510  }
511 
512 
513 
514  template <int dim, int spacedim>
515  const std::vector<unsigned int> &
517  const unsigned int) const
518  {
519  Assert(false, ExcNotImplemented());
520  return true_level_subdomain_ids_of_cells;
521  }
522  } // namespace shared
523 } // namespace parallel
524 
525 
526 #endif
527 
528 
529 
530 namespace internal
531 {
532  namespace parallel
533  {
534  namespace shared
535  {
536  template <int dim, int spacedim>
539  : shared_tria(
540  dynamic_cast<
541  const ::parallel::shared::Triangulation<dim, spacedim> *>(
542  &tria))
543  {
544  if (shared_tria && shared_tria->with_artificial_cells())
545  {
546  // Save the current set of subdomain IDs, and set subdomain IDs
547  // to the "true" owner of each cell.
548  const std::vector<types::subdomain_id> &true_subdomain_ids =
549  shared_tria->get_true_subdomain_ids_of_cells();
550 
551  saved_subdomain_ids.resize(shared_tria->n_active_cells());
552  for (const auto &cell : shared_tria->active_cell_iterators())
553  {
554  const unsigned int index = cell->active_cell_index();
555  saved_subdomain_ids[index] = cell->subdomain_id();
556  cell->set_subdomain_id(true_subdomain_ids[index]);
557  }
558  }
559  }
560 
561 
562 
563  template <int dim, int spacedim>
566  {
567  if (shared_tria && shared_tria->with_artificial_cells())
568  {
569  // Undo the subdomain modification.
570  for (const auto &cell : shared_tria->active_cell_iterators())
571  {
572  const unsigned int index = cell->active_cell_index();
573  cell->set_subdomain_id(saved_subdomain_ids[index]);
574  }
575  }
576  }
577  } // namespace shared
578  } // namespace parallel
579 } // namespace internal
580 
581 
582 /*-------------- Explicit Instantiations -------------------------------*/
583 #include "shared_tria.inst"
584 
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
const SmartPointer< const ::parallel::shared::Triangulation< dim, spacedim > > shared_tria
Definition: shared_tria.h:542
TemporarilyRestoreSubdomainIds(const Triangulation< dim, spacedim > &tria)
Definition: shared_tria.cc:538
types::subdomain_id my_subdomain
Definition: tria_base.h:310
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition: tria_base.cc:68
virtual void execute_coarsening_and_refinement() override
Definition: shared_tria.cc:355
typename ::Triangulation< dim, spacedim >::active_cell_iterator active_cell_iterator
Definition: shared_tria.h:116
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &other_tria) override
Definition: shared_tria.cc:465
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
Definition: shared_tria.cc:331
Triangulation(const MPI_Comm &mpi_communicator, const typename ::Triangulation< dim, spacedim >::MeshSmoothing=(::Triangulation< dim, spacedim >::none), const bool allow_artificial_cells=false, const Settings settings=partition_auto)
Definition: shared_tria.cc:42
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition: shared_tria.h:118
const std::vector< types::subdomain_id > & get_true_level_subdomain_ids_of_cells(const unsigned int level) const
Definition: shared_tria.cc:340
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata) override
Definition: shared_tria.cc:426
virtual bool is_multilevel_hierarchy_constructed() const override
Definition: shared_tria.cc:77
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4607
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4115
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4221
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:3874
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
T max(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13793
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13786
const types::subdomain_id artificial_subdomain_id
Definition: types.h:298
const TriangulationDescription::Settings settings
const ::Triangulation< dim, spacedim > & tria