Reference documentation for deal.II version 9.3.3
\(\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 - 2021 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>
18
20
23
25#include <deal.II/grid/tria.h>
27
29
30
31template <int structdim>
32CellData<structdim>::CellData(const unsigned int n_vertices)
34 , material_id(0)
36{}
37
38
39
40template <int structdim>
41bool
43{
44 if (vertices.size() != other.vertices.size())
45 return false;
46
47 for (unsigned int i = 0; i < vertices.size(); ++i)
48 if (vertices[i] != other.vertices[i])
49 return false;
50
51 if (material_id != other.material_id)
52 return false;
53
54 if (boundary_id != other.boundary_id)
55 return false;
56
57 if (manifold_id != other.manifold_id)
58 return false;
59
60 return true;
61}
62
63
64
65bool
66SubCellData::check_consistency(const unsigned int dim) const
67{
68 switch (dim)
69 {
70 case 1:
71 return ((boundary_lines.size() == 0) && (boundary_quads.size() == 0));
72 case 2:
73 return (boundary_quads.size() == 0);
74 }
75 return true;
76}
77
79{
80 namespace Utilities
81 {
82 namespace
83 {
87 template <int dim, int spacedim>
88 void
89 set_user_flag_and_of_its_parents(
91 {
92 cell->set_user_flag();
93 if (cell->level() != 0)
94 set_user_flag_and_of_its_parents(cell->parent());
95 }
96 } // namespace
97
98
99 template <int dim, int spacedim>
100 Description<dim, spacedim>
102 const ::Triangulation<dim, spacedim> &tria,
103 const MPI_Comm & comm,
105 const unsigned int my_rank_in)
106 {
107 if (auto tria_pdt = dynamic_cast<
109 Assert(comm == tria_pdt->get_communicator(),
110 ExcMessage("MPI communicators do not match."));
111
112 // First, figure out for what rank we are supposed to build the
113 // TriangulationDescription::Description object
114 unsigned int my_rank = my_rank_in;
117 ExcMessage("Rank has to be smaller than available processes."));
118
119 if (auto tria_pdt = dynamic_cast<
121 {
122 Assert(
126 "If parallel::distributed::Triangulation as source triangulation, my_rank has to equal global rank."));
127
129 }
130 else if (auto tria_serial =
131 dynamic_cast<const ::Triangulation<dim, spacedim> *>(
132 &tria))
133 {
134 if (my_rank == numbers::invalid_unsigned_int)
136 }
137 else
138 {
139 Assert(false,
140 ExcMessage("This type of triangulation is not supported!"));
141 }
142
143 Description<dim, spacedim> construction_data;
144
145 // store the communicator
146 construction_data.comm = comm;
147
148
149 std::map<unsigned int, std::vector<unsigned int>>
150 coinciding_vertex_groups;
151 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
152
154 coinciding_vertex_groups,
155 vertex_to_coinciding_vertex_group);
156
157 // helper function, which collects all vertices belonging to a cell
158 // (also taking into account periodicity)
159 auto add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
160 [coinciding_vertex_groups, vertex_to_coinciding_vertex_group](
162 std::vector<bool> &vertices_owned_by_locally_owned_cells) {
163 // add local vertices
164 for (const auto v : cell->vertex_indices())
165 {
166 vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] =
167 true;
168 const auto coinciding_vertex_group =
169 vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
170 if (coinciding_vertex_group !=
171 vertex_to_coinciding_vertex_group.end())
172 for (const auto &co_vertex : coinciding_vertex_groups.at(
173 coinciding_vertex_group->second))
174 vertices_owned_by_locally_owned_cells[co_vertex] = true;
175 }
176 };
177
178 construction_data.smoothing = tria.get_mesh_smoothing();
179 construction_data.settings = settings;
180
181 const bool construct_multigrid =
182 settings &
184
185 Assert(
186 !(settings &
188 (tria.get_mesh_smoothing() &
191 "Source triangulation has to be setup with limit_level_difference_at_vertices if the construction of the multigrid hierarchy is requested!"));
192
193 // 1) collect locally relevant cells (set user_flag)
194 std::vector<bool> old_user_flags;
195 tria.save_user_flags(old_user_flags);
196
197 // 1a) clear user_flags
198 const_cast<::Triangulation<dim, spacedim> &>(tria)
199 .clear_user_flags();
200
201 // 1b) loop over levels (from fine to coarse) and mark on each level
202 // the locally relevant cells
203 for (int level = tria.get_triangulation().n_global_levels() - 1;
204 level >= 0;
205 --level)
206 {
207 // collect vertices connected to a (on any level) locally owned
208 // cell
209 std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
210 tria.n_vertices());
211 for (auto cell : tria.cell_iterators_on_level(level))
212 if ((construct_multigrid &&
213 (cell->level_subdomain_id() == my_rank)) ||
214 (cell->is_active() && cell->subdomain_id() == my_rank))
215 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
216 cell, vertices_owned_by_locally_owned_cells_on_level);
217
218 for (auto cell : tria.active_cell_iterators())
219 if (cell->subdomain_id() == my_rank)
220 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
221 cell, vertices_owned_by_locally_owned_cells_on_level);
222
223 // helper function to determine if cell is locally relevant
224 // (i.e. a cell which is connected to a vertex via a locally
225 // owned cell)
226 const auto is_locally_relevant_on_level =
228 for (const auto v : cell->vertex_indices())
229 if (vertices_owned_by_locally_owned_cells_on_level
230 [cell->vertex_index(v)])
231 return true;
232 return false;
233 };
234
235 // mark all locally relevant cells
236 for (auto cell : tria.cell_iterators_on_level(level))
237 if (is_locally_relevant_on_level(cell))
238 set_user_flag_and_of_its_parents(cell);
239 }
240
241 // 2) set_up coarse-grid triangulation
242 {
243 std::map<unsigned int, unsigned int> vertices_locally_relevant;
244
245 // a) loop over all cells
246 for (auto cell : tria.cell_iterators_on_level(0))
247 {
248 if (!cell->user_flag_set())
249 continue;
250
251 // extract cell definition (with old numbering of vertices)
252 ::CellData<dim> cell_data(cell->n_vertices());
253 cell_data.material_id = cell->material_id();
254 cell_data.manifold_id = cell->manifold_id();
255 for (const auto v : cell->vertex_indices())
256 cell_data.vertices[v] = cell->vertex_index(v);
257 construction_data.coarse_cells.push_back(cell_data);
258
259 // save indices of each vertex of this cell
260 for (const auto v : cell->vertex_indices())
261 vertices_locally_relevant[cell->vertex_index(v)] =
263
264 // save translation for corase grid: lid -> gid
265 construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
266 cell->id().get_coarse_cell_id());
267 }
268
269 // b) enumerate locally relevant vertices
270 unsigned int vertex_counter = 0;
271 for (auto &vertex : vertices_locally_relevant)
272 {
273 construction_data.coarse_cell_vertices.push_back(
274 tria.get_vertices()[vertex.first]);
275 vertex.second = vertex_counter++;
276 }
277
278 // c) correct vertices of cells (make them local)
279 for (auto &cell : construction_data.coarse_cells)
280 for (unsigned int v = 0; v < cell.vertices.size(); ++v)
281 cell.vertices[v] = vertices_locally_relevant[cell.vertices[v]];
282 }
283
284
285 // 3) collect info of each cell
286 construction_data.cell_infos.resize(
287 tria.get_triangulation().n_global_levels());
288
289 // collect local vertices on active level
290 std::vector<bool> vertices_owned_by_locally_owned_active_cells(
291 tria.n_vertices());
292 for (auto cell : tria.active_cell_iterators())
293 if (cell->subdomain_id() == my_rank)
294 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
295 cell, vertices_owned_by_locally_owned_active_cells);
296
297 // helper function to determine if cell is locally relevant
298 // on active level
299 const auto is_locally_relevant_on_active_level =
301 if (cell->is_active())
302 for (const auto v : cell->vertex_indices())
303 if (vertices_owned_by_locally_owned_active_cells
304 [cell->vertex_index(v)])
305 return true;
306 return false;
307 };
308
309 for (unsigned int level = 0;
310 level < tria.get_triangulation().n_global_levels();
311 ++level)
312 {
313 // collect local vertices on level
314 std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
315 tria.n_vertices());
316 for (auto cell : tria.cell_iterators_on_level(level))
317 if ((construct_multigrid &&
318 (cell->level_subdomain_id() == my_rank)) ||
319 (cell->is_active() && cell->subdomain_id() == my_rank))
320 add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
321 cell, vertices_owned_by_locally_owned_cells_on_level);
322
323 // helper function to determine if cell is locally relevant
324 // on level
325 const auto is_locally_relevant_on_level =
327 for (const auto v : cell->vertex_indices())
328 if (vertices_owned_by_locally_owned_cells_on_level
329 [cell->vertex_index(v)])
330 return true;
331 return false;
332 };
333
334 auto &level_cell_infos = construction_data.cell_infos[level];
335 for (auto cell : tria.cell_iterators_on_level(level))
336 {
337 // check if cell is locally relevant
338 if (!(cell->user_flag_set()))
339 continue;
340
341 CellData<dim> cell_info;
342
343 // save coarse-cell id
344 cell_info.id = cell->id().template to_binary<dim>();
345
346 // save boundary_ids of each face of this cell
347 for (const auto f : cell->face_indices())
348 {
349 types::boundary_id boundary_ind =
350 cell->face(f)->boundary_id();
351 if (boundary_ind != numbers::internal_face_boundary_id)
352 cell_info.boundary_ids.emplace_back(f, boundary_ind);
353 }
354
355 // save manifold id
356 {
357 // ... of cell
358 cell_info.manifold_id = cell->manifold_id();
359
360 // ... of lines
361 if (dim >= 2)
362 for (const auto line : cell->line_indices())
363 cell_info.manifold_line_ids[line] =
364 cell->line(line)->manifold_id();
365
366 // ... of quads
367 if (dim == 3)
368 for (const auto f : cell->face_indices())
369 cell_info.manifold_quad_ids[f] =
370 cell->quad(f)->manifold_id();
371 }
372
373 // subdomain and level subdomain id
374 cell_info.subdomain_id = numbers::artificial_subdomain_id;
375 cell_info.level_subdomain_id = numbers::artificial_subdomain_id;
376
377 if (is_locally_relevant_on_active_level(cell))
378 {
379 cell_info.level_subdomain_id = cell->level_subdomain_id();
380 cell_info.subdomain_id = cell->subdomain_id();
381 }
382 else if (is_locally_relevant_on_level(cell))
383 {
384 cell_info.level_subdomain_id = cell->level_subdomain_id();
385 }
386 // else: cell is locally relevant but an artificial cell
387
388 level_cell_infos.emplace_back(cell_info);
389 }
390 }
391
392 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags(
393 old_user_flags);
394
395 return construction_data;
396 }
397
398
399
400 template <int dim, int spacedim>
401 Description<dim, spacedim>
403 const std::function<void(::Triangulation<dim, spacedim> &)>
404 & serial_grid_generator,
405 const std::function<void(::Triangulation<dim, spacedim> &,
406 const MPI_Comm &,
407 const unsigned int)> &serial_grid_partitioner,
408 const MPI_Comm & comm,
409 const int group_size,
410 const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
412 {
413#ifndef DEAL_II_WITH_MPI
414 (void)serial_grid_generator;
415 (void)serial_grid_partitioner;
416 (void)comm;
417 (void)group_size;
418 (void)smoothing;
419 (void)settings;
420
422#else
423 const unsigned int my_rank =
425 const unsigned int group_root = (my_rank / group_size) * group_size;
426
427 const int mpi_tag =
429
430 // check if process is root of the group
431 if (my_rank == group_root)
432 {
433 // Step 1: create serial triangulation
437 static_cast<
438 typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
439 smoothing |
440 Triangulation<dim,
441 spacedim>::limit_level_difference_at_vertices) :
442 smoothing);
443 serial_grid_generator(tria);
444
445 // Step 2: partition active cells and ...
446 serial_grid_partitioner(tria, comm, group_size);
447
448 // ... cells on the levels if multigrid is required
449 if (settings &
452
453 const unsigned int end_group =
454 std::min(group_root + group_size,
456
457 // 3) create Description for the other processes in group; since
458 // we expect that this function is called for huge meshes, one
459 // Description is created at a time and send away; only once the
460 // Description has been sent away, the next rank is processed.
461 for (unsigned int other_rank = group_root + 1; other_rank < end_group;
462 ++other_rank)
463 {
464 // 3a) create construction data for other ranks
465 const auto construction_data =
467 comm,
468 settings,
469 other_rank);
470 // 3b) pack
471 std::vector<char> buffer;
472 ::Utilities::pack(construction_data, buffer, false);
473
474 // 3c) send TriangulationDescription::Description
475 const auto ierr = MPI_Send(buffer.data(),
476 buffer.size(),
477 MPI_CHAR,
478 other_rank,
479 mpi_tag,
480 comm);
481 AssertThrowMPI(ierr);
482 }
483
484 // 4) create TriangulationDescription::Description for this process
485 // (root of group)
487 comm,
488 settings,
489 my_rank);
490 }
491 else
492 {
493 // 3a) recv packed TriangulationDescription::Description from
494 // group-root process
495 // (counter-part of 3c of root process)
496 MPI_Status status;
497 auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
498 AssertThrowMPI(ierr);
499
500 int len;
501 MPI_Get_count(&status, MPI_CHAR, &len);
502
503 std::vector<char> buf(len);
504 ierr = MPI_Recv(buf.data(),
505 len,
506 MPI_CHAR,
507 status.MPI_SOURCE,
508 mpi_tag,
509 comm,
510 &status);
511 AssertThrowMPI(ierr);
512
513 // 3b) unpack TriangulationDescription::Description (counter-part of
514 // 3b of root process)
515 auto construction_data =
516 ::Utilities::template unpack<Description<dim, spacedim>>(
517 buf, false);
518
519 // WARNING: serialization cannot handle the MPI communicator
520 // which is the reason why we have to set it here explicitly
521 construction_data.comm = comm;
522
523 return construction_data;
524 }
525#endif
526 }
527
528 } // namespace Utilities
529} // namespace TriangulationDescription
530
531
532
533/*-------------- Explicit Instantiations -------------------------------*/
534#include "tria_description.inst"
535
536
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4590
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1746
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:6388
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4094
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)
@ 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:128
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1218
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
static const unsigned int invalid_unsigned_int
Definition: types.h:196
const types::manifold_id flat_manifold_id
Definition: types.h:264
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned int manifold_id
Definition: types.h:141
unsigned int material_id
Definition: types.h:152
unsigned int boundary_id
Definition: types.h:129
*braid_SplitCommworld & comm
std::vector< unsigned int > vertices
bool operator==(const CellData< structdim > &other) const
types::manifold_id manifold_id
CellData(const unsigned int n_vertices=GeometryInfo< structdim >::vertices_per_cell)
types::material_id material_id
types::boundary_id boundary_id
std::vector< CellData< 2 > > boundary_quads
bool check_consistency(const unsigned int dim) const
std::vector< CellData< 1 > > boundary_lines