Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
shared_tria.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2018 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/utilities.h>
18 
19 #include <deal.II/distributed/shared_tria.h>
20 #include <deal.II/distributed/tria.h>
21 
22 #include <deal.II/grid/filtered_iterator.h>
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_tools.h>
29 
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 #ifdef DEAL_II_WITH_MPI
34 namespace parallel
35 {
36  namespace shared
37  {
38  template <int dim, int spacedim>
40  MPI_Comm mpi_communicator,
41  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
42  smooth_grid,
43  const bool allow_artificial_cells,
44  const Settings settings)
45  : ::parallel::Triangulation<dim, spacedim>(mpi_communicator,
46  smooth_grid,
47  false)
48  , settings(settings)
49  , allow_artificial_cells(allow_artificial_cells)
50  {
51  const auto partition_settings =
54  settings;
55  (void)partition_settings;
56  Assert(partition_settings == partition_auto ||
57  partition_settings == partition_metis ||
58  partition_settings == partition_zoltan ||
59  partition_settings == partition_zorder ||
60  partition_settings == partition_custom_signal,
61  ExcMessage("Settings must contain exactly one type of the active "
62  "cell partitioning scheme."));
63 
66  ExcMessage("construct_multigrid_hierarchy requires "
67  "allow_artificial_cells to be set to true."));
68  }
69 
70 
71 
72  template <int dim, int spacedim>
73  void
75  {
76 # ifdef DEBUG
77  // Check that all meshes are the same (or at least have the same
78  // total number of active cells):
79  const unsigned int max_active_cells =
80  Utilities::MPI::max(this->n_active_cells(), this->get_communicator());
81  Assert(
82  max_active_cells == this->n_active_cells(),
83  ExcMessage(
84  "A parallel::shared::Triangulation needs to be refined in the same"
85  "way on all processors, but the participating processors don't "
86  "agree on the number of active cells."));
87 # endif
88 
89  auto partition_settings = (partition_zoltan | partition_metis |
90  partition_zorder | partition_custom_signal) &
91  settings;
92  if (partition_settings == partition_auto)
93 # ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
94  partition_settings = partition_zoltan;
95 # elif defined DEAL_II_WITH_METIS
96  partition_settings = partition_metis;
97 # else
98  partition_settings = partition_zorder;
99 # endif
100 
101  if (partition_settings == partition_zoltan)
102  {
103 # ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
104  AssertThrow(false,
105  ExcMessage(
106  "Choosing 'partition_zoltan' requires the library "
107  "to be compiled with support for Zoltan! "
108  "Instead, you might use 'partition_auto' to select "
109  "a partitioning algorithm that is supported "
110  "by your current configuration."));
111 # else
113  this->n_subdomains, *this, SparsityTools::Partitioner::zoltan);
114 # endif
115  }
116  else if (partition_settings == partition_metis)
117  {
118 # ifndef DEAL_II_WITH_METIS
119  AssertThrow(false,
120  ExcMessage(
121  "Choosing 'partition_metis' requires the library "
122  "to be compiled with support for METIS! "
123  "Instead, you might use 'partition_auto' to select "
124  "a partitioning algorithm that is supported "
125  "by your current configuration."));
126 # else
127  GridTools::partition_triangulation(this->n_subdomains,
128  *this,
130 # endif
131  }
132  else if (partition_settings == partition_zorder)
133  {
134  GridTools::partition_triangulation_zorder(this->n_subdomains, *this);
135  }
136  else if (partition_settings == partition_custom_signal)
137  {
138  // User partitions mesh manually
139  }
140  else
141  {
142  AssertThrow(false, ExcInternalError())
143  }
144 
145  // do not partition multigrid levels if user is
146  // defining a custom partition
147  if ((settings & construct_multigrid_hierarchy) &&
148  !(settings & partition_custom_signal))
149  ::GridTools::partition_multigrid_levels(*this);
150 
151  true_subdomain_ids_of_cells.resize(this->n_active_cells());
152 
153  // loop over all cells and mark artificial:
154  typename parallel::shared::Triangulation<dim,
155  spacedim>::active_cell_iterator
156  cell = this->begin_active(),
157  endc = this->end();
158 
159  if (allow_artificial_cells)
160  {
161  // get active halo layer of (ghost) cells
162  // parallel::shared::Triangulation<dim>::
163  std::function<bool(
165  active_cell_iterator &)>
166  predicate = IteratorFilters::SubdomainEqualTo(this->my_subdomain);
167 
168  const std::vector<typename parallel::shared::Triangulation<
169  dim,
170  spacedim>::active_cell_iterator>
171  active_halo_layer_vector =
172  ::GridTools::compute_active_cell_halo_layer(*this,
173  predicate);
174  std::set<typename parallel::shared::Triangulation<dim, spacedim>::
175  active_cell_iterator>
176  active_halo_layer(active_halo_layer_vector.begin(),
177  active_halo_layer_vector.end());
178 
179  for (unsigned int index = 0; cell != endc; cell++, index++)
180  {
181  // store original/true subdomain ids:
182  true_subdomain_ids_of_cells[index] = cell->subdomain_id();
183 
184  if (cell->is_locally_owned() == false &&
185  active_halo_layer.find(cell) == active_halo_layer.end())
186  cell->set_subdomain_id(numbers::artificial_subdomain_id);
187  }
188 
189  // loop over all cells in multigrid hierarchy and mark artificial:
190  if (settings & construct_multigrid_hierarchy)
191  {
192  true_level_subdomain_ids_of_cells.resize(this->n_levels());
193 
194  std::function<bool(
196  cell_iterator &)>
198  for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
199  {
200  true_level_subdomain_ids_of_cells[lvl].resize(
201  this->n_cells(lvl));
202 
203  const std::vector<typename parallel::shared::Triangulation<
204  dim,
205  spacedim>::cell_iterator>
206  level_halo_layer_vector =
207  ::GridTools::compute_cell_halo_layer_on_level(
208  *this, predicate, lvl);
209  std::set<typename parallel::shared::
210  Triangulation<dim, spacedim>::cell_iterator>
211  level_halo_layer(level_halo_layer_vector.begin(),
212  level_halo_layer_vector.end());
213 
215  cell_iterator cell = this->begin(lvl),
216  endc = this->end(lvl);
217  for (unsigned int index = 0; cell != endc; cell++, index++)
218  {
219  // Store true level subdomain IDs before setting
220  // artificial
221  true_level_subdomain_ids_of_cells[lvl][index] =
222  cell->level_subdomain_id();
223 
224  // for active cells, we must have knowledge of level
225  // subdomain ids of all neighbors to our subdomain, not
226  // just neighbors on the same level. if the cells
227  // subdomain id was not set to artitficial above, we will
228  // also keep its level subdomain id since it is either
229  // owned by this processor or in the ghost layer of the
230  // active mesh.
231  if (!cell->has_children() &&
232  cell->subdomain_id() !=
234  continue;
235 
236  // we must have knowledge of our parent in the hierarchy
237  if (cell->has_children())
238  {
239  bool keep_cell = false;
240  for (unsigned int c = 0;
241  c < GeometryInfo<dim>::max_children_per_cell;
242  ++c)
243  if (cell->child(c)->level_subdomain_id() ==
244  this->my_subdomain)
245  {
246  keep_cell = true;
247  break;
248  }
249  if (keep_cell)
250  continue;
251  }
252 
253  // we must have knowledge of our neighbors on the same
254  // level
255  if (!cell->is_locally_owned_on_level() &&
256  level_halo_layer.find(cell) != level_halo_layer.end())
257  continue;
258 
259  // mark all other cells to artificial
260  cell->set_level_subdomain_id(
262  }
263  }
264  }
265  }
266  else
267  {
268  // just store true subdomain ids
269  for (unsigned int index = 0; cell != endc; cell++, index++)
270  true_subdomain_ids_of_cells[index] = cell->subdomain_id();
271  }
272 
273 # ifdef DEBUG
274  {
275  // Assert that each cell is owned by a processor
276  unsigned int n_my_cells = 0;
277  typename parallel::shared::Triangulation<dim,
278  spacedim>::active_cell_iterator
279  cell = this->begin_active(),
280  endc = this->end();
281  for (; cell != endc; ++cell)
282  if (cell->is_locally_owned())
283  n_my_cells += 1;
284 
285  const unsigned int total_cells =
286  Utilities::MPI::sum(n_my_cells, this->get_communicator());
287  Assert(total_cells == this->n_active_cells(),
288  ExcMessage("Not all cells are assigned to a processor."));
289  }
290 
291  // If running with multigrid, assert that each level
292  // cell is owned by a processor
293  if (settings & construct_multigrid_hierarchy)
294  {
295  unsigned int n_my_cells = 0;
296  typename parallel::shared::Triangulation<dim, spacedim>::cell_iterator
297  cell = this->begin(),
298  endc = this->end();
299  for (; cell != endc; ++cell)
300  if (cell->is_locally_owned_on_level())
301  n_my_cells += 1;
302 
303  const unsigned int total_cells =
304  Utilities::MPI::sum(n_my_cells, this->get_communicator());
305  Assert(total_cells == this->n_cells(),
306  ExcMessage("Not all cells are assigned to a processor."));
307  }
308 # endif
309  }
310 
311 
312 
313  template <int dim, int spacedim>
314  bool
316  {
317  return allow_artificial_cells;
318  }
319 
320 
321 
322  template <int dim, int spacedim>
323  const std::vector<types::subdomain_id> &
325  {
326  return true_subdomain_ids_of_cells;
327  }
328 
329 
330 
331  template <int dim, int spacedim>
332  const std::vector<types::subdomain_id> &
334  const unsigned int level) const
335  {
336  Assert(level < true_level_subdomain_ids_of_cells.size(),
337  ExcInternalError());
338  Assert(true_level_subdomain_ids_of_cells[level].size() ==
339  this->n_cells(level),
340  ExcInternalError());
341  return true_level_subdomain_ids_of_cells[level];
342  }
343 
344 
345 
346  template <int dim, int spacedim>
347  void
349  {
351  partition();
352  this->update_number_cache();
353  }
354 
355 
356 
357  template <int dim, int spacedim>
358  void
360  const std::vector<Point<spacedim>> &vertices,
361  const std::vector<CellData<dim>> & cells,
362  const SubCellData & subcelldata)
363  {
364  try
365  {
367  vertices, cells, subcelldata);
368  }
369  catch (
370  const typename ::Triangulation<dim, spacedim>::DistortedCellList
371  &)
372  {
373  // the underlying triangulation should not be checking for distorted
374  // cells
375  Assert(false, ExcInternalError());
376  }
377  partition();
378  this->update_number_cache();
379  }
380 
381 
382 
383  template <int dim, int spacedim>
384  void
386  const ::Triangulation<dim, spacedim> &other_tria)
387  {
388  Assert(
389  (dynamic_cast<
390  const ::parallel::distributed::Triangulation<dim, spacedim> *>(
391  &other_tria) == nullptr),
392  ExcMessage(
393  "Cannot use this function on parallel::distributed::Triangulation."));
394 
396  other_tria);
397  partition();
398  this->update_number_cache();
399  }
400 
401 
402 
403  template <int dim, int spacedim>
404  void
406  {
408 
409  if (settings & construct_multigrid_hierarchy)
411  }
412  } // namespace shared
413 } // namespace parallel
414 
415 #else
416 
417 namespace parallel
418 {
419  namespace shared
420  {
421  template <int dim, int spacedim>
422  bool
424  {
425  Assert(false, ExcNotImplemented());
426  return true;
427  }
428 
429 
430 
431  template <int dim, int spacedim>
432  const std::vector<unsigned int> &
434  {
435  Assert(false, ExcNotImplemented());
436  return true_subdomain_ids_of_cells;
437  }
438 
439 
440 
441  template <int dim, int spacedim>
442  const std::vector<unsigned int> &
444  const unsigned int) const
445  {
446  Assert(false, ExcNotImplemented());
447  return true_level_subdomain_ids_of_cells;
448  }
449  } // namespace shared
450 } // namespace parallel
451 
452 
453 #endif
454 
455 
456 /*-------------- Explicit Instantiations -------------------------------*/
457 #include "shared_tria.inst"
458 
459 DEAL_II_NAMESPACE_CLOSE
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
Definition: shared_tria.cc:324
Triangulation(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:39
virtual void execute_coarsening_and_refinement() override
Definition: shared_tria.cc:348
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
const std::vector< types::subdomain_id > & get_true_level_subdomain_ids_of_cells(const unsigned int level) const
Definition: shared_tria.cc:333
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2906
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:2679
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:13267
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition: tria_base.cc:66
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void update_number_cache() override
Definition: shared_tria.cc:405
T sum(const T &t, const MPI_Comm &mpi_communicator)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:10504
#define Assert(cond, exc)
Definition: exceptions.h:1407
virtual void update_number_cache()
Definition: tria_base.cc:163
const types::subdomain_id artificial_subdomain_id
Definition: types.h:275
__global__ void set(Number *val, const Number s, const size_type N)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata) override
Definition: shared_tria.cc:359
static ::ExceptionBase & ExcNotImplemented()
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &other_tria) override
Definition: shared_tria.cc:385
T max(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcInternalError()