Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+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\}}\)
Loading...
Searching...
No Matches
shared_tria.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/mpi.h>
16#include <deal.II/base/mpi.templates.h>
18
21
24#include <deal.II/grid/tria.h>
27
29
30#include <type_traits>
31
32
34
35#ifdef DEAL_II_WITH_MPI
36namespace parallel
37{
38 namespace shared
39 {
40 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)
51 , settings(settings)
52 , allow_artificial_cells(allow_artificial_cells)
53 {
54 const auto partition_settings =
55 (partition_zoltan | partition_metis | partition_zorder |
56 partition_custom_signal) &
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
67 if (settings & construct_multigrid_hierarchy)
68 Assert(allow_artificial_cells,
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>
77 bool Triangulation<dim, spacedim>::is_multilevel_hierarchy_constructed()
78 const
79 {
80 return (settings & construct_multigrid_hierarchy);
81 }
82
83
84
85 template <int dim, int spacedim>
87 void Triangulation<dim, spacedim>::partition()
88 {
89# ifdef DEBUG
90 // Check that all meshes are the same (or at least have the same
91 // total number of active cells):
92 const unsigned int max_active_cells =
93 Utilities::MPI::max(this->n_active_cells(), this->get_communicator());
94 Assert(
95 max_active_cells == this->n_active_cells(),
97 "A parallel::shared::Triangulation needs to be refined in the same "
98 "way on all processors, but the participating processors don't "
99 "agree on the number of active cells."));
100# endif
101
102 auto partition_settings = (partition_zoltan | partition_metis |
103 partition_zorder | partition_custom_signal) &
104 settings;
105 if (partition_settings == partition_auto)
106# ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
107 partition_settings = partition_zoltan;
108# elif defined DEAL_II_WITH_METIS
109 partition_settings = partition_metis;
110# else
111 partition_settings = partition_zorder;
112# endif
113
114 if (partition_settings == partition_zoltan)
115 {
116# ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
117 AssertThrow(false,
119 "Choosing 'partition_zoltan' requires the library "
120 "to be compiled with support for Zoltan! "
121 "Instead, you might use 'partition_auto' to select "
122 "a partitioning algorithm that is supported "
123 "by your current configuration."));
124# else
126 this->n_subdomains, *this, SparsityTools::Partitioner::zoltan);
127# endif
128 }
129 else if (partition_settings == partition_metis)
130 {
131# ifndef DEAL_II_WITH_METIS
132 AssertThrow(false,
134 "Choosing 'partition_metis' requires the library "
135 "to be compiled with support for METIS! "
136 "Instead, you might use 'partition_auto' to select "
137 "a partitioning algorithm that is supported "
138 "by your current configuration."));
139# else
140 GridTools::partition_triangulation(this->n_subdomains,
141 *this,
143# endif
144 }
145 else if (partition_settings == partition_zorder)
146 {
147 GridTools::partition_triangulation_zorder(this->n_subdomains, *this);
148 }
149 else if (partition_settings == partition_custom_signal)
150 {
151 // User partitions mesh manually
152 }
153 else
154 {
156 }
157
158 // do not partition multigrid levels if user is
159 // defining a custom partition
160 if ((settings & construct_multigrid_hierarchy) &&
161 !(settings & partition_custom_signal))
163
164 true_subdomain_ids_of_cells.resize(this->n_active_cells());
165
166 // loop over all cells and mark artificial:
168 spacedim>::active_cell_iterator
169 cell = this->begin_active(),
170 endc = this->end();
171
172 if (allow_artificial_cells)
173 {
174 // get active halo layer of (ghost) cells
175 // parallel::shared::Triangulation<dim>::
176 std::function<bool(
179 predicate = IteratorFilters::SubdomainEqualTo(this->my_subdomain);
180
181 const std::vector<typename parallel::shared::Triangulation<
182 dim,
183 spacedim>::active_cell_iterator>
184 active_halo_layer_vector =
186 predicate);
187 std::set<typename parallel::shared::Triangulation<dim, spacedim>::
188 active_cell_iterator>
189 active_halo_layer(active_halo_layer_vector.begin(),
190 active_halo_layer_vector.end());
191
192 for (unsigned int index = 0; cell != endc; cell++, index++)
193 {
194 // store original/true subdomain ids:
195 true_subdomain_ids_of_cells[index] = cell->subdomain_id();
196
197 if (cell->is_locally_owned() == false &&
198 active_halo_layer.find(cell) == active_halo_layer.end())
199 cell->set_subdomain_id(numbers::artificial_subdomain_id);
200 }
201
202 // loop over all cells in multigrid hierarchy and mark artificial:
203 if (settings & construct_multigrid_hierarchy)
204 {
205 true_level_subdomain_ids_of_cells.resize(this->n_levels());
206
207 std::function<bool(
209 cell_iterator &)>
211 for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
212 {
213 true_level_subdomain_ids_of_cells[lvl].resize(
214 this->n_cells(lvl));
215
216 const std::vector<typename parallel::shared::Triangulation<
217 dim,
218 spacedim>::cell_iterator>
219 level_halo_layer_vector =
221 *this, predicate, lvl);
222 std::set<typename parallel::shared::
223 Triangulation<dim, spacedim>::cell_iterator>
224 level_halo_layer(level_halo_layer_vector.begin(),
225 level_halo_layer_vector.end());
226
228 cell_iterator cell = this->begin(lvl),
229 endc = this->end(lvl);
230 for (unsigned int index = 0; cell != endc; cell++, index++)
231 {
232 // Store true level subdomain IDs before setting
233 // artificial
234 true_level_subdomain_ids_of_cells[lvl][index] =
235 cell->level_subdomain_id();
236
237 // for active cells, we must have knowledge of level
238 // subdomain ids of all neighbors to our subdomain, not
239 // just neighbors on the same level. if the cells
240 // subdomain id was not set to artitficial above, we will
241 // also keep its level subdomain id since it is either
242 // owned by this processor or in the ghost layer of the
243 // active mesh.
244 if (cell->is_active() &&
245 cell->subdomain_id() !=
247 continue;
248
249 // we must have knowledge of our parent in the hierarchy
250 if (cell->has_children())
251 {
252 bool keep_cell = false;
253 for (unsigned int c = 0;
254 c < GeometryInfo<dim>::max_children_per_cell;
255 ++c)
256 if (cell->child(c)->level_subdomain_id() ==
257 this->my_subdomain)
258 {
259 keep_cell = true;
260 break;
261 }
262 if (keep_cell)
263 continue;
264 }
265
266 // we must have knowledge of our neighbors on the same
267 // level
268 if (!cell->is_locally_owned_on_level() &&
269 level_halo_layer.find(cell) != level_halo_layer.end())
270 continue;
271
272 // mark all other cells to artificial
273 cell->set_level_subdomain_id(
275 }
276 }
277 }
278 }
279 else
280 {
281 // just store true subdomain ids
282 for (unsigned int index = 0; cell != endc; cell++, index++)
283 true_subdomain_ids_of_cells[index] = cell->subdomain_id();
284 }
285
286# ifdef DEBUG
287 {
288 // Assert that each cell is owned by a processor
289 const unsigned int n_my_cells = std::count_if(
290 this->begin_active(),
292 this->end()),
293 [](const auto &i) { return (i.is_locally_owned()); });
294
295 const unsigned int total_cells =
296 Utilities::MPI::sum(n_my_cells, this->get_communicator());
297 Assert(total_cells == this->n_active_cells(),
298 ExcMessage("Not all cells are assigned to a processor."));
299 }
300
301 // If running with multigrid, assert that each level
302 // cell is owned by a processor
303 if (settings & construct_multigrid_hierarchy)
304 {
305 const unsigned int n_my_cells =
306 std::count_if(this->begin(), this->end(), [](const auto &i) {
307 return (i.is_locally_owned_on_level());
308 });
309
310
311 const unsigned int total_cells =
312 Utilities::MPI::sum(n_my_cells, this->get_communicator());
313 Assert(total_cells == this->n_cells(),
314 ExcMessage("Not all cells are assigned to a processor."));
315 }
316# endif
317 }
318
319
320
321 template <int dim, int spacedim>
323 bool Triangulation<dim, spacedim>::with_artificial_cells() const
324 {
325 return allow_artificial_cells;
326 }
327
328
329
330 template <int dim, int spacedim>
332 const std::vector<types::subdomain_id>
334 {
335 return true_subdomain_ids_of_cells;
336 }
337
338
339
340 template <int dim, int spacedim>
342 const std::vector<types::subdomain_id>
344 const unsigned int level) const
345 {
346 Assert(level < true_level_subdomain_ids_of_cells.size(),
348 Assert(true_level_subdomain_ids_of_cells[level].size() ==
349 this->n_cells(level),
351 return true_level_subdomain_ids_of_cells[level];
352 }
353
354
355
356 template <int dim, int spacedim>
358 void Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
359 {
360 // make sure that all refinement/coarsening flags are the same on all
361 // processes
362 {
363 // Obtain the type used to store the different possibilities
364 // a cell can be refined. This is a bit awkward because
365 // what `cell->refine_flag_set()` returns is a struct
366 // type, RefinementCase, which internally stores a
367 // std::uint8_t, which actually holds integers of
368 // enum type RefinementPossibilities<dim>::Possibilities.
369 // In the following, use the actual name of the enum, but
370 // make sure that it is in fact a `std::uint8_t` or
371 // equally sized type.
372 using int_type = std::underlying_type_t<
374 static_assert(sizeof(int_type) == sizeof(std::uint8_t),
375 "Internal type mismatch.");
376
377 std::vector<int_type> refinement_configurations(this->n_active_cells() *
378 2,
379 int_type(0));
380 for (const auto &cell : this->active_cell_iterators())
381 if (cell->is_locally_owned())
382 {
383 refinement_configurations[cell->active_cell_index() * 2 + 0] =
384 static_cast<int_type>(cell->refine_flag_set());
385 refinement_configurations[cell->active_cell_index() * 2 + 1] =
386 static_cast<int_type>(cell->coarsen_flag_set() ? 1 : 0);
387 }
388
389 Utilities::MPI::max(refinement_configurations,
390 this->get_communicator(),
391 refinement_configurations);
392
393 for (const auto &cell : this->active_cell_iterators())
394 {
395 cell->clear_refine_flag();
396 cell->clear_coarsen_flag();
397
398 Assert(
399 (refinement_configurations[cell->active_cell_index() * 2 + 0] >
400 0 ?
401 1 :
402 0) +
403 refinement_configurations[cell->active_cell_index() * 2 +
404 1] <=
405 1,
407 "Refinement/coarsening flags of cells are not consistent in parallel!"));
408
409 if (refinement_configurations[cell->active_cell_index() * 2 + 0] !=
410 0)
411 cell->set_refine_flag(RefinementCase<dim>(
412 refinement_configurations[cell->active_cell_index() * 2 + 0]));
413
414 if (refinement_configurations[cell->active_cell_index() * 2 + 1] >
415 0)
416 cell->set_coarsen_flag();
417 }
418 }
419
421 partition();
422 this->update_number_cache();
423 }
424
425
426
427 template <int dim, int spacedim>
429 void Triangulation<dim, spacedim>::create_triangulation(
430 const std::vector<Point<spacedim>> &vertices,
431 const std::vector<CellData<dim>> &cells,
432 const SubCellData &subcelldata)
433 {
434 try
435 {
437 vertices, cells, subcelldata);
438 }
439 catch (
440 const typename ::Triangulation<dim, spacedim>::DistortedCellList
441 &)
442 {
443 // the underlying triangulation should not be checking for distorted
444 // cells
446 }
447 partition();
448 this->update_number_cache();
449 }
450
451
452
453 template <int dim, int spacedim>
455 void Triangulation<dim, spacedim>::create_triangulation(
456 const TriangulationDescription::Description<dim, spacedim>
457 &construction_data)
458 {
459 (void)construction_data;
460
462 }
463
464
465
466 template <int dim, int spacedim>
468 void Triangulation<dim, spacedim>::copy_triangulation(
469 const ::Triangulation<dim, spacedim> &other_tria)
470 {
471 Assert(
472 (dynamic_cast<
473 const ::parallel::DistributedTriangulationBase<dim, spacedim>
474 *>(&other_tria) == nullptr),
476 "Cannot use this function on parallel::distributed::Triangulation."));
477
479 other_tria);
480 partition();
481 this->update_number_cache();
482 }
483 } // namespace shared
484} // namespace parallel
485
486#else
487
488namespace parallel
489{
490 namespace shared
491 {
492 template <int dim, int spacedim>
495 {
497 return true;
498 }
499
500
501
502 template <int dim, int spacedim>
505 const
506 {
507 return false;
508 }
509
510
511
512 template <int dim, int spacedim>
514 const std::vector<unsigned int>
516 {
518 return true_subdomain_ids_of_cells;
519 }
520
521
522
523 template <int dim, int spacedim>
525 const std::vector<unsigned int>
527 const unsigned int) const
528 {
530 return true_level_subdomain_ids_of_cells;
531 }
532 } // namespace shared
533} // namespace parallel
534
535
536#endif
537
538
539
540namespace internal
541{
542 namespace parallel
543 {
544 namespace shared
545 {
546 template <int dim, int spacedim>
549 : shared_tria(
550 dynamic_cast<
551 const ::parallel::shared::Triangulation<dim, spacedim> *>(
552 &tria))
553 {
554 if (shared_tria && shared_tria->with_artificial_cells())
555 {
556 // Save the current set of subdomain IDs, and set subdomain IDs
557 // to the "true" owner of each cell.
558 const std::vector<types::subdomain_id> &true_subdomain_ids =
559 shared_tria->get_true_subdomain_ids_of_cells();
560
561 saved_subdomain_ids.resize(shared_tria->n_active_cells());
562 for (const auto &cell : shared_tria->active_cell_iterators())
563 {
564 const unsigned int index = cell->active_cell_index();
565 saved_subdomain_ids[index] = cell->subdomain_id();
566 cell->set_subdomain_id(true_subdomain_ids[index]);
567 }
568 }
569 }
570
571
572
573 template <int dim, int spacedim>
576 {
577 if (shared_tria && shared_tria->with_artificial_cells())
578 {
579 // Undo the subdomain modification.
580 for (const auto &cell : shared_tria->active_cell_iterators())
581 {
582 const unsigned int index = cell->active_cell_index();
583 cell->set_subdomain_id(saved_subdomain_ids[index]);
584 }
585 }
586 }
587 } // namespace shared
588 } // namespace parallel
589} // namespace internal
590
591
592/*-------------- Explicit Instantiations -------------------------------*/
593#include "shared_tria.inst"
594
Definition point.h:111
cell_iterator end() const
const SmartPointer< const ::parallel::shared::Triangulation< dim, spacedim > > shared_tria
TemporarilyRestoreSubdomainIds(const Triangulation< dim, spacedim > &tria)
types::subdomain_id my_subdomain
Definition tria_base.h:319
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition tria_base.cc:67
virtual void execute_coarsening_and_refinement() override
typename ::Triangulation< dim, spacedim >::active_cell_iterator active_cell_iterator
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
const std::vector< types::subdomain_id > & get_true_level_subdomain_ids_of_cells(const unsigned int level) const
virtual bool is_multilevel_hierarchy_constructed() const override
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata) override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
unsigned int level
Definition grid_out.cc:4616
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
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)
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)
T sum(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const Iterator & begin
Definition parallel.h:609
STL namespace.