Reference documentation for deal.II version GIT 90d86f7ad9 2022-05-28 14:10: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  set_user_flag_and_of_its_parents(
390  {
391  cell->set_user_flag();
392  if (cell->level() != 0)
393  set_user_flag_and_of_its_parents(cell->parent());
394  }
395 
401  template <int dim, int spacedim>
402  class CreateDescriptionFromTriangulationHelper
403  {
404  public:
405  CreateDescriptionFromTriangulationHelper(
406  const ::Triangulation<dim, spacedim> &tria,
407  const std::function<types::subdomain_id(
408  const typename ::Triangulation<dim, spacedim>::cell_iterator
409  &)> & subdomain_id_function,
410  const std::function<types::subdomain_id(
411  const typename ::Triangulation<dim, spacedim>::cell_iterator
413  const MPI_Comm & comm,
415  : tria(tria)
418  , comm(comm)
419  , settings(settings)
422  0u)
423  {
424  Assert(
428  Triangulation<dim,
429  spacedim>::limit_level_difference_at_vertices),
430  ExcMessage(
431  "Source triangulation has to be set up with "
432  "limit_level_difference_at_vertices if the construction of the "
433  "multigrid hierarchy is requested!"));
434 
437  }
438 
439  template <typename T>
440  T
441  create_description_for_rank(const unsigned int my_rank) const
442  {
443  T construction_data;
444 
445  set_additional_data(construction_data);
446 
447  // helper function, which collects all vertices belonging to a cell
448  // (also taking into account periodicity)
449  const auto
450  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
451  [this](const auto & cell,
452  std::vector<bool> &vertices_owned_by_locally_owned_cells) {
453  // add local vertices
454  for (const auto v : cell->vertex_indices())
455  {
456  vertices_owned_by_locally_owned_cells[cell->vertex_index(
457  v)] = true;
458  const auto coinciding_vertex_group =
460  cell->vertex_index(v));
461  if (coinciding_vertex_group !=
463  for (const auto &co_vertex : coinciding_vertex_groups.at(
464  coinciding_vertex_group->second))
465  vertices_owned_by_locally_owned_cells[co_vertex] = true;
466  }
467  };
468 
469  // 1) collect locally relevant cells (set user_flag)
470  std::vector<bool> old_user_flags;
471  tria.save_user_flags(old_user_flags);
472 
473  // 1a) clear user_flags
474  const_cast<::Triangulation<dim, spacedim> &>(tria)
475  .clear_user_flags();
476 
477  // 1b) loop over levels (from fine to coarse) and mark on each level
478  // the locally relevant cells
479  for (int level = tria.get_triangulation().n_global_levels() - 1;
480  level >= 0;
481  --level)
482  {
483  // collect vertices connected to a (on any level) locally owned
484  // cell
485  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
486  tria.n_vertices());
487  for (const auto &cell : tria.cell_iterators_on_level(level))
488  if (construct_multigrid &&
489  (level_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  for (const auto &cell : tria.active_cell_iterators())
494  if (subdomain_id_function(cell) == my_rank)
495  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
496  cell, vertices_owned_by_locally_owned_cells_on_level);
497 
498  // helper function to determine if cell is locally relevant
499  // (i.e. a cell which is connected to a vertex via a locally
500  // owned cell)
501  const auto is_locally_relevant_on_level = [&](const auto &cell) {
502  for (const auto v : cell->vertex_indices())
503  if (vertices_owned_by_locally_owned_cells_on_level
504  [cell->vertex_index(v)])
505  return true;
506  return false;
507  };
508 
509  // mark all locally relevant cells
510  for (const auto &cell : tria.cell_iterators_on_level(level))
511  if (is_locally_relevant_on_level(cell))
512  set_user_flag_and_of_its_parents(cell);
513  }
514 
515  // 2) set_up coarse-grid triangulation
516  {
517  std::vector<bool> vertices_locally_relevant(tria.n_vertices(),
518  false);
519 
520  // a) loop over all cells
521  for (const auto &cell : tria.cell_iterators_on_level(0))
522  {
523  if (!cell->user_flag_set())
524  continue;
525 
526  // extract cell definition (with old numbering of vertices)
527  ::CellData<dim> cell_data(cell->n_vertices());
528  cell_data.material_id = cell->material_id();
529  cell_data.manifold_id = cell->manifold_id();
530  for (const auto v : cell->vertex_indices())
531  cell_data.vertices[v] = cell->vertex_index(v);
532  construction_data.coarse_cells.push_back(cell_data);
533 
534  // save indices of each vertex of this cell
535  for (const auto v : cell->vertex_indices())
536  vertices_locally_relevant[cell->vertex_index(v)] = true;
537 
538  // save translation for corase grid: lid -> gid
539  construction_data.coarse_cell_index_to_coarse_cell_id.push_back(
540  cell->id().get_coarse_cell_id());
541  }
542 
543  add_vertices(construction_data, vertices_locally_relevant);
544  }
545 
546 
547  // 3) collect info of each cell
548  construction_data.cell_infos.resize(
550 
551  // collect local vertices on active level
552  std::vector<bool> vertices_owned_by_locally_owned_active_cells(
553  tria.n_vertices());
554  for (const auto &cell : tria.active_cell_iterators())
555  if (subdomain_id_function(cell) == my_rank)
556  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
557  cell, vertices_owned_by_locally_owned_active_cells);
558 
559  // helper function to determine if cell is locally relevant
560  // on active level
561  const auto is_locally_relevant_on_active_level =
562  [&](const auto &cell) {
563  if (cell->is_active())
564  for (const auto v : cell->vertex_indices())
565  if (vertices_owned_by_locally_owned_active_cells
566  [cell->vertex_index(v)])
567  return true;
568  return false;
569  };
570 
571  for (unsigned int level = 0;
573  ++level)
574  {
575  // collect local vertices on level
576  std::vector<bool> vertices_owned_by_locally_owned_cells_on_level(
577  tria.n_vertices());
578  for (const auto &cell : tria.cell_iterators_on_level(level))
579  if ((construct_multigrid &&
580  (level_subdomain_id_function(cell) == my_rank)) ||
581  (cell->is_active() &&
582  subdomain_id_function(cell) == my_rank))
583  add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
584  cell, vertices_owned_by_locally_owned_cells_on_level);
585 
586  // helper function to determine if cell is locally relevant
587  // on level
588  const auto is_locally_relevant_on_level = [&](const auto &cell) {
589  for (const auto v : cell->vertex_indices())
590  if (vertices_owned_by_locally_owned_cells_on_level
591  [cell->vertex_index(v)])
592  return true;
593  return false;
594  };
595 
596  auto &level_cell_infos = construction_data.cell_infos[level];
597  for (const auto &cell : tria.cell_iterators_on_level(level))
598  {
599  // check if cell is locally relevant
600  if (!(cell->user_flag_set()))
601  continue;
602 
603  CellData<dim> cell_info;
604 
605  // save coarse-cell id
606  cell_info.id = cell->id().template to_binary<dim>();
607 
608  // save boundary_ids of each face of this cell
609  for (const auto f : cell->face_indices())
610  {
611  types::boundary_id boundary_ind =
612  cell->face(f)->boundary_id();
613  if (boundary_ind != numbers::internal_face_boundary_id)
614  cell_info.boundary_ids.emplace_back(f, boundary_ind);
615  }
616 
617  // save manifold id
618  {
619  // ... of cell
620  cell_info.manifold_id = cell->manifold_id();
621 
622  // ... of lines
623  if (dim >= 2)
624  for (const auto line : cell->line_indices())
625  cell_info.manifold_line_ids[line] =
626  cell->line(line)->manifold_id();
627 
628  // ... of quads
629  if (dim == 3)
630  for (const auto f : cell->face_indices())
631  cell_info.manifold_quad_ids[f] =
632  cell->quad(f)->manifold_id();
633  }
634 
635  // subdomain and level subdomain id
636  cell_info.subdomain_id = numbers::artificial_subdomain_id;
637  cell_info.level_subdomain_id =
639 
640  if (is_locally_relevant_on_active_level(cell))
641  {
642  cell_info.subdomain_id = subdomain_id_function(cell);
643 
644  cell_info.level_subdomain_id =
646  }
647  else if (is_locally_relevant_on_level(cell))
648  {
649  cell_info.level_subdomain_id =
651  }
652  else
653  {
654  // cell is locally relevant but an artificial cell
655  }
656 
657  level_cell_infos.emplace_back(cell_info);
658  }
659  }
660 
661  const_cast<::Triangulation<dim, spacedim> &>(tria)
662  .load_user_flags(old_user_flags);
663 
664  return construction_data;
665  }
666 
667  private:
668  void
669  set_additional_data(Description<dim, spacedim> &construction_data) const
670  {
671  construction_data.comm = comm;
672  construction_data.smoothing = tria.get_mesh_smoothing();
673  construction_data.settings = settings;
674  }
675 
676  void
677  set_additional_data(DescriptionTemp<dim, spacedim> &) const
678  {
679  // nothing to do
680  }
681 
682  void
683  add_vertices(Description<dim, spacedim> &construction_data,
684  const std::vector<bool> &vertices_locally_relevant) const
685  {
686  std::vector<unsigned int> vertices_locally_relevant_indices(
687  vertices_locally_relevant.size());
688 
689  // enumerate locally relevant vertices
690  unsigned int vertex_counter = 0;
691  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
692  if (vertices_locally_relevant[i])
693  {
694  construction_data.coarse_cell_vertices.push_back(
695  tria.get_vertices()[i]);
696  vertices_locally_relevant_indices[i] = vertex_counter++;
697  }
698 
699  // correct vertices of cells (make them local)
700  for (auto &cell : construction_data.coarse_cells)
701  for (unsigned int v = 0; v < cell.vertices.size(); ++v)
702  cell.vertices[v] =
703  vertices_locally_relevant_indices[cell.vertices[v]];
704  }
705 
706  void
707  add_vertices(DescriptionTemp<dim, spacedim> &construction_data,
708  const std::vector<bool> &vertices_locally_relevant) const
709  {
710  for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
711  if (vertices_locally_relevant[i])
712  construction_data.coarse_cell_vertices.emplace_back(
713  i, tria.get_vertices()[i]);
714  }
715 
716 
717  const ::Triangulation<dim, spacedim> &tria;
718  const std::function<types::subdomain_id(
719  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
721  const std::function<types::subdomain_id(
722  const typename ::Triangulation<dim, spacedim>::cell_iterator &)>
724 
725  const MPI_Comm & comm;
728 
729  std::map<unsigned int, std::vector<unsigned int>>
731  std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
732  };
733 
734  } // namespace
735 
736 
737  template <int dim, int spacedim>
738  Description<dim, spacedim>
740  const ::Triangulation<dim, spacedim> &tria,
741  const MPI_Comm & comm,
743  const unsigned int my_rank_in)
744  {
745  if (const auto tria_pdt = dynamic_cast<
747  Assert(comm == tria_pdt->get_communicator(),
748  ExcMessage("MPI communicators do not match."));
749 
750  // First, figure out for what rank we are supposed to build the
751  // TriangulationDescription::Description object
752  unsigned int my_rank = my_rank_in;
755  ExcMessage("Rank has to be smaller than available processes."));
756 
757  if (auto tria_pdt = dynamic_cast<
759  {
762  ExcMessage(
763  "For creation from a parallel::distributed::Triangulation, "
764  "my_rank has to equal global rank."));
765 
767  }
768  else if (auto tria_serial =
769  dynamic_cast<const ::Triangulation<dim, spacedim> *>(
770  &tria))
771  {
772  if (my_rank == numbers::invalid_unsigned_int)
774  }
775  else
776  {
777  Assert(false,
778  ExcMessage("This type of triangulation is not supported!"));
779  }
780 
781  const auto subdomain_id_function = [](const auto &cell) {
782  return cell->subdomain_id();
783  };
784 
785  const auto level_subdomain_id_function = [](const auto &cell) {
786  return cell->level_subdomain_id();
787  };
788 
789  return CreateDescriptionFromTriangulationHelper<dim, spacedim>(
790  tria,
793  comm,
794  settings)
795  .template create_description_for_rank<Description<dim, spacedim>>(
796  my_rank);
797  }
798 
799 
800 
801  template <int dim, int spacedim>
804  const std::function<void(::Triangulation<dim, spacedim> &)>
805  & serial_grid_generator,
806  const std::function<void(::Triangulation<dim, spacedim> &,
807  const MPI_Comm &,
808  const unsigned int)> &serial_grid_partitioner,
809  const MPI_Comm & comm,
810  const int group_size,
811  const typename Triangulation<dim, spacedim>::MeshSmoothing smoothing,
813  {
814 #ifndef DEAL_II_WITH_MPI
815  (void)serial_grid_generator;
816  (void)serial_grid_partitioner;
817  (void)comm;
818  (void)group_size;
819  (void)smoothing;
820  (void)settings;
821 
823 #else
824  const unsigned int my_rank =
826  const unsigned int group_root = (my_rank / group_size) * group_size;
827 
828  const int mpi_tag =
830 
831  // check if process is root of the group
832  if (my_rank == group_root)
833  {
834  // Step 1: create serial triangulation
838  static_cast<
839  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
840  smoothing |
841  Triangulation<dim,
842  spacedim>::limit_level_difference_at_vertices) :
843  smoothing);
844  serial_grid_generator(tria);
845 
846  // Step 2: partition active cells and ...
847  serial_grid_partitioner(tria, comm, group_size);
848 
849  // ... cells on the levels if multigrid is required
850  if (settings &
853 
854  const unsigned int end_group =
855  std::min(group_root + group_size,
857 
858  // 3) create Description for the other processes in group; since
859  // we expect that this function is called for huge meshes, one
860  // Description is created at a time and send away; only once the
861  // Description has been sent away, the next rank is processed.
862  for (unsigned int other_rank = group_root + 1; other_rank < end_group;
863  ++other_rank)
864  {
865  // 3a) create construction data for other ranks
866  const auto construction_data =
868  comm,
869  settings,
870  other_rank);
871  // 3b) pack
872  std::vector<char> buffer;
873  ::Utilities::pack(construction_data, buffer, false);
874 
875  // 3c) send TriangulationDescription::Description
876  const auto ierr = MPI_Send(buffer.data(),
877  buffer.size(),
878  MPI_CHAR,
879  other_rank,
880  mpi_tag,
881  comm);
882  AssertThrowMPI(ierr);
883  }
884 
885  // 4) create TriangulationDescription::Description for this process
886  // (root of group)
888  comm,
889  settings,
890  my_rank);
891  }
892  else
893  {
894  // 3a) recv packed TriangulationDescription::Description from
895  // group-root process
896  // (counter-part of 3c of root process)
897  MPI_Status status;
898  auto ierr = MPI_Probe(group_root, mpi_tag, comm, &status);
899  AssertThrowMPI(ierr);
900 
901  int len;
902  MPI_Get_count(&status, MPI_CHAR, &len);
903 
904  std::vector<char> buf(len);
905  ierr = MPI_Recv(buf.data(),
906  len,
907  MPI_CHAR,
908  status.MPI_SOURCE,
909  mpi_tag,
910  comm,
911  &status);
912  AssertThrowMPI(ierr);
913 
914  // 3b) unpack TriangulationDescription::Description (counter-part of
915  // 3b of root process)
916  auto construction_data =
917  ::Utilities::template unpack<Description<dim, spacedim>>(
918  buf, false);
919 
920  // WARNING: serialization cannot handle the MPI communicator
921  // which is the reason why we have to set it here explicitly
922  construction_data.comm = comm;
923 
924  return construction_data;
925  }
926 #endif
927  }
928 
929 
930 
931  template <int dim, int spacedim>
937  {
938  const bool construct_multigrid =
939  (partition.size() > 0) &&
940  (settings &
942 
943  Assert(
944  construct_multigrid == false ||
947  ExcMessage(
948  "Source triangulation has to be set up with "
949  "limit_level_difference_at_vertices if the construction of the "
950  "multigrid hierarchy is requested!"));
951 
952  std::vector<LinearAlgebra::distributed::Vector<double>> partitions_mg;
953 
954  if (construct_multigrid) // perform first child policy
955  {
956  const auto tria_parallel =
957  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
958  &tria);
959 
960  Assert(tria_parallel, ExcInternalError());
961 
962  partition.update_ghost_values();
963 
964  partitions_mg.resize(tria.n_global_levels());
965 
966  for (unsigned int l = 0; l < tria.n_global_levels(); ++l)
967  partitions_mg[l].reinit(
968  tria_parallel->global_level_cell_index_partitioner(l).lock());
969 
970  for (int level = tria.n_global_levels() - 1; level >= 0; --level)
971  {
972  for (const auto &cell : tria.cell_iterators_on_level(level))
973  {
974  if (cell->is_locally_owned_on_level() == false)
975  continue;
976 
977  if (cell->is_active())
978  partitions_mg[level][cell->global_level_cell_index()] =
979  partition[cell->global_active_cell_index()];
980  else
981  partitions_mg[level][cell->global_level_cell_index()] =
982  partitions_mg[level + 1]
983  [cell->child(0)->global_level_cell_index()];
984  }
985 
986  partitions_mg[level].update_ghost_values();
987  }
988  }
989 
991  partition,
992  partitions_mg,
993  settings);
994  }
995 
996 
997 
998  template <int dim, int spacedim>
1004  & partitions_mg,
1005  const TriangulationDescription::Settings settings_in)
1006  {
1007 #ifdef DEAL_II_WITH_MPI
1008  if (tria.get_communicator() == MPI_COMM_NULL)
1009  AssertDimension(partition.local_size(), 0);
1010 #endif
1011 
1012  if (partition.size() == 0)
1013  {
1014  AssertDimension(partitions_mg.size(), 0);
1017  settings_in);
1018  }
1019 
1020  partition.update_ghost_values();
1021 
1022  for (const auto &partition : partitions_mg)
1023  partition.update_ghost_values();
1024 
1025  // 1) determine processes owning locally owned cells
1026  const std::vector<unsigned int> relevant_processes = [&]() {
1027  std::set<unsigned int> relevant_processes;
1028 
1029  for (unsigned int i = 0; i < partition.local_size(); ++i)
1030  relevant_processes.insert(
1031  static_cast<unsigned int>(partition.local_element(i)));
1032 
1033  for (const auto &partition : partitions_mg)
1034  for (unsigned int i = 0; i < partition.local_size(); ++i)
1035  relevant_processes.insert(
1036  static_cast<unsigned int>(partition.local_element(i)));
1037 
1038  return std::vector<unsigned int>(relevant_processes.begin(),
1039  relevant_processes.end());
1040  }();
1041 
1042  const bool construct_multigrid = partitions_mg.size() > 0;
1043 
1045 
1046  if (construct_multigrid)
1048  settings |
1050 
1051  const auto subdomain_id_function = [&partition](const auto &cell) {
1052  if ((cell->is_active() && (cell->is_artificial() == false)))
1053  return static_cast<unsigned int>(
1054  partition[cell->global_active_cell_index()]);
1055  else
1057  };
1058 
1059  const auto level_subdomain_id_function =
1060  [&construct_multigrid, &partitions_mg](const auto &cell) {
1061  if (construct_multigrid && (cell->is_artificial_on_level() == false))
1062  return static_cast<unsigned int>(
1063  partitions_mg[cell->level()][cell->global_level_cell_index()]);
1064  else
1066  };
1067 
1068  CreateDescriptionFromTriangulationHelper<dim, spacedim> helper(
1069  tria,
1073  settings);
1074 
1075  // create a description (locally owned cell and a layer of ghost cells
1076  // and all their parents)
1077  std::vector<DescriptionTemp<dim, spacedim>> descriptions_per_rank;
1078  descriptions_per_rank.reserve(relevant_processes.size());
1079 
1080  for (const auto rank : relevant_processes)
1081  descriptions_per_rank.emplace_back(
1082  helper.template create_description_for_rank<
1083  DescriptionTemp<dim, spacedim>>(rank));
1084 
1085  // collect description from all processes that used to own locally-owned
1086  // active cells of this process in a single description
1087  DescriptionTemp<dim, spacedim> description_merged;
1088  description_merged.collect(
1089  relevant_processes,
1090  descriptions_per_rank,
1091  partition.get_mpi_communicator(),
1092  dynamic_cast<
1094  &tria) == nullptr);
1095 
1096  // remove redundant entries
1097  description_merged.reduce();
1098 
1099  // convert to actual description
1100  return description_merged.convert(partition.get_mpi_communicator(),
1102  settings);
1103  }
1104 
1105  } // namespace Utilities
1106 } // namespace TriangulationDescription
1107 
1108 
1109 
1110 /*-------------- Explicit Instantiations -------------------------------*/
1111 #include "tria_description.inst"
1112 
1113 
virtual const MeshSmoothing & get_mesh_smoothing() const
virtual MPI_Comm get_communicator() const
void save_user_flags(std::ostream &out) 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:433
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:434
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4604
unsigned int level
Definition: grid_out.cc:4606
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:151
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:140
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1443
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