Reference documentation for deal.II version GIT 0b65fff18a 2023-09-27 19:30: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 - 2023 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.empty()) && (boundary_quads.empty()));
75  case 2:
76  return (boundary_quads.empty());
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 std::pair<unsigned int, Point<spacedim>> &a,
272  const std::pair<unsigned int, Point<spacedim>> &b) {
273  return a.first < b.first;
274  });
275  this->coarse_cell_vertices.erase(
276  std::unique(
277  this->coarse_cell_vertices.begin(),
278  this->coarse_cell_vertices.end(),
279  [](const std::pair<unsigned int, Point<spacedim>> &a,
280  const std::pair<unsigned int, Point<spacedim>> &b) {
281  if (a.first == b.first)
282  {
283  Assert(a.second.distance(b.second) <=
284  1e-7 *
285  std::max(a.second.norm(), b.second.norm()),
286  ExcInternalError());
287  return true;
288  }
289  return false;
290  }),
291  this->coarse_cell_vertices.end());
292  }
293 
294  // make cells unique
295  for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
296  {
297  if (this->cell_infos[i].empty())
298  continue;
299 
300  std::sort(this->cell_infos[i].begin(),
301  this->cell_infos[i].end(),
302  [](const auto &a, const auto &b) {
303  return a.id < b.id;
304  });
305 
306  std::vector<CellData<dim>> temp;
307  temp.push_back(this->cell_infos[i][0]);
308 
309  for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j)
310  if (temp.back().id == cell_infos[i][j].id)
311  {
312  temp.back().subdomain_id =
313  std::min(temp.back().subdomain_id,
314  this->cell_infos[i][j].subdomain_id);
315  temp.back().level_subdomain_id =
316  std::min(temp.back().level_subdomain_id,
317  this->cell_infos[i][j].level_subdomain_id);
318  }
319  else
320  {
321  temp.push_back(this->cell_infos[i][j]);
322  }
323 
324  this->cell_infos[i] = temp;
325  }
326  }
327 
333  Description<dim, spacedim>
334  convert(const MPI_Comm comm,
336  mesh_smoothing,
338  {
339  Description<dim, spacedim> description;
340 
341  // copy communicator
342  description.comm = comm;
343 
344  description.settings = settings;
345 
346  // use mesh smoothing from base triangulation
347  description.smoothing = mesh_smoothing;
348 
349  std::map<unsigned int, unsigned int> map;
350 
351  for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
352  {
353  description.coarse_cell_vertices.push_back(
354  this->coarse_cell_vertices[i].second);
355  map[this->coarse_cell_vertices[i].first] = i;
356  }
357 
358  description.coarse_cells = this->coarse_cells;
359 
360  for (auto &cell : description.coarse_cells)
361  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
362  cell.vertices[v] = map[cell.vertices[v]];
363 
364  description.coarse_cell_index_to_coarse_cell_id =
365  this->coarse_cell_index_to_coarse_cell_id;
366  description.cell_infos = this->cell_infos;
367 
368  return description;
369  }
370 
371  std::vector<::CellData<dim>> coarse_cells;
372 
373  std::vector<std::pair<unsigned int, Point<spacedim>>>
375 
376  std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
377 
378  std::vector<std::vector<CellData<dim>>> cell_infos;
379  };
380 
384  template <int dim, int spacedim>
385  void
386  mark_cell_and_its_parents(
388  std::vector<std::vector<bool>> &cell_marked)
389  {
390  cell_marked[cell->level()][cell->index()] = true;
391  if (cell->level() != 0)
392  mark_cell_and_its_parents(cell->parent(), cell_marked);
393  }
394 
400  template <int dim, int spacedim>
401  class CreateDescriptionFromTriangulationHelper
402  {
403  public:
404  CreateDescriptionFromTriangulationHelper(
405  const ::Triangulation<dim, spacedim> &tria,
406  const std::function<types::subdomain_id(
407  const typename ::Triangulation<dim, spacedim>::cell_iterator
409  const std::function<types::subdomain_id(
410  const typename ::Triangulation<dim, spacedim>::cell_iterator
412  const MPI_Comm comm,
414  : tria(tria)
417  , comm(comm)
418  , settings(settings)
421  0u)
422  {
423  Assert(
427  Triangulation<dim,
428  spacedim>::limit_level_difference_at_vertices),
429  ExcMessage(
430  "Source triangulation has to be set up with "
431  "limit_level_difference_at_vertices if the construction of the "
432  "multigrid hierarchy is requested!"));
433 
436  }
437 
438  template <typename T>
439  T
440  create_description_for_rank(const unsigned int my_rank) const
441  {
442  T construction_data;
443 
444  set_additional_data(construction_data);
445 
446  // helper function, which collects all vertices belonging to a cell
447  // (also taking into account periodicity)
448  const auto
449  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
450  [this](const auto &cell,
451  std::vector<bool> &vertices_owned_by_locally_owned_cells) {
452  // add local vertices
453  for (const auto v : cell->vertex_indices())
454  {
455  vertices_owned_by_locally_owned_cells[cell->vertex_index(
456  v)] = true;
457  const auto coinciding_vertex_group =
459  cell->vertex_index(v));
460  if (coinciding_vertex_group !=
462  for (const auto &co_vertex : coinciding_vertex_groups.at(
463  coinciding_vertex_group->second))
464  vertices_owned_by_locally_owned_cells[co_vertex] = true;
465  }
466  };
467 
468  // 1) loop over levels (from fine to coarse) and mark on each level
469  // the locally relevant cells
470  std::vector<std::vector<bool>> cell_marked(tria.n_levels());
471  for (unsigned int l = 0; l < tria.n_levels(); ++l)
472  cell_marked[l].resize(tria.n_raw_cells(l));
473 
474  for (int level = tria.get_triangulation().n_global_levels() - 1;
475  level >= 0;
476  --level)
477  {
478  // collect vertices connected to a (on any level) locally owned
479  // cell
480  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
481  tria.n_vertices());
482  for (const auto &cell : tria.cell_iterators_on_level(level))
483  if (construct_multigrid &&
484  (level_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  for (const auto &cell : tria.active_cell_iterators())
489  if (subdomain_id_function(cell) == my_rank)
490  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
491  cell, vertices_owned_by_locally_owned_cells_on_level);
492 
493  // helper function to determine if cell is locally relevant
494  // (i.e. a cell which is connected to a vertex via a locally
495  // owned cell)
496  const auto is_locally_relevant_on_level = [&](const auto &cell) {
497  for (const auto v : cell->vertex_indices())
498  if (vertices_owned_by_locally_owned_cells_on_level
499  [cell->vertex_index(v)])
500  return true;
501  return false;
502  };
503 
504  // mark all locally relevant cells
505  for (const auto &cell : tria.cell_iterators_on_level(level))
506  if (is_locally_relevant_on_level(cell))
507  mark_cell_and_its_parents(cell, cell_marked);
508  }
509 
510  // 2) set_up coarse-grid triangulation
511  {
512  std::vector<bool> vertices_locally_relevant(tria.n_vertices(),
513  false);
514 
515  // a) loop over all cells
516  for (const auto &cell : tria.cell_iterators_on_level(0))
517  {
518  if (!cell_marked[cell->level()][cell->index()])
519  continue;
520 
521  // extract cell definition (with old numbering of vertices)
522  ::CellData<dim> cell_data(cell->n_vertices());
523  cell_data.material_id = cell->material_id();
524  cell_data.manifold_id = cell->manifold_id();
525  for (const auto v : cell->vertex_indices())
526  cell_data.vertices[v] = cell->vertex_index(v);
527  construction_data.coarse_cells.push_back(cell_data);
528 
529  // save indices of each vertex of this cell
530  for (const auto v : cell->vertex_indices())
531  vertices_locally_relevant[cell->vertex_index(v)] = true;
532 
533  // save translation for corase grid: lid -> gid
534  construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
535  cell->id().get_coarse_cell_id());
536  }
537 
538  add_vertices(construction_data, vertices_locally_relevant);
539  }
540 
541 
542  // 3) collect info of each cell
543  construction_data.cell_infos.resize(
545 
546  // collect local vertices on active level
547  std::vector<bool> vertices_owned_by_locally_owned_active_cells(
548  tria.n_vertices());
549  for (const auto &cell : tria.active_cell_iterators())
550  if (subdomain_id_function(cell) == my_rank)
551  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
552  cell, vertices_owned_by_locally_owned_active_cells);
553 
554  // helper function to determine if cell is locally relevant
555  // on active level
556  const auto is_locally_relevant_on_active_level =
557  [&](const auto &cell) {
558  if (cell->is_active())
559  for (const auto v : cell->vertex_indices())
560  if (vertices_owned_by_locally_owned_active_cells
561  [cell->vertex_index(v)])
562  return true;
563  return false;
564  };
565 
566  for (unsigned int level = 0;
568  ++level)
569  {
570  // collect local vertices on level
571  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
572  tria.n_vertices());
573  for (const auto &cell : tria.cell_iterators_on_level(level))
574  if ((construct_multigrid &&
575  (level_subdomain_id_function(cell) == my_rank)) ||
576  (cell->is_active() &&
577  subdomain_id_function(cell) == my_rank))
578  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
579  cell, vertices_owned_by_locally_owned_cells_on_level);
580 
581  // helper function to determine if cell is locally relevant
582  // on level
583  const auto is_locally_relevant_on_level = [&](const auto &cell) {
584  for (const auto v : cell->vertex_indices())
585  if (vertices_owned_by_locally_owned_cells_on_level
586  [cell->vertex_index(v)])
587  return true;
588  return false;
589  };
590 
591  auto &level_cell_infos = construction_data.cell_infos[level];
592  for (const auto &cell : tria.cell_iterators_on_level(level))
593  {
594  // check if cell is locally relevant
595  if (!cell_marked[cell->level()][cell->index()])
596  continue;
597 
598  CellData<dim> cell_info;
599 
600  // save coarse-cell id
601  cell_info.id = cell->id().template to_binary<dim>();
602 
603  // save boundary_ids of each face of this cell
604  for (const auto f : cell->face_indices())
605  {
606  types::boundary_id boundary_ind =
607  cell->face(f)->boundary_id();
608  if (boundary_ind != numbers::internal_face_boundary_id)
609  cell_info.boundary_ids.emplace_back(f, boundary_ind);
610  }
611 
612  // save manifold id
613  {
614  // ... of cell
615  cell_info.manifold_id = cell->manifold_id();
616 
617  // ... of lines
618  if (dim >= 2)
619  for (const auto line : cell->line_indices())
620  cell_info.manifold_line_ids[line] =
621  cell->line(line)->manifold_id();
622 
623  // ... of quads
624  if (dim == 3)
625  for (const auto f : cell->face_indices())
626  cell_info.manifold_quad_ids[f] =
627  cell->quad(f)->manifold_id();
628  }
629 
630  // subdomain and level subdomain id
631  cell_info.subdomain_id = numbers::artificial_subdomain_id;
632  cell_info.level_subdomain_id =
634 
635  if (is_locally_relevant_on_active_level(cell))
636  {
637  cell_info.subdomain_id = subdomain_id_function(cell);
638 
639  cell_info.level_subdomain_id =
641  }
642  else if (is_locally_relevant_on_level(cell))
643  {
644  cell_info.level_subdomain_id =
646  }
647  else
648  {
649  // cell is locally relevant but an artificial cell
650  }
651 
652  level_cell_infos.emplace_back(cell_info);
653  }
654  }
655 
656  return construction_data;
657  }
658 
659  private:
660  void
661  set_additional_data(Description<dim, spacedim> &construction_data) const
662  {
663  construction_data.comm = comm;
664  construction_data.smoothing = tria.get_mesh_smoothing();
665  construction_data.settings = settings;
666  }
667 
668  void
669  set_additional_data(DescriptionTemp<dim, spacedim> &) const
670  {
671  // nothing to do
672  }
673 
674  void
675  add_vertices(Description<dim, spacedim> &construction_data,
676  const std::vector<bool> &vertices_locally_relevant) const
677  {
678  std::vector<unsigned int> vertices_locally_relevant_indices(
679  vertices_locally_relevant.size());
680 
681  // enumerate locally relevant vertices
682  unsigned int vertex_counter = 0;
683  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
684  if (vertices_locally_relevant[i])
685  {
686  construction_data.coarse_cell_vertices.push_back(
687  tria.get_vertices()[i]);
688  vertices_locally_relevant_indices[i] = vertex_counter++;
689  }
690 
691  // correct vertices of cells (make them local)
692  for (auto &cell : construction_data.coarse_cells)
693  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
694  cell.vertices[v] =
695  vertices_locally_relevant_indices[cell.vertices[v]];
696  }
697 
698  void
699  add_vertices(DescriptionTemp<dim, spacedim> &construction_data,
700  const std::vector<bool> &vertices_locally_relevant) const
701  {
702  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
703  if (vertices_locally_relevant[i])
704  construction_data.coarse_cell_vertices.emplace_back(
705  i, tria.get_vertices()[i]);
706  }
707 
708 
709  const ::Triangulation<dim, spacedim> &tria;
710  const std::function<types::subdomain_id(
711  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
713  const std::function<types::subdomain_id(
714  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
716 
717  const MPI_Comm comm;
720 
721  std::map<unsigned int, std::vector<unsigned int>>
723  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
724  };
725 
726  } // namespace
727 
728 
729  template <int dim, int spacedim>
730  Description<dim, spacedim>
732  const ::Triangulation<dim, spacedim> &tria,
733  const MPI_Comm comm,
735  const unsigned int my_rank_in)
736  {
737  if (const auto tria_pdt = dynamic_cast<
739  Assert(comm == tria_pdt->get_communicator(),
740  ExcMessage("MPI communicators do not match."));
741 
742  // First, figure out for what rank we are supposed to build the
743  // TriangulationDescription::Description object
744  unsigned int my_rank = my_rank_in;
747  ExcMessage("Rank has to be smaller than available processes."));
748 
749  if (auto tria_pdt = dynamic_cast<
751  {
754  ExcMessage(
755  "For creation from a parallel::distributed::Triangulation, "
756  "my_rank has to equal global rank."));
757 
759  }
760  else if (auto tria_serial =
761  dynamic_cast<const ::Triangulation<dim, spacedim> *>(
762  &tria))
763  {
764  if (my_rank == numbers::invalid_unsigned_int)
766  }
767  else
768  {
769  Assert(false,
770  ExcMessage("This type of triangulation is not supported!"));
771  }
772 
773  const auto subdomain_id_function = [](const auto &cell) {
774  return cell->subdomain_id();
775  };
776 
777  const auto level_subdomain_id_function = [](const auto &cell) {
778  return cell->level_subdomain_id();
779  };
780 
781  return CreateDescriptionFromTriangulationHelper<dim, spacedim>(
782  tria,
785  comm,
786  settings)
787  .template create_description_for_rank<Description<dim, spacedim>>(
788  my_rank);
789  }
790 
791 
792 
793  template <int dim, int spacedim>
796  const std::function<void(::Triangulation<dim, spacedim> &)>
797  &serial_grid_generator,
798  const std::function<void(::Triangulation<dim, spacedim> &,
799  const MPI_Comm,
800  const unsigned int)> &serial_grid_partitioner,
801  const MPI_Comm comm,
802  const int group_size,
803  const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
805  {
806 #ifndef DEAL_II_WITH_MPI
807  (void)serial_grid_generator;
808  (void)serial_grid_partitioner;
809  (void)comm;
810  (void)group_size;
811  (void)smoothing;
812  (void)settings;
813 
815 #else
816  const unsigned int my_rank =
818  const unsigned int group_root = (my_rank / group_size) * group_size;
819 
820  const int mpi_tag =
822 
823  // check if process is root of the group
824  if (my_rank == group_root)
825  {
826  // Step 1: create serial triangulation
830  static_cast<
831  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
832  smoothing |
833  Triangulation<dim,
834  spacedim>::limit_level_difference_at_vertices) :
835  smoothing);
836  serial_grid_generator(tria);
837 
838  // Step 2: partition active cells and ...
839  serial_grid_partitioner(tria, comm, group_size);
840 
841  // ... cells on the levels if multigrid is required
842  if (settings &
845 
846  const unsigned int end_group =
847  std::min(group_root + group_size,
849 
850  // 3) create Description for the other processes in group; since
851  // we expect that this function is called for huge meshes, one
852  // Description is created at a time and send away; only once the
853  // Description has been sent away, the next rank is processed.
854  for (unsigned int other_rank = group_root + 1; other_rank < end_group;
855  ++other_rank)
856  {
857  // 3a) create construction data for other ranks
858  const auto construction_data =
860  comm,
861  settings,
862  other_rank);
863  // 3b) pack
864  std::vector<char> buffer;
865  ::Utilities::pack(construction_data, buffer, false);
866 
867  // 3c) send TriangulationDescription::Description
868  const auto ierr = MPI_Send(buffer.data(),
869  buffer.size(),
870  MPI_CHAR,
871  other_rank,
872  mpi_tag,
873  comm);
874  AssertThrowMPI(ierr);
875  }
876 
877  // 4) create TriangulationDescription::Description for this process
878  // (root of group)
880  comm,
881  settings,
882  my_rank);
883  }
884  else
885  {
886  // 3a) recv packed TriangulationDescription::Description from
887  // group-root process
888  // (counter-part of 3c of root process)
889  MPI_Status status;
890  auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
891  AssertThrowMPI(ierr);
892 
893  int len;
894  MPI_Get_count(&status, MPI_CHAR, &len);
895 
896  std::vector<char> buf(len);
897  ierr = MPI_Recv(buf.data(),
898  len,
899  MPI_CHAR,
900  status.MPI_SOURCE,
901  mpi_tag,
902  comm,
903  &status);
904  AssertThrowMPI(ierr);
905 
906  // 3b) unpack TriangulationDescription::Description (counter-part of
907  // 3b of root process)
908  auto construction_data =
909  ::Utilities::template unpack<Description<dim, spacedim>>(
910  buf, false);
911 
912  // WARNING: serialization cannot handle the MPI communicator
913  // which is the reason why we have to set it here explicitly
914  construction_data.comm = comm;
915 
916  return construction_data;
917  }
918 #endif
919  }
920 
921 
922 
923  template <int dim, int spacedim>
929  {
930  const bool construct_multigrid =
931  (partition.size() > 0) &&
932  (settings &
934 
935  Assert(
936  construct_multigrid == false ||
939  ExcMessage(
940  "Source triangulation has to be set up with "
941  "limit_level_difference_at_vertices if the construction of the "
942  "multigrid hierarchy is requested!"));
943 
944  std::vector<LinearAlgebra::distributed::Vector<double>> partitions_mg;
945 
946  if (construct_multigrid) // perform first child policy
947  {
948  const auto tria_parallel =
949  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
950  &tria);
951 
952  Assert(tria_parallel, ExcInternalError());
953 
954  partition.update_ghost_values();
955 
956  partitions_mg.resize(tria.n_global_levels());
957 
958  for (unsigned int l = 0; l < tria.n_global_levels(); ++l)
959  partitions_mg[l].reinit(
960  tria_parallel->global_level_cell_index_partitioner(l).lock());
961 
962  for (int level = tria.n_global_levels() - 1; level >= 0; --level)
963  {
964  for (const auto &cell : tria.cell_iterators_on_level(level))
965  {
966  if (cell->is_locally_owned_on_level() == false)
967  continue;
968 
969  if (cell->is_active())
970  partitions_mg[level][cell->global_level_cell_index()] =
971  partition[cell->global_active_cell_index()];
972  else
973  partitions_mg[level][cell->global_level_cell_index()] =
974  partitions_mg[level + 1]
975  [cell->child(0)->global_level_cell_index()];
976  }
977 
978  partitions_mg[level].update_ghost_values();
979  }
980  }
981 
983  partition,
984  partitions_mg,
985  settings);
986  }
987 
988 
989 
990  template <int dim, int spacedim>
996  &partitions_mg,
997  const TriangulationDescription::Settings settings_in)
998  {
999 #ifdef DEAL_II_WITH_MPI
1000  if (tria.get_communicator() == MPI_COMM_NULL)
1001  AssertDimension(partition.locally_owned_size(), 0);
1002 #endif
1003 
1004  if (partition.size() == 0)
1005  {
1006  AssertDimension(partitions_mg.size(), 0);
1009  settings_in);
1010  }
1011 
1012  partition.update_ghost_values();
1013 
1014  for (const auto &partition : partitions_mg)
1015  partition.update_ghost_values();
1016 
1017  // 1) determine processes owning locally owned cells
1018  const std::vector<unsigned int> relevant_processes = [&]() {
1019  std::set<unsigned int> relevant_processes;
1020 
1021  for (unsigned int i = 0; i < partition.locally_owned_size(); ++i)
1022  relevant_processes.insert(
1023  static_cast<unsigned int>(partition.local_element(i)));
1024 
1025  for (const auto &partition : partitions_mg)
1026  for (unsigned int i = 0; i < partition.locally_owned_size(); ++i)
1027  relevant_processes.insert(
1028  static_cast<unsigned int>(partition.local_element(i)));
1029 
1030  return std::vector<unsigned int>(relevant_processes.begin(),
1031  relevant_processes.end());
1032  }();
1033 
1034  const bool construct_multigrid = partitions_mg.size() > 0;
1035 
1037 
1038  if (construct_multigrid)
1040  settings |
1042 
1043  const auto subdomain_id_function = [&partition](const auto &cell) {
1044  if ((cell->is_active() && (cell->is_artificial() == false)))
1045  return static_cast<unsigned int>(
1046  partition[cell->global_active_cell_index()]);
1047  else
1049  };
1050 
1051  const auto level_subdomain_id_function =
1052  [&construct_multigrid, &partitions_mg](const auto &cell) {
1053  if (construct_multigrid && (cell->is_artificial_on_level() == false))
1054  return static_cast<unsigned int>(
1055  partitions_mg[cell->level()][cell->global_level_cell_index()]);
1056  else
1058  };
1059 
1060  CreateDescriptionFromTriangulationHelper<dim, spacedim> helper(
1061  tria,
1065  settings);
1066 
1067  // create a description (locally owned cell and a layer of ghost cells
1068  // and all their parents)
1069  std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
1070  descriptions_per_rank.reserve(relevant_processes.size());
1071 
1072  for (const auto rank : relevant_processes)
1073  descriptions_per_rank.emplace_back(
1074  helper.template create_description_for_rank<
1075  DescriptionTemp<dim, spacedim>>(rank));
1076 
1077  // collect description from all processes that used to own locally-owned
1078  // active cells of this process in a single description
1079  DescriptionTemp<dim, spacedim> description_merged;
1080  description_merged.collect(
1081  relevant_processes,
1082  descriptions_per_rank,
1083  partition.get_mpi_communicator(),
1084  dynamic_cast<
1086  &tria) == nullptr);
1087 
1088  // remove redundant entries
1089  description_merged.reduce();
1090 
1091  // convert to actual description
1092  return description_merged.convert(partition.get_mpi_communicator(),
1094  settings);
1095  }
1096 
1097  } // namespace Utilities
1098 } // namespace TriangulationDescription
1099 
1100 
1101 
1102 /*-------------- Explicit Instantiations -------------------------------*/
1103 #include "tria_description.inst"
1104 
1105 
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4615
unsigned int level
Definition: grid_out.cc:4617
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:1616
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
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:7045
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4235
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 n_mpi_processes(const MPI_Comm mpi_communicator)
Definition: mpi.cc:149
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:164
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)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1343
constexpr DEAL_II_HOST 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:277
const types::subdomain_id artificial_subdomain_id
Definition: types.h:315
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const types::manifold_id flat_manifold_id
Definition: types.h:286
unsigned int manifold_id
Definition: types.h:153
global_cell_index coarse_cell_id
Definition: types.h:127
unsigned int subdomain_id
Definition: types.h:44
unsigned int material_id
Definition: types.h:164
unsigned int boundary_id
Definition: types.h:141
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