Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
17 #include <deal.II/base/mpi.h>
19 
22 
25 
27 #include <deal.II/grid/tria.h>
29 
31 
32 
33 template <int structdim>
34 CellData<structdim>::CellData(const unsigned int n_vertices)
35  : vertices(n_vertices, numbers::invalid_unsigned_int)
36  , material_id(0)
38 {}
39 
40 
41 
42 template <int structdim>
43 bool
45 {
46  if (vertices.size() != other.vertices.size())
47  return false;
48 
49  for (unsigned int i = 0; i < vertices.size(); ++i)
50  if (vertices[i] != other.vertices[i])
51  return false;
52 
53  if (material_id != other.material_id)
54  return false;
55 
56  if (boundary_id != other.boundary_id)
57  return false;
58 
59  if (manifold_id != other.manifold_id)
60  return false;
61 
62  return true;
63 }
64 
65 
66 
67 bool
68 SubCellData::check_consistency(const unsigned int dim) const
69 {
70  switch (dim)
71  {
72  case 1:
73  return ((boundary_lines.size() == 0) && (boundary_quads.size() == 0));
74  case 2:
75  return (boundary_quads.size() == 0);
76  }
77  return true;
78 }
79 
81 {
82  namespace Utilities
83  {
84  namespace
85  {
89  template <int dim, int spacedim>
90  struct DescriptionTemp
91  {
96  template <class Archive>
97  void
98  serialize(Archive &ar, const unsigned int /*version*/)
99  {
100  ar &coarse_cells;
103  ar &cell_infos;
104  }
105 
109  void
110  collect(
111  const std::vector<unsigned int> & relevant_processes,
112  const std::vector<DescriptionTemp<dim, spacedim>> &description_temp,
113  const MPI_Comm & comm,
114  const bool vertices_have_unique_ids)
115  {
116  const auto create_request = [&](const unsigned int other_rank) {
117  const auto ptr = std::find(relevant_processes.begin(),
118  relevant_processes.end(),
119  other_rank);
120 
121  Assert(ptr != relevant_processes.end(), ExcInternalError());
122 
123  const auto other_rank_index =
124  std::distance(relevant_processes.begin(), ptr);
125 
126  return description_temp[other_rank_index];
127  };
128 
129  const auto process_request =
130  [&](const unsigned int,
131  const DescriptionTemp<dim, spacedim> &request) {
132  this->merge(request, vertices_have_unique_ids);
133  };
134 
136  DescriptionTemp<dim, spacedim>>(relevant_processes,
137  create_request,
138  process_request,
139  comm);
140  }
141 
147  void
148  merge(const DescriptionTemp<dim, spacedim> &other,
149  const bool vertices_have_unique_ids)
150  {
151  this->cell_infos.resize(
152  std::max(other.cell_infos.size(), this->cell_infos.size()));
153 
154  if (vertices_have_unique_ids == false) // need to compare points
155  {
156  // comparator of points
157  const auto comp = [](const auto &a, const auto &b) {
158  const double tolerance = 1e-10;
159 
160  for (unsigned int i = 0; i < spacedim; ++i)
161  if (std::abs(a[i] - b[i]) > tolerance)
162  return a[i] < b[i];
163 
164  return false;
165  };
166 
167  // map point to local vertex index
168  std::map<Point<spacedim>, unsigned int, decltype(comp)>
169  map_point_to_local_vertex_index(comp);
170 
171  // ... initialize map with existing points
172  for (unsigned int i = 0; i < this->coarse_cell_vertices.size();
173  ++i)
174  map_point_to_local_vertex_index[coarse_cell_vertices[i]
175  .second] = i;
176 
177  // map local vertex indices within other to the new local indices
178  std::map<unsigned int, unsigned int>
179  map_old_to_new_local_vertex_index;
180 
181  // 1) renumerate vertices in other and insert into maps
182  unsigned int counter = coarse_cell_vertices.size();
183  for (const auto &p : other.coarse_cell_vertices)
184  if (map_point_to_local_vertex_index.find(p.second) ==
185  map_point_to_local_vertex_index.end())
186  {
187  this->coarse_cell_vertices.emplace_back(counter, p.second);
188  map_point_to_local_vertex_index[p.second] =
189  map_old_to_new_local_vertex_index[p.first] = counter++;
190  }
191  else
192  map_old_to_new_local_vertex_index[p.first] =
193  map_point_to_local_vertex_index[p.second];
194 
195  // 2) renumerate vertices of cells
196  auto other_coarse_cells_copy = other.coarse_cells;
197 
198  for (auto &cell : other_coarse_cells_copy)
199  for (auto &v : cell.vertices)
200  v = map_old_to_new_local_vertex_index[v];
201 
202  this->coarse_cells.insert(this->coarse_cells.end(),
203  other_coarse_cells_copy.begin(),
204  other_coarse_cells_copy.end());
205  }
206  else
207  {
208  this->coarse_cells.insert(this->coarse_cells.end(),
209  other.coarse_cells.begin(),
210  other.coarse_cells.end());
211  this->coarse_cell_vertices.insert(
212  this->coarse_cell_vertices.end(),
213  other.coarse_cell_vertices.begin(),
214  other.coarse_cell_vertices.end());
215  }
216 
219  other.coarse_cell_index_to_coarse_cell_id.begin(),
220  other.coarse_cell_index_to_coarse_cell_id.end());
221 
222  for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
223  this->cell_infos[i].insert(this->cell_infos[i].end(),
224  other.cell_infos[i].begin(),
225  other.cell_infos[i].end());
226  }
227 
231  void
232  reduce()
233  {
234  // make coarse cells unique
235  {
236  std::vector<std::tuple<types::coarse_cell_id,
238  unsigned int>>
239  temp;
240 
241  for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
242  temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
243  this->coarse_cells[i],
244  i);
245 
246  std::sort(temp.begin(),
247  temp.end(),
248  [](const auto &a, const auto &b) {
249  return std::get<0>(a) < std::get<0>(b);
250  });
251  temp.erase(std::unique(temp.begin(),
252  temp.end(),
253  [](const auto &a, const auto &b) {
254  return std::get<0>(a) == std::get<0>(b);
255  }),
256  temp.end());
257  std::sort(temp.begin(),
258  temp.end(),
259  [](const auto &a, const auto &b) {
260  return std::get<2>(a) < std::get<2>(b);
261  });
262 
263  this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
264  this->coarse_cells.resize(temp.size());
265 
266  for (unsigned int i = 0; i < temp.size(); ++i)
267  {
269  std::get<0>(temp[i]);
270  this->coarse_cells[i] = std::get<1>(temp[i]);
271  }
272  }
273 
274  // make coarse cell vertices unique
275  {
276  std::sort(this->coarse_cell_vertices.begin(),
277  this->coarse_cell_vertices.end(),
278  [](const auto &a, const auto &b) {
279  return a.first < b.first;
280  });
281  this->coarse_cell_vertices.erase(
282  std::unique(this->coarse_cell_vertices.begin(),
283  this->coarse_cell_vertices.end(),
284  [](const auto &a, const auto &b) {
285  if (a.first == b.first)
286  {
287  Assert(a.second.distance(b.second) < 10e-8,
288  ExcInternalError());
289  return true;
290  }
291  return false;
292  }),
293  this->coarse_cell_vertices.end());
294  }
295 
296  // make cells unique
297  for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
298  {
299  if (this->cell_infos[i].size() == 0)
300  continue;
301 
302  std::sort(this->cell_infos[i].begin(),
303  this->cell_infos[i].end(),
304  [](const auto &a, const auto &b) {
305  return a.id < b.id;
306  });
307 
308  std::vector<CellData<dim>> temp;
309  temp.push_back(this->cell_infos[i][0]);
310 
311  for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j)
312  if (temp.back().id == cell_infos[i][j].id)
313  {
314  temp.back().subdomain_id =
315  std::min(temp.back().subdomain_id,
316  this->cell_infos[i][j].subdomain_id);
317  temp.back().level_subdomain_id =
318  std::min(temp.back().level_subdomain_id,
319  this->cell_infos[i][j].level_subdomain_id);
320  }
321  else
322  {
323  temp.push_back(this->cell_infos[i][j]);
324  }
325 
326  this->cell_infos[i] = temp;
327  }
328  }
329 
335  Description<dim, spacedim>
336  convert(const MPI_Comm comm,
338  mesh_smoothing,
340  {
341  Description<dim, spacedim> description;
342 
343  // copy communicator
344  description.comm = comm;
345 
346  description.settings = settings;
347 
348  // use mesh smoothing from base triangulation
349  description.smoothing = mesh_smoothing;
350 
351  std::map<unsigned int, unsigned int> map;
352 
353  for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
354  {
355  description.coarse_cell_vertices.push_back(
356  this->coarse_cell_vertices[i].second);
357  map[this->coarse_cell_vertices[i].first] = i;
358  }
359 
360  description.coarse_cells = this->coarse_cells;
361 
362  for (auto &cell : description.coarse_cells)
363  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
364  cell.vertices[v] = map[cell.vertices[v]];
365 
366  description.coarse_cell_index_to_coarse_cell_id =
367  this->coarse_cell_index_to_coarse_cell_id;
368  description.cell_infos = this->cell_infos;
369 
370  return description;
371  }
372 
373  std::vector<::CellData<dim>> coarse_cells;
374 
375  std::vector<std::pair<unsigned int, Point<spacedim>>>
377 
378  std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
379 
380  std::vector<std::vector<CellData<dim>>> cell_infos;
381  };
382 
386  template <int dim, int spacedim>
387  void
388  mark_cell_and_its_parents(
390  std::vector<std::vector<bool>> & cell_marked)
391  {
392  cell_marked[cell->level()][cell->index()] = true;
393  if (cell->level() != 0)
394  mark_cell_and_its_parents(cell->parent(), cell_marked);
395  }
396 
402  template <int dim, int spacedim>
403  class CreateDescriptionFromTriangulationHelper
404  {
405  public:
406  CreateDescriptionFromTriangulationHelper(
407  const ::Triangulation<dim, spacedim> &tria,
408  const std::function<types::subdomain_id(
409  const typename ::Triangulation<dim, spacedim>::cell_iterator
410  &)> & subdomain_id_function,
411  const std::function<types::subdomain_id(
412  const typename ::Triangulation<dim, spacedim>::cell_iterator
414  const MPI_Comm & comm,
416  : tria(tria)
419  , comm(comm)
420  , settings(settings)
423  0u)
424  {
425  Assert(
429  Triangulation<dim,
430  spacedim>::limit_level_difference_at_vertices),
431  ExcMessage(
432  "Source triangulation has to be set up with "
433  "limit_level_difference_at_vertices if the construction of the "
434  "multigrid hierarchy is requested!"));
435 
438  }
439 
440  template <typename T>
441  T
442  create_description_for_rank(const unsigned int my_rank) const
443  {
444  T construction_data;
445 
446  set_additional_data(construction_data);
447 
448  // helper function, which collects all vertices belonging to a cell
449  // (also taking into account periodicity)
450  const auto
451  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
452  [this](const auto & cell,
453  std::vector<bool> &vertices_owned_by_locally_owned_cells) {
454  // add local vertices
455  for (const auto v : cell->vertex_indices())
456  {
457  vertices_owned_by_locally_owned_cells[cell->vertex_index(
458  v)] = true;
459  const auto coinciding_vertex_group =
461  cell->vertex_index(v));
462  if (coinciding_vertex_group !=
464  for (const auto &co_vertex : coinciding_vertex_groups.at(
465  coinciding_vertex_group->second))
466  vertices_owned_by_locally_owned_cells[co_vertex] = true;
467  }
468  };
469 
470  // 1) loop over levels (from fine to coarse) and mark on each level
471  // the locally relevant cells
472  std::vector<std::vector<bool>> cell_marked(tria.n_levels());
473  for (unsigned int l = 0; l < tria.n_levels(); ++l)
474  cell_marked[l].resize(tria.n_raw_cells(l));
475 
476  for (int level = tria.get_triangulation().n_global_levels() - 1;
477  level >= 0;
478  --level)
479  {
480  // collect vertices connected to a (on any level) locally owned
481  // cell
482  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
483  tria.n_vertices());
484  for (const auto &cell : tria.cell_iterators_on_level(level))
485  if (construct_multigrid &&
486  (level_subdomain_id_function(cell) == my_rank))
487  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
488  cell, vertices_owned_by_locally_owned_cells_on_level);
489 
490  for (const auto &cell : tria.active_cell_iterators())
491  if (subdomain_id_function(cell) == my_rank)
492  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
493  cell, vertices_owned_by_locally_owned_cells_on_level);
494 
495  // helper function to determine if cell is locally relevant
496  // (i.e. a cell which is connected to a vertex via a locally
497  // owned cell)
498  const auto is_locally_relevant_on_level = [&](const auto &cell) {
499  for (const auto v : cell->vertex_indices())
500  if (vertices_owned_by_locally_owned_cells_on_level
501  [cell->vertex_index(v)])
502  return true;
503  return false;
504  };
505 
506  // mark all locally relevant cells
507  for (const auto &cell : tria.cell_iterators_on_level(level))
508  if (is_locally_relevant_on_level(cell))
509  mark_cell_and_its_parents(cell, cell_marked);
510  }
511 
512  // 2) set_up coarse-grid triangulation
513  {
514  std::vector<bool> vertices_locally_relevant(tria.n_vertices(),
515  false);
516 
517  // a) loop over all cells
518  for (const auto &cell : tria.cell_iterators_on_level(0))
519  {
520  if (!cell_marked[cell->level()][cell->index()])
521  continue;
522 
523  // extract cell definition (with old numbering of vertices)
524  ::CellData<dim> cell_data(cell->n_vertices());
525  cell_data.material_id = cell->material_id();
526  cell_data.manifold_id = cell->manifold_id();
527  for (const auto v : cell->vertex_indices())
528  cell_data.vertices[v] = cell->vertex_index(v);
529  construction_data.coarse_cells.push_back(cell_data);
530 
531  // save indices of each vertex of this cell
532  for (const auto v : cell->vertex_indices())
533  vertices_locally_relevant[cell->vertex_index(v)] = true;
534 
535  // save translation for corase grid: lid -> gid
536  construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
537  cell->id().get_coarse_cell_id());
538  }
539 
540  add_vertices(construction_data, vertices_locally_relevant);
541  }
542 
543 
544  // 3) collect info of each cell
545  construction_data.cell_infos.resize(
547 
548  // collect local vertices on active level
549  std::vector<bool> vertices_owned_by_locally_owned_active_cells(
550  tria.n_vertices());
551  for (const auto &cell : tria.active_cell_iterators())
552  if (subdomain_id_function(cell) == my_rank)
553  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
554  cell, vertices_owned_by_locally_owned_active_cells);
555 
556  // helper function to determine if cell is locally relevant
557  // on active level
558  const auto is_locally_relevant_on_active_level =
559  [&](const auto &cell) {
560  if (cell->is_active())
561  for (const auto v : cell->vertex_indices())
562  if (vertices_owned_by_locally_owned_active_cells
563  [cell->vertex_index(v)])
564  return true;
565  return false;
566  };
567 
568  for (unsigned int level = 0;
570  ++level)
571  {
572  // collect local vertices on level
573  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
574  tria.n_vertices());
575  for (const auto &cell : tria.cell_iterators_on_level(level))
576  if ((construct_multigrid &&
577  (level_subdomain_id_function(cell) == my_rank)) ||
578  (cell->is_active() &&
579  subdomain_id_function(cell) == my_rank))
580  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
581  cell, vertices_owned_by_locally_owned_cells_on_level);
582 
583  // helper function to determine if cell is locally relevant
584  // on level
585  const auto is_locally_relevant_on_level = [&](const auto &cell) {
586  for (const auto v : cell->vertex_indices())
587  if (vertices_owned_by_locally_owned_cells_on_level
588  [cell->vertex_index(v)])
589  return true;
590  return false;
591  };
592 
593  auto &level_cell_infos = construction_data.cell_infos[level];
594  for (const auto &cell : tria.cell_iterators_on_level(level))
595  {
596  // check if cell is locally relevant
597  if (!cell_marked[cell->level()][cell->index()])
598  continue;
599 
600  CellData<dim> cell_info;
601 
602  // save coarse-cell id
603  cell_info.id = cell->id().template to_binary<dim>();
604 
605  // save boundary_ids of each face of this cell
606  for (const auto f : cell->face_indices())
607  {
608  types::boundary_id boundary_ind =
609  cell->face(f)->boundary_id();
610  if (boundary_ind != numbers::internal_face_boundary_id)
611  cell_info.boundary_ids.emplace_back(f, boundary_ind);
612  }
613 
614  // save manifold id
615  {
616  // ... of cell
617  cell_info.manifold_id = cell->manifold_id();
618 
619  // ... of lines
620  if (dim >= 2)
621  for (const auto line : cell->line_indices())
622  cell_info.manifold_line_ids[line] =
623  cell->line(line)->manifold_id();
624 
625  // ... of quads
626  if (dim == 3)
627  for (const auto f : cell->face_indices())
628  cell_info.manifold_quad_ids[f] =
629  cell->quad(f)->manifold_id();
630  }
631 
632  // subdomain and level subdomain id
633  cell_info.subdomain_id = numbers::artificial_subdomain_id;
634  cell_info.level_subdomain_id =
636 
637  if (is_locally_relevant_on_active_level(cell))
638  {
639  cell_info.subdomain_id = subdomain_id_function(cell);
640 
641  cell_info.level_subdomain_id =
643  }
644  else if (is_locally_relevant_on_level(cell))
645  {
646  cell_info.level_subdomain_id =
648  }
649  else
650  {
651  // cell is locally relevant but an artificial cell
652  }
653 
654  level_cell_infos.emplace_back(cell_info);
655  }
656  }
657 
658  return construction_data;
659  }
660 
661  private:
662  void
663  set_additional_data(Description<dim, spacedim> &construction_data) const
664  {
665  construction_data.comm = comm;
666  construction_data.smoothing = tria.get_mesh_smoothing();
667  construction_data.settings = settings;
668  }
669 
670  void
671  set_additional_data(DescriptionTemp<dim, spacedim> &) const
672  {
673  // nothing to do
674  }
675 
676  void
677  add_vertices(Description<dim, spacedim> &construction_data,
678  const std::vector<bool> &vertices_locally_relevant) const
679  {
680  std::vector<unsigned int> vertices_locally_relevant_indices(
681  vertices_locally_relevant.size());
682 
683  // enumerate locally relevant vertices
684  unsigned int vertex_counter = 0;
685  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
686  if (vertices_locally_relevant[i])
687  {
688  construction_data.coarse_cell_vertices.push_back(
689  tria.get_vertices()[i]);
690  vertices_locally_relevant_indices[i] = vertex_counter++;
691  }
692 
693  // correct vertices of cells (make them local)
694  for (auto &cell : construction_data.coarse_cells)
695  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
696  cell.vertices[v] =
697  vertices_locally_relevant_indices[cell.vertices[v]];
698  }
699 
700  void
701  add_vertices(DescriptionTemp<dim, spacedim> &construction_data,
702  const std::vector<bool> &vertices_locally_relevant) const
703  {
704  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
705  if (vertices_locally_relevant[i])
706  construction_data.coarse_cell_vertices.emplace_back(
707  i, tria.get_vertices()[i]);
708  }
709 
710 
711  const ::Triangulation<dim, spacedim> &tria;
712  const std::function<types::subdomain_id(
713  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
715  const std::function<types::subdomain_id(
716  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
718 
719  const MPI_Comm & comm;
722 
723  std::map<unsigned int, std::vector<unsigned int>>
725  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
726  };
727 
728  } // namespace
729 
730 
731  template <int dim, int spacedim>
732  Description<dim, spacedim>
734  const ::Triangulation<dim, spacedim> &tria,
735  const MPI_Comm & comm,
737  const unsigned int my_rank_in)
738  {
739  if (const auto tria_pdt = dynamic_cast<
741  Assert(comm == tria_pdt->get_communicator(),
742  ExcMessage("MPI communicators do not match."));
743 
744  // First, figure out for what rank we are supposed to build the
745  // TriangulationDescription::Description object
746  unsigned int my_rank = my_rank_in;
749  ExcMessage("Rank has to be smaller than available processes."));
750 
751  if (auto tria_pdt = dynamic_cast<
753  {
756  ExcMessage(
757  "For creation from a parallel::distributed::Triangulation, "
758  "my_rank has to equal global rank."));
759 
761  }
762  else if (auto tria_serial =
763  dynamic_cast<const ::Triangulation<dim, spacedim> *>(
764  &tria))
765  {
766  if (my_rank == numbers::invalid_unsigned_int)
768  }
769  else
770  {
771  Assert(false,
772  ExcMessage("This type of triangulation is not supported!"));
773  }
774 
775  const auto subdomain_id_function = [](const auto &cell) {
776  return cell->subdomain_id();
777  };
778 
779  const auto level_subdomain_id_function = [](const auto &cell) {
780  return cell->level_subdomain_id();
781  };
782 
783  return CreateDescriptionFromTriangulationHelper<dim, spacedim>(
784  tria,
787  comm,
788  settings)
789  .template create_description_for_rank<Description<dim, spacedim>>(
790  my_rank);
791  }
792 
793 
794 
795  template <int dim, int spacedim>
798  const std::function<void(::Triangulation<dim, spacedim> &)>
799  & serial_grid_generator,
800  const std::function<void(::Triangulation<dim, spacedim> &,
801  const MPI_Comm &,
802  const unsigned int)> &serial_grid_partitioner,
803  const MPI_Comm & comm,
804  const int group_size,
805  const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
807  {
808 #ifndef DEAL_II_WITH_MPI
809  (void)serial_grid_generator;
810  (void)serial_grid_partitioner;
811  (void)comm;
812  (void)group_size;
813  (void)smoothing;
814  (void)settings;
815 
817 #else
818  const unsigned int my_rank =
820  const unsigned int group_root = (my_rank / group_size) * group_size;
821 
822  const int mpi_tag =
824 
825  // check if process is root of the group
826  if (my_rank == group_root)
827  {
828  // Step 1: create serial triangulation
832  static_cast<
833  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
834  smoothing |
835  Triangulation<dim,
836  spacedim>::limit_level_difference_at_vertices) :
837  smoothing);
838  serial_grid_generator(tria);
839 
840  // Step 2: partition active cells and ...
841  serial_grid_partitioner(tria, comm, group_size);
842 
843  // ... cells on the levels if multigrid is required
844  if (settings &
847 
848  const unsigned int end_group =
849  std::min(group_root + group_size,
851 
852  // 3) create Description for the other processes in group; since
853  // we expect that this function is called for huge meshes, one
854  // Description is created at a time and send away; only once the
855  // Description has been sent away, the next rank is processed.
856  for (unsigned int other_rank = group_root + 1; other_rank < end_group;
857  ++other_rank)
858  {
859  // 3a) create construction data for other ranks
860  const auto construction_data =
862  comm,
863  settings,
864  other_rank);
865  // 3b) pack
866  std::vector<char> buffer;
867  ::Utilities::pack(construction_data, buffer, false);
868 
869  // 3c) send TriangulationDescription::Description
870  const auto ierr = MPI_Send(buffer.data(),
871  buffer.size(),
872  MPI_CHAR,
873  other_rank,
874  mpi_tag,
875  comm);
876  AssertThrowMPI(ierr);
877  }
878 
879  // 4) create TriangulationDescription::Description for this process
880  // (root of group)
882  comm,
883  settings,
884  my_rank);
885  }
886  else
887  {
888  // 3a) recv packed TriangulationDescription::Description from
889  // group-root process
890  // (counter-part of 3c of root process)
891  MPI_Status status;
892  auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
893  AssertThrowMPI(ierr);
894 
895  int len;
896  MPI_Get_count(&status, MPI_CHAR, &len);
897 
898  std::vector<char> buf(len);
899  ierr = MPI_Recv(buf.data(),
900  len,
901  MPI_CHAR,
902  status.MPI_SOURCE,
903  mpi_tag,
904  comm,
905  &status);
906  AssertThrowMPI(ierr);
907 
908  // 3b) unpack TriangulationDescription::Description (counter-part of
909  // 3b of root process)
910  auto construction_data =
911  ::Utilities::template unpack<Description<dim, spacedim>>(
912  buf, false);
913 
914  // WARNING: serialization cannot handle the MPI communicator
915  // which is the reason why we have to set it here explicitly
916  construction_data.comm = comm;
917 
918  return construction_data;
919  }
920 #endif
921  }
922 
923 
924 
925  template <int dim, int spacedim>
931  {
932  const bool construct_multigrid =
933  (partition.size() > 0) &&
934  (settings &
936 
937  Assert(
938  construct_multigrid == false ||
941  ExcMessage(
942  "Source triangulation has to be set up with "
943  "limit_level_difference_at_vertices if the construction of the "
944  "multigrid hierarchy is requested!"));
945 
946  std::vector<LinearAlgebra::distributed::Vector<double>> partitions_mg;
947 
948  if (construct_multigrid) // perform first child policy
949  {
950  const auto tria_parallel =
951  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
952  &tria);
953 
954  Assert(tria_parallel, ExcInternalError());
955 
956  partition.update_ghost_values();
957 
958  partitions_mg.resize(tria.n_global_levels());
959 
960  for (unsigned int l = 0; l < tria.n_global_levels(); ++l)
961  partitions_mg[l].reinit(
962  tria_parallel->global_level_cell_index_partitioner(l).lock());
963 
964  for (int level = tria.n_global_levels() - 1; level >= 0; --level)
965  {
966  for (const auto &cell : tria.cell_iterators_on_level(level))
967  {
968  if (cell->is_locally_owned_on_level() == false)
969  continue;
970 
971  if (cell->is_active())
972  partitions_mg[level][cell->global_level_cell_index()] =
973  partition[cell->global_active_cell_index()];
974  else
975  partitions_mg[level][cell->global_level_cell_index()] =
976  partitions_mg[level + 1]
977  [cell->child(0)->global_level_cell_index()];
978  }
979 
980  partitions_mg[level].update_ghost_values();
981  }
982  }
983 
985  partition,
986  partitions_mg,
987  settings);
988  }
989 
990 
991 
992  template <int dim, int spacedim>
998  & partitions_mg,
999  const TriangulationDescription::Settings settings_in)
1000  {
1001 #ifdef DEAL_II_WITH_MPI
1002  if (tria.get_communicator() == MPI_COMM_NULL)
1003  AssertDimension(partition.local_size(), 0);
1004 #endif
1005 
1006  if (partition.size() == 0)
1007  {
1008  AssertDimension(partitions_mg.size(), 0);
1011  settings_in);
1012  }
1013 
1014  partition.update_ghost_values();
1015 
1016  for (const auto &partition : partitions_mg)
1017  partition.update_ghost_values();
1018 
1019  // 1) determine processes owning locally owned cells
1020  const std::vector<unsigned int> relevant_processes = [&]() {
1021  std::set<unsigned int> relevant_processes;
1022 
1023  for (unsigned int i = 0; i < partition.local_size(); ++i)
1024  relevant_processes.insert(
1025  static_cast<unsigned int>(partition.local_element(i)));
1026 
1027  for (const auto &partition : partitions_mg)
1028  for (unsigned int i = 0; i < partition.local_size(); ++i)
1029  relevant_processes.insert(
1030  static_cast<unsigned int>(partition.local_element(i)));
1031 
1032  return std::vector<unsigned int>(relevant_processes.begin(),
1033  relevant_processes.end());
1034  }();
1035 
1036  const bool construct_multigrid = partitions_mg.size() > 0;
1037 
1039 
1040  if (construct_multigrid)
1042  settings |
1044 
1045  const auto subdomain_id_function = [&partition](const auto &cell) {
1046  if ((cell->is_active() && (cell->is_artificial() == false)))
1047  return static_cast<unsigned int>(
1048  partition[cell->global_active_cell_index()]);
1049  else
1051  };
1052 
1053  const auto level_subdomain_id_function =
1054  [&construct_multigrid, &partitions_mg](const auto &cell) {
1055  if (construct_multigrid && (cell->is_artificial_on_level() == false))
1056  return static_cast<unsigned int>(
1057  partitions_mg[cell->level()][cell->global_level_cell_index()]);
1058  else
1060  };
1061 
1062  CreateDescriptionFromTriangulationHelper<dim, spacedim> helper(
1063  tria,
1067  settings);
1068 
1069  // create a description (locally owned cell and a layer of ghost cells
1070  // and all their parents)
1071  std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
1072  descriptions_per_rank.reserve(relevant_processes.size());
1073 
1074  for (const auto rank : relevant_processes)
1075  descriptions_per_rank.emplace_back(
1076  helper.template create_description_for_rank<
1077  DescriptionTemp<dim, spacedim>>(rank));
1078 
1079  // collect description from all processes that used to own locally-owned
1080  // active cells of this process in a single description
1081  DescriptionTemp<dim, spacedim> description_merged;
1082  description_merged.collect(
1083  relevant_processes,
1084  descriptions_per_rank,
1085  partition.get_mpi_communicator(),
1086  dynamic_cast<
1088  &tria) == nullptr);
1089 
1090  // remove redundant entries
1091  description_merged.reduce();
1092 
1093  // convert to actual description
1094  return description_merged.convert(partition.get_mpi_communicator(),
1096  settings);
1097  }
1098 
1099  } // namespace Utilities
1100 } // namespace TriangulationDescription
1101 
1102 
1103 
1104 /*-------------- Explicit Instantiations -------------------------------*/
1105 #include "tria_description.inst"
1106 
1107 
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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4605
unsigned int level
Definition: grid_out.cc:4607
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:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
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:6398
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4178
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_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)
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:156
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:145
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1483
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:260
const types::subdomain_id artificial_subdomain_id
Definition: types.h:298
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::manifold_id flat_manifold_id
Definition: types.h:269
unsigned int manifold_id
Definition: types.h:141
global_cell_index coarse_cell_id
Definition: types.h:114
unsigned int subdomain_id
Definition: types.h:43
unsigned int material_id
Definition: types.h:152
unsigned int boundary_id
Definition: types.h:129
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