Reference documentation for deal.II version 9.4.0
\(\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\}}\)
tria_description.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2019 - 2022 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
17#include <deal.II/base/mpi.h>
19
22
25
27#include <deal.II/grid/tria.h>
29
31
32
33template <int structdim>
34CellData<structdim>::CellData(const unsigned int n_vertices)
36 , material_id(0)
38{}
39
40
41
42template <int structdim>
43bool
45{
46 if (vertices.size() != other.vertices.size())
47 return false;
48
49 for (unsigned int i = 0; i < vertices.size(); ++i)
50 if (vertices[i] != other.vertices[i])
51 return false;
52
53 if (material_id != other.material_id)
54 return false;
55
56 if (boundary_id != other.boundary_id)
57 return false;
58
59 if (manifold_id != other.manifold_id)
60 return false;
61
62 return true;
63}
64
65
66
67bool
68SubCellData::check_consistency(const unsigned int dim) const
69{
70 switch (dim)
71 {
72 case 1:
73 return ((boundary_lines.size() == 0) && (boundary_quads.size() == 0));
74 case 2:
75 return (boundary_quads.size() == 0);
76 }
77 return true;
78}
79
81{
82 namespace Utilities
83 {
84 namespace
85 {
89 template <int dim, int spacedim>
90 struct DescriptionTemp
91 {
96 template <class Archive>
97 void
98 serialize(Archive &ar, const unsigned int /*version*/)
99 {
100 ar &coarse_cells;
103 ar &cell_infos;
104 }
105
109 void
110 collect(
111 const std::vector<unsigned int> & relevant_processes,
112 const std::vector<DescriptionTemp<dim, spacedim>> &description_temp,
113 const MPI_Comm & comm,
114 const bool vertices_have_unique_ids)
115 {
116 const auto create_request = [&](const unsigned int other_rank) {
117 const auto ptr = std::find(relevant_processes.begin(),
118 relevant_processes.end(),
119 other_rank);
120
121 Assert(ptr != relevant_processes.end(), ExcInternalError());
122
123 const auto other_rank_index =
124 std::distance(relevant_processes.begin(), ptr);
125
126 return description_temp[other_rank_index];
127 };
128
129 const auto process_request =
130 [&](const unsigned int,
131 const DescriptionTemp<dim, spacedim> &request) {
132 this->merge(request, vertices_have_unique_ids);
133 };
134
136 DescriptionTemp<dim, spacedim>>(relevant_processes,
137 create_request,
138 process_request,
139 comm);
140 }
141
147 void
148 merge(const DescriptionTemp<dim, spacedim> &other,
149 const bool vertices_have_unique_ids)
150 {
151 this->cell_infos.resize(
152 std::max(other.cell_infos.size(), this->cell_infos.size()));
153
154 if (vertices_have_unique_ids == false) // need to compare points
155 {
156 // comparator of points
157 const auto comp = [](const auto &a, const auto &b) {
158 const double tolerance = 1e-10;
159
160 for (unsigned int i = 0; i < spacedim; ++i)
161 if (std::abs(a[i] - b[i]) > tolerance)
162 return a[i] < b[i];
163
164 return false;
165 };
166
167 // map point to local vertex index
168 std::map<Point<spacedim>, unsigned int, decltype(comp)>
169 map_point_to_local_vertex_index(comp);
170
171 // ... initialize map with existing points
172 for (unsigned int i = 0; i < this->coarse_cell_vertices.size();
173 ++i)
174 map_point_to_local_vertex_index[coarse_cell_vertices[i]
175 .second] = i;
176
177 // map local vertex indices within other to the new local indices
178 std::map<unsigned int, unsigned int>
179 map_old_to_new_local_vertex_index;
180
181 // 1) renumerate vertices in other and insert into maps
182 unsigned int counter = coarse_cell_vertices.size();
183 for (const auto &p : other.coarse_cell_vertices)
184 if (map_point_to_local_vertex_index.find(p.second) ==
185 map_point_to_local_vertex_index.end())
186 {
187 this->coarse_cell_vertices.emplace_back(counter, p.second);
188 map_point_to_local_vertex_index[p.second] =
189 map_old_to_new_local_vertex_index[p.first] = counter++;
190 }
191 else
192 map_old_to_new_local_vertex_index[p.first] =
193 map_point_to_local_vertex_index[p.second];
194
195 // 2) renumerate vertices of cells
196 auto other_coarse_cells_copy = other.coarse_cells;
197
198 for (auto &cell : other_coarse_cells_copy)
199 for (auto &v : cell.vertices)
200 v = map_old_to_new_local_vertex_index[v];
201
202 this->coarse_cells.insert(this->coarse_cells.end(),
203 other_coarse_cells_copy.begin(),
204 other_coarse_cells_copy.end());
205 }
206 else
207 {
208 this->coarse_cells.insert(this->coarse_cells.end(),
209 other.coarse_cells.begin(),
210 other.coarse_cells.end());
211 this->coarse_cell_vertices.insert(
212 this->coarse_cell_vertices.end(),
213 other.coarse_cell_vertices.begin(),
214 other.coarse_cell_vertices.end());
215 }
216
219 other.coarse_cell_index_to_coarse_cell_id.begin(),
220 other.coarse_cell_index_to_coarse_cell_id.end());
221
222 for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
223 this->cell_infos[i].insert(this->cell_infos[i].end(),
224 other.cell_infos[i].begin(),
225 other.cell_infos[i].end());
226 }
227
231 void
232 reduce()
233 {
234 // make coarse cells unique
235 {
236 std::vector<std::tuple<types::coarse_cell_id,
238 unsigned int>>
239 temp;
240
241 for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
242 temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
243 this->coarse_cells[i],
244 i);
245
246 std::sort(temp.begin(),
247 temp.end(),
248 [](const auto &a, const auto &b) {
249 return std::get<0>(a) < std::get<0>(b);
250 });
251 temp.erase(std::unique(temp.begin(),
252 temp.end(),
253 [](const auto &a, const auto &b) {
254 return std::get<0>(a) == std::get<0>(b);
255 }),
256 temp.end());
257 std::sort(temp.begin(),
258 temp.end(),
259 [](const auto &a, const auto &b) {
260 return std::get<2>(a) < std::get<2>(b);
261 });
262
263 this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
264 this->coarse_cells.resize(temp.size());
265
266 for (unsigned int i = 0; i < temp.size(); ++i)
267 {
269 std::get<0>(temp[i]);
270 this->coarse_cells[i] = std::get<1>(temp[i]);
271 }
272 }
273
274 // make coarse cell vertices unique
275 {
276 std::sort(this->coarse_cell_vertices.begin(),
277 this->coarse_cell_vertices.end(),
278 [](const auto &a, const auto &b) {
279 return a.first < b.first;
280 });
281 this->coarse_cell_vertices.erase(
282 std::unique(this->coarse_cell_vertices.begin(),
283 this->coarse_cell_vertices.end(),
284 [](const auto &a, const auto &b) {
285 if (a.first == b.first)
286 {
287 Assert(a.second.distance(b.second) < 10e-8,
288 ExcInternalError());
289 return true;
290 }
291 return false;
292 }),
293 this->coarse_cell_vertices.end());
294 }
295
296 // make cells unique
297 for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
298 {
299 if (this->cell_infos[i].size() == 0)
300 continue;
301
302 std::sort(this->cell_infos[i].begin(),
303 this->cell_infos[i].end(),
304 [](const auto &a, const auto &b) {
305 return a.id < b.id;
306 });
307
308 std::vector<CellData<dim>> temp;
309 temp.push_back(this->cell_infos[i][0]);
310
311 for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j)
312 if (temp.back().id == cell_infos[i][j].id)
313 {
314 temp.back().subdomain_id =
315 std::min(temp.back().subdomain_id,
316 this->cell_infos[i][j].subdomain_id);
317 temp.back().level_subdomain_id =
318 std::min(temp.back().level_subdomain_id,
319 this->cell_infos[i][j].level_subdomain_id);
320 }
321 else
322 {
323 temp.push_back(this->cell_infos[i][j]);
324 }
325
326 this->cell_infos[i] = temp;
327 }
328 }
329
335 Description<dim, spacedim>
336 convert(const MPI_Comm comm,
338 mesh_smoothing,
340 {
341 Description<dim, spacedim> description;
342
343 // copy communicator
344 description.comm = comm;
345
346 description.settings = settings;
347
348 // use mesh smoothing from base triangulation
349 description.smoothing = mesh_smoothing;
350
351 std::map<unsigned int, unsigned int> map;
352
353 for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
354 {
355 description.coarse_cell_vertices.push_back(
357 map[this->coarse_cell_vertices[i].first] = i;
358 }
359
360 description.coarse_cells = this->coarse_cells;
361
362 for (auto &cell : description.coarse_cells)
363 for (unsigned int v = 0; v < cell.vertices.size(); ++v)
364 cell.vertices[v] = map[cell.vertices[v]];
365
366 description.coarse_cell_index_to_coarse_cell_id =
367 this->coarse_cell_index_to_coarse_cell_id;
368 description.cell_infos = this->cell_infos;
369
370 return description;
371 }
372
373 std::vector<::CellData<dim>> coarse_cells;
374
375 std::vector<std::pair<unsigned int, Point<spacedim>>>
377
378 std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
379
380 std::vector<std::vector<CellData<dim>>> cell_infos;
381 };
382
386 template <int dim, int spacedim>
387 void
388 set_user_flag_and_of_its_parents(
390 {
391 cell->set_user_flag();
392 if (cell->level() != 0)
393 set_user_flag_and_of_its_parents(cell->parent());
394 }
395
401 template <int dim, int spacedim>
402 class CreateDescriptionFromTriangulationHelper
403 {
404 public:
405 CreateDescriptionFromTriangulationHelper(
406 const ::Triangulation<dim, spacedim> &tria,
407 const std::function<types::subdomain_id(
408 const typename ::Triangulation<dim, spacedim>::cell_iterator
410 const std::function<types::subdomain_id(
411 const typename ::Triangulation<dim, spacedim>::cell_iterator
413 const MPI_Comm & comm,
415 : tria(tria)
418 , comm(comm)
422 0u)
423 {
424 Assert(
428 Triangulation<dim,
429 spacedim>::limit_level_difference_at_vertices),
431 "Source triangulation has to be set up with "
432 "limit_level_difference_at_vertices if the construction of the "
433 "multigrid hierarchy is requested!"));
434
437 }
438
439 template <typename T>
440 T
441 create_description_for_rank(const unsigned int my_rank) const
442 {
443 T construction_data;
444
445 set_additional_data(construction_data);
446
447 // helper function, which collects all vertices belonging to a cell
448 // (also taking into account periodicity)
449 const auto
450 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
451 [this](const auto & cell,
452 std::vector<bool> &vertices_owned_by_locally_owned_cells) {
453 // add local vertices
454 for (const auto v : cell->vertex_indices())
455 {
456 vertices_owned_by_locally_owned_cells[cell->vertex_index(
457 v)] = true;
458 const auto coinciding_vertex_group =
460 cell->vertex_index(v));
461 if (coinciding_vertex_group !=
463 for (const auto &co_vertex : coinciding_vertex_groups.at(
464 coinciding_vertex_group->second))
465 vertices_owned_by_locally_owned_cells[co_vertex] = true;
466 }
467 };
468
469 // 1) collect locally relevant cells (set user_flag)
470 std::vector<bool> old_user_flags;
471 tria.save_user_flags(old_user_flags);
472
473 // 1a) clear user_flags
475 .clear_user_flags();
476
477 // 1b) loop over levels (from fine to coarse) and mark on each level
478 // the locally relevant cells
479 for (int level = tria.get_triangulation().n_global_levels() - 1;
480 level >= 0;
481 --level)
482 {
483 // collect vertices connected to a (on any level) locally owned
484 // cell
485 std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
486 tria.n_vertices());
487 for (const auto &cell : tria.cell_iterators_on_level(level))
489 (level_subdomain_id_function(cell) == my_rank))
490 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
491 cell, vertices_owned_by_locally_owned_cells_on_level);
492
493 for (const auto &cell : tria.active_cell_iterators())
494 if (subdomain_id_function(cell) == my_rank)
495 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
496 cell, vertices_owned_by_locally_owned_cells_on_level);
497
498 // helper function to determine if cell is locally relevant
499 // (i.e. a cell which is connected to a vertex via a locally
500 // owned cell)
501 const auto is_locally_relevant_on_level = [&](const auto &cell) {
502 for (const auto v : cell->vertex_indices())
503 if (vertices_owned_by_locally_owned_cells_on_level
504 [cell->vertex_index(v)])
505 return true;
506 return false;
507 };
508
509 // mark all locally relevant cells
510 for (const auto &cell : tria.cell_iterators_on_level(level))
511 if (is_locally_relevant_on_level(cell))
512 set_user_flag_and_of_its_parents(cell);
513 }
514
515 // 2) set_up coarse-grid triangulation
516 {
517 std::vector<bool> vertices_locally_relevant(tria.n_vertices(),
518 false);
519
520 // a) loop over all cells
521 for (const auto &cell : tria.cell_iterators_on_level(0))
522 {
523 if (!cell->user_flag_set())
524 continue;
525
526 // extract cell definition (with old numbering of vertices)
527 ::CellData<dim> cell_data(cell->n_vertices());
528 cell_data.material_id = cell->material_id();
529 cell_data.manifold_id = cell->manifold_id();
530 for (const auto v : cell->vertex_indices())
531 cell_data.vertices[v] = cell->vertex_index(v);
532 construction_data.coarse_cells.push_back(cell_data);
533
534 // save indices of each vertex of this cell
535 for (const auto v : cell->vertex_indices())
536 vertices_locally_relevant[cell->vertex_index(v)] = true;
537
538 // save translation for corase grid: lid -> gid
539 construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
540 cell->id().get_coarse_cell_id());
541 }
542
543 add_vertices(construction_data, vertices_locally_relevant);
544 }
545
546
547 // 3) collect info of each cell
548 construction_data.cell_infos.resize(
550
551 // collect local vertices on active level
552 std::vector<bool> vertices_owned_by_locally_owned_active_cells(
553 tria.n_vertices());
554 for (const auto &cell : tria.active_cell_iterators())
555 if (subdomain_id_function(cell) == my_rank)
556 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
557 cell, vertices_owned_by_locally_owned_active_cells);
558
559 // helper function to determine if cell is locally relevant
560 // on active level
561 const auto is_locally_relevant_on_active_level =
562 [&](const auto &cell) {
563 if (cell->is_active())
564 for (const auto v : cell->vertex_indices())
565 if (vertices_owned_by_locally_owned_active_cells
566 [cell->vertex_index(v)])
567 return true;
568 return false;
569 };
570
571 for (unsigned int level = 0;
573 ++level)
574 {
575 // collect local vertices on level
576 std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
577 tria.n_vertices());
578 for (const auto &cell : tria.cell_iterators_on_level(level))
579 if ((construct_multigrid &&
580 (level_subdomain_id_function(cell) == my_rank)) ||
581 (cell->is_active() &&
582 subdomain_id_function(cell) == my_rank))
583 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
584 cell, vertices_owned_by_locally_owned_cells_on_level);
585
586 // helper function to determine if cell is locally relevant
587 // on level
588 const auto is_locally_relevant_on_level = [&](const auto &cell) {
589 for (const auto v : cell->vertex_indices())
590 if (vertices_owned_by_locally_owned_cells_on_level
591 [cell->vertex_index(v)])
592 return true;
593 return false;
594 };
595
596 auto &level_cell_infos = construction_data.cell_infos[level];
597 for (const auto &cell : tria.cell_iterators_on_level(level))
598 {
599 // check if cell is locally relevant
600 if (!(cell->user_flag_set()))
601 continue;
602
603 CellData<dim> cell_info;
604
605 // save coarse-cell id
606 cell_info.id = cell->id().template to_binary<dim>();
607
608 // save boundary_ids of each face of this cell
609 for (const auto f : cell->face_indices())
610 {
611 types::boundary_id boundary_ind =
612 cell->face(f)->boundary_id();
613 if (boundary_ind != numbers::internal_face_boundary_id)
614 cell_info.boundary_ids.emplace_back(f, boundary_ind);
615 }
616
617 // save manifold id
618 {
619 // ... of cell
620 cell_info.manifold_id = cell->manifold_id();
621
622 // ... of lines
623 if (dim >= 2)
624 for (const auto line : cell->line_indices())
625 cell_info.manifold_line_ids[line] =
626 cell->line(line)->manifold_id();
627
628 // ... of quads
629 if (dim == 3)
630 for (const auto f : cell->face_indices())
631 cell_info.manifold_quad_ids[f] =
632 cell->quad(f)->manifold_id();
633 }
634
635 // subdomain and level subdomain id
636 cell_info.subdomain_id = numbers::artificial_subdomain_id;
637 cell_info.level_subdomain_id =
639
640 if (is_locally_relevant_on_active_level(cell))
641 {
642 cell_info.subdomain_id = subdomain_id_function(cell);
643
644 cell_info.level_subdomain_id =
646 }
647 else if (is_locally_relevant_on_level(cell))
648 {
649 cell_info.level_subdomain_id =
651 }
652 else
653 {
654 // cell is locally relevant but an artificial cell
655 }
656
657 level_cell_infos.emplace_back(cell_info);
658 }
659 }
660
662 .load_user_flags(old_user_flags);
663
664 return construction_data;
665 }
666
667 private:
668 void
669 set_additional_data(Description<dim, spacedim> &construction_data) const
670 {
671 construction_data.comm = comm;
672 construction_data.smoothing = tria.get_mesh_smoothing();
673 construction_data.settings = settings;
674 }
675
676 void
677 set_additional_data(DescriptionTemp<dim, spacedim> &) const
678 {
679 // nothing to do
680 }
681
682 void
683 add_vertices(Description<dim, spacedim> &construction_data,
684 const std::vector<bool> &vertices_locally_relevant) const
685 {
686 std::vector<unsigned int> vertices_locally_relevant_indices(
687 vertices_locally_relevant.size());
688
689 // enumerate locally relevant vertices
690 unsigned int vertex_counter = 0;
691 for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
692 if (vertices_locally_relevant[i])
693 {
694 construction_data.coarse_cell_vertices.push_back(
695 tria.get_vertices()[i]);
696 vertices_locally_relevant_indices[i] = vertex_counter++;
697 }
698
699 // correct vertices of cells (make them local)
700 for (auto &cell : construction_data.coarse_cells)
701 for (unsigned int v = 0; v < cell.vertices.size(); ++v)
702 cell.vertices[v] =
703 vertices_locally_relevant_indices[cell.vertices[v]];
704 }
705
706 void
707 add_vertices(DescriptionTemp<dim, spacedim> &construction_data,
708 const std::vector<bool> &vertices_locally_relevant) const
709 {
710 for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
711 if (vertices_locally_relevant[i])
712 construction_data.coarse_cell_vertices.emplace_back(
713 i, tria.get_vertices()[i]);
714 }
715
716
717 const ::Triangulation<dim, spacedim> &tria;
718 const std::function<types::subdomain_id(
719 const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
721 const std::function<types::subdomain_id(
722 const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
724
725 const MPI_Comm & comm;
728
729 std::map<unsigned int, std::vector<unsigned int>>
731 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
732 };
733
734 } // namespace
735
736
737 template <int dim, int spacedim>
738 Description<dim, spacedim>
740 const ::Triangulation<dim, spacedim> &tria,
741 const MPI_Comm & comm,
743 const unsigned int my_rank_in)
744 {
745 if (const auto tria_pdt = dynamic_cast<
747 Assert(comm == tria_pdt->get_communicator(),
748 ExcMessage("MPI communicators do not match."));
749
750 // First, figure out for what rank we are supposed to build the
751 // TriangulationDescription::Description object
752 unsigned int my_rank = my_rank_in;
755 ExcMessage("Rank has to be smaller than available processes."));
756
757 if (auto tria_pdt = dynamic_cast<
759 {
763 "For creation from a parallel::distributed::Triangulation, "
764 "my_rank has to equal global rank."));
765
767 }
768 else if (auto tria_serial =
769 dynamic_cast<const ::Triangulation<dim, spacedim> *>(
770 &tria))
771 {
772 if (my_rank == numbers::invalid_unsigned_int)
774 }
775 else
776 {
777 Assert(false,
778 ExcMessage("This type of triangulation is not supported!"));
779 }
780
781 const auto subdomain_id_function = [](const auto &cell) {
782 return cell->subdomain_id();
783 };
784
785 const auto level_subdomain_id_function = [](const auto &cell) {
786 return cell->level_subdomain_id();
787 };
788
789 return CreateDescriptionFromTriangulationHelper<dim, spacedim>(
790 tria,
793 comm,
794 settings)
795 .template create_description_for_rank<Description<dim, spacedim>>(
796 my_rank);
797 }
798
799
800
801 template <int dim, int spacedim>
804 const std::function<void(::Triangulation<dim, spacedim> &)>
805 & serial_grid_generator,
806 const std::function<void(::Triangulation<dim, spacedim> &,
807 const MPI_Comm &,
808 const unsigned int)> &serial_grid_partitioner,
809 const MPI_Comm & comm,
810 const int group_size,
811 const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
813 {
814#ifndef DEAL_II_WITH_MPI
815 (void)serial_grid_generator;
816 (void)serial_grid_partitioner;
817 (void)comm;
818 (void)group_size;
819 (void)smoothing;
820 (void)settings;
821
823#else
824 const unsigned int my_rank =
826 const unsigned int group_root = (my_rank / group_size) * group_size;
827
828 const int mpi_tag =
830
831 // check if process is root of the group
832 if (my_rank == group_root)
833 {
834 // Step 1: create serial triangulation
838 static_cast<
839 typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
840 smoothing |
841 Triangulation<dim,
842 spacedim>::limit_level_difference_at_vertices) :
843 smoothing);
844 serial_grid_generator(tria);
845
846 // Step 2: partition active cells and ...
847 serial_grid_partitioner(tria, comm, group_size);
848
849 // ... cells on the levels if multigrid is required
850 if (settings &
853
854 const unsigned int end_group =
855 std::min(group_root + group_size,
857
858 // 3) create Description for the other processes in group; since
859 // we expect that this function is called for huge meshes, one
860 // Description is created at a time and send away; only once the
861 // Description has been sent away, the next rank is processed.
862 for (unsigned int other_rank = group_root + 1; other_rank < end_group;
863 ++other_rank)
864 {
865 // 3a) create construction data for other ranks
866 const auto construction_data =
868 comm,
869 settings,
870 other_rank);
871 // 3b) pack
872 std::vector<char> buffer;
873 ::Utilities::pack(construction_data, buffer, false);
874
875 // 3c) send TriangulationDescription::Description
876 const auto ierr = MPI_Send(buffer.data(),
877 buffer.size(),
878 MPI_CHAR,
879 other_rank,
880 mpi_tag,
881 comm);
882 AssertThrowMPI(ierr);
883 }
884
885 // 4) create TriangulationDescription::Description for this process
886 // (root of group)
888 comm,
889 settings,
890 my_rank);
891 }
892 else
893 {
894 // 3a) recv packed TriangulationDescription::Description from
895 // group-root process
896 // (counter-part of 3c of root process)
897 MPI_Status status;
898 auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
899 AssertThrowMPI(ierr);
900
901 int len;
902 MPI_Get_count(&status, MPI_CHAR, &len);
903
904 std::vector<char> buf(len);
905 ierr = MPI_Recv(buf.data(),
906 len,
907 MPI_CHAR,
908 status.MPI_SOURCE,
909 mpi_tag,
910 comm,
911 &status);
912 AssertThrowMPI(ierr);
913
914 // 3b) unpack TriangulationDescription::Description (counter-part of
915 // 3b of root process)
916 auto construction_data =
917 ::Utilities::template unpack<Description<dim, spacedim>>(
918 buf, false);
919
920 // WARNING: serialization cannot handle the MPI communicator
921 // which is the reason why we have to set it here explicitly
922 construction_data.comm = comm;
923
924 return construction_data;
925 }
926#endif
927 }
928
929
930
931 template <int dim, int spacedim>
937 {
938 const bool construct_multigrid =
939 (partition.size() > 0) &&
940 (settings &
942
943 Assert(
944 construct_multigrid == false ||
948 "Source triangulation has to be set up with "
949 "limit_level_difference_at_vertices if the construction of the "
950 "multigrid hierarchy is requested!"));
951
952 std::vector<LinearAlgebra::distributed::Vector<double>> partitions_mg;
953
954 if (construct_multigrid) // perform first child policy
955 {
956 const auto tria_parallel =
957 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
958 &tria);
959
960 Assert(tria_parallel, ExcInternalError());
961
962 partition.update_ghost_values();
963
964 partitions_mg.resize(tria.n_global_levels());
965
966 for (unsigned int l = 0; l < tria.n_global_levels(); ++l)
967 partitions_mg[l].reinit(
968 tria_parallel->global_level_cell_index_partitioner(l).lock());
969
970 for (int level = tria.n_global_levels() - 1; level >= 0; --level)
971 {
972 for (const auto &cell : tria.cell_iterators_on_level(level))
973 {
974 if (cell->is_locally_owned_on_level() == false)
975 continue;
976
977 if (cell->is_active())
978 partitions_mg[level][cell->global_level_cell_index()] =
979 partition[cell->global_active_cell_index()];
980 else
981 partitions_mg[level][cell->global_level_cell_index()] =
982 partitions_mg[level + 1]
983 [cell->child(0)->global_level_cell_index()];
984 }
985
986 partitions_mg[level].update_ghost_values();
987 }
988 }
989
991 partition,
992 partitions_mg,
993 settings);
994 }
995
996
997
998 template <int dim, int spacedim>
1004 & partitions_mg,
1005 const TriangulationDescription::Settings settings_in)
1006 {
1007#ifdef DEAL_II_WITH_MPI
1008 if (tria.get_communicator() == MPI_COMM_NULL)
1009 AssertDimension(partition.local_size(), 0);
1010#endif
1011
1012 if (partition.size() == 0)
1013 {
1014 AssertDimension(partitions_mg.size(), 0);
1017 settings_in);
1018 }
1019
1020 partition.update_ghost_values();
1021
1022 for (const auto &partition : partitions_mg)
1023 partition.update_ghost_values();
1024
1025 // 1) determine processes owning locally owned cells
1026 const std::vector<unsigned int> relevant_processes = [&]() {
1027 std::set<unsigned int> relevant_processes;
1028
1029 for (unsigned int i = 0; i < partition.local_size(); ++i)
1030 relevant_processes.insert(
1031 static_cast<unsigned int>(partition.local_element(i)));
1032
1033 for (const auto &partition : partitions_mg)
1034 for (unsigned int i = 0; i < partition.local_size(); ++i)
1035 relevant_processes.insert(
1036 static_cast<unsigned int>(partition.local_element(i)));
1037
1038 return std::vector<unsigned int>(relevant_processes.begin(),
1039 relevant_processes.end());
1040 }();
1041
1042 const bool construct_multigrid = partitions_mg.size() > 0;
1043
1045
1048 settings |
1050
1051 const auto subdomain_id_function = [&partition](const auto &cell) {
1052 if ((cell->is_active() && (cell->is_artificial() == false)))
1053 return static_cast<unsigned int>(
1054 partition[cell->global_active_cell_index()]);
1055 else
1057 };
1058
1059 const auto level_subdomain_id_function =
1060 [&construct_multigrid, &partitions_mg](const auto &cell) {
1061 if (construct_multigrid && (cell->is_artificial_on_level() == false))
1062 return static_cast<unsigned int>(
1063 partitions_mg[cell->level()][cell->global_level_cell_index()]);
1064 else
1066 };
1067
1068 CreateDescriptionFromTriangulationHelper<dim, spacedim> helper(
1069 tria,
1073 settings);
1074
1075 // create a description (locally owned cell and a layer of ghost cells
1076 // and all their parents)
1077 std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
1078 descriptions_per_rank.reserve(relevant_processes.size());
1079
1080 for (const auto rank : relevant_processes)
1081 descriptions_per_rank.emplace_back(
1082 helper.template create_description_for_rank<
1083 DescriptionTemp<dim, spacedim>>(rank));
1084
1085 // collect description from all processes that used to own locally-owned
1086 // active cells of this process in a single description
1087 DescriptionTemp<dim, spacedim> description_merged;
1088 description_merged.collect(
1089 relevant_processes,
1090 descriptions_per_rank,
1091 partition.get_mpi_communicator(),
1092 dynamic_cast<
1094 &tria) == nullptr);
1095
1096 // remove redundant entries
1097 description_merged.reduce();
1098
1099 // convert to actual description
1100 return description_merged.convert(partition.get_mpi_communicator(),
1102 settings);
1103 }
1104
1105 } // namespace Utilities
1106} // namespace TriangulationDescription
1107
1108
1109
1110/*-------------- Explicit Instantiations -------------------------------*/
1111#include "tria_description.inst"
1112
1113
virtual MPI_Comm get_communicator() const
virtual const MeshSmoothing & get_mesh_smoothing() const
void save_user_flags(std::ostream &out) const
const std::vector< Point< spacedim > > & get_vertices() const
virtual unsigned int n_global_levels() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4604
unsigned int level
Definition: grid_out.cc:4606
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
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:6398
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4178
static const char T
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
Description< dim, spacedim > create_description_from_triangulation_in_groups(const std::function< void(::Triangulation< dim, spacedim > &)> &serial_grid_generator, const std::function< void(::Triangulation< dim, spacedim > &, const MPI_Comm &, const unsigned int)> &serial_grid_partitioner, const MPI_Comm &comm, const int group_size=1, const typename Triangulation< dim, spacedim >::MeshSmoothing smoothing=::Triangulation< dim, spacedim >::none, const TriangulationDescription::Settings setting=TriangulationDescription::Settings::default_setting)
Description< dim, spacedim > create_description_from_triangulation(const ::Triangulation< dim, spacedim > &tria, const MPI_Comm &comm, const TriangulationDescription::Settings settings=TriangulationDescription::Settings::default_setting, const unsigned int my_rank_in=numbers::invalid_unsigned_int)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::vector< unsigned int > selector(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm &comm)
@ fully_distributed_create
TriangulationDescription::Utilities::create_description_from_triangulation()
Definition: mpi_tags.h:99
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:151
T reduce(const T &local_value, const MPI_Comm &comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:140
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1483
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
const types::boundary_id internal_face_boundary_id
Definition: types.h:260
const types::subdomain_id artificial_subdomain_id
Definition: types.h:298
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::manifold_id flat_manifold_id
Definition: types.h:269
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int manifold_id
Definition: types.h:141
global_cell_index coarse_cell_id
Definition: types.h:114
unsigned int subdomain_id
Definition: types.h:43
unsigned int material_id
Definition: types.h:152
unsigned int boundary_id
Definition: types.h:129
std::vector< unsigned int > vertices
bool operator==(const CellData< structdim > &other) const
types::manifold_id manifold_id
types::material_id material_id
types::boundary_id boundary_id
CellData(const unsigned int n_vertices=ReferenceCells::get_hypercube< structdim >().n_vertices())
std::vector< CellData< 2 > > boundary_quads
bool check_consistency(const unsigned int dim) const
std::vector< CellData< 1 > > boundary_lines
std::vector< std::pair< unsigned int, Point< spacedim > > > coarse_cell_vertices
const TriangulationDescription::Settings settings
const bool construct_multigrid
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
const std::function< types::subdomain_id(const typename ::Triangulation< dim, spacedim >::cell_iterator &)> subdomain_id_function
const std::function< types::subdomain_id(const typename ::Triangulation< dim, spacedim >::cell_iterator &)> level_subdomain_id_function
std::vector< types::coarse_cell_id > coarse_cell_index_to_coarse_cell_id
const MPI_Comm & comm
std::vector< std::vector< CellData< dim > > > cell_infos
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
const ::Triangulation< dim, spacedim > & tria
std::vector<::CellData< dim > > coarse_cells