Reference documentation for deal.II version Git 32ab9f15fd 2020-11-24 23:04:10 -0500
\(\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 - 2020 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 
31 template <int structdim>
32 CellData<structdim>::CellData(const unsigned int n_vertices)
33  : vertices(n_vertices, numbers::invalid_unsigned_int)
34  , material_id(0)
36 {}
37 
38 
39 
40 template <int structdim>
41 bool
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 
65 bool
66 SubCellData::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,
104  const TriangulationDescription::Settings settings,
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;
116  my_rank < ::Utilities::MPI::n_mpi_processes(comm),
117  ExcMessage("Rank has to be smaller than available processes."));
118 
119  if (auto tria_pdt = dynamic_cast<
121  {
122  Assert(
123  my_rank == numbers::invalid_unsigned_int ||
124  my_rank == ::Utilities::MPI::this_mpi_process(comm),
125  ExcMessage(
126  "If parallel::distributed::Triangulation as source triangulation, my_rank has to equal global rank."));
127 
128  my_rank = ::Utilities::MPI::this_mpi_process(comm);
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)
135  my_rank = ::Utilities::MPI::this_mpi_process(comm);
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 &
187  TriangulationDescription::Settings::construct_multigrid_hierarchy) ||
188  (tria.get_mesh_smoothing() &
190  ExcMessage(
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->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->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->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
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>
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,
411  const TriangulationDescription::Settings settings)
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 
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
unsigned int manifold_id
Definition: types.h:141
std::vector<::CellData< dim > > coarse_cells
TriangulationDescription::Utilities::create_description_from_triangulation()
Definition: mpi_tags.h:99
std::vector< std::pair< unsigned int, types::boundary_id > > boundary_ids
CellData(const unsigned int n_vertices=GeometryInfo< structdim >::vertices_per_cell)
bool operator==(const CellData< structdim > &other) const
std::vector< unsigned int > vertices
std::vector< types::coarse_cell_id > coarse_cell_index_to_coarse_cell_id
unsigned int material_id
Definition: types.h:152
std::array< types::manifold_id, GeometryInfo< dim >::lines_per_cell > manifold_line_ids
types::boundary_id boundary_id
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3013
Triangulation< dim, spacedim >::MeshSmoothing smoothing
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)
static ::ExceptionBase & ExcMessage(std::string arg1)
bool check_consistency(const unsigned int dim) const
std::array< types::manifold_id, dim==1 ? 1 :GeometryInfo< 3 >::quads_per_cell > manifold_quad_ids
#define Assert(cond, exc)
Definition: exceptions.h:1466
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)
types::material_id material_id
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
unsigned int level
Definition: grid_out.cc:4343
Point< 3 > vertices[4]
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:5542
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1182
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
types::manifold_id manifold_id
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
Definition: cuda.h:32
*braid_SplitCommworld & comm
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1747
std::vector< std::vector< CellData< dim > > > cell_infos
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
T min(const T &t, const MPI_Comm &mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
std::vector< Point< spacedim > > coarse_cell_vertices