Reference documentation for deal.II version Git 4e68a80cad 2021-10-22 15:50:12 -0600
\(\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>
19 
21 
24 
26 #include <deal.II/grid/tria.h>
28 
30 
31 
32 template <int structdim>
33 CellData<structdim>::CellData(const unsigned int n_vertices)
34  : vertices(n_vertices, numbers::invalid_unsigned_int)
35  , material_id(0)
37 {}
38 
39 
40 
41 template <int structdim>
42 bool
44 {
45  if (vertices.size() != other.vertices.size())
46  return false;
47 
48  for (unsigned int i = 0; i < vertices.size(); ++i)
49  if (vertices[i] != other.vertices[i])
50  return false;
51 
52  if (material_id != other.material_id)
53  return false;
54 
55  if (boundary_id != other.boundary_id)
56  return false;
57 
58  if (manifold_id != other.manifold_id)
59  return false;
60 
61  return true;
62 }
63 
64 
65 
66 bool
67 SubCellData::check_consistency(const unsigned int dim) const
68 {
69  switch (dim)
70  {
71  case 1:
72  return ((boundary_lines.size() == 0) && (boundary_quads.size() == 0));
73  case 2:
74  return (boundary_quads.size() == 0);
75  }
76  return true;
77 }
78 
80 {
81  namespace Utilities
82  {
83  namespace
84  {
88  template <int dim, int spacedim>
89  struct DescriptionTemp
90  {
95  template <class Archive>
96  void
97  serialize(Archive &ar, const unsigned int /*version*/)
98  {
99  ar &coarse_cells;
102  ar &cell_infos;
103  }
104 
108  void
109  collect(
110  const std::vector<unsigned int> & relevant_processes,
111  const std::vector<DescriptionTemp<dim, spacedim>> &description_temp,
112  const MPI_Comm & comm)
113  {
115  char>
116  process(
117  [&]() { return relevant_processes; },
118  [&](const unsigned int other_rank,
119  std::vector<char> &send_buffer) {
120  const auto ptr = std::find(relevant_processes.begin(),
121  relevant_processes.end(),
122  other_rank);
123 
124  Assert(ptr != relevant_processes.end(), ExcInternalError());
125 
126  const auto other_rank_index =
127  std::distance(relevant_processes.begin(), ptr);
128 
129  send_buffer =
130  ::Utilities::pack(description_temp[other_rank_index],
131  false);
132  },
133  [&](const unsigned int &,
134  const std::vector<char> &recv_buffer,
135  std::vector<char> &) {
136  this->merge(
137  ::Utilities::unpack<DescriptionTemp<dim, spacedim>>(
138  recv_buffer, false));
139  });
140 
142  process, comm)
143  .run();
144  }
145 
151  void
152  merge(const DescriptionTemp<dim, spacedim> &other)
153  {
154  this->cell_infos.resize(
155  std::max(other.cell_infos.size(), this->cell_infos.size()));
156 
157  this->coarse_cells.insert(this->coarse_cells.end(),
158  other.coarse_cells.begin(),
159  other.coarse_cells.end());
160  this->coarse_cell_vertices.insert(this->coarse_cell_vertices.end(),
161  other.coarse_cell_vertices.begin(),
162  other.coarse_cell_vertices.end());
165  other.coarse_cell_index_to_coarse_cell_id.begin(),
166  other.coarse_cell_index_to_coarse_cell_id.end());
167 
168  for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
169  this->cell_infos[i].insert(this->cell_infos[i].end(),
170  other.cell_infos[i].begin(),
171  other.cell_infos[i].end());
172  }
173 
177  void
178  reduce()
179  {
180  // make coarse cells unique
181  {
182  std::vector<std::tuple<types::coarse_cell_id,
184  unsigned int>>
185  temp;
186 
187  for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
188  temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
189  this->coarse_cells[i],
190  i);
191 
192  std::sort(temp.begin(),
193  temp.end(),
194  [](const auto &a, const auto &b) {
195  return std::get<0>(a) < std::get<0>(b);
196  });
197  temp.erase(std::unique(temp.begin(),
198  temp.end(),
199  [](const auto &a, const auto &b) {
200  return std::get<0>(a) == std::get<0>(b);
201  }),
202  temp.end());
203  std::sort(temp.begin(),
204  temp.end(),
205  [](const auto &a, const auto &b) {
206  return std::get<2>(a) < std::get<2>(b);
207  });
208 
209  this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
210  this->coarse_cells.resize(temp.size());
211 
212  for (unsigned int i = 0; i < temp.size(); ++i)
213  {
215  std::get<0>(temp[i]);
216  this->coarse_cells[i] = std::get<1>(temp[i]);
217  }
218  }
219 
220  // make coarse cell vertices unique
221  {
222  std::sort(this->coarse_cell_vertices.begin(),
223  this->coarse_cell_vertices.end(),
224  [](const auto &a, const auto &b) {
225  return a.first < b.first;
226  });
227  this->coarse_cell_vertices.erase(
228  std::unique(this->coarse_cell_vertices.begin(),
229  this->coarse_cell_vertices.end(),
230  [](const auto &a, const auto &b) {
231  return a.first == b.first;
232  }),
233  this->coarse_cell_vertices.end());
234  }
235 
236  // make cells unique
237  for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
238  {
239  if (this->cell_infos[i].size() == 0)
240  continue;
241 
242  std::sort(this->cell_infos[i].begin(),
243  this->cell_infos[i].end(),
244  [](const auto &a, const auto &b) {
245  return a.id < b.id;
246  });
247 
248  std::vector<CellData<dim>> temp;
249  temp.push_back(this->cell_infos[i][0]);
250 
251  for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j)
252  if (temp.back().id == cell_infos[i][j].id)
253  {
254  temp.back().subdomain_id =
255  std::min(temp.back().subdomain_id,
256  this->cell_infos[i][j].subdomain_id);
257  temp.back().level_subdomain_id =
258  std::min(temp.back().level_subdomain_id,
259  this->cell_infos[i][j].level_subdomain_id);
260  }
261  else
262  {
263  temp.push_back(this->cell_infos[i][j]);
264  }
265 
266  this->cell_infos[i] = temp;
267  }
268  }
269 
275  Description<dim, spacedim>
276  convert(const MPI_Comm comm,
278  mesh_smoothing,
280  {
281  Description<dim, spacedim> description;
282 
283  // copy communicator
284  description.comm = comm;
285 
286  description.settings = settings;
287 
288  // use mesh smoothing from base triangulation
289  description.smoothing = mesh_smoothing;
290 
291  std::map<unsigned int, unsigned int> map;
292 
293  for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
294  {
295  description.coarse_cell_vertices.push_back(
296  this->coarse_cell_vertices[i].second);
297  map[this->coarse_cell_vertices[i].first] = i;
298  }
299 
300  description.coarse_cells = this->coarse_cells;
301 
302  for (auto &cell : description.coarse_cells)
303  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
304  cell.vertices[v] = map[cell.vertices[v]];
305 
306  description.coarse_cell_index_to_coarse_cell_id =
307  this->coarse_cell_index_to_coarse_cell_id;
308  description.cell_infos = this->cell_infos;
309 
310  return description;
311  }
312 
313  std::vector<::CellData<dim>> coarse_cells;
314 
315  std::vector<std::pair<unsigned int, Point<spacedim>>>
317 
318  std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
319 
320  std::vector<std::vector<CellData<dim>>> cell_infos;
321  };
322 
326  template <int dim, int spacedim>
327  void
328  set_user_flag_and_of_its_parents(
330  {
331  cell->set_user_flag();
332  if (cell->level() != 0)
333  set_user_flag_and_of_its_parents(cell->parent());
334  }
335 
341  template <int dim, int spacedim>
342  class CreateDescriptionFromTriangulationHelper
343  {
344  public:
345  CreateDescriptionFromTriangulationHelper(
346  const ::Triangulation<dim, spacedim> &tria,
347  const std::function<types::subdomain_id(
348  const typename ::Triangulation<dim, spacedim>::cell_iterator
349  &)> & subdomain_id_function,
350  const std::function<types::subdomain_id(
351  const typename ::Triangulation<dim, spacedim>::cell_iterator
353  const MPI_Comm & comm,
355  : tria(tria)
358  , comm(comm)
359  , settings(settings)
361  settings &
363  {
364  Assert(
367  (tria.get_mesh_smoothing() &
368  Triangulation<dim,
369  spacedim>::limit_level_difference_at_vertices),
370  ExcMessage(
371  "Source triangulation has to be set up with "
372  "limit_level_difference_at_vertices if the construction of the "
373  "multigrid hierarchy is requested!"));
374 
377  }
378 
379  template <typename T>
380  T
381  create_description_for_rank(const unsigned int my_rank) const
382  {
383  T construction_data;
384 
385  set_additional_data(construction_data);
386 
387  // helper function, which collects all vertices belonging to a cell
388  // (also taking into account periodicity)
389  const auto
390  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
391  [this](const auto & cell,
392  std::vector<bool> &vertices_owned_by_locally_owned_cells) {
393  // add local vertices
394  for (const auto v : cell->vertex_indices())
395  {
396  vertices_owned_by_locally_owned_cells[cell->vertex_index(
397  v)] = true;
398  const auto coinciding_vertex_group =
400  cell->vertex_index(v));
401  if (coinciding_vertex_group !=
403  for (const auto &co_vertex : coinciding_vertex_groups.at(
404  coinciding_vertex_group->second))
405  vertices_owned_by_locally_owned_cells[co_vertex] = true;
406  }
407  };
408 
409  // 1) collect locally relevant cells (set user_flag)
410  std::vector<bool> old_user_flags;
411  tria.save_user_flags(old_user_flags);
412 
413  // 1a) clear user_flags
414  const_cast<::Triangulation<dim, spacedim> &>(tria)
415  .clear_user_flags();
416 
417  // 1b) loop over levels (from fine to coarse) and mark on each level
418  // the locally relevant cells
419  for (int level = tria.get_triangulation().n_global_levels() - 1;
420  level >= 0;
421  --level)
422  {
423  // collect vertices connected to a (on any level) locally owned
424  // cell
425  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
426  tria.n_vertices());
427  for (const auto &cell : tria.cell_iterators_on_level(level))
428  if (construct_multigrid &&
429  (level_subdomain_id_function(cell) == my_rank))
430  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
431  cell, vertices_owned_by_locally_owned_cells_on_level);
432 
433  for (const auto &cell : tria.active_cell_iterators())
434  if (subdomain_id_function(cell) == my_rank)
435  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
436  cell, vertices_owned_by_locally_owned_cells_on_level);
437 
438  // helper function to determine if cell is locally relevant
439  // (i.e. a cell which is connected to a vertex via a locally
440  // owned cell)
441  const auto is_locally_relevant_on_level = [&](const auto &cell) {
442  for (const auto v : cell->vertex_indices())
443  if (vertices_owned_by_locally_owned_cells_on_level
444  [cell->vertex_index(v)])
445  return true;
446  return false;
447  };
448 
449  // mark all locally relevant cells
450  for (const auto &cell : tria.cell_iterators_on_level(level))
451  if (is_locally_relevant_on_level(cell))
452  set_user_flag_and_of_its_parents(cell);
453  }
454 
455  // 2) set_up coarse-grid triangulation
456  {
457  std::vector<bool> vertices_locally_relevant(tria.n_vertices(),
458  false);
459 
460  // a) loop over all cells
461  for (const auto &cell : tria.cell_iterators_on_level(0))
462  {
463  if (!cell->user_flag_set())
464  continue;
465 
466  // extract cell definition (with old numbering of vertices)
467  ::CellData<dim> cell_data(cell->n_vertices());
468  cell_data.material_id = cell->material_id();
469  cell_data.manifold_id = cell->manifold_id();
470  for (const auto v : cell->vertex_indices())
471  cell_data.vertices[v] = cell->vertex_index(v);
472  construction_data.coarse_cells.push_back(cell_data);
473 
474  // save indices of each vertex of this cell
475  for (const auto v : cell->vertex_indices())
476  vertices_locally_relevant[cell->vertex_index(v)] = true;
477 
478  // save translation for corase grid: lid -> gid
479  construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
480  cell->id().get_coarse_cell_id());
481  }
482 
483  add_vertices(construction_data, vertices_locally_relevant);
484  }
485 
486 
487  // 3) collect info of each cell
488  construction_data.cell_infos.resize(
489  tria.get_triangulation().n_global_levels());
490 
491  // collect local vertices on active level
492  std::vector<bool> vertices_owned_by_locally_owned_active_cells(
493  tria.n_vertices());
494  for (const auto &cell : tria.active_cell_iterators())
495  if (subdomain_id_function(cell) == my_rank)
496  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
497  cell, vertices_owned_by_locally_owned_active_cells);
498 
499  // helper function to determine if cell is locally relevant
500  // on active level
501  const auto is_locally_relevant_on_active_level =
502  [&](const auto &cell) {
503  if (cell->is_active())
504  for (const auto v : cell->vertex_indices())
505  if (vertices_owned_by_locally_owned_active_cells
506  [cell->vertex_index(v)])
507  return true;
508  return false;
509  };
510 
511  for (unsigned int level = 0;
512  level < tria.get_triangulation().n_global_levels();
513  ++level)
514  {
515  // collect local vertices on level
516  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
517  tria.n_vertices());
518  for (const auto &cell : tria.cell_iterators_on_level(level))
519  if ((construct_multigrid &&
520  (level_subdomain_id_function(cell) == my_rank)) ||
521  (cell->is_active() &&
522  subdomain_id_function(cell) == my_rank))
523  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
524  cell, vertices_owned_by_locally_owned_cells_on_level);
525 
526  // helper function to determine if cell is locally relevant
527  // on level
528  const auto is_locally_relevant_on_level = [&](const auto &cell) {
529  for (const auto v : cell->vertex_indices())
530  if (vertices_owned_by_locally_owned_cells_on_level
531  [cell->vertex_index(v)])
532  return true;
533  return false;
534  };
535 
536  auto &level_cell_infos = construction_data.cell_infos[level];
537  for (const auto &cell : tria.cell_iterators_on_level(level))
538  {
539  // check if cell is locally relevant
540  if (!(cell->user_flag_set()))
541  continue;
542 
543  CellData<dim> cell_info;
544 
545  // save coarse-cell id
546  cell_info.id = cell->id().template to_binary<dim>();
547 
548  // save boundary_ids of each face of this cell
549  for (const auto f : cell->face_indices())
550  {
551  types::boundary_id boundary_ind =
552  cell->face(f)->boundary_id();
553  if (boundary_ind != numbers::internal_face_boundary_id)
554  cell_info.boundary_ids.emplace_back(f, boundary_ind);
555  }
556 
557  // save manifold id
558  {
559  // ... of cell
560  cell_info.manifold_id = cell->manifold_id();
561 
562  // ... of lines
563  if (dim >= 2)
564  for (const auto line : cell->line_indices())
565  cell_info.manifold_line_ids[line] =
566  cell->line(line)->manifold_id();
567 
568  // ... of quads
569  if (dim == 3)
570  for (const auto f : cell->face_indices())
571  cell_info.manifold_quad_ids[f] =
572  cell->quad(f)->manifold_id();
573  }
574 
575  // subdomain and level subdomain id
576  cell_info.subdomain_id = numbers::artificial_subdomain_id;
577  cell_info.level_subdomain_id =
579 
580  if (is_locally_relevant_on_active_level(cell))
581  {
582  cell_info.subdomain_id = subdomain_id_function(cell);
583 
584  cell_info.level_subdomain_id =
586  }
587  else if (is_locally_relevant_on_level(cell))
588  {
589  cell_info.level_subdomain_id =
591  }
592  else
593  {
594  // cell is locally relevant but an artificial cell
595  }
596 
597  level_cell_infos.emplace_back(cell_info);
598  }
599  }
600 
601  const_cast<::Triangulation<dim, spacedim> &>(tria)
602  .load_user_flags(old_user_flags);
603 
604  return construction_data;
605  }
606 
607  private:
608  void
609  set_additional_data(Description<dim, spacedim> &construction_data) const
610  {
611  construction_data.comm = comm;
612  construction_data.smoothing = tria.get_mesh_smoothing();
613  construction_data.settings = settings;
614  }
615 
616  void
617  set_additional_data(DescriptionTemp<dim, spacedim> &) const
618  {
619  // nothing to do
620  }
621 
622  void
623  add_vertices(Description<dim, spacedim> &construction_data,
624  const std::vector<bool> &vertices_locally_relevant) const
625  {
626  std::vector<unsigned int> vertices_locally_relevant_indices(
627  vertices_locally_relevant.size());
628 
629  // enumerate locally relevant vertices
630  unsigned int vertex_counter = 0;
631  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
632  if (vertices_locally_relevant[i])
633  {
634  construction_data.coarse_cell_vertices.push_back(
635  tria.get_vertices()[i]);
636  vertices_locally_relevant_indices[i] = vertex_counter++;
637  }
638 
639  // correct vertices of cells (make them local)
640  for (auto &cell : construction_data.coarse_cells)
641  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
642  cell.vertices[v] =
643  vertices_locally_relevant_indices[cell.vertices[v]];
644  }
645 
646  void
647  add_vertices(DescriptionTemp<dim, spacedim> &construction_data,
648  const std::vector<bool> &vertices_locally_relevant) const
649  {
650  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
651  if (vertices_locally_relevant[i])
652  construction_data.coarse_cell_vertices.emplace_back(
653  i, tria.get_vertices()[i]);
654  }
655 
656 
657  const ::Triangulation<dim, spacedim> &tria;
658  const std::function<types::subdomain_id(
659  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
661  const std::function<types::subdomain_id(
662  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
664 
665  const MPI_Comm & comm;
668 
669  std::map<unsigned int, std::vector<unsigned int>>
671  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
672  };
673 
674  } // namespace
675 
676 
677  template <int dim, int spacedim>
678  Description<dim, spacedim>
680  const ::Triangulation<dim, spacedim> &tria,
681  const MPI_Comm & comm,
683  const unsigned int my_rank_in)
684  {
685  if (const auto tria_pdt = dynamic_cast<
687  Assert(comm == tria_pdt->get_communicator(),
688  ExcMessage("MPI communicators do not match."));
689 
690  // First, figure out for what rank we are supposed to build the
691  // TriangulationDescription::Description object
692  unsigned int my_rank = my_rank_in;
694  my_rank < ::Utilities::MPI::n_mpi_processes(comm),
695  ExcMessage("Rank has to be smaller than available processes."));
696 
697  if (auto tria_pdt = dynamic_cast<
699  {
701  my_rank == ::Utilities::MPI::this_mpi_process(comm),
702  ExcMessage(
703  "For creation from a parallel::distributed::Triangulation, "
704  "my_rank has to equal global rank."));
705 
706  my_rank = ::Utilities::MPI::this_mpi_process(comm);
707  }
708  else if (auto tria_serial =
709  dynamic_cast<const ::Triangulation<dim, spacedim> *>(
710  &tria))
711  {
712  if (my_rank == numbers::invalid_unsigned_int)
713  my_rank = ::Utilities::MPI::this_mpi_process(comm);
714  }
715  else
716  {
717  Assert(false,
718  ExcMessage("This type of triangulation is not supported!"));
719  }
720 
721  const auto subdomain_id_function = [](const auto &cell) {
722  return cell->subdomain_id();
723  };
724 
725  const auto level_subdomain_id_function = [](const auto &cell) {
726  return cell->level_subdomain_id();
727  };
728 
729  return CreateDescriptionFromTriangulationHelper<dim, spacedim>(
730  tria,
733  comm,
734  settings)
735  .template create_description_for_rank<Description<dim, spacedim>>(
736  my_rank);
737  }
738 
739 
740 
741  template <int dim, int spacedim>
744  const std::function<void(::Triangulation<dim, spacedim> &)>
745  & serial_grid_generator,
746  const std::function<void(::Triangulation<dim, spacedim> &,
747  const MPI_Comm &,
748  const unsigned int)> &serial_grid_partitioner,
749  const MPI_Comm & comm,
750  const int group_size,
751  const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
753  {
754 #ifndef DEAL_II_WITH_MPI
755  (void)serial_grid_generator;
756  (void)serial_grid_partitioner;
757  (void)comm;
758  (void)group_size;
759  (void)smoothing;
760  (void)settings;
761 
763 #else
764  const unsigned int my_rank =
766  const unsigned int group_root = (my_rank / group_size) * group_size;
767 
768  const int mpi_tag =
770 
771  // check if process is root of the group
772  if (my_rank == group_root)
773  {
774  // Step 1: create serial triangulation
778  static_cast<
779  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
780  smoothing |
781  Triangulation<dim,
782  spacedim>::limit_level_difference_at_vertices) :
783  smoothing);
784  serial_grid_generator(tria);
785 
786  // Step 2: partition active cells and ...
787  serial_grid_partitioner(tria, comm, group_size);
788 
789  // ... cells on the levels if multigrid is required
790  if (settings &
793 
794  const unsigned int end_group =
795  std::min(group_root + group_size,
797 
798  // 3) create Description for the other processes in group; since
799  // we expect that this function is called for huge meshes, one
800  // Description is created at a time and send away; only once the
801  // Description has been sent away, the next rank is processed.
802  for (unsigned int other_rank = group_root + 1; other_rank < end_group;
803  ++other_rank)
804  {
805  // 3a) create construction data for other ranks
806  const auto construction_data =
808  comm,
809  settings,
810  other_rank);
811  // 3b) pack
812  std::vector<char> buffer;
813  ::Utilities::pack(construction_data, buffer, false);
814 
815  // 3c) send TriangulationDescription::Description
816  const auto ierr = MPI_Send(buffer.data(),
817  buffer.size(),
818  MPI_CHAR,
819  other_rank,
820  mpi_tag,
821  comm);
822  AssertThrowMPI(ierr);
823  }
824 
825  // 4) create TriangulationDescription::Description for this process
826  // (root of group)
828  comm,
829  settings,
830  my_rank);
831  }
832  else
833  {
834  // 3a) recv packed TriangulationDescription::Description from
835  // group-root process
836  // (counter-part of 3c of root process)
837  MPI_Status status;
838  auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
839  AssertThrowMPI(ierr);
840 
841  int len;
842  MPI_Get_count(&status, MPI_CHAR, &len);
843 
844  std::vector<char> buf(len);
845  ierr = MPI_Recv(buf.data(),
846  len,
847  MPI_CHAR,
848  status.MPI_SOURCE,
849  mpi_tag,
850  comm,
851  &status);
852  AssertThrowMPI(ierr);
853 
854  // 3b) unpack TriangulationDescription::Description (counter-part of
855  // 3b of root process)
856  auto construction_data =
857  ::Utilities::template unpack<Description<dim, spacedim>>(
858  buf, false);
859 
860  // WARNING: serialization cannot handle the MPI communicator
861  // which is the reason why we have to set it here explicitly
862  construction_data.comm = comm;
863 
864  return construction_data;
865  }
866 #endif
867  }
868 
869 
870 
871  template <int dim, int spacedim>
877  {
878  const bool construct_multigrid =
879  (partition.size() > 0) &&
880  (settings &
882 
883  Assert(
884  construct_multigrid == false ||
885  (tria.get_mesh_smoothing() &
887  ExcMessage(
888  "Source triangulation has to be set up with "
889  "limit_level_difference_at_vertices if the construction of the "
890  "multigrid hierarchy is requested!"));
891 
892  std::vector<LinearAlgebra::distributed::Vector<double>> partitions_mg;
893 
894  if (construct_multigrid) // perform first child policy
895  {
896  const auto tria_parallel =
897  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
898  &tria);
899 
900  Assert(tria_parallel, ExcInternalError());
901 
902  partition.update_ghost_values();
903 
904  partitions_mg.resize(tria.n_global_levels());
905 
906  for (unsigned int l = 0; l < tria.n_global_levels(); ++l)
907  partitions_mg[l].reinit(
908  tria_parallel->global_level_cell_index_partitioner(l).lock());
909 
910  for (int level = tria.n_global_levels() - 1; level >= 0; --level)
911  {
912  for (const auto &cell : tria.cell_iterators_on_level(level))
913  {
914  if (cell->is_locally_owned_on_level() == false)
915  continue;
916 
917  if (cell->is_active())
918  partitions_mg[level][cell->global_level_cell_index()] =
919  partition[cell->global_active_cell_index()];
920  else
921  partitions_mg[level][cell->global_level_cell_index()] =
922  partitions_mg[level + 1]
923  [cell->child(0)->global_level_cell_index()];
924  }
925 
926  partitions_mg[level].update_ghost_values();
927  }
928  }
929 
931  partition,
932  partitions_mg,
933  settings);
934  }
935 
936 
937 
938  template <int dim, int spacedim>
944  & partitions_mg,
945  const TriangulationDescription::Settings settings_in)
946  {
947 #ifdef DEAL_II_WITH_MPI
948  if (tria.get_communicator() == MPI_COMM_NULL)
949  AssertDimension(partition.local_size(), 0);
950 #endif
951 
952  if (partition.size() == 0)
953  {
954  AssertDimension(partitions_mg.size(), 0);
956  tria.get_communicator(),
957  settings_in);
958  }
959 
960  partition.update_ghost_values();
961 
962  for (const auto &partition : partitions_mg)
963  partition.update_ghost_values();
964 
965  // 1) determine processes owning locally owned cells
966  const std::vector<unsigned int> relevant_processes = [&]() {
967  std::set<unsigned int> relevant_processes;
968 
969  for (unsigned int i = 0; i < partition.local_size(); ++i)
970  relevant_processes.insert(
971  static_cast<unsigned int>(partition.local_element(i)));
972 
973  for (const auto &partition : partitions_mg)
974  for (unsigned int i = 0; i < partition.local_size(); ++i)
975  relevant_processes.insert(
976  static_cast<unsigned int>(partition.local_element(i)));
977 
978  return std::vector<unsigned int>(relevant_processes.begin(),
979  relevant_processes.end());
980  }();
981 
982  const bool construct_multigrid = partitions_mg.size() > 0;
983 
985 
986  if (construct_multigrid)
987  settings = static_cast<TriangulationDescription::Settings>(
988  settings |
990 
991  const auto subdomain_id_function = [&partition](const auto &cell) {
992  if ((cell->is_active() && (cell->is_artificial() == false)))
993  return static_cast<unsigned int>(
994  partition[cell->global_active_cell_index()]);
995  else
997  };
998 
999  const auto level_subdomain_id_function =
1000  [&construct_multigrid, &partitions_mg](const auto &cell) {
1001  if (construct_multigrid && (cell->is_artificial_on_level() == false))
1002  return static_cast<unsigned int>(
1003  partitions_mg[cell->level()][cell->global_level_cell_index()]);
1004  else
1006  };
1007 
1008  CreateDescriptionFromTriangulationHelper<dim, spacedim> helper(
1009  tria,
1012  tria.get_communicator(),
1013  settings);
1014 
1015  // create a description (locally owned cell and a layer of ghost cells
1016  // and all their parents)
1017  std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
1018  descriptions_per_rank.reserve(relevant_processes.size());
1019 
1020  for (const auto rank : relevant_processes)
1021  descriptions_per_rank.emplace_back(
1022  helper.template create_description_for_rank<
1023  DescriptionTemp<dim, spacedim>>(rank));
1024 
1025  // collect description from all processes that used to own locally-owned
1026  // active cells of this process in a single description
1027  DescriptionTemp<dim, spacedim> description_merged;
1028  description_merged.collect(relevant_processes,
1029  descriptions_per_rank,
1030  partition.get_mpi_communicator());
1031 
1032  // remove redundant entries
1033  description_merged.reduce();
1034 
1035  // convert to actual description
1036  return description_merged.convert(partition.get_mpi_communicator(),
1037  tria.get_mesh_smoothing(),
1038  settings);
1039  }
1040 
1041  } // namespace Utilities
1042 } // namespace TriangulationDescription
1043 
1044 
1045 
1046 /*-------------- Explicit Instantiations -------------------------------*/
1047 #include "tria_description.inst"
1048 
1049 
const bool construct_multigrid
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
unsigned int manifold_id
Definition: types.h:141
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
virtual MPI_Comm get_communicator() const
Definition: tria.cc:11233
TriangulationDescription::Utilities::create_description_from_triangulation()
Definition: mpi_tags.h:99
const ::Triangulation< dim, spacedim > & tria
CellData(const unsigned int n_vertices=GeometryInfo< structdim >::vertices_per_cell)
bool operator==(const CellData< structdim > &other) const
std::vector< unsigned int > vertices
unsigned int material_id
Definition: types.h:152
Point< 2 > second
Definition: grid_out.cc:4587
std::vector< types::coarse_cell_id > coarse_cell_index_to_coarse_cell_id
types::boundary_id boundary_id
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4126
const std::function< types::subdomain_id(const typename ::Triangulation< dim, spacedim >::cell_iterator &)> subdomain_id_function
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)
Number local_element(const size_type local_index) const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition: tria.cc:13275
static ::ExceptionBase & ExcMessage(std::string arg1)
bool check_consistency(const unsigned int dim) const
unsigned int subdomain_id
Definition: types.h:43
static const char T
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
#define Assert(cond, exc)
Definition: exceptions.h:1461
const std::function< types::subdomain_id(const typename ::Triangulation< dim, spacedim >::cell_iterator &)> level_subdomain_id_function
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::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
types::material_id material_id
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
VectorType::value_type * end(VectorType &V)
Point< 3 > vertices[4]
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
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:6342
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1219
std::vector< std::vector< CellData< dim > > > cell_infos
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
types::manifold_id manifold_id
std::vector<::CellData< dim > > coarse_cells
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
virtual size_type size() const override
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1778
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
const TriangulationDescription::Settings settings
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
virtual unsigned int n_global_levels() const
VectorType::value_type * begin(VectorType &V)
const MPI_Comm & comm
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1322
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
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)
const MPI_Comm & get_mpi_communicator() const
global_cell_index coarse_cell_id
Definition: types.h:114
std::vector< std::pair< unsigned int, Point< spacedim > > > coarse_cell_vertices
virtual const MeshSmoothing & get_mesh_smoothing() const
Definition: tria.cc:11254
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:691
void serialize(Archive &ar, const unsigned int version)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()