Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3087-ga35b476a78 2025-04-19 20:30: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 if constexpr (running_in_debug_mode())
90 {
91 // Check that all meshes are the same (or at least have the same
92 // total number of active cells):
93 const unsigned int max_active_cells =
94 Utilities::MPI::max(this->n_active_cells(),
95 this->get_mpi_communicator());
96 Assert(
97 max_active_cells == this->n_active_cells(),
99 "A parallel::shared::Triangulation needs to be refined in the same "
100 "way on all processors, but the participating processors don't "
101 "agree on the number of active cells."));
102 }
103
104 auto partition_settings = (partition_zoltan | partition_metis |
105 partition_zorder | partition_custom_signal) &
106 settings;
107 if (partition_settings == partition_auto)
108# ifdef DEAL_II_TRILINOS_WITH_ZOLTAN
109 partition_settings = partition_zoltan;
110# elif defined DEAL_II_WITH_METIS
111 partition_settings = partition_metis;
112# else
113 partition_settings = partition_zorder;
114# endif
115
116 if (partition_settings == partition_zoltan)
117 {
118# ifndef DEAL_II_TRILINOS_WITH_ZOLTAN
119 AssertThrow(false,
121 "Choosing 'partition_zoltan' requires the library "
122 "to be compiled with support for Zoltan! "
123 "Instead, you might use 'partition_auto' to select "
124 "a partitioning algorithm that is supported "
125 "by your current configuration."));
126# else
128 this->n_subdomains, *this, SparsityTools::Partitioner::zoltan);
129# endif
130 }
131 else if (partition_settings == partition_metis)
132 {
133# ifndef DEAL_II_WITH_METIS
134 AssertThrow(false,
136 "Choosing 'partition_metis' requires the library "
137 "to be compiled with support for METIS! "
138 "Instead, you might use 'partition_auto' to select "
139 "a partitioning algorithm that is supported "
140 "by your current configuration."));
141# else
142 GridTools::partition_triangulation(this->n_subdomains,
143 *this,
145# endif
146 }
147 else if (partition_settings == partition_zorder)
148 {
149 GridTools::partition_triangulation_zorder(this->n_subdomains, *this);
150 }
151 else if (partition_settings == partition_custom_signal)
152 {
153 // User partitions mesh manually
154 }
155 else
156 {
158 }
159
160 // do not partition multigrid levels if user is
161 // defining a custom partition
162 if ((settings & construct_multigrid_hierarchy) &&
163 !(settings & partition_custom_signal))
165
166 true_subdomain_ids_of_cells.resize(this->n_active_cells());
167
168 // loop over all cells and mark artificial:
170 spacedim>::active_cell_iterator
171 cell = this->begin_active(),
172 endc = this->end();
173
174 if (allow_artificial_cells)
175 {
176 // get active halo layer of (ghost) cells
177 // parallel::shared::Triangulation<dim>::
178 std::function<bool(
181 predicate = IteratorFilters::SubdomainEqualTo(this->my_subdomain);
182
183 const std::vector<typename parallel::shared::Triangulation<
184 dim,
185 spacedim>::active_cell_iterator>
186 active_halo_layer_vector =
188 predicate);
189 std::set<typename parallel::shared::Triangulation<dim, spacedim>::
190 active_cell_iterator>
191 active_halo_layer(active_halo_layer_vector.begin(),
192 active_halo_layer_vector.end());
193
194 for (unsigned int index = 0; cell != endc; cell++, index++)
195 {
196 // store original/true subdomain ids:
197 true_subdomain_ids_of_cells[index] = cell->subdomain_id();
198
199 if (cell->is_locally_owned() == false &&
200 active_halo_layer.find(cell) == active_halo_layer.end())
201 cell->set_subdomain_id(numbers::artificial_subdomain_id);
202 }
203
204 // loop over all cells in multigrid hierarchy and mark artificial:
205 if (settings & construct_multigrid_hierarchy)
206 {
207 true_level_subdomain_ids_of_cells.resize(this->n_levels());
208
209 std::function<bool(
211 cell_iterator &)>
213 for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
214 {
215 true_level_subdomain_ids_of_cells[lvl].resize(
216 this->n_cells(lvl));
217
218 const std::vector<typename parallel::shared::Triangulation<
219 dim,
220 spacedim>::cell_iterator>
221 level_halo_layer_vector =
223 *this, predicate, lvl);
224 std::set<typename parallel::shared::
225 Triangulation<dim, spacedim>::cell_iterator>
226 level_halo_layer(level_halo_layer_vector.begin(),
227 level_halo_layer_vector.end());
228
230 cell_iterator cell = this->begin(lvl),
231 endc = this->end(lvl);
232 for (unsigned int index = 0; cell != endc; cell++, index++)
233 {
234 // Store true level subdomain IDs before setting
235 // artificial
236 true_level_subdomain_ids_of_cells[lvl][index] =
237 cell->level_subdomain_id();
238
239 // for active cells, we must have knowledge of level
240 // subdomain ids of all neighbors to our subdomain, not
241 // just neighbors on the same level. if the cells
242 // subdomain id was not set to artitficial above, we will
243 // also keep its level subdomain id since it is either
244 // owned by this processor or in the ghost layer of the
245 // active mesh.
246 if (cell->is_active() &&
247 cell->subdomain_id() !=
249 continue;
250
251 // we must have knowledge of our parent in the hierarchy
252 if (cell->has_children())
253 {
254 bool keep_cell = false;
255 for (unsigned int c = 0;
256 c < GeometryInfo<dim>::max_children_per_cell;
257 ++c)
258 if (cell->child(c)->level_subdomain_id() ==
259 this->my_subdomain)
260 {
261 keep_cell = true;
262 break;
263 }
264 if (keep_cell)
265 continue;
266 }
267
268 // we must have knowledge of our neighbors on the same
269 // level
270 if (!cell->is_locally_owned_on_level() &&
271 level_halo_layer.find(cell) != level_halo_layer.end())
272 continue;
273
274 // mark all other cells to artificial
275 cell->set_level_subdomain_id(
277 }
278 }
279 }
280 }
281 else
282 {
283 // just store true subdomain ids
284 for (unsigned int index = 0; cell != endc; cell++, index++)
285 true_subdomain_ids_of_cells[index] = cell->subdomain_id();
286 }
287
288 if constexpr (running_in_debug_mode())
289 {
290 {
291 // Assert that each cell is owned by a processor
292 const unsigned int n_my_cells = std::count_if(
293 this->begin_active(),
295 this->end()),
296 [](const auto &i) { return (i.is_locally_owned()); });
297
298 const unsigned int total_cells =
299 Utilities::MPI::sum(n_my_cells, this->get_mpi_communicator());
300 Assert(total_cells == this->n_active_cells(),
301 ExcMessage("Not all cells are assigned to a processor."));
302 }
303
304 // If running with multigrid, assert that each level
305 // cell is owned by a processor
306 if (settings & construct_multigrid_hierarchy)
307 {
308 const unsigned int n_my_cells =
309 std::count_if(this->begin(), this->end(), [](const auto &i) {
310 return (i.is_locally_owned_on_level());
311 });
312
313
314 const unsigned int total_cells =
315 Utilities::MPI::sum(n_my_cells, this->get_mpi_communicator());
316 Assert(total_cells == this->n_cells(),
317 ExcMessage("Not all cells are assigned to a processor."));
318 }
319 }
320 }
321
322
323
324 template <int dim, int spacedim>
326 bool Triangulation<dim, spacedim>::with_artificial_cells() const
327 {
328 return allow_artificial_cells;
329 }
330
331
332
333 template <int dim, int spacedim>
335 const std::vector<types::subdomain_id>
337 {
338 return true_subdomain_ids_of_cells;
339 }
340
341
342
343 template <int dim, int spacedim>
345 const std::vector<types::subdomain_id>
347 const unsigned int level) const
348 {
349 Assert(level < true_level_subdomain_ids_of_cells.size(),
351 Assert(true_level_subdomain_ids_of_cells[level].size() ==
352 this->n_cells(level),
354 return true_level_subdomain_ids_of_cells[level];
355 }
356
357
358
359 template <int dim, int spacedim>
361 void Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
362 {
363 // make sure that all refinement/coarsening flags are the same on all
364 // processes
365 {
366 // Obtain the type used to store the different possibilities
367 // a cell can be refined. This is a bit awkward because
368 // what `cell->refine_flag_set()` returns is a struct
369 // type, RefinementCase, which internally stores a
370 // std::uint8_t, which actually holds integers of
371 // enum type RefinementPossibilities<dim>::Possibilities.
372 // In the following, use the actual name of the enum, but
373 // make sure that it is in fact a `std::uint8_t` or
374 // equally sized type.
375 using int_type = std::underlying_type_t<
377 static_assert(sizeof(int_type) == sizeof(std::uint8_t),
378 "Internal type mismatch.");
379
380 std::vector<int_type> refinement_configurations(this->n_active_cells() *
381 2,
382 int_type(0));
383 for (const auto &cell : this->active_cell_iterators())
384 if (cell->is_locally_owned())
385 {
386 refinement_configurations[cell->active_cell_index() * 2 + 0] =
387 static_cast<int_type>(cell->refine_flag_set());
388 refinement_configurations[cell->active_cell_index() * 2 + 1] =
389 static_cast<int_type>(cell->coarsen_flag_set() ? 1 : 0);
390 }
391
392 Utilities::MPI::max(refinement_configurations,
393 this->get_mpi_communicator(),
394 refinement_configurations);
395
396 for (const auto &cell : this->active_cell_iterators())
397 {
398 cell->clear_refine_flag();
399 cell->clear_coarsen_flag();
400
401 Assert(
402 (refinement_configurations[cell->active_cell_index() * 2 + 0] >
403 0 ?
404 1 :
405 0) +
406 refinement_configurations[cell->active_cell_index() * 2 +
407 1] <=
408 1,
410 "Refinement/coarsening flags of cells are not consistent in parallel!"));
411
412 if (refinement_configurations[cell->active_cell_index() * 2 + 0] !=
413 0)
414 cell->set_refine_flag(RefinementCase<dim>(
415 refinement_configurations[cell->active_cell_index() * 2 + 0]));
416
417 if (refinement_configurations[cell->active_cell_index() * 2 + 1] >
418 0)
419 cell->set_coarsen_flag();
420 }
421 }
422
424 partition();
425 this->update_number_cache();
426 }
427
428
429
430 template <int dim, int spacedim>
432 void Triangulation<dim, spacedim>::create_triangulation(
433 const std::vector<Point<spacedim>> &vertices,
434 const std::vector<CellData<dim>> &cells,
435 const SubCellData &subcelldata)
436 {
437 try
438 {
440 vertices, cells, subcelldata);
441 }
442 catch (
443 const typename ::Triangulation<dim, spacedim>::DistortedCellList
444 &)
445 {
446 // the underlying triangulation should not be checking for distorted
447 // cells
449 }
450 partition();
451 this->update_number_cache();
452 }
453
454
455
456 template <int dim, int spacedim>
458 void Triangulation<dim, spacedim>::create_triangulation(
459 const TriangulationDescription::Description<dim, spacedim>
460 &construction_data)
461 {
462 (void)construction_data;
463
465 }
466
467
468
469 template <int dim, int spacedim>
471 void Triangulation<dim, spacedim>::copy_triangulation(
472 const ::Triangulation<dim, spacedim> &other_tria)
473 {
474 Assert(
475 (dynamic_cast<
476 const ::parallel::DistributedTriangulationBase<dim, spacedim>
477 *>(&other_tria) == nullptr),
479 "Cannot use this function on parallel::distributed::Triangulation."));
480
482 other_tria);
483 partition();
484 this->update_number_cache();
485 }
486 } // namespace shared
487} // namespace parallel
488
489#else
490
491namespace parallel
492{
493 namespace shared
494 {
495 template <int dim, int spacedim>
498 {
500 return true;
501 }
502
503
504
505 template <int dim, int spacedim>
508 const
509 {
510 return false;
511 }
512
513
514
515 template <int dim, int spacedim>
517 const std::vector<unsigned int>
519 {
521 return true_subdomain_ids_of_cells;
522 }
523
524
525
526 template <int dim, int spacedim>
528 const std::vector<unsigned int>
530 const unsigned int) const
531 {
533 return true_level_subdomain_ids_of_cells;
534 }
535 } // namespace shared
536} // namespace parallel
537
538
539#endif
540
541
542
543namespace internal
544{
545 namespace parallel
546 {
547 namespace shared
548 {
549 template <int dim, int spacedim>
552 : shared_tria(
553 dynamic_cast<
554 const ::parallel::shared::Triangulation<dim, spacedim> *>(
555 &tria))
556 {
557 if (shared_tria && shared_tria->with_artificial_cells())
558 {
559 // Save the current set of subdomain IDs, and set subdomain IDs
560 // to the "true" owner of each cell.
561 const std::vector<types::subdomain_id> &true_subdomain_ids =
562 shared_tria->get_true_subdomain_ids_of_cells();
563
564 saved_subdomain_ids.resize(shared_tria->n_active_cells());
565 for (const auto &cell : shared_tria->active_cell_iterators())
566 {
567 const unsigned int index = cell->active_cell_index();
568 saved_subdomain_ids[index] = cell->subdomain_id();
569 cell->set_subdomain_id(true_subdomain_ids[index]);
570 }
571 }
572 }
573
574
575
576 template <int dim, int spacedim>
579 {
580 if (shared_tria && shared_tria->with_artificial_cells())
581 {
582 // Undo the subdomain modification.
583 for (const auto &cell : shared_tria->active_cell_iterators())
584 {
585 const unsigned int index = cell->active_cell_index();
586 cell->set_subdomain_id(saved_subdomain_ids[index]);
587 }
588 }
589 }
590 } // namespace shared
591 } // namespace parallel
592} // namespace internal
593
594
595/*-------------- Explicit Instantiations -------------------------------*/
596#include "distributed/shared_tria.inst"
597
Definition point.h:113
cell_iterator end() const
const ObserverPointer< const ::parallel::shared::Triangulation< dim, spacedim > > shared_tria
TemporarilyRestoreSubdomainIds(const Triangulation< dim, spacedim > &tria)
types::subdomain_id my_subdomain
Definition tria_base.h:326
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:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:242
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int level
Definition grid_out.cc:4635
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::size_t size
Definition mpi.cc:745
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)
constexpr types::subdomain_id artificial_subdomain_id
Definition types.h:402
STL namespace.