Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+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
tria_description.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 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
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)
35 : vertices(n_vertices, numbers::invalid_unsigned_int)
36 , material_id(0)
37 , manifold_id(numbers::flat_manifold_id)
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.empty()) && (boundary_quads.empty()));
74 case 2:
75 return (boundary_quads.empty());
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> &future_owners_of_locally_owned_cells,
112 const std::vector<DescriptionTemp<dim, spacedim>> &description_temp,
113 const MPI_Comm comm,
114 const bool vertices_have_unique_ids)
115 {
116 // Use the some-to-some version of the consensus algorithm framework
117 // whereby we send requests to other processes that then deal with
118 // them but do not send anything back.
119 //
120 // Note that the input (description_temp) *may* contain an entry for
121 // the current process. As documented, the consensus algorithm will
122 // simply copy that into the output queue, i.e., it will call
123 // process_request() on it as well, and the data will simply come
124 // back out on the local process.
125 const auto create_request = [&](const unsigned int other_rank) {
126 const auto ptr =
127 std::find(future_owners_of_locally_owned_cells.begin(),
128 future_owners_of_locally_owned_cells.end(),
129 other_rank);
130
131 Assert(ptr != future_owners_of_locally_owned_cells.end(),
133
134 const auto other_rank_index =
135 std::distance(future_owners_of_locally_owned_cells.begin(), ptr);
136
137 return description_temp[other_rank_index];
138 };
139
140 const auto process_request =
141 [&](const unsigned int,
142 const DescriptionTemp<dim, spacedim> &request) -> void {
143 this->merge(request, vertices_have_unique_ids);
144 };
145
147 DescriptionTemp<dim, spacedim>>(
148 future_owners_of_locally_owned_cells,
149 create_request,
150 process_request,
151 comm);
152 }
153
159 void
160 merge(const DescriptionTemp<dim, spacedim> &other,
161 const bool vertices_have_unique_ids)
162 {
163 this->cell_infos.resize(
164 std::max(other.cell_infos.size(), this->cell_infos.size()));
165
166 if (vertices_have_unique_ids == false) // need to compare points
167 {
168 // map point to local vertex index
169 std::map<Point<spacedim>,
170 unsigned int,
172 map_point_to_local_vertex_index(
174
175 // ... initialize map with existing points
176 for (unsigned int i = 0; i < this->coarse_cell_vertices.size();
177 ++i)
178 map_point_to_local_vertex_index[coarse_cell_vertices[i]
179 .second] = i;
180
181 // map local vertex indices within other to the new local indices
182 std::map<unsigned int, unsigned int>
183 map_old_to_new_local_vertex_index;
184
185 // 1) re-enumerate vertices in other and insert into maps
186 unsigned int counter = coarse_cell_vertices.size();
187 for (const auto &p : other.coarse_cell_vertices)
188 if (map_point_to_local_vertex_index.find(p.second) ==
189 map_point_to_local_vertex_index.end())
190 {
191 this->coarse_cell_vertices.emplace_back(counter, p.second);
192 map_point_to_local_vertex_index[p.second] =
193 map_old_to_new_local_vertex_index[p.first] = counter++;
194 }
195 else
196 map_old_to_new_local_vertex_index[p.first] =
197 map_point_to_local_vertex_index[p.second];
198
199 // 2) re-enumerate vertices of cells
200 auto other_coarse_cells_copy = other.coarse_cells;
201
202 for (auto &cell : other_coarse_cells_copy)
203 for (auto &v : cell.vertices)
204 v = map_old_to_new_local_vertex_index[v];
205
206 this->coarse_cells.insert(this->coarse_cells.end(),
207 other_coarse_cells_copy.begin(),
208 other_coarse_cells_copy.end());
209 }
210 else
211 {
212 this->coarse_cells.insert(this->coarse_cells.end(),
213 other.coarse_cells.begin(),
214 other.coarse_cells.end());
215 this->coarse_cell_vertices.insert(
216 this->coarse_cell_vertices.end(),
217 other.coarse_cell_vertices.begin(),
218 other.coarse_cell_vertices.end());
219 }
220
223 other.coarse_cell_index_to_coarse_cell_id.begin(),
224 other.coarse_cell_index_to_coarse_cell_id.end());
225
226 for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
227 this->cell_infos[i].insert(this->cell_infos[i].end(),
228 other.cell_infos[i].begin(),
229 other.cell_infos[i].end());
230 }
231
235 void
236 reduce()
237 {
238 // make coarse cells unique
239 {
240 std::vector<std::tuple<types::coarse_cell_id,
242 unsigned int>>
243 temp;
244
245 for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
246 temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
247 this->coarse_cells[i],
248 i);
249
250 std::sort(temp.begin(),
251 temp.end(),
252 [](const auto &a, const auto &b) {
253 return std::get<0>(a) < std::get<0>(b);
254 });
255 temp.erase(std::unique(temp.begin(),
256 temp.end(),
257 [](const auto &a, const auto &b) {
258 return std::get<0>(a) == std::get<0>(b);
259 }),
260 temp.end());
261 std::sort(temp.begin(),
262 temp.end(),
263 [](const auto &a, const auto &b) {
264 return std::get<2>(a) < std::get<2>(b);
265 });
266
267 this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
268 this->coarse_cells.resize(temp.size());
269
270 for (unsigned int i = 0; i < temp.size(); ++i)
271 {
273 std::get<0>(temp[i]);
274 this->coarse_cells[i] = std::get<1>(temp[i]);
275 }
276 }
277
278 // make coarse cell vertices unique
279 {
280 std::sort(this->coarse_cell_vertices.begin(),
281 this->coarse_cell_vertices.end(),
282 [](const std::pair<unsigned int, Point<spacedim>> &a,
283 const std::pair<unsigned int, Point<spacedim>> &b) {
284 return a.first < b.first;
285 });
286 this->coarse_cell_vertices.erase(
287 std::unique(
288 this->coarse_cell_vertices.begin(),
289 this->coarse_cell_vertices.end(),
290 [](const std::pair<unsigned int, Point<spacedim>> &a,
291 const std::pair<unsigned int, Point<spacedim>> &b) {
292 if (a.first == b.first)
293 {
294 Assert(a.second.distance(b.second) <=
295 1e-7 *
296 std::max(a.second.norm(), b.second.norm()),
297 ExcMessage(
298 "In the process of merging the vertices of "
299 "the coarse meshes used on different processes, "
300 "there were two processes that used the same "
301 "vertex index for points that are not the same. "
302 "This suggests that you are using different "
303 "coarse meshes on different processes. This "
304 "should not happen."));
305 return true;
306 }
307 return false;
308 }),
309 this->coarse_cell_vertices.end());
310 }
311
312 // make cells unique
313 for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
314 {
315 if (this->cell_infos[i].empty())
316 continue;
317
318 std::sort(this->cell_infos[i].begin(),
319 this->cell_infos[i].end(),
320 [](const auto &a, const auto &b) {
321 return a.id < b.id;
322 });
323
324 std::vector<CellData<dim>> temp;
325 temp.push_back(this->cell_infos[i][0]);
326
327 for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j)
328 if (temp.back().id == cell_infos[i][j].id)
329 {
330 temp.back().subdomain_id =
331 std::min(temp.back().subdomain_id,
332 this->cell_infos[i][j].subdomain_id);
333 temp.back().level_subdomain_id =
334 std::min(temp.back().level_subdomain_id,
335 this->cell_infos[i][j].level_subdomain_id);
336 }
337 else
338 {
339 temp.push_back(this->cell_infos[i][j]);
340 }
341
342 this->cell_infos[i] = temp;
343 }
344 }
345
351 Description<dim, spacedim>
352 convert(const MPI_Comm comm,
354 mesh_smoothing,
356 {
357 Description<dim, spacedim> description;
358
359 // copy communicator
360 description.comm = comm;
361
362 description.settings = settings;
363
364 // use mesh smoothing from base triangulation
365 description.smoothing = mesh_smoothing;
366
367 std::map<unsigned int, unsigned int> map;
368
369 for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
370 {
371 description.coarse_cell_vertices.push_back(
373 map[this->coarse_cell_vertices[i].first] = i;
374 }
375
376 description.coarse_cells = this->coarse_cells;
377
378 for (auto &cell : description.coarse_cells)
379 for (unsigned int v = 0; v < cell.vertices.size(); ++v)
380 cell.vertices[v] = map[cell.vertices[v]];
381
382 description.coarse_cell_index_to_coarse_cell_id =
384 description.cell_infos = this->cell_infos;
385
386 return description;
387 }
388
389 std::vector<::CellData<dim>> coarse_cells;
390
391 std::vector<std::pair<unsigned int, Point<spacedim>>>
393
394 std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
395
396 std::vector<std::vector<CellData<dim>>> cell_infos;
397 };
398
402 template <int dim, int spacedim>
403 void
404 mark_cell_and_its_parents(
406 std::vector<std::vector<bool>> &cell_marked)
407 {
408 cell_marked[cell->level()][cell->index()] = true;
409 if (cell->level() != 0)
410 mark_cell_and_its_parents(cell->parent(), cell_marked);
411 }
412
418 template <typename DescriptionType, int dim, int spacedim>
419 DescriptionType
420 create_description_for_rank(
421 const ::Triangulation<dim, spacedim> &tria,
422 const std::function<types::subdomain_id(
423 const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
424 &subdomain_id_function,
425 const std::function<types::subdomain_id(
426 const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
427 &level_subdomain_id_function,
428 const std::map<unsigned int, std::vector<unsigned int>>
429 &coinciding_vertex_groups,
430 const std::map<unsigned int, unsigned int>
431 &vertex_to_coinciding_vertex_group,
432 const MPI_Comm comm,
433 const unsigned int my_rank,
435 {
436 static_assert(
437 std::is_same_v<DescriptionType, Description<dim, spacedim>> ||
438 std::is_same_v<DescriptionType, DescriptionTemp<dim, spacedim>>,
439 "Wrong template type.");
440 Assert(
443 (tria.get_mesh_smoothing() &
444 Triangulation<dim, spacedim>::limit_level_difference_at_vertices),
446 "Source triangulation has to be set up with "
447 "limit_level_difference_at_vertices if the construction of the "
448 "multigrid hierarchy is requested!"));
449
450 const bool construct_multigrid =
453
454 DescriptionType construction_data;
455 if constexpr (std::is_same_v<DescriptionType,
456 Description<dim, spacedim>>)
457 {
458 construction_data.comm = comm;
459 construction_data.smoothing = tria.get_mesh_smoothing();
460 construction_data.settings = settings;
461 }
462 else
463 (void)comm;
464
465 // A helper function that marks the indices of all vertices belonging
466 // to a cell (also taking into account their periodic breathren) in
467 // the bit vector passed as second argument.
468 const auto
469 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
470 [&coinciding_vertex_groups, &vertex_to_coinciding_vertex_group](
471 const typename ::Triangulation<dim, spacedim>::cell_iterator
472 &cell,
473 std::vector<bool> &vertices_on_locally_owned_cells) {
474 for (const unsigned int v : cell->vertex_indices())
475 {
476 const auto global_vertex_index = cell->vertex_index(v);
477 vertices_on_locally_owned_cells[global_vertex_index] = true;
478
479 const auto coinciding_vertex_group =
480 vertex_to_coinciding_vertex_group.find(global_vertex_index);
481 if (coinciding_vertex_group !=
482 vertex_to_coinciding_vertex_group.end())
483 for (const auto &co_vertex : coinciding_vertex_groups.at(
484 coinciding_vertex_group->second))
485 vertices_on_locally_owned_cells[co_vertex] = true;
486 }
487 };
488
489 const auto add_vertices =
490 [&tria](const std::vector<bool> &vertices_locally_relevant,
491 DescriptionType &construction_data) {
492 if constexpr (std::is_same_v<DescriptionType,
493 Description<dim, spacedim>>)
494 {
495 std::vector<unsigned int> vertices_locally_relevant_indices(
496 vertices_locally_relevant.size());
497
498 // enumerate locally relevant vertices
499 unsigned int vertex_counter = 0;
500 for (unsigned int i = 0; i < vertices_locally_relevant.size();
501 ++i)
502 if (vertices_locally_relevant[i])
503 {
504 construction_data.coarse_cell_vertices.push_back(
505 tria.get_vertices()[i]);
506 vertices_locally_relevant_indices[i] = vertex_counter++;
507 }
508
509 // correct vertices of cells (make them local)
510 for (auto &cell : construction_data.coarse_cells)
511 for (unsigned int v = 0; v < cell.vertices.size(); ++v)
512 cell.vertices[v] =
513 vertices_locally_relevant_indices[cell.vertices[v]];
514 }
515 else
516 {
517 for (unsigned int i = 0; i < vertices_locally_relevant.size();
518 ++i)
519 if (vertices_locally_relevant[i])
520 construction_data.coarse_cell_vertices.emplace_back(
521 i, tria.get_vertices()[i]);
522 }
523 };
524
525
526 // 1) loop over levels (from fine to coarse) and mark on each level
527 // the locally relevant cells
528 std::vector<std::vector<bool>> cell_marked(tria.n_levels());
529 for (unsigned int l = 0; l < tria.n_levels(); ++l)
530 cell_marked[l].resize(tria.n_raw_cells(l));
531
532 for (int level = tria.get_triangulation().n_global_levels() - 1;
533 level >= 0;
534 --level)
535 {
536 // collect vertices connected to a (on any level) locally owned
537 // cell
538 std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
539 tria.n_vertices());
540 for (const auto &cell : tria.cell_iterators_on_level(level))
541 if (construct_multigrid &&
542 (level_subdomain_id_function(cell) == my_rank))
543 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
544 cell, vertices_owned_by_locally_owned_cells_on_level);
545
546 for (const auto &cell : tria.active_cell_iterators())
547 if (subdomain_id_function(cell) == my_rank)
548 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
549 cell, vertices_owned_by_locally_owned_cells_on_level);
550
551 // helper function to determine if cell is locally relevant
552 // (i.e. a cell which is connected to a vertex via a locally
553 // owned cell)
554 const auto is_locally_relevant_on_level = [&](const auto &cell) {
555 for (const auto v : cell->vertex_indices())
556 if (vertices_owned_by_locally_owned_cells_on_level
557 [cell->vertex_index(v)])
558 return true;
559 return false;
560 };
561
562 // mark all locally relevant cells
563 for (const auto &cell : tria.cell_iterators_on_level(level))
564 if (is_locally_relevant_on_level(cell))
565 mark_cell_and_its_parents(cell, cell_marked);
566 }
567
568 // 2) set_up coarse-grid triangulation
569 {
570 std::vector<bool> vertices_locally_relevant(tria.n_vertices(), false);
571
572 // a) loop over all cells
573 for (const auto &cell : tria.cell_iterators_on_level(0))
574 {
575 if (!cell_marked[cell->level()][cell->index()])
576 continue;
577
578 // extract cell definition (with old numbering of vertices)
579 ::CellData<dim> cell_data(cell->n_vertices());
580 cell_data.material_id = cell->material_id();
581 cell_data.manifold_id = cell->manifold_id();
582 for (const auto v : cell->vertex_indices())
583 cell_data.vertices[v] = cell->vertex_index(v);
584 construction_data.coarse_cells.push_back(cell_data);
585
586 // save indices of each vertex of this cell
587 for (const auto v : cell->vertex_indices())
588 vertices_locally_relevant[cell->vertex_index(v)] = true;
589
590 // save translation for corase grid: lid -> gid
591 construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
592 cell->id().get_coarse_cell_id());
593 }
594
595 add_vertices(vertices_locally_relevant, construction_data);
596 }
597
598
599 // 3) collect info of each cell
600 construction_data.cell_infos.resize(
602
603 // collect local vertices on active level
604 std::vector<bool> vertices_owned_by_locally_owned_active_cells(
605 tria.n_vertices());
606 for (const auto &cell : tria.active_cell_iterators())
607 if (subdomain_id_function(cell) == my_rank)
608 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
609 cell, vertices_owned_by_locally_owned_active_cells);
610
611 // helper function to determine if cell is locally relevant
612 // on active level
613 const auto is_locally_relevant_on_active_level = [&](const auto &cell) {
614 if (cell->is_active())
615 for (const auto v : cell->vertex_indices())
616 if (vertices_owned_by_locally_owned_active_cells
617 [cell->vertex_index(v)])
618 return true;
619 return false;
620 };
621
622 for (unsigned int level = 0;
623 level < tria.get_triangulation().n_global_levels();
624 ++level)
625 {
626 // collect local vertices on level
627 std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
628 tria.n_vertices());
629 for (const auto &cell : tria.cell_iterators_on_level(level))
630 if ((construct_multigrid &&
631 (level_subdomain_id_function(cell) == my_rank)) ||
632 (cell->is_active() && subdomain_id_function(cell) == my_rank))
633 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
634 cell, vertices_owned_by_locally_owned_cells_on_level);
635
636 // helper function to determine if cell is locally relevant
637 // on level
638 const auto is_locally_relevant_on_level = [&](const auto &cell) {
639 for (const auto v : cell->vertex_indices())
640 if (vertices_owned_by_locally_owned_cells_on_level
641 [cell->vertex_index(v)])
642 return true;
643 return false;
644 };
645
646 auto &level_cell_infos = construction_data.cell_infos[level];
647 for (const auto &cell : tria.cell_iterators_on_level(level))
648 {
649 // check if cell is locally relevant
650 if (!cell_marked[cell->level()][cell->index()])
651 continue;
652
653 CellData<dim> cell_info;
654
655 // save coarse-cell id
656 cell_info.id = cell->id().template to_binary<dim>();
657
658 // save boundary_ids of each face of this cell
659 for (const auto f : cell->face_indices())
660 {
661 types::boundary_id boundary_ind =
662 cell->face(f)->boundary_id();
663 if (boundary_ind != numbers::internal_face_boundary_id)
664 cell_info.boundary_ids.emplace_back(f, boundary_ind);
665 }
666
667 // save manifold id
668 {
669 // ... of cell
670 cell_info.manifold_id = cell->manifold_id();
671
672 // ... of lines
673 if (dim >= 2)
674 for (const auto line : cell->line_indices())
675 cell_info.manifold_line_ids[line] =
676 cell->line(line)->manifold_id();
677
678 // ... of quads
679 if (dim == 3)
680 for (const auto f : cell->face_indices())
681 cell_info.manifold_quad_ids[f] =
682 cell->quad(f)->manifold_id();
683 }
684
685 // subdomain and level subdomain id
686 cell_info.subdomain_id = numbers::artificial_subdomain_id;
687 cell_info.level_subdomain_id = numbers::artificial_subdomain_id;
688
689 if (is_locally_relevant_on_active_level(cell))
690 {
691 cell_info.subdomain_id = subdomain_id_function(cell);
692
693 cell_info.level_subdomain_id =
694 level_subdomain_id_function(cell);
695 }
696 else if (is_locally_relevant_on_level(cell))
697 {
698 cell_info.level_subdomain_id =
699 level_subdomain_id_function(cell);
700 }
701 else
702 {
703 // cell is locally relevant but an artificial cell
704 }
705
706 level_cell_infos.emplace_back(cell_info);
707 }
708 }
709
710 return construction_data;
711 }
712 } // namespace
713
714
715 template <int dim, int spacedim>
716 Description<dim, spacedim>
718 const ::Triangulation<dim, spacedim> &tria,
719 const MPI_Comm comm,
721 const unsigned int my_rank_in)
722 {
723 if (const auto ptria =
724 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
725 &tria))
726 {
727 Assert(comm == ptria->get_communicator(),
728 ExcMessage("MPI communicators do not match."));
732 "For creation from a parallel::Triangulation, "
733 "my_rank has to equal the rank of the current process "
734 "in the given communicator."));
735 }
736
737 // If we are dealing with a sequential triangulation, then someone
738 // will have needed to set the subdomain_ids by hand. Make sure that
739 // all ids we see are less than the number of processes we are
740 // supposed to split the triangulation into.
741 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
742 &tria) == nullptr)
743 {
744#if DEBUG
745 const unsigned int n_mpi_processes =
747 for (const auto &cell : tria.active_cell_iterators())
748 Assert(cell->subdomain_id() < n_mpi_processes,
749 ExcMessage("You can't have a cell with subdomain_id of " +
750 std::to_string(cell->subdomain_id()) +
751 " when splitting the triangulation using an MPI "
752 " communicator with only " +
753 std::to_string(n_mpi_processes) + " processes."));
754#endif
755 }
756
757 // First, figure out for what rank we are supposed to build the
758 // TriangulationDescription::Description object
759 const unsigned int my_rank =
760 (my_rank_in == numbers::invalid_unsigned_int ?
762 my_rank_in);
763
764 const auto subdomain_id_function = [](const auto &cell) {
765 return cell->subdomain_id();
766 };
767
768 const auto level_subdomain_id_function = [](const auto &cell) {
769 return cell->level_subdomain_id();
770 };
771
772 std::map<unsigned int, std::vector<unsigned int>>
773 coinciding_vertex_groups;
774 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
776 coinciding_vertex_groups,
777 vertex_to_coinciding_vertex_group);
778
779 return create_description_for_rank<Description<dim, spacedim>>(
780 tria,
781 subdomain_id_function,
782 level_subdomain_id_function,
783 coinciding_vertex_groups,
784 vertex_to_coinciding_vertex_group,
785 comm,
786 my_rank,
787 settings);
788 }
789
790
791
792 template <int dim, int spacedim>
795 const std::function<void(::Triangulation<dim, spacedim> &)>
796 &serial_grid_generator,
797 const std::function<void(::Triangulation<dim, spacedim> &,
798 const MPI_Comm,
799 const unsigned int)> &serial_grid_partitioner,
800 const MPI_Comm comm,
801 const int group_size,
802 const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
804 {
805#ifndef DEAL_II_WITH_MPI
806 (void)serial_grid_generator;
807 (void)serial_grid_partitioner;
808 (void)comm;
809 (void)group_size;
810 (void)smoothing;
811 (void)settings;
812
814#else
815 const unsigned int my_rank =
817 const unsigned int group_root = (my_rank / group_size) * group_size;
818
819 const int mpi_tag =
821
822 // check if process is root of the group
823 if (my_rank == group_root)
824 {
825 // Step 1: create serial triangulation
829 static_cast<
830 typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
831 smoothing |
832 Triangulation<dim,
833 spacedim>::limit_level_difference_at_vertices) :
834 smoothing);
835 serial_grid_generator(tria);
836
837 // Step 2: partition active cells and ...
838 serial_grid_partitioner(tria, comm, group_size);
839
840 // ... cells on the levels if multigrid is required
841 if (settings &
844
845 const unsigned int end_group =
846 std::min(group_root + group_size,
848
849 // 3) create Description for the other processes in group; since
850 // we expect that this function is called for huge meshes, one
851 // Description is created at a time and send away; only once the
852 // Description has been sent away, the next rank is processed.
853 for (unsigned int other_rank = group_root + 1; other_rank < end_group;
854 ++other_rank)
855 {
856 // 3a) create construction data for other ranks
857 const auto construction_data =
859 comm,
860 settings,
861 other_rank);
862 // 3b) pack
863 std::vector<char> buffer;
864 ::Utilities::pack(construction_data, buffer, false);
865
866 // 3c) send TriangulationDescription::Description
867 const auto ierr = MPI_Send(buffer.data(),
868 buffer.size(),
869 MPI_CHAR,
870 other_rank,
871 mpi_tag,
872 comm);
873 AssertThrowMPI(ierr);
874 }
875
876 // 4) create TriangulationDescription::Description for this process
877 // (root of group)
879 comm,
880 settings,
881 my_rank);
882 }
883 else
884 {
885 // 3a) recv packed TriangulationDescription::Description from
886 // group-root process
887 // (counter-part of 3c of root process)
888 MPI_Status status;
889 auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
890 AssertThrowMPI(ierr);
891
892 int len;
893 MPI_Get_count(&status, MPI_CHAR, &len);
894
895 std::vector<char> buf(len);
896 ierr = MPI_Recv(buf.data(),
897 len,
898 MPI_CHAR,
899 status.MPI_SOURCE,
900 mpi_tag,
901 comm,
902 &status);
903 AssertThrowMPI(ierr);
904
905 // 3b) unpack TriangulationDescription::Description (counter-part of
906 // 3b of root process)
907 auto construction_data =
908 ::Utilities::template unpack<Description<dim, spacedim>>(
909 buf, false);
910
911 // WARNING: serialization cannot handle the MPI communicator
912 // which is the reason why we have to set it here explicitly
913 construction_data.comm = comm;
914
915 return construction_data;
916 }
917#endif
918 }
919
920
921
922 template <int dim, int spacedim>
928 {
929 const bool construct_multigrid =
930 (partition.size() > 0) &&
931 (settings &
933
934 Assert(
935 construct_multigrid == false ||
936 (tria.get_mesh_smoothing() &
939 "Source triangulation has to be set up with "
940 "limit_level_difference_at_vertices if the construction of the "
941 "multigrid hierarchy is requested!"));
942
943 std::vector<LinearAlgebra::distributed::Vector<double>> partitions_mg;
944
945 // If desired, also create a multigrid hierarchy. For this, we have to
946 // build a hierarchy of partitions (one for each level of the
947 // triangulation) in which each cell is assigned to the same process
948 // as its first child (if not active) or to the same process that already
949 // owns the cell (for an active level-cell).
950 if (construct_multigrid)
951 {
952 const auto tria_parallel =
953 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
954 &tria);
955 Assert(tria_parallel, ExcNotImplemented());
956
957 // Give the level partitioners the right size:
958 partitions_mg.resize(tria.n_global_levels());
959 for (unsigned int l = 0; l < tria.n_global_levels(); ++l)
960 partitions_mg[l].reinit(
961 tria_parallel->global_level_cell_index_partitioner(l).lock());
962
963 // Make sure we know about all of the owners of the active cells,
964 // whether locally owned or not. Then we traverse the triangulation
965 // from the finest level to the coarsest level.
966 //
967 // On each level, traverse the cell. If the cell is not locally
968 // owned, we don't care about it. If it is active, we copy the
969 // owner process from the cell's non-level owner. Otherwise,
970 // use the owner of the first cell.
971 partition.update_ghost_values();
972 for (int level = tria.n_global_levels() - 1; level >= 0; --level)
973 {
974 for (const auto &cell : tria.cell_iterators_on_level(level))
975 {
976 if (cell->is_locally_owned_on_level())
977 {
978 if (cell->is_active())
979 partitions_mg[level][cell->global_level_cell_index()] =
980 partition[cell->global_active_cell_index()];
981 else
982 partitions_mg[level][cell->global_level_cell_index()] =
983 partitions_mg[level + 1]
984 [cell->child(0)
985 ->global_level_cell_index()];
986 }
987 }
988
989 // Having touched all of the locally owned cells on the
990 // current level, exchange information with the other processes
991 // about the cells that are ghosts so that on the next coarser
992 // level we can access information about children again:
993 partitions_mg[level].update_ghost_values();
994 }
995 }
996
997 // Forward to the other function.
999 partition,
1000 partitions_mg,
1001 settings);
1002 }
1003
1004
1005
1006 template <int dim, int spacedim>
1009 const Triangulation<dim, spacedim> &tria,
1012 &partitions_mg,
1013 const TriangulationDescription::Settings settings_in)
1014 {
1015#ifdef DEAL_II_WITH_MPI
1016 if (tria.get_communicator() == MPI_COMM_NULL)
1017 AssertDimension(partition.locally_owned_size(), 0);
1018#endif
1019
1020 if (partition.size() == 0)
1021 {
1022 AssertDimension(partitions_mg.size(), 0);
1024 tria.get_communicator(),
1025 settings_in);
1026 }
1027
1028 // Update partitioner ghost elements because we will later want
1029 // to ask also about the future owners of ghost cells.
1030 partition.update_ghost_values();
1031 for (const auto &partition : partitions_mg)
1032 partition.update_ghost_values();
1033
1034 // 1) Determine process ids that appear on locally owned cells. Create
1035 // a sorted vector by first creating a std::set and then copying
1036 // the result. (Note that we get only locally *owned* cells in
1037 // the output because we only loop over the locally *owned*
1038 // entries of the partitioning vector, even though
1039 // 'partition.local_element(i)' could also return locally relevant
1040 // elements if 'i' were to exceed the number of locally owned
1041 // elements.)
1042 const std::vector<unsigned int> future_owners_of_locally_owned_cells =
1043 [&partition, &partitions_mg]() {
1044 std::set<unsigned int> relevant_process_set;
1045
1046 const unsigned int n_mpi_ranks =
1048 partition.get_mpi_communicator());
1049 (void)n_mpi_ranks;
1050
1051 for (unsigned int i = 0; i < partition.locally_owned_size(); ++i)
1052 {
1053 Assert(static_cast<unsigned int>(partition.local_element(i)) ==
1054 partition.local_element(i),
1055 ExcMessage(
1056 "The elements of a partition vector must be integers."));
1057 Assert(
1058 partition.local_element(i) < n_mpi_ranks,
1059 ExcMessage(
1060 "The elements of a partition vector must be between zero "
1061 "and the number of processes in the communicator "
1062 "to be used for partitioning the triangulation."));
1063 relevant_process_set.insert(
1064 static_cast<unsigned int>(partition.local_element(i)));
1065 }
1066
1067 for (const auto &partition : partitions_mg)
1068 for (unsigned int i = 0; i < partition.locally_owned_size(); ++i)
1069 {
1070 Assert(
1071 static_cast<unsigned int>(partition.local_element(i)) ==
1072 partition.local_element(i),
1073 ExcMessage(
1074 "The elements of a partition vector must be integers."));
1075 Assert(
1076 partition.local_element(i) < n_mpi_ranks,
1077 ExcMessage(
1078 "The elements of a partition vector must be between zero "
1079 "and the number of processes in the communicator "
1080 "to be used for partitioning the triangulation."));
1081 relevant_process_set.insert(
1082 static_cast<unsigned int>(partition.local_element(i)));
1083 }
1084
1085 return std::vector<unsigned int>(relevant_process_set.begin(),
1086 relevant_process_set.end());
1087 }();
1088
1089 const bool construct_multigrid = (partitions_mg.size() > 0);
1090
1091 const TriangulationDescription::Settings settings =
1092 (construct_multigrid ?
1096 settings_in);
1097
1098
1099 // Set up a function that returns the future owner rank for a cell.
1100 // Same then also for the level owner. These functions work for
1101 // locally owned and ghost cells.
1102 const auto cell_to_future_owner =
1103 [&partition](const auto &cell) -> types::subdomain_id {
1104 if ((cell->is_active() && (cell->is_artificial() == false)))
1105 return static_cast<types::subdomain_id>(
1106 partition[cell->global_active_cell_index()]);
1107 else
1109 };
1110
1111 const auto mg_cell_to_future_owner =
1112 [&construct_multigrid,
1113 &partitions_mg](const auto &cell) -> types::subdomain_id {
1114 if (construct_multigrid && (cell->is_artificial_on_level() == false))
1115 return static_cast<types::subdomain_id>(
1116 partitions_mg[cell->level()][cell->global_level_cell_index()]);
1117 else
1119 };
1120
1121 // Create a description (locally owned cell and a layer of ghost cells
1122 // and all their parents). We first create a description in the
1123 // 'temporary' format (using class DescriptionTemp), which we will
1124 // later convert to its final form.
1125 std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
1126 descriptions_per_rank.reserve(
1127 future_owners_of_locally_owned_cells.size());
1128
1129 std::map<unsigned int, std::vector<unsigned int>>
1130 coinciding_vertex_groups;
1131 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
1133 coinciding_vertex_groups,
1134 vertex_to_coinciding_vertex_group);
1135
1136 for (const auto rank : future_owners_of_locally_owned_cells)
1137 descriptions_per_rank.emplace_back(
1138 create_description_for_rank<DescriptionTemp<dim, spacedim>>(
1139 tria,
1140 cell_to_future_owner,
1141 mg_cell_to_future_owner,
1142 coinciding_vertex_groups,
1143 vertex_to_coinciding_vertex_group,
1144 tria.get_communicator(),
1145 rank,
1146 settings));
1147
1148 // Collect description from all processes that used to own locally-owned
1149 // active cells of this process in a single description
1150 DescriptionTemp<dim, spacedim> description_merged;
1151 description_merged.collect(
1152 future_owners_of_locally_owned_cells,
1153 descriptions_per_rank,
1154 partition.get_mpi_communicator(),
1155 dynamic_cast<
1157 &tria) == nullptr);
1158
1159 // remove redundant entries
1160 description_merged.reduce();
1161
1162 // convert to actual description
1163 return description_merged.convert(partition.get_mpi_communicator(),
1164 tria.get_mesh_smoothing(),
1165 settings);
1166 }
1167
1168 } // namespace Utilities
1169} // namespace TriangulationDescription
1170
1171
1172
1173/*-------------- Explicit Instantiations -------------------------------*/
1174#include "tria_description.inst"
1175
1176
Definition point.h:111
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4624
unsigned int level
Definition grid_out.cc:4626
unsigned int vertex_indices[2]
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
const Triangulation< dim, spacedim > & tria
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)
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
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)
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:98
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107
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:1381
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:312
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
static const unsigned int invalid_unsigned_int
Definition types.h:220
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const Iterator & begin
Definition parallel.h:609
::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:156
std::uint64_t global_vertex_index
Definition types.h:48
global_cell_index coarse_cell_id
Definition types.h:130
*braid_SplitCommworld & comm
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
std::vector< types::coarse_cell_id > coarse_cell_index_to_coarse_cell_id
std::vector< std::vector< CellData< dim > > > cell_infos
std::vector<::CellData< dim > > coarse_cells