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