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