Reference documentation for deal.II version Git a285ab78cc 2020-08-04 19:09:24 +0200
\(\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.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 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 
16 
17 #include <deal.II/base/logstream.h>
19 #include <deal.II/base/utilities.h>
20 
23 
25 #include <deal.II/grid/tria.h>
28 
31 
32 #include <algorithm>
33 #include <fstream>
34 #include <iostream>
35 #include <numeric>
36 
37 
39 
40 
41 #ifdef DEAL_II_WITH_P4EST
42 
43 namespace
44 {
45  template <int dim, int spacedim>
46  void
47  get_vertex_to_cell_mappings(
49  std::vector<unsigned int> & vertex_touch_count,
50  std::vector<std::list<
52  unsigned int>>> & vertex_to_cell)
53  {
54  vertex_touch_count.resize(triangulation.n_vertices());
55  vertex_to_cell.resize(triangulation.n_vertices());
56 
57  for (const auto &cell : triangulation.active_cell_iterators())
58  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
59  {
60  ++vertex_touch_count[cell->vertex_index(v)];
61  vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
62  }
63  }
64 
65 
66 
67  template <int dim, int spacedim>
68  void
69  get_edge_to_cell_mappings(
70  const Triangulation<dim, spacedim> &triangulation,
71  std::vector<unsigned int> & edge_touch_count,
72  std::vector<std::list<
74  unsigned int>>> & edge_to_cell)
75  {
76  Assert(triangulation.n_levels() == 1, ExcInternalError());
77 
78  edge_touch_count.resize(triangulation.n_active_lines());
79  edge_to_cell.resize(triangulation.n_active_lines());
80 
81  for (const auto &cell : triangulation.active_cell_iterators())
82  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
83  {
84  ++edge_touch_count[cell->line(l)->index()];
85  edge_to_cell[cell->line(l)->index()].emplace_back(cell, l);
86  }
87  }
88 
89 
90 
95  template <int dim, int spacedim>
96  void
97  set_vertex_and_cell_info(
98  const Triangulation<dim, spacedim> &triangulation,
99  const std::vector<unsigned int> & vertex_touch_count,
100  const std::vector<std::list<
102  unsigned int>>> & vertex_to_cell,
103  const std::vector<types::global_dof_index>
104  & coarse_cell_to_p4est_tree_permutation,
105  const bool set_vertex_info,
106  typename internal::p4est::types<dim>::connectivity *connectivity)
107  {
108  // copy the vertices into the connectivity structure. the triangulation
109  // exports the array of vertices, but some of the entries are sometimes
110  // unused; this shouldn't be the case for a newly created triangulation,
111  // but make sure
112  //
113  // note that p4est stores coordinates as a triplet of values even in 2d
114  Assert(triangulation.get_used_vertices().size() ==
115  triangulation.get_vertices().size(),
116  ExcInternalError());
117  Assert(std::find(triangulation.get_used_vertices().begin(),
118  triangulation.get_used_vertices().end(),
119  false) == triangulation.get_used_vertices().end(),
120  ExcInternalError());
121  if (set_vertex_info == true)
122  for (unsigned int v = 0; v < triangulation.n_vertices(); ++v)
123  {
124  connectivity->vertices[3 * v] = triangulation.get_vertices()[v][0];
125  connectivity->vertices[3 * v + 1] =
126  triangulation.get_vertices()[v][1];
127  connectivity->vertices[3 * v + 2] =
128  (spacedim == 2 ? 0 : triangulation.get_vertices()[v][2]);
129  }
130 
131  // next store the tree_to_vertex indices (each tree is here only a single
132  // cell in the coarse mesh). p4est requires vertex numbering in clockwise
133  // orientation
134  //
135  // while we're at it, also copy the neighborship information between cells
137  cell = triangulation.begin_active(),
138  endc = triangulation.end();
139  for (; cell != endc; ++cell)
140  {
141  const unsigned int index =
142  coarse_cell_to_p4est_tree_permutation[cell->index()];
143 
144  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
145  {
146  if (set_vertex_info == true)
147  connectivity
148  ->tree_to_vertex[index * GeometryInfo<dim>::vertices_per_cell +
149  v] = cell->vertex_index(v);
150  connectivity
151  ->tree_to_corner[index * GeometryInfo<dim>::vertices_per_cell +
152  v] = cell->vertex_index(v);
153  }
154 
155  // neighborship information. if a cell is at a boundary, then enter
156  // the index of the cell itself here
157  for (auto f : GeometryInfo<dim>::face_indices())
158  if (cell->face(f)->at_boundary() == false)
159  connectivity
160  ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
161  coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->index()];
162  else
163  connectivity
164  ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
165  coarse_cell_to_p4est_tree_permutation[cell->index()];
166 
167  // fill tree_to_face, which is essentially neighbor_to_neighbor;
168  // however, we have to remap the resulting face number as well
169  for (auto f : GeometryInfo<dim>::face_indices())
170  if (cell->face(f)->at_boundary() == false)
171  {
172  switch (dim)
173  {
174  case 2:
175  {
176  connectivity->tree_to_face
177  [index * GeometryInfo<dim>::faces_per_cell + f] =
178  cell->neighbor_of_neighbor(f);
179  break;
180  }
181 
182  case 3:
183  {
184  /*
185  * The values for tree_to_face are in 0..23 where ttf % 6
186  * gives the face number and ttf / 4 the face orientation
187  * code. The orientation is determined as follows. Let
188  * my_face and other_face be the two face numbers of the
189  * connecting trees in 0..5. Then the first face vertex
190  * of the lower of my_face and other_face connects to a
191  * face vertex numbered 0..3 in the higher of my_face and
192  * other_face. The face orientation is defined as this
193  * number. If my_face == other_face, treating either of
194  * both faces as the lower one leads to the same result.
195  */
196 
197  connectivity->tree_to_face[index * 6 + f] =
198  cell->neighbor_of_neighbor(f);
199 
200  unsigned int face_idx_list[2] = {
201  f, cell->neighbor_of_neighbor(f)};
203  cell_list[2] = {cell, cell->neighbor(f)};
204  unsigned int smaller_idx = 0;
205 
206  if (f > cell->neighbor_of_neighbor(f))
207  smaller_idx = 1;
208 
209  unsigned int larger_idx = (smaller_idx + 1) % 2;
210  // smaller = *_list[smaller_idx]
211  // larger = *_list[larger_idx]
212 
213  unsigned int v = 0;
214 
215  // global vertex index of vertex 0 on face of cell with
216  // smaller local face index
217  unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
219  face_idx_list[smaller_idx],
220  0,
221  cell_list[smaller_idx]->face_orientation(
222  face_idx_list[smaller_idx]),
223  cell_list[smaller_idx]->face_flip(
224  face_idx_list[smaller_idx]),
225  cell_list[smaller_idx]->face_rotation(
226  face_idx_list[smaller_idx])));
227 
228  // loop over vertices on face from other cell and compare
229  // global vertex numbers
230  for (unsigned int i = 0;
231  i < GeometryInfo<dim>::vertices_per_face;
232  ++i)
233  {
234  unsigned int idx =
235  cell_list[larger_idx]->vertex_index(
237  face_idx_list[larger_idx], i));
238 
239  if (idx == g_idx)
240  {
241  v = i;
242  break;
243  }
244  }
245 
246  connectivity->tree_to_face[index * 6 + f] += 6 * v;
247  break;
248  }
249 
250  default:
251  Assert(false, ExcNotImplemented());
252  }
253  }
254  else
255  connectivity
256  ->tree_to_face[index * GeometryInfo<dim>::faces_per_cell + f] = f;
257  }
258 
259  // now fill the vertex information
260  connectivity->ctt_offset[0] = 0;
261  std::partial_sum(vertex_touch_count.begin(),
262  vertex_touch_count.end(),
263  &connectivity->ctt_offset[1]);
264 
265  const typename internal::p4est::types<dim>::locidx num_vtt =
266  std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
267  (void)num_vtt;
268  Assert(connectivity->ctt_offset[triangulation.n_vertices()] == num_vtt,
269  ExcInternalError());
270 
271  for (unsigned int v = 0; v < triangulation.n_vertices(); ++v)
272  {
273  Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
274  ExcInternalError());
275 
276  typename std::list<
277  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
278  unsigned int>>::const_iterator p =
279  vertex_to_cell[v].begin();
280  for (unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
281  {
282  connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
283  coarse_cell_to_p4est_tree_permutation[p->first->index()];
284  connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
285  p->second;
286  }
287  }
288  }
289 
290 
291 
292  template <int dim, int spacedim>
293  bool
295  const typename internal::p4est::types<dim>::forest *parallel_forest,
296  const typename internal::p4est::types<dim>::topidx coarse_grid_cell)
297  {
298  Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
299  ExcInternalError());
300  return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
301  (coarse_grid_cell <= parallel_forest->last_local_tree));
302  }
303 
304 
305  template <int dim, int spacedim>
306  void
307  delete_all_children_and_self(
308  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
309  {
310  if (cell->has_children())
311  for (unsigned int c = 0; c < cell->n_children(); ++c)
312  delete_all_children_and_self<dim, spacedim>(cell->child(c));
313  else
314  cell->set_coarsen_flag();
315  }
316 
317 
318 
319  template <int dim, int spacedim>
320  void
321  delete_all_children(
322  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
323  {
324  if (cell->has_children())
325  for (unsigned int c = 0; c < cell->n_children(); ++c)
326  delete_all_children_and_self<dim, spacedim>(cell->child(c));
327  }
328 
329 
330  template <int dim, int spacedim>
331  void
332  determine_level_subdomain_id_recursively(
333  const typename internal::p4est::types<dim>::tree & tree,
335  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
336  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
337  typename internal::p4est::types<dim>::forest & forest,
338  const types::subdomain_id my_subdomain,
339  const std::vector<std::vector<bool>> & marked_vertices)
340  {
341  if (dealii_cell->level_subdomain_id() == numbers::artificial_subdomain_id)
342  {
343  // important: only assign the level_subdomain_id if it is a ghost cell
344  // even though we could fill in all.
345  bool used = false;
346  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
347  {
348  if (marked_vertices[dealii_cell->level()]
349  [dealii_cell->vertex_index(v)])
350  {
351  used = true;
352  break;
353  }
354  }
355 
356  // Special case: if this cell is active we might be a ghost neighbor
357  // to a locally owned cell across a vertex that is finer.
358  // Example (M= my, O=dealii_cell, owned by somebody else):
359  // *------*
360  // | |
361  // | O |
362  // | |
363  // *---*---*------*
364  // | M | M |
365  // *---*---*
366  // | | M |
367  // *---*---*
368  if (!used && dealii_cell->is_active() &&
369  dealii_cell->is_artificial() == false &&
370  dealii_cell->level() + 1 < static_cast<int>(marked_vertices.size()))
371  {
372  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
373  {
374  if (marked_vertices[dealii_cell->level() + 1]
375  [dealii_cell->vertex_index(v)])
376  {
377  used = true;
378  break;
379  }
380  }
381  }
382 
383  // Like above, but now the other way around
384  if (!used && dealii_cell->is_active() &&
385  dealii_cell->is_artificial() == false && dealii_cell->level() > 0)
386  {
387  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
388  {
389  if (marked_vertices[dealii_cell->level() - 1]
390  [dealii_cell->vertex_index(v)])
391  {
392  used = true;
393  break;
394  }
395  }
396  }
397 
398  if (used)
399  {
401  &forest, tree_index, &p4est_cell, my_subdomain);
402  Assert((owner != -2) && (owner != -1),
403  ExcMessage("p4est should know the owner."));
404  dealii_cell->set_level_subdomain_id(owner);
405  }
406  }
407 
408  if (dealii_cell->has_children())
409  {
412  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
413  ++c)
414  switch (dim)
415  {
416  case 2:
417  P4EST_QUADRANT_INIT(&p4est_child[c]);
418  break;
419  case 3:
420  P8EST_QUADRANT_INIT(&p4est_child[c]);
421  break;
422  default:
423  Assert(false, ExcNotImplemented());
424  }
425 
426 
428  p4est_child);
429 
430  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
431  ++c)
432  {
433  determine_level_subdomain_id_recursively<dim, spacedim>(
434  tree,
435  tree_index,
436  dealii_cell->child(c),
437  p4est_child[c],
438  forest,
439  my_subdomain,
440  marked_vertices);
441  }
442  }
443  }
444 
445 
446  template <int dim, int spacedim>
447  void
448  match_tree_recursively(
449  const typename internal::p4est::types<dim>::tree & tree,
450  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
451  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
452  const typename internal::p4est::types<dim>::forest & forest,
453  const types::subdomain_id my_subdomain)
454  {
455  // check if this cell exists in the local p4est cell
456  if (sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
457  &p4est_cell,
459  -1)
460  {
461  // yes, cell found in local part of p4est
462  delete_all_children<dim, spacedim>(dealii_cell);
463  if (dealii_cell->is_active())
464  dealii_cell->set_subdomain_id(my_subdomain);
465  }
466  else
467  {
468  // no, cell not found in local part of p4est. this means that the
469  // local part is more refined than the current cell. if this cell has
470  // no children of its own, we need to refine it, and if it does
471  // already have children then loop over all children and see if they
472  // are locally available as well
473  if (dealii_cell->is_active())
474  dealii_cell->set_refine_flag();
475  else
476  {
479  for (unsigned int c = 0;
480  c < GeometryInfo<dim>::max_children_per_cell;
481  ++c)
482  switch (dim)
483  {
484  case 2:
485  P4EST_QUADRANT_INIT(&p4est_child[c]);
486  break;
487  case 3:
488  P8EST_QUADRANT_INIT(&p4est_child[c]);
489  break;
490  default:
491  Assert(false, ExcNotImplemented());
492  }
493 
494 
496  p4est_child);
497 
498  for (unsigned int c = 0;
499  c < GeometryInfo<dim>::max_children_per_cell;
500  ++c)
502  const_cast<typename internal::p4est::types<dim>::tree *>(
503  &tree),
504  &p4est_child[c]) == false)
505  {
506  // no, this child is locally not available in the p4est.
507  // delete all its children but, because this may not be
508  // successful, make sure to mark all children recursively
509  // as not local.
510  delete_all_children<dim, spacedim>(dealii_cell->child(c));
511  dealii_cell->child(c)->recursively_set_subdomain_id(
513  }
514  else
515  {
516  // at least some part of the tree rooted in this child is
517  // locally available
518  match_tree_recursively<dim, spacedim>(tree,
519  dealii_cell->child(c),
520  p4est_child[c],
521  forest,
522  my_subdomain);
523  }
524  }
525  }
526  }
527 
528 
529  template <int dim, int spacedim>
530  void
531  match_quadrant(
532  const ::Triangulation<dim, spacedim> * tria,
533  unsigned int dealii_index,
534  const typename internal::p4est::types<dim>::quadrant &ghost_quadrant,
535  types::subdomain_id ghost_owner)
536  {
537  const int l = ghost_quadrant.level;
538 
539  for (int i = 0; i < l; ++i)
540  {
542  i,
543  dealii_index);
544  if (cell->is_active())
545  {
546  cell->clear_coarsen_flag();
547  cell->set_refine_flag();
548  return;
549  }
550 
551  const int child_id =
553  i + 1);
554  dealii_index = cell->child_index(child_id);
555  }
556 
558  l,
559  dealii_index);
560  if (cell->has_children())
561  delete_all_children<dim, spacedim>(cell);
562  else
563  {
564  cell->clear_coarsen_flag();
565  cell->set_subdomain_id(ghost_owner);
566  }
567  }
568 
569 
570 
576  template <int dim, int spacedim>
577  class RefineAndCoarsenList
578  {
579  public:
580  RefineAndCoarsenList(const Triangulation<dim, spacedim> &triangulation,
581  const std::vector<types::global_dof_index>
582  &p4est_tree_to_coarse_cell_permutation,
583  const types::subdomain_id my_subdomain);
584 
593  static int
594  refine_callback(
595  typename internal::p4est::types<dim>::forest * forest,
596  typename internal::p4est::types<dim>::topidx coarse_cell_index,
597  typename internal::p4est::types<dim>::quadrant *quadrant);
598 
603  static int
604  coarsen_callback(
605  typename internal::p4est::types<dim>::forest * forest,
606  typename internal::p4est::types<dim>::topidx coarse_cell_index,
607  typename internal::p4est::types<dim>::quadrant *children[]);
608 
609  bool
610  pointers_are_at_end() const;
611 
612  private:
613  std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
614  typename std::vector<typename internal::p4est::types<dim>::quadrant>::
615  const_iterator current_refine_pointer;
616 
617  std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
618  typename std::vector<typename internal::p4est::types<dim>::quadrant>::
619  const_iterator current_coarsen_pointer;
620 
621  void
622  build_lists(
623  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
624  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
625  const types::subdomain_id myid);
626  };
627 
628 
629 
630  template <int dim, int spacedim>
631  bool
632  RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end() const
633  {
634  return ((current_refine_pointer == refine_list.end()) &&
635  (current_coarsen_pointer == coarsen_list.end()));
636  }
637 
638 
639 
640  template <int dim, int spacedim>
641  RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
642  const Triangulation<dim, spacedim> &triangulation,
643  const std::vector<types::global_dof_index>
644  & p4est_tree_to_coarse_cell_permutation,
645  const types::subdomain_id my_subdomain)
646  {
647  // count how many flags are set and allocate that much memory
648  unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
649  for (const auto &cell : triangulation.active_cell_iterators())
650  {
651  // skip cells that are not local
652  if (cell->subdomain_id() != my_subdomain)
653  continue;
654 
655  if (cell->refine_flag_set())
656  ++n_refine_flags;
657  else if (cell->coarsen_flag_set())
658  ++n_coarsen_flags;
659  }
660 
661  refine_list.reserve(n_refine_flags);
662  coarsen_list.reserve(n_coarsen_flags);
663 
664 
665  // now build the lists of cells that are flagged. note that p4est will
666  // traverse its cells in the order in which trees appear in the
667  // forest. this order is not the same as the order of coarse cells in the
668  // deal.II Triangulation because we have translated everything by the
669  // coarse_cell_to_p4est_tree_permutation permutation. in order to make
670  // sure that the output array is already in the correct order, traverse
671  // our coarse cells in the same order in which p4est will:
672  for (unsigned int c = 0; c < triangulation.n_cells(0); ++c)
673  {
674  unsigned int coarse_cell_index =
675  p4est_tree_to_coarse_cell_permutation[c];
676 
678  &triangulation, 0, coarse_cell_index);
679 
680  typename internal::p4est::types<dim>::quadrant p4est_cell;
682  /*level=*/0,
683  /*index=*/0);
684  p4est_cell.p.which_tree = c;
685  build_lists(cell, p4est_cell, my_subdomain);
686  }
687 
688 
689  Assert(refine_list.size() == n_refine_flags, ExcInternalError());
690  Assert(coarsen_list.size() == n_coarsen_flags, ExcInternalError());
691 
692  // make sure that our ordering in fact worked
693  for (unsigned int i = 1; i < refine_list.size(); ++i)
694  Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
695  ExcInternalError());
696  for (unsigned int i = 1; i < coarsen_list.size(); ++i)
697  Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
698  ExcInternalError());
699 
700  current_refine_pointer = refine_list.begin();
701  current_coarsen_pointer = coarsen_list.begin();
702  }
703 
704 
705 
706  template <int dim, int spacedim>
707  void
708  RefineAndCoarsenList<dim, spacedim>::build_lists(
709  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
710  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
711  const types::subdomain_id my_subdomain)
712  {
713  if (cell->is_active())
714  {
715  if (cell->subdomain_id() == my_subdomain)
716  {
717  if (cell->refine_flag_set())
718  refine_list.push_back(p4est_cell);
719  else if (cell->coarsen_flag_set())
720  coarsen_list.push_back(p4est_cell);
721  }
722  }
723  else
724  {
727  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
728  ++c)
729  switch (dim)
730  {
731  case 2:
732  P4EST_QUADRANT_INIT(&p4est_child[c]);
733  break;
734  case 3:
735  P8EST_QUADRANT_INIT(&p4est_child[c]);
736  break;
737  default:
738  Assert(false, ExcNotImplemented());
739  }
741  p4est_child);
742  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
743  ++c)
744  {
745  p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
746  build_lists(cell->child(c), p4est_child[c], my_subdomain);
747  }
748  }
749  }
750 
751 
752  template <int dim, int spacedim>
753  int
754  RefineAndCoarsenList<dim, spacedim>::refine_callback(
755  typename internal::p4est::types<dim>::forest * forest,
756  typename internal::p4est::types<dim>::topidx coarse_cell_index,
757  typename internal::p4est::types<dim>::quadrant *quadrant)
758  {
759  RefineAndCoarsenList<dim, spacedim> *this_object =
760  reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *>(
761  forest->user_pointer);
762 
763  // if there are no more cells in our list the current cell can't be
764  // flagged for refinement
765  if (this_object->current_refine_pointer == this_object->refine_list.end())
766  return false;
767 
768  Assert(coarse_cell_index <=
769  this_object->current_refine_pointer->p.which_tree,
770  ExcInternalError());
771 
772  // if p4est hasn't yet reached the tree of the next flagged cell the
773  // current cell can't be flagged for refinement
774  if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
775  return false;
776 
777  // now we're in the right tree in the forest
778  Assert(coarse_cell_index <=
779  this_object->current_refine_pointer->p.which_tree,
780  ExcInternalError());
781 
782  // make sure that the p4est loop over cells hasn't gotten ahead of our own
783  // pointer
785  quadrant, &*this_object->current_refine_pointer) <= 0,
786  ExcInternalError());
787 
788  // now, if the p4est cell is one in the list, it is supposed to be refined
790  quadrant, &*this_object->current_refine_pointer))
791  {
792  ++this_object->current_refine_pointer;
793  return true;
794  }
795 
796  // p4est cell is not in list
797  return false;
798  }
799 
800 
801 
802  template <int dim, int spacedim>
803  int
804  RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
805  typename internal::p4est::types<dim>::forest * forest,
806  typename internal::p4est::types<dim>::topidx coarse_cell_index,
807  typename internal::p4est::types<dim>::quadrant *children[])
808  {
809  RefineAndCoarsenList<dim, spacedim> *this_object =
810  reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *>(
811  forest->user_pointer);
812 
813  // if there are no more cells in our list the current cell can't be
814  // flagged for coarsening
815  if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
816  return false;
817 
818  Assert(coarse_cell_index <=
819  this_object->current_coarsen_pointer->p.which_tree,
820  ExcInternalError());
821 
822  // if p4est hasn't yet reached the tree of the next flagged cell the
823  // current cell can't be flagged for coarsening
824  if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
825  return false;
826 
827  // now we're in the right tree in the forest
828  Assert(coarse_cell_index <=
829  this_object->current_coarsen_pointer->p.which_tree,
830  ExcInternalError());
831 
832  // make sure that the p4est loop over cells hasn't gotten ahead of our own
833  // pointer
835  children[0], &*this_object->current_coarsen_pointer) <= 0,
836  ExcInternalError());
837 
838  // now, if the p4est cell is one in the list, it is supposed to be
839  // coarsened
841  children[0], &*this_object->current_coarsen_pointer))
842  {
843  // move current pointer one up
844  ++this_object->current_coarsen_pointer;
845 
846  // note that the next 3 cells in our list need to correspond to the
847  // other siblings of the cell we have just found
848  for (unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
849  ++c)
850  {
852  children[c], &*this_object->current_coarsen_pointer),
853  ExcInternalError());
854  ++this_object->current_coarsen_pointer;
855  }
856 
857  return true;
858  }
859 
860  // p4est cell is not in list
861  return false;
862  }
863 
864 
865 
872  template <int dim, int spacedim>
873  class PartitionWeights
874  {
875  public:
881  explicit PartitionWeights(const std::vector<unsigned int> &cell_weights);
882 
890  static int
891  cell_weight(typename internal::p4est::types<dim>::forest *forest,
892  typename internal::p4est::types<dim>::topidx coarse_cell_index,
893  typename internal::p4est::types<dim>::quadrant *quadrant);
894 
895  private:
896  std::vector<unsigned int> cell_weights_list;
897  std::vector<unsigned int>::const_iterator current_pointer;
898  };
899 
900 
901  template <int dim, int spacedim>
902  PartitionWeights<dim, spacedim>::PartitionWeights(
903  const std::vector<unsigned int> &cell_weights)
904  : cell_weights_list(cell_weights)
905  {
906  // set the current pointer to the first element of the list, given that
907  // we will walk through it sequentially
908  current_pointer = cell_weights_list.begin();
909  }
910 
911 
912  template <int dim, int spacedim>
913  int
914  PartitionWeights<dim, spacedim>::cell_weight(
915  typename internal::p4est::types<dim>::forest *forest,
918  {
919  // the function gets two additional arguments, but we don't need them
920  // since we know in which order p4est will walk through the cells
921  // and have already built our weight lists in this order
922 
923  PartitionWeights<dim, spacedim> *this_object =
924  reinterpret_cast<PartitionWeights<dim, spacedim> *>(forest->user_pointer);
925 
926  Assert(this_object->current_pointer >=
927  this_object->cell_weights_list.begin(),
928  ExcInternalError());
929  Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
930  ExcInternalError());
931 
932  // get the weight, increment the pointer, and return the weight
933  return *this_object->current_pointer++;
934  }
935 
936 
937 
938  template <int dim, int spacedim>
939  using quadrant_cell_relation_t = typename std::tuple<
940  typename ::internal::p4est::types<dim>::quadrant *,
941  typename ::Triangulation<dim, spacedim>::CellStatus,
942  typename ::Triangulation<dim, spacedim>::cell_iterator>;
943 
944 
945 
955  template <int dim, int spacedim>
956  inline void
957  add_single_quadrant_cell_relation(
958  std::vector<quadrant_cell_relation_t<dim, spacedim>> & quad_cell_rel,
959  const typename ::internal::p4est::types<dim>::tree & tree,
960  const unsigned int idx,
961  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
962  const typename Triangulation<dim, spacedim>::CellStatus status)
963  {
964  const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
965 
966  const auto q =
967  static_cast<typename ::internal::p4est::types<dim>::quadrant *>(
968  sc_array_index(const_cast<sc_array_t *>(&tree.quadrants), idx));
969 
970  // check if we will be writing into valid memory
971  Assert(local_quadrant_index < quad_cell_rel.size(), ExcInternalError());
972 
973  // store relation
974  quad_cell_rel[local_quadrant_index] =
975  std::make_tuple(q, status, dealii_cell);
976  }
977 
978 
979 
989  template <int dim, int spacedim>
990  void
991  update_quadrant_cell_relations_recursively(
992  std::vector<quadrant_cell_relation_t<dim, spacedim>> & quad_cell_rel,
993  const typename ::internal::p4est::types<dim>::tree & tree,
994  const typename Triangulation<dim, spacedim>::cell_iterator & dealii_cell,
995  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
996  {
997  // find index of p4est_cell in the quadrants array of the corresponding tree
998  const int idx = sc_array_bsearch(
999  const_cast<sc_array_t *>(&tree.quadrants),
1000  &p4est_cell,
1002  if (idx == -1 &&
1004  const_cast<typename ::internal::p4est::types<dim>::tree *>(
1005  &tree),
1006  &p4est_cell) == false))
1007  // this quadrant and none of its children belong to us.
1008  return;
1009 
1010  // recurse further if both p4est and dealii still have children
1011  const bool p4est_has_children = (idx == -1);
1012  if (p4est_has_children && dealii_cell->has_children())
1013  {
1014  // recurse further
1015  typename ::internal::p4est::types<dim>::quadrant
1017 
1018  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1019  ++c)
1020  switch (dim)
1021  {
1022  case 2:
1023  P4EST_QUADRANT_INIT(&p4est_child[c]);
1024  break;
1025  case 3:
1026  P8EST_QUADRANT_INIT(&p4est_child[c]);
1027  break;
1028  default:
1029  Assert(false, ExcNotImplemented());
1030  }
1031 
1033  &p4est_cell, p4est_child);
1034 
1035  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1036  ++c)
1037  {
1038  update_quadrant_cell_relations_recursively<dim, spacedim>(
1039  quad_cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1040  }
1041  }
1042  else if (!p4est_has_children && !dealii_cell->has_children())
1043  {
1044  // this active cell didn't change
1045  // save tuple into corresponding position
1046  add_single_quadrant_cell_relation<dim, spacedim>(
1047  quad_cell_rel,
1048  tree,
1049  idx,
1050  dealii_cell,
1052  }
1053  else if (p4est_has_children) // based on the conditions above, we know that
1054  // dealii_cell has no children
1055  {
1056  // this cell got refined in p4est, but the dealii_cell has not yet been
1057  // refined
1058 
1059  // this quadrant is not active
1060  // generate its children, and store information in those
1061  typename ::internal::p4est::types<dim>::quadrant
1063  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1064  ++c)
1065  switch (dim)
1066  {
1067  case 2:
1068  P4EST_QUADRANT_INIT(&p4est_child[c]);
1069  break;
1070  case 3:
1071  P8EST_QUADRANT_INIT(&p4est_child[c]);
1072  break;
1073  default:
1074  Assert(false, ExcNotImplemented());
1075  }
1076 
1078  &p4est_cell, p4est_child);
1079 
1080  // mark first child with CELL_REFINE and the remaining children with
1081  // CELL_INVALID, but associate them all with the parent cell unpack
1082  // algorithm will be called only on CELL_REFINE flagged quadrant
1083  int child_idx;
1084  typename Triangulation<dim, spacedim>::CellStatus cell_status;
1085  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1086  ++i)
1087  {
1088  child_idx = sc_array_bsearch(
1089  const_cast<sc_array_t *>(&tree.quadrants),
1090  &p4est_child[i],
1092 
1093  cell_status = (i == 0) ? Triangulation<dim, spacedim>::CELL_REFINE :
1095 
1096  add_single_quadrant_cell_relation<dim, spacedim>(
1097  quad_cell_rel, tree, child_idx, dealii_cell, cell_status);
1098  }
1099  }
1100  else // based on the conditions above, we know that p4est_cell has no
1101  // children, and the dealii_cell does
1102  {
1103  // its children got coarsened into this cell in p4est,
1104  // but the dealii_cell still has its children
1105  add_single_quadrant_cell_relation<dim, spacedim>(
1106  quad_cell_rel,
1107  tree,
1108  idx,
1109  dealii_cell,
1111  }
1112  }
1113 } // namespace
1114 
1115 
1116 
1117 namespace parallel
1118 {
1119  namespace distributed
1120  {
1121  /* ------------------ class DataTransfer<dim,spacedim> ----------------- */
1122 
1123 
1124  template <int dim, int spacedim>
1126  MPI_Comm mpi_communicator)
1127  : mpi_communicator(mpi_communicator)
1128  , variable_size_data_stored(false)
1129  {}
1130 
1131 
1132 
1133  template <int dim, int spacedim>
1134  void
1136  const std::vector<quadrant_cell_relation_t> &quad_cell_relations,
1137  const std::vector<typename CellAttachedData::pack_callback_t>
1138  &pack_callbacks_fixed,
1139  const std::vector<typename CellAttachedData::pack_callback_t>
1140  &pack_callbacks_variable)
1141  {
1142  Assert(src_data_fixed.size() == 0,
1143  ExcMessage("Previously packed data has not been released yet!"));
1144  Assert(src_sizes_variable.size() == 0, ExcInternalError());
1145 
1146  const unsigned int n_callbacks_fixed = pack_callbacks_fixed.size();
1147  const unsigned int n_callbacks_variable = pack_callbacks_variable.size();
1148 
1149  // Store information that we packed variable size data in
1150  // a member variable for later.
1151  variable_size_data_stored = (n_callbacks_variable > 0);
1152 
1153  // If variable transfer is scheduled, we will store the data size that
1154  // each variable size callback function writes in this auxiliary
1155  // container. The information will be stored by each cell in this vector
1156  // temporarily.
1157  std::vector<unsigned int> cell_sizes_variable_cumulative(
1158  n_callbacks_variable);
1159 
1160  // Prepare the buffer structure, in which each callback function will
1161  // store its data for each active cell.
1162  // The outmost shell in this container construct corresponds to the
1163  // data packed per cell. The next layer resembles the data that
1164  // each callback function packs on the corresponding cell. These
1165  // buffers are chains of chars stored in an std::vector<char>.
1166  // A visualisation of the data structure:
1167  /* clang-format off */
1168  // | cell_1 | | cell_2 | ...
1169  // || callback_1 || callback_2 |...| || callback_1 || callback_2 |...| ...
1170  // |||char|char|...|||char|char|...|...| |||char|char|...|||char|char|...|...| ...
1171  /* clang-format on */
1172  std::vector<std::vector<std::vector<char>>> packed_fixed_size_data(
1173  quad_cell_relations.size());
1174  std::vector<std::vector<std::vector<char>>> packed_variable_size_data(
1175  variable_size_data_stored ? quad_cell_relations.size() : 0);
1176 
1177  //
1178  // --------- Pack data for fixed and variable size transfer ---------
1179  //
1180  // Iterate over all cells, call all callback functions on each cell,
1181  // and store their data in the corresponding buffer scope.
1182  {
1183  auto quad_cell_rel_it = quad_cell_relations.cbegin();
1184  auto data_cell_fixed_it = packed_fixed_size_data.begin();
1185  auto data_cell_variable_it = packed_variable_size_data.begin();
1186  for (; quad_cell_rel_it != quad_cell_relations.cend();
1187  ++quad_cell_rel_it, ++data_cell_fixed_it)
1188  {
1189  const auto &cell_status = std::get<1>(*quad_cell_rel_it);
1190  const auto &dealii_cell = std::get<2>(*quad_cell_rel_it);
1191 
1192  // Assertions about the tree structure.
1193  switch (cell_status)
1194  {
1198  CELL_REFINE:
1199  // double check the condition that we will only ever attach
1200  // data to active cells when we get here
1201  Assert(dealii_cell->is_active(), ExcInternalError());
1202  break;
1203 
1206  // double check the condition that we will only ever attach
1207  // data to cells with children when we get here. however, we
1208  // can only tolerate one level of coarsening at a time, so
1209  // check that the children are all active
1210  Assert(dealii_cell->is_active() == false, ExcInternalError());
1211  for (unsigned int c = 0;
1212  c < GeometryInfo<dim>::max_children_per_cell;
1213  ++c)
1214  Assert(dealii_cell->child(c)->is_active(),
1215  ExcInternalError());
1216  break;
1217 
1220  // do nothing on invalid cells
1221  break;
1222 
1223  default:
1224  Assert(false, ExcInternalError());
1225  break;
1226  }
1227 
1228  // Reserve memory corresponding to the number of callback
1229  // functions that will be called.
1230  // If variable size transfer is scheduled, we need to leave
1231  // room for an array that holds information about how many
1232  // bytes each of the variable size callback functions will
1233  // write.
1234  // On cells flagged with CELL_INVALID, only its CellStatus
1235  // will be stored.
1236  const unsigned int n_fixed_size_data_sets_on_cell =
1237  1 +
1238  ((cell_status ==
1240  spacedim>::CELL_INVALID) ?
1241  0 :
1242  ((variable_size_data_stored ? 1 : 0) + n_callbacks_fixed));
1243  data_cell_fixed_it->resize(n_fixed_size_data_sets_on_cell);
1244 
1245  // We continue with packing all data on this specific cell.
1246  auto data_fixed_it = data_cell_fixed_it->begin();
1247 
1248  // First, we pack the CellStatus information.
1249  // to get consistent data sizes on each cell for the fixed size
1250  // transfer, we won't allow compression
1251  *data_fixed_it =
1252  Utilities::pack(cell_status, /*allow_compression=*/false);
1253  ++data_fixed_it;
1254 
1255  // Proceed with all registered callback functions.
1256  // Skip cells with the CELL_INVALID flag.
1257  if (cell_status !=
1259  spacedim>::CELL_INVALID)
1260  {
1261  // Pack fixed size data.
1262  for (auto callback_it = pack_callbacks_fixed.cbegin();
1263  callback_it != pack_callbacks_fixed.cend();
1264  ++callback_it, ++data_fixed_it)
1265  {
1266  *data_fixed_it = (*callback_it)(dealii_cell, cell_status);
1267  }
1268 
1269  // Pack variable size data.
1270  // If we store variable size data, we need to transfer
1271  // the sizes of each corresponding callback function
1272  // via fixed size transfer as well.
1274  {
1275  const unsigned int n_variable_size_data_sets_on_cell =
1276  ((cell_status ==
1278  CELL_INVALID) ?
1279  0 :
1280  n_callbacks_variable);
1281  data_cell_variable_it->resize(
1282  n_variable_size_data_sets_on_cell);
1283 
1284  auto callback_it = pack_callbacks_variable.cbegin();
1285  auto data_variable_it = data_cell_variable_it->begin();
1286  auto sizes_variable_it =
1287  cell_sizes_variable_cumulative.begin();
1288  for (; callback_it != pack_callbacks_variable.cend();
1289  ++callback_it, ++data_variable_it, ++sizes_variable_it)
1290  {
1291  *data_variable_it =
1292  (*callback_it)(dealii_cell, cell_status);
1293 
1294  // Store data sizes for each callback function first.
1295  // Make it cumulative below.
1296  *sizes_variable_it = data_variable_it->size();
1297  }
1298 
1299  // Turn size vector into its cumulative representation.
1300  std::partial_sum(cell_sizes_variable_cumulative.begin(),
1301  cell_sizes_variable_cumulative.end(),
1302  cell_sizes_variable_cumulative.begin());
1303 
1304  // Serialize cumulative variable size vector value-by-value.
1305  // This way we can circumvent the overhead of storing the
1306  // container object as a whole, since we know its size by
1307  // the number of registered callback functions.
1308  data_fixed_it->resize(n_callbacks_variable *
1309  sizeof(unsigned int));
1310  for (unsigned int i = 0; i < n_callbacks_variable; ++i)
1311  std::memcpy(&(data_fixed_it->at(i *
1312  sizeof(unsigned int))),
1313  &(cell_sizes_variable_cumulative.at(i)),
1314  sizeof(unsigned int));
1315 
1316  ++data_fixed_it;
1317  }
1318 
1319  // Double check that we packed everything we wanted
1320  // in the fixed size buffers.
1321  Assert(data_fixed_it == data_cell_fixed_it->end(),
1322  ExcInternalError());
1323  }
1324 
1325  // Increment the variable size data iterator
1326  // only if we actually pack this kind of data
1327  // to avoid getting out of bounds.
1329  ++data_cell_variable_it;
1330  } // loop over quad_cell_relations
1331  }
1332 
1333  //
1334  // ----------- Gather data sizes for fixed size transfer ------------
1335  //
1336  // Generate a vector which stores the sizes of each callback function,
1337  // including the packed CellStatus transfer.
1338  // Find the very first cell that we wrote to with all callback
1339  // functions (i.e. a cell that was not flagged with CELL_INVALID)
1340  // and store the sizes of each buffer.
1341  //
1342  // To deal with the case that at least one of the processors does not own
1343  // any cell at all, we will exchange the information about the data sizes
1344  // among them later. The code in between is still well-defined, since the
1345  // following loops will be skipped.
1346  std::vector<unsigned int> local_sizes_fixed(
1347  1 + n_callbacks_fixed + (variable_size_data_stored ? 1 : 0));
1348  for (const auto &data_cell : packed_fixed_size_data)
1349  {
1350  if (data_cell.size() == local_sizes_fixed.size())
1351  {
1352  auto sizes_fixed_it = local_sizes_fixed.begin();
1353  auto data_fixed_it = data_cell.cbegin();
1354  for (; data_fixed_it != data_cell.cend();
1355  ++data_fixed_it, ++sizes_fixed_it)
1356  {
1357  *sizes_fixed_it = data_fixed_it->size();
1358  }
1359 
1360  break;
1361  }
1362  }
1363 
1364  // Check if all cells have valid sizes.
1365  for (auto data_cell_fixed_it = packed_fixed_size_data.cbegin();
1366  data_cell_fixed_it != packed_fixed_size_data.cend();
1367  ++data_cell_fixed_it)
1368  {
1369  Assert((data_cell_fixed_it->size() == 1) ||
1370  (data_cell_fixed_it->size() == local_sizes_fixed.size()),
1371  ExcInternalError());
1372  }
1373 
1374  // Share information about the packed data sizes
1375  // of all callback functions across all processors, in case one
1376  // of them does not own any cells at all.
1377  std::vector<unsigned int> global_sizes_fixed(local_sizes_fixed.size());
1378  Utilities::MPI::max(local_sizes_fixed,
1379  this->mpi_communicator,
1380  global_sizes_fixed);
1381 
1382  // Construct cumulative sizes, since this is the only information
1383  // we need from now on.
1384  sizes_fixed_cumulative.resize(global_sizes_fixed.size());
1385  std::partial_sum(global_sizes_fixed.begin(),
1386  global_sizes_fixed.end(),
1387  sizes_fixed_cumulative.begin());
1388 
1389  //
1390  // ---------- Gather data sizes for variable size transfer ----------
1391  //
1393  {
1394  src_sizes_variable.reserve(packed_variable_size_data.size());
1395  for (const auto &data_cell : packed_variable_size_data)
1396  {
1397  int variable_data_size_on_cell = 0;
1398 
1399  for (const auto &data : data_cell)
1400  variable_data_size_on_cell += data.size();
1401 
1402  src_sizes_variable.push_back(variable_data_size_on_cell);
1403  }
1404  }
1405 
1406  //
1407  // ------------------------ Build buffers ---------------------------
1408  //
1409  const unsigned int expected_size_fixed =
1410  quad_cell_relations.size() * sizes_fixed_cumulative.back();
1411  const unsigned int expected_size_variable =
1412  std::accumulate(src_sizes_variable.begin(),
1413  src_sizes_variable.end(),
1415 
1416  // Move every piece of packed fixed size data into the consecutive buffer.
1417  src_data_fixed.reserve(expected_size_fixed);
1418  for (const auto &data_cell_fixed : packed_fixed_size_data)
1419  {
1420  // Move every fraction of packed data into the buffer
1421  // reserved for this particular cell.
1422  for (const auto &data_fixed : data_cell_fixed)
1423  std::move(data_fixed.begin(),
1424  data_fixed.end(),
1425  std::back_inserter(src_data_fixed));
1426 
1427  // If we only packed the CellStatus information
1428  // (i.e. encountered a cell flagged CELL_INVALID),
1429  // fill the remaining space with invalid entries.
1430  // We can skip this if there is nothing else to pack.
1431  if ((data_cell_fixed.size() == 1) &&
1432  (sizes_fixed_cumulative.size() > 1))
1433  {
1434  const std::size_t bytes_skipped =
1436 
1437  src_data_fixed.insert(src_data_fixed.end(),
1438  bytes_skipped,
1439  static_cast<char>(-1)); // invalid_char
1440  }
1441  }
1442 
1443  // Move every piece of packed variable size data into the consecutive
1444  // buffer.
1446  {
1447  src_data_variable.reserve(expected_size_variable);
1448  for (const auto &data_cell : packed_variable_size_data)
1449  {
1450  // Move every fraction of packed data into the buffer
1451  // reserved for this particular cell.
1452  for (const auto &data : data_cell)
1453  std::move(data.begin(),
1454  data.end(),
1455  std::back_inserter(src_data_variable));
1456  }
1457  }
1458 
1459  // Double check that we packed everything correctly.
1460  Assert(src_data_fixed.size() == expected_size_fixed, ExcInternalError());
1461  Assert(src_data_variable.size() == expected_size_variable,
1462  ExcInternalError());
1463  }
1464 
1465 
1466 
1467  template <int dim, int spacedim>
1468  void
1470  const typename ::internal::p4est::types<dim>::forest
1471  *parallel_forest,
1472  const typename ::internal::p4est::types<dim>::gloidx
1473  *previous_global_first_quadrant)
1474  {
1475  Assert(sizes_fixed_cumulative.size() > 0,
1476  ExcMessage("No data has been packed!"));
1477 
1478  // Resize memory according to the data that we will receive.
1479  dest_data_fixed.resize(parallel_forest->local_num_quadrants *
1480  sizes_fixed_cumulative.back());
1481 
1482  // Execute non-blocking fixed size transfer.
1483  typename ::internal::p4est::types<dim>::transfer_context
1484  *tf_context;
1485  tf_context =
1487  parallel_forest->global_first_quadrant,
1488  previous_global_first_quadrant,
1489  parallel_forest->mpicomm,
1490  0,
1491  dest_data_fixed.data(),
1492  src_data_fixed.data(),
1493  sizes_fixed_cumulative.back());
1494 
1496  {
1497  // Resize memory according to the data that we will receive.
1498  dest_sizes_variable.resize(parallel_forest->local_num_quadrants);
1499 
1500  // Execute fixed size transfer of data sizes for variable size
1501  // transfer.
1503  parallel_forest->global_first_quadrant,
1504  previous_global_first_quadrant,
1505  parallel_forest->mpicomm,
1506  1,
1507  dest_sizes_variable.data(),
1508  src_sizes_variable.data(),
1509  sizeof(int));
1510  }
1511 
1513 
1514  // Release memory of previously packed data.
1515  src_data_fixed.clear();
1516  src_data_fixed.shrink_to_fit();
1517 
1519  {
1520  // Resize memory according to the data that we will receive.
1521  dest_data_variable.resize(
1522  std::accumulate(dest_sizes_variable.begin(),
1523  dest_sizes_variable.end(),
1525 
1526 # if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1527 # else
1528  // ----- WORKAROUND -----
1529  // An assertion in p4est prevents us from sending/receiving no data
1530  // at all, which is mandatory if one of our processes does not own
1531  // any quadrant. This bypasses the assertion from being triggered.
1532  // - see: https://github.com/cburstedde/p4est/issues/48
1533  if (src_sizes_variable.size() == 0)
1534  src_sizes_variable.resize(1);
1535  if (dest_sizes_variable.size() == 0)
1536  dest_sizes_variable.resize(1);
1537 # endif
1538 
1539  // Execute variable size transfer.
1541  parallel_forest->global_first_quadrant,
1542  previous_global_first_quadrant,
1543  parallel_forest->mpicomm,
1544  1,
1545  dest_data_variable.data(),
1546  dest_sizes_variable.data(),
1547  src_data_variable.data(),
1548  src_sizes_variable.data());
1549 
1550  // Release memory of previously packed data.
1551  src_sizes_variable.clear();
1552  src_sizes_variable.shrink_to_fit();
1553  src_data_variable.clear();
1554  src_data_variable.shrink_to_fit();
1555  }
1556  }
1557 
1558 
1559 
1560  template <int dim, int spacedim>
1561  void
1563  std::vector<quadrant_cell_relation_t> &quad_cell_relations) const
1564  {
1565  Assert(sizes_fixed_cumulative.size() > 0,
1566  ExcMessage("No data has been packed!"));
1567  if (quad_cell_relations.size() > 0)
1568  {
1569  Assert(dest_data_fixed.size() > 0,
1570  ExcMessage("No data has been received!"));
1571  }
1572 
1573  // Size of CellStatus object that will be unpacked on each cell.
1574  const unsigned int size = sizes_fixed_cumulative.front();
1575 
1576  // Iterate over all cells and overwrite the CellStatus
1577  // information from the transferred data.
1578  // Proceed buffer iterator position to next cell after
1579  // each iteration.
1580  auto quad_cell_rel_it = quad_cell_relations.begin();
1581  auto dest_fixed_it = dest_data_fixed.cbegin();
1582  for (; quad_cell_rel_it != quad_cell_relations.end();
1583  ++quad_cell_rel_it, dest_fixed_it += sizes_fixed_cumulative.back())
1584  {
1585  std::get<1>(*quad_cell_rel_it) = // cell_status
1588  dest_fixed_it,
1589  dest_fixed_it + size,
1590  /*allow_compression=*/false);
1591  }
1592  }
1593 
1594 
1595 
1596  template <int dim, int spacedim>
1597  void
1599  const std::vector<quadrant_cell_relation_t> &quad_cell_relations,
1600  const unsigned int handle,
1601  const std::function<void(
1602  const typename ::Triangulation<dim, spacedim>::cell_iterator &,
1603  const typename ::Triangulation<dim, spacedim>::CellStatus &,
1604  const boost::iterator_range<std::vector<char>::const_iterator> &)>
1605  &unpack_callback) const
1606  {
1607  // We decode the handle returned by register_data_attach() back into
1608  // a format we can use. All even handles belong to those callback
1609  // functions which write/read variable size data, all odd handles interact
1610  // with fixed size buffers.
1611  const bool callback_variable_transfer = (handle % 2 == 0);
1612  const unsigned int callback_index = handle / 2;
1613 
1614  // Cells will always receive fixed size data (i.e., CellStatus
1615  // information), but not necessarily variable size data (e.g., with a
1616  // ParticleHandler a cell might not contain any particle at all).
1617  // Thus it is sufficient to check if fixed size data has been received.
1618  Assert(sizes_fixed_cumulative.size() > 0,
1619  ExcMessage("No data has been packed!"));
1620  if (quad_cell_relations.size() > 0)
1621  {
1622  Assert(dest_data_fixed.size() > 0,
1623  ExcMessage("No data has been received!"));
1624  }
1625 
1626  std::vector<char>::const_iterator dest_data_it;
1627  std::vector<char>::const_iterator dest_sizes_cell_it;
1628 
1629  // Depending on whether our callback function unpacks fixed or
1630  // variable size data, we have to pursue different approaches
1631  // to localize the correct fraction of the buffer from which
1632  // we are allowed to read.
1633  unsigned int offset = numbers::invalid_unsigned_int;
1634  unsigned int size = numbers::invalid_unsigned_int;
1635  unsigned int data_increment = numbers::invalid_unsigned_int;
1636 
1637  if (callback_variable_transfer)
1638  {
1639  // For the variable size data, we need to extract the
1640  // data size from the fixed size buffer on each cell.
1641  //
1642  // We packed this information last, so the last packed
1643  // object in the fixed size buffer corresponds to the
1644  // variable data sizes.
1645  //
1646  // The last entry of sizes_fixed_cumulative corresponds
1647  // to the size of all fixed size data packed on the cell.
1648  // To get the offset for the last packed object, we need
1649  // to get the next-to-last entry.
1650  const unsigned int offset_variable_data_sizes =
1652 
1653  // This iterator points to the data size that the
1654  // callback_function packed for each specific cell.
1655  // Adjust buffer iterator to the offset of the callback
1656  // function so that we only have to advance its position
1657  // to the next cell after each iteration.
1658  dest_sizes_cell_it = dest_data_fixed.cbegin() +
1659  offset_variable_data_sizes +
1660  callback_index * sizeof(unsigned int);
1661 
1662  // Let the data iterator point to the correct buffer.
1663  dest_data_it = dest_data_variable.cbegin();
1664  }
1665  else
1666  {
1667  // For the fixed size data, we can get the information about
1668  // the buffer location on each cell directly from the
1669  // sizes_fixed_cumulative vector.
1670  offset = sizes_fixed_cumulative[callback_index];
1671  size = sizes_fixed_cumulative[callback_index + 1] - offset;
1672  data_increment = sizes_fixed_cumulative.back();
1673 
1674  // Let the data iterator point to the correct buffer.
1675  // Adjust buffer iterator to the offset of the callback
1676  // function so that we only have to advance its position
1677  // to the next cell after each iteration.
1678  dest_data_it = dest_data_fixed.cbegin() + offset;
1679  }
1680 
1681  // Iterate over all cells and unpack the transferred data.
1682  auto quad_cell_rel_it = quad_cell_relations.begin();
1683  auto dest_sizes_it = dest_sizes_variable.cbegin();
1684  for (; quad_cell_rel_it != quad_cell_relations.end();
1685  ++quad_cell_rel_it, dest_data_it += data_increment)
1686  {
1687  const auto &cell_status = std::get<1>(*quad_cell_rel_it);
1688  const auto &dealii_cell = std::get<2>(*quad_cell_rel_it);
1689 
1690  if (callback_variable_transfer)
1691  {
1692  // Update the increment according to the whole data size
1693  // of the current cell.
1694  data_increment = *dest_sizes_it;
1695 
1696  if (cell_status !=
1698  spacedim>::CELL_INVALID)
1699  {
1700  // Extract the corresponding values for offset and size from
1701  // the cumulative sizes array stored in the fixed size buffer.
1702  if (callback_index == 0)
1703  offset = 0;
1704  else
1705  std::memcpy(&offset,
1706  &(*(dest_sizes_cell_it - sizeof(unsigned int))),
1707  sizeof(unsigned int));
1708 
1709  std::memcpy(&size,
1710  &(*dest_sizes_cell_it),
1711  sizeof(unsigned int));
1712 
1713  size -= offset;
1714 
1715  // Move the data iterator to the corresponding position
1716  // of the callback function and adjust the increment
1717  // accordingly.
1718  dest_data_it += offset;
1719  data_increment -= offset;
1720  }
1721 
1722  // Advance data size iterators to the next cell.
1723  dest_sizes_cell_it += sizes_fixed_cumulative.back();
1724  ++dest_sizes_it;
1725  }
1726 
1727  switch (cell_status)
1728  {
1730  spacedim>::CELL_PERSIST:
1732  spacedim>::CELL_COARSEN:
1733  unpack_callback(dealii_cell,
1734  cell_status,
1735  boost::make_iterator_range(dest_data_it,
1736  dest_data_it +
1737  size));
1738  break;
1739 
1741  spacedim>::CELL_REFINE:
1742  unpack_callback(dealii_cell->parent(),
1743  cell_status,
1744  boost::make_iterator_range(dest_data_it,
1745  dest_data_it +
1746  size));
1747  break;
1748 
1750  spacedim>::CELL_INVALID:
1751  // Skip this cell.
1752  break;
1753 
1754  default:
1755  Assert(false, ExcInternalError());
1756  break;
1757  }
1758  }
1759  }
1760 
1761 
1762 
1763  template <int dim, int spacedim>
1764  void
1766  const typename ::internal::p4est::types<dim>::forest
1767  * parallel_forest,
1768  const std::string &filename) const
1769  {
1770  // Large fractions of this function have been copied from
1771  // DataOutInterface::write_vtu_in_parallel.
1772  // TODO: Write general MPIIO interface.
1773 
1774  Assert(sizes_fixed_cumulative.size() > 0,
1775  ExcMessage("No data has been packed!"));
1776 
1778 
1779  //
1780  // ---------- Fixed size data ----------
1781  //
1782  {
1783  const std::string fname_fixed = std::string(filename) + "_fixed.data";
1784 
1785  MPI_Info info;
1786  int ierr = MPI_Info_create(&info);
1787  AssertThrowMPI(ierr);
1788 
1789  MPI_File fh;
1790  ierr = MPI_File_open(mpi_communicator,
1791  DEAL_II_MPI_CONST_CAST(fname_fixed.c_str()),
1792  MPI_MODE_CREATE | MPI_MODE_WRONLY,
1793  info,
1794  &fh);
1795  AssertThrowMPI(ierr);
1796 
1797  ierr = MPI_File_set_size(fh, 0); // delete the file contents
1798  AssertThrowMPI(ierr);
1799  // this barrier is necessary, because otherwise others might already
1800  // write while one core is still setting the size to zero.
1801  ierr = MPI_Barrier(mpi_communicator);
1802  AssertThrowMPI(ierr);
1803  ierr = MPI_Info_free(&info);
1804  AssertThrowMPI(ierr);
1805  // ------------------
1806 
1807  // Check if number of processors is lined up with p4est partitioning.
1808  Assert(myrank < parallel_forest->mpisize, ExcInternalError());
1809 
1810  // Write cumulative sizes to file.
1811  // Since each processor owns the same information about the data sizes,
1812  // it is sufficient to let only the first processor perform this task.
1813  if (myrank == 0)
1814  {
1815  const unsigned int *data = sizes_fixed_cumulative.data();
1816 
1817  ierr = MPI_File_write_at(fh,
1818  0,
1819  DEAL_II_MPI_CONST_CAST(data),
1820  sizes_fixed_cumulative.size(),
1821  MPI_UNSIGNED,
1822  MPI_STATUS_IGNORE);
1823  AssertThrowMPI(ierr);
1824  }
1825 
1826  // Write packed data to file simultaneously.
1827  const unsigned int offset_fixed =
1828  sizes_fixed_cumulative.size() * sizeof(unsigned int);
1829 
1830  const char *data = src_data_fixed.data();
1831 
1832  ierr = MPI_File_write_at(
1833  fh,
1834  offset_fixed +
1835  parallel_forest->global_first_quadrant[myrank] *
1836  sizes_fixed_cumulative.back(), // global position in file
1837  DEAL_II_MPI_CONST_CAST(data),
1838  src_data_fixed.size(), // local buffer
1839  MPI_CHAR,
1840  MPI_STATUS_IGNORE);
1841  AssertThrowMPI(ierr);
1842 
1843  ierr = MPI_File_close(&fh);
1844  AssertThrowMPI(ierr);
1845  }
1846 
1847  //
1848  // ---------- Variable size data ----------
1849  //
1851  {
1852  const std::string fname_variable =
1853  std::string(filename) + "_variable.data";
1854 
1855  MPI_Info info;
1856  int ierr = MPI_Info_create(&info);
1857  AssertThrowMPI(ierr);
1858 
1859  MPI_File fh;
1860  ierr = MPI_File_open(mpi_communicator,
1861  DEAL_II_MPI_CONST_CAST(fname_variable.c_str()),
1862  MPI_MODE_CREATE | MPI_MODE_WRONLY,
1863  info,
1864  &fh);
1865  AssertThrowMPI(ierr);
1866 
1867  ierr = MPI_File_set_size(fh, 0); // delete the file contents
1868  AssertThrowMPI(ierr);
1869  // this barrier is necessary, because otherwise others might already
1870  // write while one core is still setting the size to zero.
1871  ierr = MPI_Barrier(mpi_communicator);
1872  AssertThrowMPI(ierr);
1873  ierr = MPI_Info_free(&info);
1874  AssertThrowMPI(ierr);
1875 
1876  // Write sizes of each cell into file simultaneously.
1877  {
1878  const int *data = src_sizes_variable.data();
1879  ierr =
1880  MPI_File_write_at(fh,
1881  parallel_forest->global_first_quadrant[myrank] *
1882  sizeof(int), // global position in file
1883  DEAL_II_MPI_CONST_CAST(data),
1884  src_sizes_variable.size(), // local buffer
1885  MPI_INT,
1886  MPI_STATUS_IGNORE);
1887  AssertThrowMPI(ierr);
1888  }
1889 
1890 
1891  const unsigned int offset_variable =
1892  parallel_forest->global_num_quadrants * sizeof(int);
1893 
1894  // Gather size of data in bytes we want to store from this processor.
1895  const unsigned int size_on_proc = src_data_variable.size();
1896 
1897  // Compute prefix sum
1898  unsigned int prefix_sum = 0;
1899  ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST(&size_on_proc),
1900  &prefix_sum,
1901  1,
1902  MPI_UNSIGNED,
1903  MPI_SUM,
1905  AssertThrowMPI(ierr);
1906 
1907  const char *data = src_data_variable.data();
1908 
1909  // Write data consecutively into file.
1910  ierr = MPI_File_write_at(fh,
1911  offset_variable +
1912  prefix_sum, // global position in file
1913  DEAL_II_MPI_CONST_CAST(data),
1914  src_data_variable.size(), // local buffer
1915  MPI_CHAR,
1916  MPI_STATUS_IGNORE);
1917  AssertThrowMPI(ierr);
1918 
1919  ierr = MPI_File_close(&fh);
1920  AssertThrowMPI(ierr);
1921  }
1922  }
1923 
1924 
1925 
1926  template <int dim, int spacedim>
1927  void
1929  const typename ::internal::p4est::types<dim>::forest
1930  * parallel_forest,
1931  const std::string &filename,
1932  const unsigned int n_attached_deserialize_fixed,
1933  const unsigned int n_attached_deserialize_variable)
1934  {
1935  // Large fractions of this function have been copied from
1936  // DataOutInterface::write_vtu_in_parallel.
1937  // TODO: Write general MPIIO interface.
1938 
1939  Assert(dest_data_fixed.size() == 0,
1940  ExcMessage("Previously loaded data has not been released yet!"));
1941 
1942  variable_size_data_stored = (n_attached_deserialize_variable > 0);
1943 
1945 
1946  //
1947  // ---------- Fixed size data ----------
1948  //
1949  {
1950  const std::string fname_fixed = std::string(filename) + "_fixed.data";
1951 
1952  MPI_Info info;
1953  int ierr = MPI_Info_create(&info);
1954  AssertThrowMPI(ierr);
1955 
1956  MPI_File fh;
1957  ierr = MPI_File_open(mpi_communicator,
1958  DEAL_II_MPI_CONST_CAST(fname_fixed.c_str()),
1959  MPI_MODE_RDONLY,
1960  info,
1961  &fh);
1962  AssertThrowMPI(ierr);
1963 
1964  ierr = MPI_Info_free(&info);
1965  AssertThrowMPI(ierr);
1966 
1967  // Check if number of processors is lined up with p4est partitioning.
1968  Assert(myrank < parallel_forest->mpisize, ExcInternalError());
1969 
1970  // Read cumulative sizes from file.
1971  // Since all processors need the same information about the data sizes,
1972  // let each of them retrieve it by reading from the same location in
1973  // the file.
1974  sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
1975  (variable_size_data_stored ? 1 : 0));
1976  ierr = MPI_File_read_at(fh,
1977  0,
1978  sizes_fixed_cumulative.data(),
1979  sizes_fixed_cumulative.size(),
1980  MPI_UNSIGNED,
1981  MPI_STATUS_IGNORE);
1982  AssertThrowMPI(ierr);
1983 
1984  // Allocate sufficient memory.
1985  dest_data_fixed.resize(parallel_forest->local_num_quadrants *
1986  sizes_fixed_cumulative.back());
1987 
1988  // Read packed data from file simultaneously.
1989  const unsigned int offset =
1990  sizes_fixed_cumulative.size() * sizeof(unsigned int);
1991 
1992  ierr = MPI_File_read_at(
1993  fh,
1994  offset + parallel_forest->global_first_quadrant[myrank] *
1995  sizes_fixed_cumulative.back(), // global position in file
1996  dest_data_fixed.data(),
1997  dest_data_fixed.size(), // local buffer
1998  MPI_CHAR,
1999  MPI_STATUS_IGNORE);
2000  AssertThrowMPI(ierr);
2001 
2002  ierr = MPI_File_close(&fh);
2003  AssertThrowMPI(ierr);
2004  }
2005 
2006  //
2007  // ---------- Variable size data ----------
2008  //
2009  if (variable_size_data_stored)
2010  {
2011  const std::string fname_variable =
2012  std::string(filename) + "_variable.data";
2013 
2014  MPI_Info info;
2015  int ierr = MPI_Info_create(&info);
2016  AssertThrowMPI(ierr);
2017 
2018  MPI_File fh;
2019  ierr = MPI_File_open(mpi_communicator,
2020  DEAL_II_MPI_CONST_CAST(fname_variable.c_str()),
2021  MPI_MODE_RDONLY,
2022  info,
2023  &fh);
2024  AssertThrowMPI(ierr);
2025 
2026  ierr = MPI_Info_free(&info);
2027  AssertThrowMPI(ierr);
2028 
2029  // Read sizes of all locally owned cells.
2030  dest_sizes_variable.resize(parallel_forest->local_num_quadrants);
2031  ierr =
2032  MPI_File_read_at(fh,
2033  parallel_forest->global_first_quadrant[myrank] *
2034  sizeof(int),
2035  dest_sizes_variable.data(),
2036  dest_sizes_variable.size(),
2037  MPI_INT,
2038  MPI_STATUS_IGNORE);
2039  AssertThrowMPI(ierr);
2040 
2041  const unsigned int offset =
2042  parallel_forest->global_num_quadrants * sizeof(int);
2043 
2044  const unsigned int size_on_proc =
2045  std::accumulate(dest_sizes_variable.begin(),
2046  dest_sizes_variable.end(),
2047  0);
2048 
2049  // share information among all processors by prefix sum
2050  unsigned int prefix_sum = 0;
2051  ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST(&size_on_proc),
2052  &prefix_sum,
2053  1,
2054  MPI_UNSIGNED,
2055  MPI_SUM,
2057  AssertThrowMPI(ierr);
2058 
2059  dest_data_variable.resize(size_on_proc);
2060  ierr = MPI_File_read_at(fh,
2061  offset + prefix_sum,
2062  dest_data_variable.data(),
2063  dest_data_variable.size(),
2064  MPI_CHAR,
2065  MPI_STATUS_IGNORE);
2066  AssertThrowMPI(ierr);
2067 
2068  ierr = MPI_File_close(&fh);
2069  AssertThrowMPI(ierr);
2070  }
2071  }
2072 
2073 
2074 
2075  template <int dim, int spacedim>
2076  void
2078  {
2079  variable_size_data_stored = false;
2080 
2081  // free information about data sizes
2082  sizes_fixed_cumulative.clear();
2083  sizes_fixed_cumulative.shrink_to_fit();
2084 
2085  // free fixed size transfer data
2086  src_data_fixed.clear();
2087  src_data_fixed.shrink_to_fit();
2088 
2089  dest_data_fixed.clear();
2090  dest_data_fixed.shrink_to_fit();
2091 
2092  // free variable size transfer data
2093  src_sizes_variable.clear();
2094  src_sizes_variable.shrink_to_fit();
2095 
2096  src_data_variable.clear();
2097  src_data_variable.shrink_to_fit();
2098 
2099  dest_sizes_variable.clear();
2100  dest_sizes_variable.shrink_to_fit();
2101 
2102  dest_data_variable.clear();
2103  dest_data_variable.shrink_to_fit();
2104  }
2105 
2106 
2107 
2108  /* ----------------- class Triangulation<dim,spacedim> ----------------- */
2109 
2110 
2111  template <int dim, int spacedim>
2113  MPI_Comm mpi_communicator,
2114  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
2115  smooth_grid,
2116  const Settings settings)
2117  : // Do not check for distorted cells.
2118  // For multigrid, we need limit_level_difference_at_vertices
2119  // to make sure the transfer operators only need to consider two levels.
2120  ::parallel::DistributedTriangulationBase<dim, spacedim>(
2121  mpi_communicator,
2122  (settings & construct_multigrid_hierarchy) ?
2123  static_cast<
2124  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
2125  smooth_grid |
2127  smooth_grid,
2128  false)
2129  , settings(settings)
2130  , triangulation_has_content(false)
2131  , connectivity(nullptr)
2132  , parallel_forest(nullptr)
2133  , cell_attached_data({0, 0, {}, {}})
2134  , data_transfer(mpi_communicator)
2135  {
2136  parallel_ghost = nullptr;
2137  }
2138 
2139 
2140 
2141  template <int dim, int spacedim>
2143  {
2144  // virtual functions called in constructors and destructors never use the
2145  // override in a derived class
2146  // for clarity be explicit on which function is called
2147  try
2148  {
2150  }
2151  catch (...)
2152  {}
2153 
2154  AssertNothrow(triangulation_has_content == false, ExcInternalError());
2155  AssertNothrow(connectivity == nullptr, ExcInternalError());
2156  AssertNothrow(parallel_forest == nullptr, ExcInternalError());
2157  }
2158 
2159 
2160 
2161  template <int dim, int spacedim>
2162  void
2164  const std::vector<Point<spacedim>> &vertices,
2165  const std::vector<CellData<dim>> & cells,
2166  const SubCellData & subcelldata)
2167  {
2168  try
2169  {
2171  vertices, cells, subcelldata);
2172  }
2173  catch (
2174  const typename ::Triangulation<dim, spacedim>::DistortedCellList
2175  &)
2176  {
2177  // the underlying triangulation should not be checking for distorted
2178  // cells
2179  Assert(false, ExcInternalError());
2180  }
2181 
2182  // note that now we have some content in the p4est objects and call the
2183  // functions that do the actual work (which are dimension dependent, so
2184  // separate)
2185  triangulation_has_content = true;
2186 
2187  setup_coarse_cell_to_p4est_tree_permutation();
2188 
2189  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
2190 
2191  try
2192  {
2193  copy_local_forest_to_triangulation();
2194  }
2195  catch (const typename Triangulation<dim>::DistortedCellList &)
2196  {
2197  // the underlying triangulation should not be checking for distorted
2198  // cells
2199  Assert(false, ExcInternalError());
2200  }
2201 
2202  this->update_periodic_face_map();
2203  this->update_number_cache();
2204  }
2205 
2206 
2207 
2208  template <int dim, int spacedim>
2209  void
2212  &construction_data)
2213  {
2214  (void)construction_data;
2215 
2216  Assert(false, ExcInternalError());
2217  }
2218 
2219 
2220 
2221  // This anonymous namespace contains utility for
2222  // the function Triangulation::communicate_locally_moved_vertices
2223  namespace CommunicateLocallyMovedVertices
2224  {
2225  namespace
2226  {
2232  template <int dim, int spacedim>
2233  struct CellInfo
2234  {
2235  // store all the tree_indices we send/receive consecutively (n_cells
2236  // entries)
2237  std::vector<unsigned int> tree_index;
2238  // store all the quadrants we send/receive consecutively (n_cells
2239  // entries)
2240  std::vector<typename ::internal::p4est::types<dim>::quadrant>
2242  // store for each cell the number of vertices we send/receive
2243  // and then the vertex indices (for each cell: n_vertices+1 entries)
2244  std::vector<unsigned int> vertex_indices;
2245  // store for each cell the vertices we send/receive
2246  // (for each cell n_vertices entries)
2247  std::vector<::Point<spacedim>> vertices;
2248  // for receiving and unpacking data we need to store pointers to the
2249  // first vertex and vertex_index on each cell additionally
2250  // both vectors have as many entries as there are cells
2251  std::vector<unsigned int *> first_vertex_indices;
2252  std::vector<::Point<spacedim> *> first_vertices;
2253 
2254  unsigned int
2255  bytes_for_buffer() const
2256  {
2257  return sizeof(unsigned int) +
2258  tree_index.size() * sizeof(unsigned int) +
2259  quadrants.size() *
2260  sizeof(
2261  typename ::internal::p4est::types<dim>::quadrant) +
2262  vertex_indices.size() * sizeof(unsigned int) +
2263  vertices.size() * sizeof(::Point<spacedim>);
2264  }
2265 
2266  void
2267  pack_data(std::vector<char> &buffer) const
2268  {
2269  buffer.resize(bytes_for_buffer());
2270 
2271  char *ptr = buffer.data();
2272 
2273  const unsigned int num_cells = tree_index.size();
2274  std::memcpy(ptr, &num_cells, sizeof(unsigned int));
2275  ptr += sizeof(unsigned int);
2276 
2277  std::memcpy(ptr,
2278  tree_index.data(),
2279  num_cells * sizeof(unsigned int));
2280  ptr += num_cells * sizeof(unsigned int);
2281 
2282  std::memcpy(
2283  ptr,
2284  quadrants.data(),
2285  num_cells *
2286  sizeof(typename ::internal::p4est::types<dim>::quadrant));
2287  ptr +=
2288  num_cells *
2289  sizeof(typename ::internal::p4est::types<dim>::quadrant);
2290 
2291  std::memcpy(ptr,
2292  vertex_indices.data(),
2293  vertex_indices.size() * sizeof(unsigned int));
2294  ptr += vertex_indices.size() * sizeof(unsigned int);
2295  std::memcpy(ptr,
2296  vertices.data(),
2297  vertices.size() * sizeof(::Point<spacedim>));
2298  ptr += vertices.size() * sizeof(::Point<spacedim>);
2299 
2300  Assert(ptr == buffer.data() + buffer.size(), ExcInternalError());
2301  }
2302 
2303  void
2304  unpack_data(const std::vector<char> &buffer)
2305  {
2306  const char * ptr = buffer.data();
2307  unsigned int cells;
2308  memcpy(&cells, ptr, sizeof(unsigned int));
2309  ptr += sizeof(unsigned int);
2310 
2311  tree_index.resize(cells);
2312  memcpy(tree_index.data(), ptr, sizeof(unsigned int) * cells);
2313  ptr += sizeof(unsigned int) * cells;
2314 
2315  quadrants.resize(cells);
2316  memcpy(quadrants.data(),
2317  ptr,
2318  sizeof(
2319  typename ::internal::p4est::types<dim>::quadrant) *
2320  cells);
2321  ptr +=
2322  sizeof(typename ::internal::p4est::types<dim>::quadrant) *
2323  cells;
2324 
2325  vertex_indices.clear();
2326  first_vertex_indices.resize(cells);
2327  std::vector<unsigned int> n_vertices_on_cell(cells);
2328  std::vector<unsigned int> first_indices(cells);
2329  for (unsigned int c = 0; c < cells; ++c)
2330  {
2331  // The first 'vertex index' is the number of vertices.
2332  // Additionally, we need to store the pointer to this
2333  // vertex index with respect to the std::vector
2334  const unsigned int *const vertex_index =
2335  reinterpret_cast<const unsigned int *>(ptr);
2336  first_indices[c] = vertex_indices.size();
2337  vertex_indices.push_back(*vertex_index);
2338  n_vertices_on_cell[c] = *vertex_index;
2339  ptr += sizeof(unsigned int);
2340  // Now copy all the 'real' vertex_indices
2341  vertex_indices.resize(vertex_indices.size() +
2342  n_vertices_on_cell[c]);
2343  memcpy(&vertex_indices[vertex_indices.size() -
2344  n_vertices_on_cell[c]],
2345  ptr,
2346  n_vertices_on_cell[c] * sizeof(unsigned int));
2347  ptr += n_vertices_on_cell[c] * sizeof(unsigned int);
2348  }
2349  for (unsigned int c = 0; c < cells; ++c)
2350  first_vertex_indices[c] = &vertex_indices[first_indices[c]];
2351 
2352  vertices.clear();
2353  first_vertices.resize(cells);
2354  for (unsigned int c = 0; c < cells; ++c)
2355  {
2356  first_indices[c] = vertices.size();
2357  vertices.resize(vertices.size() + n_vertices_on_cell[c]);
2358  memcpy(&vertices[vertices.size() - n_vertices_on_cell[c]],
2359  ptr,
2360  n_vertices_on_cell[c] * sizeof(::Point<spacedim>));
2361  ptr += n_vertices_on_cell[c] * sizeof(::Point<spacedim>);
2362  }
2363  for (unsigned int c = 0; c < cells; ++c)
2364  first_vertices[c] = &vertices[first_indices[c]];
2365 
2366  Assert(ptr == buffer.data() + buffer.size(), ExcInternalError());
2367  }
2368  };
2369 
2370 
2371 
2372  // This function is responsible for gathering the information
2373  // we want to send to each process.
2374  // For each dealii cell on the coarsest level the corresponding
2375  // p4est_cell has to be provided when calling this function.
2376  // By recursing through all children we consider each active cell.
2377  // vertices_with_ghost_neighbors tells us which vertices
2378  // are in the ghost layer and for which processes they might
2379  // be interesting.
2380  // Whether a vertex has actually been updated locally is
2381  // stored in vertex_locally_moved. Only those are considered
2382  // for sending.
2383  // The gathered information is saved into needs_to_get_cell.
2384  template <int dim, int spacedim>
2385  void
2386  fill_vertices_recursively(
2388  & tria,
2389  const unsigned int tree_index,
2391  &dealii_cell,
2392  const typename ::internal::p4est::types<dim>::quadrant
2393  &p4est_cell,
2394  const std::map<unsigned int, std::set<::types::subdomain_id>>
2396  const std::vector<bool> &vertex_locally_moved,
2397  std::map<::types::subdomain_id, CellInfo<dim, spacedim>>
2398  &needs_to_get_cell)
2399  {
2400  // see if we have to
2401  // recurse...
2402  if (dealii_cell->has_children())
2403  {
2404  typename ::internal::p4est::types<dim>::quadrant
2406  ::internal::p4est::init_quadrant_children<dim>(p4est_cell,
2407  p4est_child);
2408 
2409 
2410  for (unsigned int c = 0;
2411  c < GeometryInfo<dim>::max_children_per_cell;
2412  ++c)
2413  fill_vertices_recursively<dim, spacedim>(
2414  tria,
2415  tree_index,
2416  dealii_cell->child(c),
2417  p4est_child[c],
2419  vertex_locally_moved,
2420  needs_to_get_cell);
2421  return;
2422  }
2423 
2424  // We're at a leaf cell. If the cell is locally owned, we may
2425  // have to send its vertices to other processors if any of
2426  // its vertices is adjacent to a ghost cell and has been moved
2427  //
2428  // If one of the vertices of the cell is interesting,
2429  // send all moved vertices of the cell to all processors
2430  // adjacent to all cells adjacent to this vertex
2431  if (dealii_cell->is_locally_owned())
2432  {
2433  std::set<::types::subdomain_id> send_to;
2434  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2435  {
2436  const std::map<unsigned int,
2437  std::set<::types::subdomain_id>>::
2438  const_iterator neighbor_subdomains_of_vertex =
2440  dealii_cell->vertex_index(v));
2441 
2442  if (neighbor_subdomains_of_vertex !=
2444  {
2445  Assert(neighbor_subdomains_of_vertex->second.size() != 0,
2446  ExcInternalError());
2447  send_to.insert(
2448  neighbor_subdomains_of_vertex->second.begin(),
2449  neighbor_subdomains_of_vertex->second.end());
2450  }
2451  }
2452 
2453  if (send_to.size() > 0)
2454  {
2455  std::vector<unsigned int> vertex_indices;
2456  std::vector<::Point<spacedim>> local_vertices;
2457  for (const unsigned int v :
2459  if (vertex_locally_moved[dealii_cell->vertex_index(v)])
2460  {
2461  vertex_indices.push_back(v);
2462  local_vertices.push_back(dealii_cell->vertex(v));
2463  }
2464 
2465  if (vertex_indices.size() > 0)
2466  for (const auto subdomain : send_to)
2467  {
2468  // get an iterator to what needs to be sent to that
2469  // subdomain (if already exists), or create such an
2470  // object
2471  const typename std::map<
2473  CellInfo<dim, spacedim>>::iterator p =
2474  needs_to_get_cell
2475  .insert(std::make_pair(subdomain,
2476  CellInfo<dim, spacedim>()))
2477  .first;
2478 
2479  p->second.tree_index.push_back(tree_index);
2480  p->second.quadrants.push_back(p4est_cell);
2481 
2482  p->second.vertex_indices.push_back(
2483  vertex_indices.size());
2484  p->second.vertex_indices.insert(
2485  p->second.vertex_indices.end(),
2486  vertex_indices.begin(),
2487  vertex_indices.end());
2488 
2489  p->second.vertices.insert(p->second.vertices.end(),
2490  local_vertices.begin(),
2491  local_vertices.end());
2492  }
2493  }
2494  }
2495  }
2496 
2497 
2498 
2499  // After the cell data has been received this function is responsible
2500  // for moving the vertices in the corresponding ghost layer locally.
2501  // As in fill_vertices_recursively for each dealii cell on the
2502  // coarsest level the corresponding p4est_cell has to be provided
2503  // when calling this function. By recursing through through all
2504  // children we consider each active cell.
2505  // Additionally, we need to give a pointer to the first vertex indices
2506  // and vertices. Since the first information saved in vertex_indices
2507  // is the number of vertices this all the information we need.
2508  template <int dim, int spacedim>
2509  void
2510  set_vertices_recursively(
2512  const typename ::internal::p4est::types<dim>::quadrant
2513  &p4est_cell,
2515  &dealii_cell,
2516  const typename ::internal::p4est::types<dim>::quadrant
2517  & quadrant,
2518  const ::Point<spacedim> *const vertices,
2519  const unsigned int *const vertex_indices)
2520  {
2521  if (::internal::p4est::quadrant_is_equal<dim>(p4est_cell,
2522  quadrant))
2523  {
2524  Assert(!dealii_cell->is_artificial(), ExcInternalError());
2525  Assert(dealii_cell->is_active(), ExcInternalError());
2526  Assert(!dealii_cell->is_locally_owned(), ExcInternalError());
2527 
2528  const unsigned int n_vertices = vertex_indices[0];
2529 
2530  // update dof indices of cell
2531  for (unsigned int i = 0; i < n_vertices; ++i)
2532  dealii_cell->vertex(vertex_indices[i + 1]) = vertices[i];
2533 
2534  return;
2535  }
2536 
2537  if (dealii_cell->is_active())
2538  return;
2539 
2540  if (!::internal::p4est::quadrant_is_ancestor<dim>(p4est_cell,
2541  quadrant))
2542  return;
2543 
2544  typename ::internal::p4est::types<dim>::quadrant
2546  ::internal::p4est::init_quadrant_children<dim>(p4est_cell,
2547  p4est_child);
2548 
2549  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
2550  ++c)
2551  set_vertices_recursively<dim, spacedim>(tria,
2552  p4est_child[c],
2553  dealii_cell->child(c),
2554  quadrant,
2555  vertices,
2556  vertex_indices);
2557  }
2558  } // namespace
2559  } // namespace CommunicateLocallyMovedVertices
2560 
2561 
2562 
2563  template <int dim, int spacedim>
2564  void
2566  {
2567  triangulation_has_content = false;
2568 
2569  cell_attached_data = {0, 0, {}, {}};
2570  data_transfer.clear();
2571 
2572  if (parallel_ghost != nullptr)
2573  {
2575  parallel_ghost);
2576  parallel_ghost = nullptr;
2577  }
2578 
2579  if (parallel_forest != nullptr)
2580  {
2582  parallel_forest = nullptr;
2583  }
2584 
2585  if (connectivity != nullptr)
2586  {
2588  connectivity);
2589  connectivity = nullptr;
2590  }
2591 
2592  coarse_cell_to_p4est_tree_permutation.resize(0);
2593  p4est_tree_to_coarse_cell_permutation.resize(0);
2594 
2596 
2597  this->update_number_cache();
2598  }
2599 
2600 
2601 
2602  template <int dim, int spacedim>
2603  bool
2605  {
2606  return settings &
2608  }
2609 
2610 
2611 
2612  template <int dim, int spacedim>
2613  bool
2615  {
2616  if (this->n_global_levels() <= 1)
2617  return false; // can not have hanging nodes without refined cells
2618 
2619  // if there are any active cells with level less than n_global_levels()-1,
2620  // then there is obviously also one with level n_global_levels()-1, and
2621  // consequently there must be a hanging node somewhere.
2622  //
2623  // The problem is that we cannot just ask for the first active cell, but
2624  // instead need to filter over locally owned cells.
2625  const bool have_coarser_cell =
2626  std::any_of(this->begin_active(this->n_global_levels() - 2),
2627  this->end_active(this->n_global_levels() - 2),
2628  [](const CellAccessor<dim, spacedim> &cell) {
2629  return cell.is_locally_owned();
2630  });
2631 
2632  // return true if at least one process has a coarser cell
2633  return 0 < Utilities::MPI::max(have_coarser_cell ? 1 : 0,
2634  this->mpi_communicator);
2635  }
2636 
2637 
2638 
2639  template <int dim, int spacedim>
2640  void
2642  {
2643  DynamicSparsityPattern cell_connectivity;
2645  cell_connectivity);
2646  coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
2648  cell_connectivity, coarse_cell_to_p4est_tree_permutation);
2649 
2650  p4est_tree_to_coarse_cell_permutation =
2651  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
2652  }
2653 
2654 
2655 
2656  template <int dim, int spacedim>
2657  void
2659  const std::string &file_basename) const
2660  {
2661  Assert(parallel_forest != nullptr,
2662  ExcMessage("Can't produce output when no forest is created yet."));
2664  parallel_forest, nullptr, file_basename.c_str());
2665  }
2666 
2667 
2668 
2669  template <int dim, int spacedim>
2670  void
2671  Triangulation<dim, spacedim>::save(const std::string &filename) const
2672  {
2673  Assert(
2674  cell_attached_data.n_attached_deserialize == 0,
2675  ExcMessage(
2676  "not all SolutionTransfer's got deserialized after the last load()"));
2677  Assert(this->n_cells() > 0,
2678  ExcMessage("Can not save() an empty Triangulation."));
2679 
2680  // signal that serialization is going to happen
2681  this->signals.pre_distributed_save();
2682 
2683  if (this->my_subdomain == 0)
2684  {
2685  std::string fname = std::string(filename) + ".info";
2686  std::ofstream f(fname.c_str());
2687  f << "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2688  << std::endl
2689  << 4 << " "
2690  << Utilities::MPI::n_mpi_processes(this->mpi_communicator) << " "
2691  << cell_attached_data.pack_callbacks_fixed.size() << " "
2692  << cell_attached_data.pack_callbacks_variable.size() << " "
2693  << this->n_cells(0) << std::endl;
2694  }
2695 
2696  // each cell should have been flagged `CELL_PERSIST`
2697  for (const auto &quad_cell_rel : local_quadrant_cell_relations)
2698  {
2699  (void)quad_cell_rel;
2700  Assert(
2701  (std::get<1>(quad_cell_rel) == // cell_status
2703  ExcInternalError());
2704  }
2705 
2706  if (cell_attached_data.n_attached_data_sets > 0)
2707  {
2708  // cast away constness
2709  auto tria = const_cast<
2711  this);
2712 
2713  // pack attached data first
2714  tria->data_transfer.pack_data(
2715  local_quadrant_cell_relations,
2716  cell_attached_data.pack_callbacks_fixed,
2717  cell_attached_data.pack_callbacks_variable);
2718 
2719  // then store buffers in file
2720  tria->data_transfer.save(parallel_forest, filename);
2721 
2722  // and release the memory afterwards
2723  tria->data_transfer.clear();
2724  }
2725 
2726  ::internal::p4est::functions<dim>::save(filename.c_str(),
2727  parallel_forest,
2728  false);
2729 
2730  // clear all of the callback data, as explained in the documentation of
2731  // register_data_attach()
2732  {
2734  const_cast<
2736  this);
2737 
2738  tria->cell_attached_data.n_attached_data_sets = 0;
2739  tria->cell_attached_data.pack_callbacks_fixed.clear();
2740  tria->cell_attached_data.pack_callbacks_variable.clear();
2741  }
2742 
2743  // signal that serialization has finished
2745  }
2746 
2747 
2748 
2749  template <int dim, int spacedim>
2750  void
2751  Triangulation<dim, spacedim>::load(const std::string &filename,
2752  const bool autopartition)
2753  {
2754  Assert(
2755  this->n_cells() > 0,
2756  ExcMessage(
2757  "load() only works if the Triangulation already contains a coarse mesh!"));
2758  Assert(
2759  this->n_levels() == 1,
2760  ExcMessage(
2761  "Triangulation may only contain coarse cells when calling load()."));
2762 
2763  // signal that de-serialization is going to happen
2764  this->signals.pre_distributed_load();
2765 
2766  if (parallel_ghost != nullptr)
2767  {
2769  parallel_ghost);
2770  parallel_ghost = nullptr;
2771  }
2773  parallel_forest = nullptr;
2775  connectivity);
2776  connectivity = nullptr;
2777 
2778  unsigned int version, numcpus, attached_count_fixed,
2779  attached_count_variable, n_coarse_cells;
2780  {
2781  std::string fname = std::string(filename) + ".info";
2782  std::ifstream f(fname.c_str());
2783  AssertThrow(f, ExcIO());
2784  std::string firstline;
2785  getline(f, firstline); // skip first line
2786  f >> version >> numcpus >> attached_count_fixed >>
2787  attached_count_variable >> n_coarse_cells;
2788  }
2789 
2790  AssertThrow(version == 4,
2791  ExcMessage("Incompatible version found in .info file."));
2792  Assert(this->n_cells(0) == n_coarse_cells,
2793  ExcMessage("Number of coarse cells differ!"));
2794 
2795  // clear all of the callback data, as explained in the documentation of
2796  // register_data_attach()
2797  cell_attached_data.n_attached_data_sets = 0;
2798  cell_attached_data.n_attached_deserialize =
2799  attached_count_fixed + attached_count_variable;
2800 
2802  filename.c_str(),
2803  this->mpi_communicator,
2804  0,
2805  false,
2806  autopartition,
2807  0,
2808  this,
2809  &connectivity);
2810 
2811  if (numcpus != Utilities::MPI::n_mpi_processes(this->mpi_communicator))
2812  {
2813  // We are changing the number of CPUs so we need to repartition.
2814  // Note that p4est actually distributes the cells between the changed
2815  // number of CPUs and so everything works without this call, but
2816  // this command changes the distribution for some reason, so we
2817  // will leave it in here.
2818  if (this->signals.cell_weight.num_slots() == 0)
2819  {
2820  // no cell weights given -- call p4est's 'partition' without a
2821  // callback for cell weights
2823  parallel_forest,
2824  /* prepare coarsening */ 1,
2825  /* weight_callback */ nullptr);
2826  }
2827  else
2828  {
2829  // get cell weights for a weighted repartitioning.
2830  const std::vector<unsigned int> cell_weights = get_cell_weights();
2831 
2832  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
2833 
2834  // attach (temporarily) a pointer to the cell weights through
2835  // p4est's user_pointer object
2836  Assert(parallel_forest->user_pointer == this, ExcInternalError());
2837  parallel_forest->user_pointer = &partition_weights;
2838 
2840  parallel_forest,
2841  /* prepare coarsening */ 1,
2842  /* weight_callback */
2843  &PartitionWeights<dim, spacedim>::cell_weight);
2844 
2845  // reset the user pointer to its previous state
2846  parallel_forest->user_pointer = this;
2847  }
2848  }
2849 
2850  try
2851  {
2852  copy_local_forest_to_triangulation();
2853  }
2854  catch (const typename Triangulation<dim>::DistortedCellList &)
2855  {
2856  // the underlying
2857  // triangulation should not
2858  // be checking for
2859  // distorted cells
2860  Assert(false, ExcInternalError());
2861  }
2862 
2863  // load saved data, if any was stored
2864  if (cell_attached_data.n_attached_deserialize > 0)
2865  {
2866  data_transfer.load(parallel_forest,
2867  filename,
2868  attached_count_fixed,
2869  attached_count_variable);
2870 
2871  data_transfer.unpack_cell_status(local_quadrant_cell_relations);
2872 
2873  // the CellStatus of all stored cells should always be CELL_PERSIST.
2874  for (const auto &quad_cell_rel : local_quadrant_cell_relations)
2875  {
2876  (void)quad_cell_rel;
2877  Assert(
2878  (std::get<1>(quad_cell_rel) == // cell_status
2880  spacedim>::CELL_PERSIST),
2881  ExcInternalError());
2882  }
2883  }
2884 
2885  this->update_periodic_face_map();
2886  this->update_number_cache();
2887 
2888  // signal that de-serialization is finished
2890  }
2891 
2892 
2893 
2894  template <int dim, int spacedim>
2895  unsigned int
2897  {
2898  Assert(parallel_forest != nullptr,
2899  ExcMessage(
2900  "Can't produce a check sum when no forest is created yet."));
2901  return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2902  }
2903 
2904 
2905 
2906  template <int dim, int spacedim>
2907  const typename ::internal::p4est::types<dim>::forest *
2909  {
2910  Assert(parallel_forest != nullptr,
2911  ExcMessage("The forest has not been allocated yet."));
2912  return parallel_forest;
2913  }
2914 
2915 
2916 
2917  template <int dim, int spacedim>
2918  typename ::internal::p4est::types<dim>::tree *
2920  const int dealii_coarse_cell_index) const
2921  {
2922  const unsigned int tree_index =
2923  coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2924  typename ::internal::p4est::types<dim>::tree *tree =
2925  static_cast<typename ::internal::p4est::types<dim>::tree *>(
2926  sc_array_index(parallel_forest->trees, tree_index));
2927 
2928  return tree;
2929  }
2930 
2931 
2932 
2933  // Note: this has been added here to prevent that these functions
2934  // appear in the Doxygen documentation of ::Triangulation
2935 # ifndef DOXYGEN
2936 
2937  template <>
2938  void
2940  std::integral_constant<int, 2>)
2941  {
2942  const unsigned int dim = 2, spacedim = 2;
2943  Assert(this->n_cells(0) > 0, ExcInternalError());
2944  Assert(this->n_levels() == 1, ExcInternalError());
2945 
2946  // data structures that counts how many cells touch each vertex
2947  // (vertex_touch_count), and which cells touch a given vertex (together
2948  // with the local numbering of that vertex within the cells that touch
2949  // it)
2950  std::vector<unsigned int> vertex_touch_count;
2951  std::vector<
2952  std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2953  unsigned int>>>
2954  vertex_to_cell;
2955  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
2956  const ::internal::p4est::types<2>::locidx num_vtt =
2957  std::accumulate(vertex_touch_count.begin(),
2958  vertex_touch_count.end(),
2959  0u);
2960 
2961  // now create a connectivity object with the right sizes for all
2962  // arrays. set vertex information only in debug mode (saves a few bytes
2963  // in optimized mode)
2964  const bool set_vertex_info
2965 # ifdef DEBUG
2966  = true
2967 # else
2968  = false
2969 # endif
2970  ;
2971 
2973  (set_vertex_info == true ? this->n_vertices() : 0),
2974  this->n_cells(0),
2975  this->n_vertices(),
2976  num_vtt);
2977 
2978  set_vertex_and_cell_info(*this,
2979  vertex_touch_count,
2980  vertex_to_cell,
2981  coarse_cell_to_p4est_tree_permutation,
2982  set_vertex_info,
2983  connectivity);
2984 
2985  Assert(p4est_connectivity_is_valid(connectivity) == 1,
2986  ExcInternalError());
2987 
2988  // now create a forest out of the connectivity data structure
2990  this->mpi_communicator,
2991  connectivity,
2992  /* minimum initial number of quadrants per tree */ 0,
2993  /* minimum level of upfront refinement */ 0,
2994  /* use uniform upfront refinement */ 1,
2995  /* user_data_size = */ 0,
2996  /* user_data_constructor = */ nullptr,
2997  /* user_pointer */ this);
2998  }
2999 
3000 
3001 
3002  // TODO: This is a verbatim copy of the 2,2 case. However, we can't just
3003  // specialize the dim template argument, but let spacedim open
3004  template <>
3005  void
3007  std::integral_constant<int, 2>)
3008  {
3009  const unsigned int dim = 2, spacedim = 3;
3010  Assert(this->n_cells(0) > 0, ExcInternalError());
3011  Assert(this->n_levels() == 1, ExcInternalError());
3012 
3013  // data structures that counts how many cells touch each vertex
3014  // (vertex_touch_count), and which cells touch a given vertex (together
3015  // with the local numbering of that vertex within the cells that touch
3016  // it)
3017  std::vector<unsigned int> vertex_touch_count;
3018  std::vector<
3019  std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
3020  unsigned int>>>
3021  vertex_to_cell;
3022  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
3023  const ::internal::p4est::types<2>::locidx num_vtt =
3024  std::accumulate(vertex_touch_count.begin(),
3025  vertex_touch_count.end(),
3026  0u);
3027 
3028  // now create a connectivity object with the right sizes for all
3029  // arrays. set vertex information only in debug mode (saves a few bytes
3030  // in optimized mode)
3031  const bool set_vertex_info
3032 # ifdef DEBUG
3033  = true
3034 # else
3035  = false
3036 # endif
3037  ;
3038 
3040  (set_vertex_info == true ? this->n_vertices() : 0),
3041  this->n_cells(0),
3042  this->n_vertices(),
3043  num_vtt);
3044 
3045  set_vertex_and_cell_info(*this,
3046  vertex_touch_count,
3047  vertex_to_cell,
3048  coarse_cell_to_p4est_tree_permutation,
3049  set_vertex_info,
3050  connectivity);
3051 
3052  Assert(p4est_connectivity_is_valid(connectivity) == 1,
3053  ExcInternalError());
3054 
3055  // now create a forest out of the connectivity data structure
3057  this->mpi_communicator,
3058  connectivity,
3059  /* minimum initial number of quadrants per tree */ 0,
3060  /* minimum level of upfront refinement */ 0,
3061  /* use uniform upfront refinement */ 1,
3062  /* user_data_size = */ 0,
3063  /* user_data_constructor = */ nullptr,
3064  /* user_pointer */ this);
3065  }
3066 
3067 
3068 
3069  template <>
3070  void
3072  std::integral_constant<int, 3>)
3073  {
3074  const int dim = 3, spacedim = 3;
3075  Assert(this->n_cells(0) > 0, ExcInternalError());
3076  Assert(this->n_levels() == 1, ExcInternalError());
3077 
3078  // data structures that counts how many cells touch each vertex
3079  // (vertex_touch_count), and which cells touch a given vertex (together
3080  // with the local numbering of that vertex within the cells that touch
3081  // it)
3082  std::vector<unsigned int> vertex_touch_count;
3083  std::vector<std::list<
3084  std::pair<Triangulation<3>::active_cell_iterator, unsigned int>>>
3085  vertex_to_cell;
3086  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
3087  const ::internal::p4est::types<2>::locidx num_vtt =
3088  std::accumulate(vertex_touch_count.begin(),
3089  vertex_touch_count.end(),
3090  0u);
3091 
3092  std::vector<unsigned int> edge_touch_count;
3093  std::vector<std::list<
3094  std::pair<Triangulation<3>::active_cell_iterator, unsigned int>>>
3095  edge_to_cell;
3096  get_edge_to_cell_mappings(*this, edge_touch_count, edge_to_cell);
3097  const ::internal::p4est::types<2>::locidx num_ett =
3098  std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
3099 
3100  // now create a connectivity object with the right sizes for all arrays
3101  const bool set_vertex_info
3102 # ifdef DEBUG
3103  = true
3104 # else
3105  = false
3106 # endif
3107  ;
3108 
3110  (set_vertex_info == true ? this->n_vertices() : 0),
3111  this->n_cells(0),
3112  this->n_active_lines(),
3113  num_ett,
3114  this->n_vertices(),
3115  num_vtt);
3116 
3117  set_vertex_and_cell_info(*this,
3118  vertex_touch_count,
3119  vertex_to_cell,
3120  coarse_cell_to_p4est_tree_permutation,
3121  set_vertex_info,
3122  connectivity);
3123 
3124  // next to tree-to-edge
3125  // data. note that in p4est lines
3126  // are ordered as follows
3127  // *---3---* *---3---*
3128  // /| | / /|
3129  // 6 | 11 6 7 11
3130  // / 10 | / / |
3131  // * | | *---2---* |
3132  // | *---1---* | | *
3133  // | / / | 9 /
3134  // 8 4 5 8 | 5
3135  // |/ / | |/
3136  // *---0---* *---0---*
3137  // whereas in deal.II they are like this:
3138  // *---7---* *---7---*
3139  // /| | / /|
3140  // 4 | 11 4 5 11
3141  // / 10 | / / |
3142  // * | | *---6---* |
3143  // | *---3---* | | *
3144  // | / / | 9 /
3145  // 8 0 1 8 | 1
3146  // |/ / | |/
3147  // *---2---* *---2---*
3148 
3149  const unsigned int deal_to_p4est_line_index[12] = {
3150  4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
3151 
3153  this->begin_active();
3154  cell != this->end();
3155  ++cell)
3156  {
3157  const unsigned int index =
3158  coarse_cell_to_p4est_tree_permutation[cell->index()];
3159  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++e)
3160  connectivity->tree_to_edge[index * GeometryInfo<3>::lines_per_cell +
3161  deal_to_p4est_line_index[e]] =
3162  cell->line(e)->index();
3163  }
3164 
3165  // now also set edge-to-tree
3166  // information
3167  connectivity->ett_offset[0] = 0;
3168  std::partial_sum(edge_touch_count.begin(),
3169  edge_touch_count.end(),
3170  &connectivity->ett_offset[1]);
3171 
3172  Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
3173  ExcInternalError());
3174 
3175  for (unsigned int v = 0; v < this->n_active_lines(); ++v)
3176  {
3177  Assert(edge_to_cell[v].size() == edge_touch_count[v],
3178  ExcInternalError());
3179 
3180  std::list<
3181  std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
3182  unsigned int>>::const_iterator p =
3183  edge_to_cell[v].begin();
3184  for (unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
3185  {
3186  connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
3187  coarse_cell_to_p4est_tree_permutation[p->first->index()];
3188  connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
3189  deal_to_p4est_line_index[p->second];
3190  }
3191  }
3192 
3193  Assert(p8est_connectivity_is_valid(connectivity) == 1,
3194  ExcInternalError());
3195 
3196  // now create a forest out of the connectivity data structure
3198  this->mpi_communicator,
3199  connectivity,
3200  /* minimum initial number of quadrants per tree */ 0,
3201  /* minimum level of upfront refinement */ 0,
3202  /* use uniform upfront refinement */ 1,
3203  /* user_data_size = */ 0,
3204  /* user_data_constructor = */ nullptr,
3205  /* user_pointer */ this);
3206  }
3207 # endif
3208 
3209 
3210 
3211  namespace
3212  {
3213  // ensures the 2:1 mesh balance for periodic boundary conditions in the
3214  // artificial cell layer (the active cells are taken care of by p4est)
3215  template <int dim, int spacedim>
3216  bool
3217  enforce_mesh_balance_over_periodic_boundaries(
3219  {
3220  if (tria.get_periodic_face_map().size() == 0)
3221  return false;
3222 
3223  std::vector<bool> flags_before[2];
3224  tria.save_coarsen_flags(flags_before[0]);
3225  tria.save_refine_flags(flags_before[1]);
3226 
3227  std::vector<unsigned int> topological_vertex_numbering(
3228  tria.n_vertices());
3229  for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
3230  topological_vertex_numbering[i] = i;
3231  // combine vertices that have different locations (and thus, different
3232  // vertex_index) but represent the same topological entity over periodic
3233  // boundaries. The vector topological_vertex_numbering contains a linear
3234  // map from 0 to n_vertices at input and at output relates periodic
3235  // vertices with only one vertex index. The output is used to always
3236  // identify the same vertex according to the periodicity, e.g. when
3237  // finding the maximum cell level around a vertex.
3238  //
3239  // Example: On a 3D cell with vertices numbered from 0 to 7 and periodic
3240  // boundary conditions in x direction, the vector
3241  // topological_vertex_numbering will contain the numbers
3242  // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5},
3243  // {6,7} belong together, respectively). If periodicity is set in x and
3244  // z direction, the output is {0,0,2,2,0,0,2,2}, and if periodicity is
3245  // in all directions, the output is simply {0,0,0,0,0,0,0,0}.
3246  using cell_iterator =
3248  typename std::map<std::pair<cell_iterator, unsigned int>,
3249  std::pair<std::pair<cell_iterator, unsigned int>,
3250  std::bitset<3>>>::const_iterator it;
3251  for (it = tria.get_periodic_face_map().begin();
3252  it != tria.get_periodic_face_map().end();
3253  ++it)
3254  {
3255  const cell_iterator &cell_1 = it->first.first;
3256  const unsigned int face_no_1 = it->first.second;
3257  const cell_iterator &cell_2 = it->second.first.first;
3258  const unsigned int face_no_2 = it->second.first.second;
3259  const std::bitset<3> face_orientation = it->second.second;
3260 
3261  if (cell_1->level() == cell_2->level())
3262  {
3263  for (unsigned int v = 0;
3264  v < GeometryInfo<dim - 1>::vertices_per_cell;
3265  ++v)
3266  {
3267  // take possible non-standard orientation of face on cell[0]
3268  // into account
3269  const unsigned int vface0 =
3271  v,
3272  face_orientation[0],
3273  face_orientation[1],
3274  face_orientation[2]);
3275  const unsigned int vi0 =
3276  topological_vertex_numbering[cell_1->face(face_no_1)
3277  ->vertex_index(vface0)];
3278  const unsigned int vi1 =
3279  topological_vertex_numbering[cell_2->face(face_no_2)
3280  ->vertex_index(v)];
3281  const unsigned int min_index = std::min(vi0, vi1);
3282  topological_vertex_numbering[cell_1->face(face_no_1)
3283  ->vertex_index(vface0)] =
3284  topological_vertex_numbering[cell_2->face(face_no_2)
3285  ->vertex_index(v)] =
3286  min_index;
3287  }
3288  }
3289  }
3290 
3291  // There must not be any chains!
3292  for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
3293  {
3294  const unsigned int j = topological_vertex_numbering[i];
3295  if (j != i)
3296  Assert(topological_vertex_numbering[j] == j, ExcInternalError());
3297  }
3298 
3299 
3300  // this code is replicated from grid/tria.cc but using an indirection
3301  // for periodic boundary conditions
3302  bool continue_iterating = true;
3303  std::vector<int> vertex_level(tria.n_vertices());
3304  while (continue_iterating)
3305  {
3306  // store highest level one of the cells adjacent to a vertex
3307  // belongs to
3308  std::fill(vertex_level.begin(), vertex_level.end(), 0);
3310  cell = tria.begin_active(),
3311  endc = tria.end();
3312  for (; cell != endc; ++cell)
3313  {
3314  if (cell->refine_flag_set())
3315  for (const unsigned int vertex :
3317  vertex_level[topological_vertex_numbering
3318  [cell->vertex_index(vertex)]] =
3319  std::max(vertex_level[topological_vertex_numbering
3320  [cell->vertex_index(vertex)]],
3321  cell->level() + 1);
3322  else if (!cell->coarsen_flag_set())
3323  for (const unsigned int vertex :
3325  vertex_level[topological_vertex_numbering
3326  [cell->vertex_index(vertex)]] =
3327  std::max(vertex_level[topological_vertex_numbering
3328  [cell->vertex_index(vertex)]],
3329  cell->level());
3330  else
3331  {
3332  // if coarsen flag is set then tentatively assume
3333  // that the cell will be coarsened. this isn't
3334  // always true (the coarsen flag could be removed
3335  // again) and so we may make an error here. we try
3336  // to correct this by iterating over the entire
3337  // process until we are converged
3338  Assert(cell->coarsen_flag_set(), ExcInternalError());
3339  for (const unsigned int vertex :
3341  vertex_level[topological_vertex_numbering
3342  [cell->vertex_index(vertex)]] =
3343  std::max(vertex_level[topological_vertex_numbering
3344  [cell->vertex_index(vertex)]],
3345  cell->level() - 1);
3346  }
3347  }
3348 
3349  continue_iterating = false;
3350 
3351  // loop over all cells in reverse order. do so because we
3352  // can then update the vertex levels on the adjacent
3353  // vertices and maybe already flag additional cells in this
3354  // loop
3355  //
3356  // note that not only may we have to add additional
3357  // refinement flags, but we will also have to remove
3358  // coarsening flags on cells adjacent to vertices that will
3359  // see refinement
3360  for (cell = tria.last_active(); cell != endc; --cell)
3361  if (cell->refine_flag_set() == false)
3362  {
3363  for (const unsigned int vertex :
3365  if (vertex_level[topological_vertex_numbering
3366  [cell->vertex_index(vertex)]] >=
3367  cell->level() + 1)
3368  {
3369  // remove coarsen flag...
3370  cell->clear_coarsen_flag();
3371 
3372  // ...and if necessary also refine the current
3373  // cell, at the same time updating the level
3374  // information about vertices
3375  if (vertex_level[topological_vertex_numbering
3376  [cell->vertex_index(vertex)]] >
3377  cell->level() + 1)
3378  {
3379  cell->set_refine_flag();
3380  continue_iterating = true;
3381 
3382  for (const unsigned int v :
3384  vertex_level[topological_vertex_numbering
3385  [cell->vertex_index(v)]] =
3386  std::max(
3387  vertex_level[topological_vertex_numbering
3388  [cell->vertex_index(v)]],
3389  cell->level() + 1);
3390  }
3391 
3392  // continue and see whether we may, for example,
3393  // go into the inner 'if' above based on a
3394  // different vertex
3395  }
3396  }
3397 
3398  // clear coarsen flag if not all children were marked
3399  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3400  tria.begin();
3401  cell != tria.end();
3402  ++cell)
3403  {
3404  // nothing to do if we are already on the finest level
3405  if (cell->is_active())
3406  continue;
3407 
3408  const unsigned int n_children = cell->n_children();
3409  unsigned int flagged_children = 0;
3410  for (unsigned int child = 0; child < n_children; ++child)
3411  if (cell->child(child)->is_active() &&
3412  cell->child(child)->coarsen_flag_set())
3413  ++flagged_children;
3414 
3415  // if not all children were flagged for coarsening, remove
3416  // coarsen flags
3417  if (flagged_children < n_children)
3418  for (unsigned int child = 0; child < n_children; ++child)
3419  if (cell->child(child)->is_active())
3420  cell->child(child)->clear_coarsen_flag();
3421  }
3422  }
3423  std::vector<bool> flags_after[2];
3424  tria.save_coarsen_flags(flags_after[0]);
3425  tria.save_refine_flags(flags_after[1]);
3426  return ((flags_before[0] != flags_after[0]) ||
3427  (flags_before[1] != flags_after[1]));
3428  }
3429  } // namespace
3430 
3431 
3432 
3433  template <int dim, int spacedim>
3434  bool
3436  {
3437  std::vector<bool> flags_before[2];
3438  this->save_coarsen_flags(flags_before[0]);
3439  this->save_refine_flags(flags_before[1]);
3440 
3441  bool mesh_changed = false;
3442  unsigned int loop_counter = 0;
3443  do
3444  {
3447  this->update_periodic_face_map();
3448  // enforce 2:1 mesh balance over periodic boundaries
3449  mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this);
3450 
3451  // We can't be sure that we won't run into a situation where we can
3452  // not reconcile mesh smoothing and balancing of periodic faces. As we
3453  // don't know what else to do, at least abort with an error message.
3454  ++loop_counter;
3455  AssertThrow(
3456  loop_counter < 32,
3457  ExcMessage(
3458  "Infinite loop in "
3459  "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
3460  "for periodic boundaries detected. Aborting."));
3461  }
3462  while (mesh_changed);
3463 
3464  // check if any of the refinement flags were changed during this
3465  // function and return that value
3466  std::vector<bool> flags_after[2];
3467  this->save_coarsen_flags(flags_after[0]);
3468  this->save_refine_flags(flags_after[1]);
3469  return ((flags_before[0] != flags_after[0]) ||
3470  (flags_before[1] != flags_after[1]));
3471  }
3472 
3473 
3474 
3475  template <int dim, int spacedim>
3476  void
3478  {
3479  // disable mesh smoothing for recreating the deal.II triangulation,
3480  // otherwise we might not be able to reproduce the p4est mesh
3481  // exactly. We restore the original smoothing at the end of this
3482  // function. Note that the smoothing flag is used in the normal
3483  // refinement process.
3484  typename Triangulation<dim, spacedim>::MeshSmoothing save_smooth =
3485  this->smooth_grid;
3486 
3487  // We will refine manually to match the p4est further down, which
3488  // obeys a level difference of 2 at each vertex (see the balance call
3489  // to p4est). We can disable this here so we store fewer artificial
3490  // cells (in some cases).
3491  // For geometric multigrid it turns out that
3492  // we will miss level cells at shared vertices if we ignore this.
3493  // See tests/mpi/mg_06. In particular, the flag is still necessary
3494  // even though we force it for the original smooth_grid in the
3495  // constructor.
3496  if (settings & construct_multigrid_hierarchy)
3497  this->smooth_grid =
3498  ::Triangulation<dim,
3500  else
3502 
3503  bool mesh_changed = false;
3504 
3505  // remove all deal.II refinements. Note that we could skip this and
3506  // start from our current state, because the algorithm later coarsens as
3507  // necessary. This has the advantage of being faster when large parts
3508  // of the local partition changes (likely) and gives a deterministic
3509  // ordering of the cells (useful for snapshot/resume).
3510  // TODO: is there a more efficient way to do this?
3511  if (settings & mesh_reconstruction_after_repartitioning)
3512  while (this->begin_active()->level() > 0)
3513  {
3514  for (const auto &cell : this->active_cell_iterators())
3515  {
3516  cell->set_coarsen_flag();
3517  }
3518 
3520  try
3521  {
3524  }
3525  catch (
3527  {
3528  // the underlying triangulation should not be checking for
3529  // distorted cells
3530  Assert(false, ExcInternalError());
3531  }
3532  }
3533 
3534 
3535  // query p4est for the ghost cells
3536  if (parallel_ghost != nullptr)
3537  {
3539  parallel_ghost);
3540  parallel_ghost = nullptr;
3541  }
3543  parallel_forest,
3544  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3545  P4EST_CONNECT_CORNER) :
3546  typename ::internal::p4est::types<dim>::balance_type(
3547  P8EST_CONNECT_CORNER)));
3548 
3549  Assert(parallel_ghost, ExcInternalError());
3550 
3551 
3552  // set all cells to artificial. we will later set it to the correct
3553  // subdomain in match_tree_recursively
3554  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3555  this->begin(0);
3556  cell != this->end(0);
3557  ++cell)
3558  cell->recursively_set_subdomain_id(numbers::artificial_subdomain_id);
3559 
3560  do
3561  {
3562  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3563  this->begin(0);
3564  cell != this->end(0);
3565  ++cell)
3566  {
3567  // if this processor stores no part of the forest that comes out
3568  // of this coarse grid cell, then we need to delete all children
3569  // of this cell (the coarse grid cell remains)
3570  if (tree_exists_locally<dim, spacedim>(
3571  parallel_forest,
3572  coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
3573  false)
3574  {
3575  delete_all_children<dim, spacedim>(cell);
3576  if (cell->is_active())
3577  cell->set_subdomain_id(numbers::artificial_subdomain_id);
3578  }
3579 
3580  else
3581  {
3582  // this processor stores at least a part of the tree that
3583  // comes out of this cell.
3584 
3585  typename ::internal::p4est::types<dim>::quadrant
3586  p4est_coarse_cell;
3587  typename ::internal::p4est::types<dim>::tree *tree =
3588  init_tree(cell->index());
3589 
3590  ::internal::p4est::init_coarse_quadrant<dim>(
3591  p4est_coarse_cell);
3592 
3593  match_tree_recursively<dim, spacedim>(*tree,
3594  cell,
3595  p4est_coarse_cell,
3596  *parallel_forest,
3597  this->my_subdomain);
3598  }
3599  }
3600 
3601  // check mesh for ghost cells, refine as necessary. iterate over
3602  // every ghostquadrant, find corresponding deal coarsecell and
3603  // recurse.
3604  typename ::internal::p4est::types<dim>::quadrant *quadr;
3605  types::subdomain_id ghost_owner = 0;
3606  typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
3607 
3608  for (unsigned int g_idx = 0;
3609  g_idx < parallel_ghost->ghosts.elem_count;
3610  ++g_idx)
3611  {
3612  while (g_idx >= static_cast<unsigned int>(
3613  parallel_ghost->proc_offsets[ghost_owner + 1]))
3614  ++ghost_owner;
3615  while (g_idx >= static_cast<unsigned int>(
3616  parallel_ghost->tree_offsets[ghost_tree + 1]))
3617  ++ghost_tree;
3618 
3619  quadr = static_cast<
3620  typename ::internal::p4est::types<dim>::quadrant *>(
3621  sc_array_index(&parallel_ghost->ghosts, g_idx));
3622 
3623  unsigned int coarse_cell_index =
3624  p4est_tree_to_coarse_cell_permutation[ghost_tree];
3625 
3626  match_quadrant<dim, spacedim>(this,
3627  coarse_cell_index,
3628  *quadr,
3629  ghost_owner);
3630  }
3631 
3632  // fix all the flags to make sure we have a consistent mesh
3634 
3635  // see if any flags are still set
3636  mesh_changed =
3637  std::any_of(this->begin_active(),
3638  active_cell_iterator{this->end()},
3639  [](const CellAccessor<dim, spacedim> &cell) {
3640  return cell.refine_flag_set() ||
3641  cell.coarsen_flag_set();
3642  });
3643 
3644  // actually do the refinement to change the local mesh by
3645  // calling the base class refinement function directly
3646  try
3647  {
3650  }
3651  catch (
3653  {
3654  // the underlying triangulation should not be checking for
3655  // distorted cells
3656  Assert(false, ExcInternalError());
3657  }
3658  }
3659  while (mesh_changed);
3660 
3661 # ifdef DEBUG
3662  // check if correct number of ghosts is created
3663  unsigned int num_ghosts = 0;
3664 
3665  for (const auto &cell : this->active_cell_iterators())
3666  {
3667  if (cell->subdomain_id() != this->my_subdomain &&
3668  cell->subdomain_id() != numbers::artificial_subdomain_id)
3669  ++num_ghosts;
3670  }
3671 
3672  Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
3673  ExcInternalError());
3674 # endif
3675 
3676 
3677 
3678  // fill level_subdomain_ids for geometric multigrid
3679  // the level ownership of a cell is defined as the owner if the cell is
3680  // active or as the owner of child(0) we need this information for all our
3681  // ancestors and the same-level neighbors of our own cells (=level ghosts)
3682  if (settings & construct_multigrid_hierarchy)
3683  {
3684  // step 1: We set our own ids all the way down and all the others to
3685  // -1. Note that we do not fill other cells we could figure out the
3686  // same way, because we might accidentally set an id for a cell that
3687  // is not a ghost cell.
3688  for (unsigned int lvl = this->n_levels(); lvl > 0;)
3689  {
3690  --lvl;
3692  endc = this->end(lvl);
3693  for (cell = this->begin(lvl); cell != endc; ++cell)
3694  {
3695  if ((cell->is_active() &&
3696  cell->subdomain_id() ==
3697  this->locally_owned_subdomain()) ||
3698  (cell->has_children() &&
3699  cell->child(0)->level_subdomain_id() ==
3700  this->locally_owned_subdomain()))
3701  cell->set_level_subdomain_id(
3702  this->locally_owned_subdomain());
3703  else
3704  {
3705  // not our cell
3706  cell->set_level_subdomain_id(
3708  }
3709  }
3710  }
3711 
3712  // step 2: make sure all the neighbors to our level_cells exist. Need
3713  // to look up in p4est...
3714  std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3715  for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3716  marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3717 
3718  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3719  this->begin(0);
3720  cell != this->end(0);
3721  ++cell)
3722  {
3723  typename ::internal::p4est::types<dim>::quadrant
3724  p4est_coarse_cell;
3725  const unsigned int tree_index =
3726  coarse_cell_to_p4est_tree_permutation[cell->index()];
3727  typename ::internal::p4est::types<dim>::tree *tree =
3728  init_tree(cell->index());
3729 
3730  ::internal::p4est::init_coarse_quadrant<dim>(
3731  p4est_coarse_cell);
3732 
3733  determine_level_subdomain_id_recursively<dim, spacedim>(
3734  *tree,
3735  tree_index,
3736  cell,
3737  p4est_coarse_cell,
3738  *parallel_forest,
3739  this->my_subdomain,
3740  marked_vertices);
3741  }
3742 
3743  // step 3: make sure we have the parent of our level cells
3744  for (unsigned int lvl = this->n_levels(); lvl > 0;)
3745  {
3746  --lvl;
3748  endc = this->end(lvl);
3749  for (cell = this->begin(lvl); cell != endc; ++cell)
3750  {
3751  if (cell->has_children())
3752  for (unsigned int c = 0;
3753  c < GeometryInfo<dim>::max_children_per_cell;
3754  ++c)
3755  {
3756  if (cell->child(c)->level_subdomain_id() ==
3757  this->locally_owned_subdomain())
3758  {
3759  // at least one of the children belongs to us, so
3760  // make sure we set the level subdomain id
3761  const types::subdomain_id mark =
3762  cell->child(0)->level_subdomain_id();
3764  ExcInternalError()); // we should know the
3765  // child(0)
3766  cell->set_level_subdomain_id(mark);
3767  break;
3768  }
3769  }
3770  }
3771  }
3772  }
3773 
3774 
3775 
3776  // check that our local copy has exactly as many cells as the p4est
3777  // original (at least if we are on only one processor); for parallel
3778  // computations, we want to check that we have at least as many as p4est
3779  // stores locally (in the future we should check that we have exactly as
3780  // many non-artificial cells as parallel_forest->local_num_quadrants)
3781  {
3782  const unsigned int total_local_cells = this->n_active_cells();
3783  (void)total_local_cells;
3784 
3785  if (Utilities::MPI::n_mpi_processes(this->mpi_communicator) == 1)
3786  {
3787  Assert(static_cast<unsigned int>(
3788  parallel_forest->local_num_quadrants) == total_local_cells,
3789  ExcInternalError());
3790  }
3791  else
3792  {
3793  Assert(static_cast<unsigned int>(
3794  parallel_forest->local_num_quadrants) <= total_local_cells,
3795  ExcInternalError());
3796  }
3797 
3798  // count the number of owned, active cells and compare with p4est.
3799  unsigned int n_owned = 0;
3800  for (const auto &cell : this->active_cell_iterators())
3801  {
3802  if (cell->subdomain_id() == this->my_subdomain)
3803  ++n_owned;
3804  }
3805 
3806  Assert(static_cast<unsigned int>(
3807  parallel_forest->local_num_quadrants) == n_owned,
3808  ExcInternalError());
3809  }
3810 
3811  this->smooth_grid = save_smooth;
3812 
3813  // finally, after syncing the parallel_forest with the triangulation,
3814  // also update the quadrant_cell_relations, which will be used for
3815  // repartitioning, further refinement/coarsening, and unpacking
3816  // of stored or transferred data.
3817  update_quadrant_cell_relations();
3818  }
3819 
3820 
3821 
3822  template <int dim, int spacedim>
3823  void
3825  {
3826  // do not allow anisotropic refinement
3827 # ifdef DEBUG
3828  for (const auto &cell : this->active_cell_iterators())
3829  if (cell->is_locally_owned() && cell->refine_flag_set())
3830  Assert(cell->refine_flag_set() ==
3832  ExcMessage(
3833  "This class does not support anisotropic refinement"));
3834 # endif
3835 
3836 
3837  // safety check: p4est has an upper limit on the level of a cell
3838  if (this->n_levels() ==
3840  {
3842  cell = this->begin_active(
3844  cell !=
3846  1);
3847  ++cell)
3848  {
3849  AssertThrow(
3850  !(cell->refine_flag_set()),
3851  ExcMessage(
3852  "Fatal Error: maximum refinement level of p4est reached."));
3853  }
3854  }
3855 
3857 
3858  // signal that refinement is going to happen
3860 
3861  // now do the work we're supposed to do when we are in charge
3862  // make sure all flags are cleared on cells we don't own, since nothing
3863  // good can come of that if they are still around
3864  for (const auto &cell : this->active_cell_iterators())
3865  if (cell->is_ghost() || cell->is_artificial())
3866  {
3867  cell->clear_refine_flag();
3868  cell->clear_coarsen_flag();
3869  }
3870 
3871 
3872  // count how many cells will be refined and coarsened, and allocate that
3873  // much memory
3874  RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3875  *this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3876 
3877  // copy refine and coarsen flags into p4est and execute the refinement
3878  // and coarsening. this uses the refine_and_coarsen_list just built,
3879  // which is communicated to the callback functions through
3880  // p4est's user_pointer object
3881  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3882  parallel_forest->user_pointer = &refine_and_coarsen_list;
3883 
3884  if (parallel_ghost != nullptr)
3885  {
3887  parallel_ghost);
3888  parallel_ghost = nullptr;
3889  }
3891  parallel_forest,
3892  /* refine_recursive */ false,
3893  &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3894  /*init_callback=*/nullptr);
3896  parallel_forest,
3897  /* coarsen_recursive */ false,
3898  &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3899  /*init_callback=*/nullptr);
3900 
3901  // make sure all cells in the lists have been consumed
3902  Assert(refine_and_coarsen_list.pointers_are_at_end(), ExcInternalError());
3903 
3904  // reset the pointer
3905  parallel_forest->user_pointer = this;
3906 
3907  // enforce 2:1 hanging node condition
3909  parallel_forest,
3910  /* face and corner balance */
3911  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3912  P4EST_CONNECT_FULL) :
3913  typename ::internal::p4est::types<dim>::balance_type(
3914  P8EST_CONNECT_FULL)),
3915  /*init_callback=*/nullptr);
3916 
3917  // since refinement and/or coarsening on the parallel forest
3918  // has happened, we need to update the quadrant cell relations
3919  update_quadrant_cell_relations();
3920 
3921  // before repartitioning the mesh, store the current distribution
3922  // of the p4est quadrants and let others attach mesh related info
3923  // (such as SolutionTransfer data)
3924  std::vector<typename ::internal::p4est::types<dim>::gloidx>
3925  previous_global_first_quadrant;
3926 
3927  // pack data only if anything has been attached
3928  if (cell_attached_data.n_attached_data_sets > 0)
3929  {
3930  data_transfer.pack_data(local_quadrant_cell_relations,
3931  cell_attached_data.pack_callbacks_fixed,
3932  cell_attached_data.pack_callbacks_variable);
3933 
3934  // before repartitioning the p4est object, save a copy of the
3935  // positions of the global first quadrants for data transfer later
3936  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3937  std::memcpy(previous_global_first_quadrant.data(),
3938  parallel_forest->global_first_quadrant,
3939  sizeof(
3940  typename ::internal::p4est::types<dim>::gloidx) *
3941  (parallel_forest->mpisize + 1));
3942  }
3943 
3944  if (!(settings & no_automatic_repartitioning))
3945  {
3946  // partition the new mesh between all processors. If cell weights have
3947  // not been given balance the number of cells.
3948  if (this->signals.cell_weight.num_slots() == 0)
3950  parallel_forest,
3951  /* prepare coarsening */ 1,
3952  /* weight_callback */ nullptr);
3953  else
3954  {
3955  // get cell weights for a weighted repartitioning.
3956  const std::vector<unsigned int> cell_weights = get_cell_weights();
3957 
3958  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3959 
3960  // attach (temporarily) a pointer to the cell weights through
3961  // p4est's user_pointer object
3962  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3963  parallel_forest->user_pointer = &partition_weights;
3964 
3966  parallel_forest,
3967  /* prepare coarsening */ 1,
3968  /* weight_callback */
3969  &PartitionWeights<dim, spacedim>::cell_weight);
3970 
3971  // release data
3973  parallel_forest, 0, nullptr, nullptr);
3974  // reset the user pointer to its previous state
3975  parallel_forest->user_pointer = this;
3976  }
3977  }
3978 
3979  // finally copy back from local part of tree to deal.II
3980  // triangulation. before doing so, make sure there are no refine or
3981  // coarsen flags pending
3982  for (const auto &cell : this->active_cell_iterators())
3983  {
3984  cell->clear_refine_flag();
3985  cell->clear_coarsen_flag();
3986  }
3987 
3988  try
3989  {
3990  copy_local_forest_to_triangulation();
3991  }
3992  catch (const typename Triangulation<dim>::DistortedCellList &)
3993  {
3994  // the underlying triangulation should not be checking for distorted
3995  // cells
3996  Assert(false, ExcInternalError());
3997  }
3998 
3999  // transfer data
4000  // only if anything has been attached
4001  if (cell_attached_data.n_attached_data_sets > 0)
4002  {
4003  // execute transfer after triangulation got updated
4004  data_transfer.execute_transfer(parallel_forest,
4005  previous_global_first_quadrant.data());
4006 
4007  // also update the CellStatus information on the new mesh
4008  data_transfer.unpack_cell_status(local_quadrant_cell_relations);
4009  }
4010 
4011 # ifdef DEBUG
4012  // Check that we know the level subdomain ids of all our neighbors. This
4013  // also involves coarser cells that share a vertex if they are active.
4014  //
4015  // Example (M= my, O=other):
4016  // *------*
4017  // | |
4018  // | O |
4019  // | |
4020  // *---*---*------*
4021  // | M | M |
4022  // *---*---*
4023  // | | M |
4024  // *---*---*
4025  // ^- the parent can be owned by somebody else, so O is not a neighbor
4026  // one level coarser
4027  if (settings & construct_multigrid_hierarchy)
4028  {
4029  for (unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
4030  {
4031  std::vector<bool> active_verts =
4032  this->mark_locally_active_vertices_on_level(lvl);
4033 
4034  const unsigned int maybe_coarser_lvl =
4035  (lvl > 0) ? (lvl - 1) : lvl;
4037  cell = this->begin(maybe_coarser_lvl),
4038  endc = this->end(lvl);
4039  for (; cell != endc; ++cell)
4040  if (cell->level() == static_cast<int>(lvl) || cell->is_active())
4041  {
4042  const bool is_level_artificial =
4043  (cell->level_subdomain_id() ==
4045  bool need_to_know = false;
4046  for (const unsigned int vertex :
4048  if (active_verts[cell->vertex_index(vertex)])
4049  {
4050  need_to_know = true;
4051  break;
4052  }
4053 
4054  Assert(
4055  !need_to_know || !is_level_artificial,
4056  ExcMessage(
4057  "Internal error: the owner of cell" +
4058  cell->id().to_string() +
4059  " is unknown even though it is needed for geometric multigrid."));
4060  }
4061  }
4062  }
4063 # endif
4064 
4065  this->update_periodic_face_map();
4066  this->update_number_cache();
4067 
4068  // signal that refinement is finished
4070  }
4071 
4072 
4073 
4074  template <int dim, int spacedim>
4075  void
4077  {
4078 # ifdef DEBUG
4079  for (const auto &cell : this->active_cell_iterators())
4080  if (cell->is_locally_owned())
4081  Assert(
4082  !cell->refine_flag_set() && !cell->coarsen_flag_set(),
4083  ExcMessage(
4084  "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
4085 # endif
4086 
4087  // signal that repartitioning is going to happen
4089 
4090  // before repartitioning the mesh let others attach mesh related info
4091  // (such as SolutionTransfer data) to the p4est
4092  std::vector<typename ::internal::p4est::types<dim>::gloidx>
4093  previous_global_first_quadrant;
4094 
4095  // pack data only if anything has been attached
4096  if (cell_attached_data.n_attached_data_sets > 0)
4097  {
4098  data_transfer.pack_data(local_quadrant_cell_relations,
4099  cell_attached_data.pack_callbacks_fixed,
4100  cell_attached_data.pack_callbacks_variable);
4101 
4102  // before repartitioning the p4est object, save a copy of the
4103  // positions of quadrant for data transfer later
4104  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
4105  std::memcpy(previous_global_first_quadrant.data(),
4106  parallel_forest->global_first_quadrant,
4107  sizeof(
4108  typename ::internal::p4est::types<dim>::gloidx) *
4109  (parallel_forest->mpisize + 1));
4110  }
4111 
4112  if (this->signals.cell_weight.num_slots() == 0)
4113  {
4114  // no cell weights given -- call p4est's 'partition' without a
4115  // callback for cell weights
4117  parallel_forest,
4118  /* prepare coarsening */ 1,
4119  /* weight_callback */ nullptr);
4120  }
4121  else
4122  {
4123  // get cell weights for a weighted repartitioning.
4124  const std::vector<unsigned int> cell_weights = get_cell_weights();
4125 
4126  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
4127 
4128  // attach (temporarily) a pointer to the cell weights through p4est's
4129  // user_pointer object
4130  Assert(parallel_forest->user_pointer == this, ExcInternalError());
4131  parallel_forest->user_pointer = &partition_weights;
4132 
4134  parallel_forest,
4135  /* prepare coarsening */ 1,
4136  /* weight_callback */
4137  &PartitionWeights<dim, spacedim>::cell_weight);
4138 
4139  // reset the user pointer to its previous state
4140  parallel_forest->user_pointer = this;
4141  }
4142 
4143  try
4144  {
4145  copy_local_forest_to_triangulation();
4146  }
4147  catch (const typename Triangulation<dim>::DistortedCellList &)
4148  {
4149  // the underlying triangulation should not be checking for distorted
4150  // cells
4151  Assert(false, ExcInternalError());
4152  }
4153 
4154  // transfer data
4155  // only if anything has been attached
4156  if (cell_attached_data.n_attached_data_sets > 0)
4157  {
4158  // execute transfer after triangulation got updated
4159  data_transfer.execute_transfer(parallel_forest,
4160  previous_global_first_quadrant.data());
4161  }
4162 
4163  this->update_periodic_face_map();
4164 
4165  // update how many cells, edges, etc, we store locally
4166  this->update_number_cache();
4167 
4168  // signal that repartitioning is finished
4170  }
4171 
4172 
4173 
4174  template <int dim, int spacedim>
4175  void
4177  const std::vector<bool> &vertex_locally_moved)
4178  {
4179  Assert(vertex_locally_moved.size() == this->n_vertices(),
4180  ExcDimensionMismatch(vertex_locally_moved.size(),
4181  this->n_vertices()));
4182 # ifdef DEBUG
4183  {
4184  const std::vector<bool> locally_owned_vertices =
4186  for (unsigned int i = 0; i < locally_owned_vertices.size(); ++i)
4187  Assert((vertex_locally_moved[i] == false) ||
4188  (locally_owned_vertices[i] == true),
4189  ExcMessage("The vertex_locally_moved argument must not "
4190  "contain vertices that are not locally owned"));
4191  }
4192 # endif
4193 
4194  // First find out which process should receive which vertices.
4195  // These are specifically the ones that are located on cells at the
4196  // boundary of the subdomain this process owns and the receiving
4197  // process taking periodic faces into account.
4198  // Here, it is sufficient to collect all vertices that are located
4199  // at that boundary.
4200  const std::map<unsigned int, std::set<::types::subdomain_id>>
4203 
4204  // now collect cells and their vertices
4205  // for the interested neighbors
4206  using cellmap_t =
4207  std::map<::types::subdomain_id,
4208  CommunicateLocallyMovedVertices::CellInfo<dim, spacedim>>;
4209  cellmap_t needs_to_get_cells;
4210 
4211  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
4212  this->begin(0);
4213  cell != this->end(0);
4214  ++cell)
4215  {
4216  typename ::internal::p4est::types<dim>::quadrant
4217  p4est_coarse_cell;
4218  ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4219 
4220  CommunicateLocallyMovedVertices::fill_vertices_recursively<dim,
4221  spacedim>(
4222  *this,
4223  this->get_coarse_cell_to_p4est_tree_permutation()[cell->index()],
4224  cell,
4225  p4est_coarse_cell,
4227  vertex_locally_moved,
4228  needs_to_get_cells);
4229  }
4230 
4231  // Send information
4232 
4233  // We need to protect this communication below using a mutex:
4236  mutex, this->get_communicator());
4237 
4238  const int mpi_tag = Utilities::MPI::internal::Tags::
4240 
4241  std::vector<std::vector<char>> sendbuffers(needs_to_get_cells.size());
4242  std::vector<std::vector<char>>::iterator buffer = sendbuffers.begin();
4243  std::vector<MPI_Request> requests(needs_to_get_cells.size());
4244  std::vector<unsigned int> destinations;
4245 
4246  unsigned int idx = 0;
4247 
4248  for (typename cellmap_t::iterator it = needs_to_get_cells.begin();
4249  it != needs_to_get_cells.end();
4250  ++it, ++buffer, ++idx)
4251  {
4252  const unsigned int num_cells = it->second.tree_index.size();
4253  (void)num_cells;
4254  destinations.push_back(it->first);
4255 
4256  Assert(num_cells == it->second.quadrants.size(), ExcInternalError());
4257  Assert(num_cells > 0, ExcInternalError());
4258 
4259  // pack all the data into
4260  // the buffer for this
4261  // recipient and send
4262  // it. keep data around
4263  // till we can make sure
4264  // that the packet has been
4265  // received
4266  it->second.pack_data(*buffer);
4267  const int ierr = MPI_Isend(buffer->data(),
4268  buffer->size(),
4269  MPI_BYTE,
4270  it->first,
4271  mpi_tag,
4272  this->get_communicator(),
4273  &requests[idx]);
4274  AssertThrowMPI(ierr);
4275  }
4276 
4277  Assert(destinations.size() == needs_to_get_cells.size(),
4278  ExcInternalError());
4279 
4280  // collect the neighbors
4281  // that are going to send stuff to us
4282  const unsigned int n_senders =
4284  this->get_communicator(), destinations);
4285 
4286  // receive ghostcelldata
4287  std::vector<char> receive;
4288  CommunicateLocallyMovedVertices::CellInfo<dim, spacedim> cellinfo;
4289  for (unsigned int i = 0; i < n_senders; ++i)
4290  {
4291  MPI_Status status;
4292  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4293  mpi_tag,
4294  this->get_communicator(),
4295  &status);
4296  AssertThrowMPI(ierr);
4297 
4298  int len;
4299  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4300  AssertThrowMPI(ierr);
4301  receive.resize(len);
4302 
4303  char *ptr = receive.data();
4304  ierr = MPI_Recv(ptr,
4305  len,
4306  MPI_BYTE,
4307  status.MPI_SOURCE,
4308  status.MPI_TAG,
4309  this->get_communicator(),
4310  &status);
4311  AssertThrowMPI(ierr);
4312 
4313  cellinfo.unpack_data(receive);
4314  const unsigned int cells = cellinfo.tree_index.size();
4315  for (unsigned int c = 0; c < cells; ++c)
4316  {
4317  typename ::parallel::distributed::
4318  Triangulation<dim, spacedim>::cell_iterator cell(
4319  this,
4320  0,
4321  this->get_p4est_tree_to_coarse_cell_permutation()
4322  [cellinfo.tree_index[c]]);
4323 
4324  typename ::internal::p4est::types<dim>::quadrant
4325  p4est_coarse_cell;
4326  ::internal::p4est::init_coarse_quadrant<dim>(
4327  p4est_coarse_cell);
4328 
4329  CommunicateLocallyMovedVertices::set_vertices_recursively<
4330  dim,
4331  spacedim>(*this,
4332  p4est_coarse_cell,
4333  cell,
4334  cellinfo.quadrants[c],
4335  cellinfo.first_vertices[c],
4336  cellinfo.first_vertex_indices[c]);
4337  }
4338  }
4339 
4340  // complete all sends, so that we can
4341  // safely destroy the buffers.
4342  if (requests.size() > 0)
4343  {
4344  const int ierr =
4345  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4346  AssertThrowMPI(ierr);
4347  }
4348 
4349  // check all msgs got sent and received
4350  Assert(Utilities::MPI::sum(needs_to_get_cells.size(),
4351  this->get_communicator()) ==
4352  Utilities::MPI::sum(n_senders, this->get_communicator()),
4353  ExcInternalError());
4354  }
4355 
4356 
4357 
4358  template <int dim, int spacedim>
4359  unsigned int
4361  const std::function<std::vector<char>(const cell_iterator &,
4362  const CellStatus)> &pack_callback,
4363  const bool returns_variable_size_data)
4364  {
4365  unsigned int handle = numbers::invalid_unsigned_int;
4366 
4367  // Add new callback function to the corresponding register.
4368  // Encode handles according to returns_variable_size_data.
4369  if (returns_variable_size_data)
4370  {
4371  handle = 2 * cell_attached_data.pack_callbacks_variable.size();
4372  cell_attached_data.pack_callbacks_variable.push_back(pack_callback);
4373  }
4374  else
4375  {
4376  handle = 2 * cell_attached_data.pack_callbacks_fixed.size() + 1;
4377  cell_attached_data.pack_callbacks_fixed.push_back(pack_callback);
4378  }
4379 
4380  // Increase overall counter.
4381  ++cell_attached_data.n_attached_data_sets;
4382 
4383  return handle;
4384  }
4385 
4386 
4387 
4388  template <int dim, int spacedim>
4389  void
4391  const unsigned int handle,
4392  const std::function<
4393  void(const cell_iterator &,
4394  const CellStatus,
4395  const boost::iterator_range<std::vector<char>::const_iterator> &)>
4396  &unpack_callback)
4397  {
4398  Assert(cell_attached_data.n_attached_data_sets > 0,
4399  ExcMessage("The notify_ready_to_unpack() has already been called "
4400  "once for each registered callback."));
4401 
4402  // check if local_quadrant_cell_relations have been previously gathered
4403  // correctly
4404  Assert(local_quadrant_cell_relations.size() ==
4405  static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4406  ExcInternalError());
4407 
4408 # ifdef DEBUG
4409  // check validity of handle and deregister pack_callback function.
4410  // first reset with invalid entries to preserve ambiguity of
4411  // handles, then free memory when all were unpacked (see below).
4412  const unsigned int callback_index = handle / 2;
4413  if (handle % 2 == 0)
4414  {
4415  Assert(callback_index <
4416  cell_attached_data.pack_callbacks_variable.size(),
4417  ExcMessage("Invalid handle."));
4418 
4419  Assert(cell_attached_data.pack_callbacks_variable[callback_index] !=
4420  nullptr,
4421  ExcInternalError());
4422  cell_attached_data.pack_callbacks_variable[callback_index] = nullptr;
4423  }
4424  else
4425  {
4426  Assert(callback_index <
4427  cell_attached_data.pack_callbacks_fixed.size(),
4428  ExcMessage("Invalid handle."));
4429 
4430  Assert(cell_attached_data.pack_callbacks_fixed[callback_index] !=
4431  nullptr,
4432  ExcInternalError());
4433  cell_attached_data.pack_callbacks_fixed[callback_index] = nullptr;
4434  }
4435 # endif
4436 
4437  // perform unpacking
4438  data_transfer.unpack_data(local_quadrant_cell_relations,
4439  handle,
4440  unpack_callback);
4441 
4442  // decrease counters
4443  --cell_attached_data.n_attached_data_sets;
4444  if (cell_attached_data.n_attached_deserialize > 0)
4445  --cell_attached_data.n_attached_deserialize;
4446 
4447  // important: only remove data if we are not in the deserialization
4448  // process. There, each SolutionTransfer registers and unpacks before
4449  // the next one does this, so n_attached_data_sets is only 1 here. This
4450  // would destroy the saved data before the second SolutionTransfer can
4451  // get it. This created a bug that is documented in
4452  // tests/mpi/p4est_save_03 with more than one SolutionTransfer.
4453  if (cell_attached_data.n_attached_data_sets == 0 &&
4454  cell_attached_data.n_attached_deserialize == 0)
4455  {
4456  // everybody got their data, time for cleanup!
4457  cell_attached_data.pack_callbacks_fixed.clear();
4458  cell_attached_data.pack_callbacks_variable.clear();
4459  data_transfer.clear();
4460 
4461  // reset all cell_status entries after coarsening/refinement
4462  for (auto &quad_cell_rel : local_quadrant_cell_relations)
4463  std::get<1>(quad_cell_rel) =
4465  }
4466  }
4467 
4468 
4469 
4470  template <int dim, int spacedim>
4471  const std::vector<types::global_dof_index> &
4473  const
4474  {
4475  return p4est_tree_to_coarse_cell_permutation;
4476  }
4477 
4478 
4479 
4480  template <int dim, int spacedim>
4481  const std::vector<types::global_dof_index> &
4483  const
4484  {
4485  return coarse_cell_to_p4est_tree_permutation;
4486  }
4487 
4488 
4489 
4490  template <int dim, int spacedim>
4491  std::vector<bool>
4493  const int level) const
4494  {
4495  Assert(dim > 1, ExcNotImplemented());
4496 
4497  std::vector<bool> marked_vertices(this->n_vertices(), false);
4498  cell_iterator cell = this->begin(level), endc = this->end(level);
4499  for (; cell != endc; ++cell)
4500  if (cell->level_subdomain_id() == this->locally_owned_subdomain())
4501  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
4502  marked_vertices[cell->vertex_index(v)] = true;
4503 
4509  typename std::map<std::pair<cell_iterator, unsigned int>,
4510  std::pair<std::pair<cell_iterator, unsigned int>,
4511  std::bitset<3>>>::const_iterator it;
4512 
4513  // When a connectivity in the code below is detected, the assignment
4514  // 'marked_vertices[v1] = marked_vertices[v2] = true' makes sure that
4515  // the information about the periodicity propagates back to vertices on
4516  // cells that are not owned locally. However, in the worst case we want
4517  // to connect to a vertex that is 'dim' hops away from the locally owned
4518  // cell. Depending on the order of the periodic face map, we might
4519  // connect to that point by chance or miss it. However, after looping
4520  // through all the periodic directions (which are at most as many as
4521  // the number of space dimensions) we can be sure that all connections
4522  // to vertices have been created.
4523  for (unsigned int repetition = 0; repetition < dim; ++repetition)
4524  for (it = this->get_periodic_face_map().begin();
4525  it != this->get_periodic_face_map().end();
4526  ++it)
4527  {
4528  const cell_iterator & cell_1 = it->first.first;
4529  const unsigned int face_no_1 = it->first.second;
4530  const cell_iterator & cell_2 = it->second.first.first;
4531  const unsigned int face_no_2 = it->second.first.second;
4532  const std::bitset<3> &face_orientation = it->second.second;
4533 
4534  if (cell_1->level() == level && cell_2->level() == level)
4535  {
4536  for (unsigned int v = 0;
4537  v < GeometryInfo<dim - 1>::vertices_per_cell;
4538  ++v)
4539  {
4540  // take possible non-standard orientation of faces into
4541  // account
4542  const unsigned int vface0 =
4544  v,
4545  face_orientation[0],
4546  face_orientation[1],
4547  face_orientation[2]);
4548  if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
4549  vface0)] ||
4550  marked_vertices[cell_2->face(face_no_2)->vertex_index(
4551  v)])
4552  marked_vertices[cell_1->face(face_no_1)->vertex_index(
4553  vface0)] =
4554  marked_vertices[cell_2->face(face_no_2)->vertex_index(
4555  v)] = true;
4556  }
4557  }
4558  }
4559 
4560  return marked_vertices;
4561  }
4562 
4563 
4564 
4565  template <int dim, int spacedim>
4566  unsigned int
4569  {
4570  return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
4571  }
4572 
4573 
4574 
4575  template <int dim, int spacedim>
4578  const unsigned int coarse_cell_index) const
4579  {
4580  return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
4581  }
4582 
4583 
4584 
4585  template <int dim, int spacedim>
4586  void
4588  const std::vector<::GridTools::PeriodicFacePair<cell_iterator>>
4589  &periodicity_vector)
4590  {
4591  Assert(triangulation_has_content == true,
4592  ExcMessage("The triangulation is empty!"));
4593  Assert(this->n_levels() == 1,
4594  ExcMessage("The triangulation is refined!"));
4595 
4596  // call the base class for storing the periodicity information; we must
4597  // do this before going to p4est and rebuilding the triangulation to get
4598  // the level subdomain ids correct in the multigrid case
4600 
4601  for (const auto &face_pair : periodicity_vector)
4602  {
4603  const cell_iterator first_cell = face_pair.cell[0];
4604  const cell_iterator second_cell = face_pair.cell[1];
4605  const unsigned int face_left = face_pair.face_idx[0];
4606  const unsigned int face_right = face_pair.face_idx[1];
4607 
4608  // respective cells of the matching faces in p4est
4609  const unsigned int tree_left =
4610  coarse_cell_to_p4est_tree_permutation[first_cell->index()];
4611  const unsigned int tree_right =
4612  coarse_cell_to_p4est_tree_permutation[second_cell->index()];
4613 
4614  // p4est wants to know which corner the first corner on
4615  // the face with the lower id is mapped to on the face with
4616  // with the higher id. For d==2 there are only two possibilities
4617  // that are determined by it->orientation[1].
4618  // For d==3 we have to use GridTools::OrientationLookupTable.
4619  // The result is given below.
4620 
4621  unsigned int p4est_orientation = 0;
4622  if (dim == 2)
4623  p4est_orientation = face_pair.orientation[1];
4624  else
4625  {
4626  const unsigned int face_idx_list[] = {face_left, face_right};
4627  const cell_iterator cell_list[] = {first_cell, second_cell};
4628  unsigned int lower_idx, higher_idx;
4629  if (face_left <= face_right)
4630  {
4631  higher_idx = 1;
4632  lower_idx = 0;
4633  }
4634  else
4635  {
4636  higher_idx = 0;
4637  lower_idx = 1;
4638  }
4639 
4640  // get the cell index of the first index on the face with the
4641  // lower id
4642  unsigned int first_p4est_idx_on_cell =
4643  p8est_face_corners[face_idx_list[lower_idx]][0];
4644  unsigned int first_dealii_idx_on_face =
4646  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
4647  ++i)
4648  {
4649  const unsigned int first_dealii_idx_on_cell =
4651  face_idx_list[lower_idx],
4652  i,
4653  cell_list[lower_idx]->face_orientation(
4654  face_idx_list[lower_idx]),
4655  cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
4656  cell_list[lower_idx]->face_rotation(
4657  face_idx_list[lower_idx]));
4658  if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
4659  {
4660  first_dealii_idx_on_face = i;
4661  break;
4662  }
4663  }
4664  Assert(first_dealii_idx_on_face != numbers::invalid_unsigned_int,
4665  ExcInternalError());
4666  // Now map dealii_idx_on_face according to the orientation
4667  constexpr unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
4668  {0, 1, 2, 3},
4669  {3, 1, 2, 0},
4670  {3, 2, 1, 0},
4671  {2, 3, 0, 1},
4672  {1, 3, 0, 2},
4673  {1, 0, 3, 2},
4674  {2, 0, 3, 1}};
4675  constexpr unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
4676  {0, 1, 2, 3},
4677  {3, 1, 2, 0},
4678  {3, 2, 1, 0},
4679  {2, 3, 0, 1},
4680  {2, 0, 3, 1},
4681  {1, 0, 3, 2},
4682  {1, 3, 0, 2}};
4683  const unsigned int second_dealii_idx_on_face =
4684  lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
4685  [first_dealii_idx_on_face] :
4686  right_to_left[face_pair.orientation.to_ulong()]
4687  [first_dealii_idx_on_face];
4688  const unsigned int second_dealii_idx_on_cell =
4690  face_idx_list[higher_idx],
4691  second_dealii_idx_on_face,
4692  cell_list[higher_idx]->face_orientation(
4693  face_idx_list[higher_idx]),
4694  cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
4695  cell_list[higher_idx]->face_rotation(
4696  face_idx_list[higher_idx]));
4697  // map back to p4est
4698  const unsigned int second_p4est_idx_on_face =
4699  p8est_corner_face_corners[second_dealii_idx_on_cell]
4700  [face_idx_list[higher_idx]];
4701  p4est_orientation = second_p4est_idx_on_face;
4702  }
4703 
4705  connectivity,
4706  tree_left,
4707  tree_right,
4708  face_left,
4709  face_right,
4710  p4est_orientation);
4711  }
4712 
4713 
4715  connectivity) == 1,
4716  ExcInternalError());
4717 
4718  // now create a forest out of the connectivity data structure
4721  this->mpi_communicator,
4722  connectivity,
4723  /* minimum initial number of quadrants per tree */ 0,
4724  /* minimum level of upfront refinement */ 0,
4725  /* use uniform upfront refinement */ 1,
4726  /* user_data_size = */ 0,
4727  /* user_data_constructor = */ nullptr,
4728  /* user_pointer */ this);
4729 
4730  try
4731  {
4732  copy_local_forest_to_triangulation();
4733  }
4734  catch (const typename Triangulation<dim>::DistortedCellList &)
4735  {
4736  // the underlying triangulation should not be checking for distorted
4737  // cells
4738  Assert(false, ExcInternalError());
4739  }
4740 
4741  // The range of ghost_owners might have changed so update that information
4742  this->update_number_cache();
4743  }
4744 
4745 
4746 
4747  template <int dim, int spacedim>
4748  std::size_t
4750  {
4751  std::size_t mem =
4754  MemoryConsumption::memory_consumption(triangulation_has_content) +
4756  MemoryConsumption::memory_consumption(parallel_forest) +
4758  cell_attached_data.n_attached_data_sets) +
4759  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_fixed)
4760  // +
4761  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_variable)
4762  // +
4763  // TODO[TH]: how?
4765  coarse_cell_to_p4est_tree_permutation) +
4767  p4est_tree_to_coarse_cell_permutation) +
4768  memory_consumption_p4est();
4769 
4770  return mem;
4771  }
4772 
4773 
4774 
4775  template <int dim, int spacedim>
4776  std::size_t
4778  {
4779  return ::internal::p4est::functions<dim>::forest_memory_used(
4780  parallel_forest) +
4782  connectivity);
4783  }
4784 
4785 
4786 
4787  template <int dim, int spacedim>
4788  void
4790  const ::Triangulation<dim, spacedim> &other_tria)
4791  {
4792  try
4793  {
4795  copy_triangulation(other_tria);
4796  }
4797  catch (
4798  const typename ::Triangulation<dim, spacedim>::DistortedCellList
4799  &)
4800  {
4801  // the underlying triangulation should not be checking for distorted
4802  // cells
4803  Assert(false, ExcInternalError());
4804  }
4805 
4806  // note that now we have some content in the p4est objects and call the
4807  // functions that do the actual work (which are dimension dependent, so
4808  // separate)
4809  triangulation_has_content = true;
4810 
4811  Assert(other_tria.n_levels() == 1,
4812  ExcMessage(
4813  "Parallel distributed triangulations can only be copied, "
4814  "if they are not refined!"));
4815 
4816  if (const ::parallel::distributed::Triangulation<dim, spacedim>
4817  *other_tria_x =
4818  dynamic_cast<const ::parallel::distributed::
4819  Triangulation<dim, spacedim> *>(&other_tria))
4820  {
4821  coarse_cell_to_p4est_tree_permutation =
4822  other_tria_x->coarse_cell_to_p4est_tree_permutation;
4823  p4est_tree_to_coarse_cell_permutation =
4824  other_tria_x->p4est_tree_to_coarse_cell_permutation;
4825  cell_attached_data = other_tria_x->cell_attached_data;
4826  data_transfer = other_tria_x->data_transfer;
4827 
4828  settings = other_tria_x->settings;
4829  }
4830  else
4831  {
4832  setup_coarse_cell_to_p4est_tree_permutation();
4833  }
4834 
4835  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
4836 
4837  try
4838  {
4839  copy_local_forest_to_triangulation();
4840  }
4841  catch (const typename Triangulation<dim>::DistortedCellList &)
4842  {
4843  // the underlying triangulation should not be checking for distorted
4844  // cells
4845  Assert(false, ExcInternalError());
4846  }
4847 
4848  this->update_periodic_face_map();
4849  this->update_number_cache();
4850  }
4851 
4852 
4853 
4854  template <int dim, int spacedim>
4855  void
4857  {
4858  // reorganize memory for local_quadrant_cell_relations
4859  local_quadrant_cell_relations.resize(
4860  parallel_forest->local_num_quadrants);
4861  local_quadrant_cell_relations.shrink_to_fit();
4862 
4863  // recurse over p4est
4864  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
4865  this->begin(0);
4866  cell != this->end(0);
4867  ++cell)
4868  {
4869  // skip coarse cells that are not ours
4870  if (tree_exists_locally<dim, spacedim>(
4871  parallel_forest,
4872  coarse_cell_to_p4est_tree_permutation[cell->index()]) == false)
4873  continue;
4874 
4875  // initialize auxiliary top level p4est quadrant
4876  typename ::internal::p4est::types<dim>::quadrant
4877  p4est_coarse_cell;
4878  ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4879 
4880  // determine tree to start recursion on
4881  typename ::internal::p4est::types<dim>::tree *tree =
4882  init_tree(cell->index());
4883 
4884  update_quadrant_cell_relations_recursively<dim, spacedim>(
4885  local_quadrant_cell_relations, *tree, cell, p4est_coarse_cell);
4886  }
4887  }
4888 
4889 
4890 
4891  template <int dim, int spacedim>
4892  std::vector<unsigned int>
4894  {
4895  // check if local_quadrant_cell_relations have been previously gathered
4896  // correctly
4897  Assert(local_quadrant_cell_relations.size() ==
4898  static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4899  ExcInternalError());
4900 
4901  // Allocate the space for the weights. In fact we do not know yet, how
4902  // many cells we own after the refinement (only p4est knows that
4903  // at this point). We simply reserve n_active_cells space and if many
4904  // more cells are refined than coarsened than additional reallocation
4905  // will be done inside get_cell_weights_recursively.
4906  std::vector<unsigned int> weights;
4907  weights.reserve(this->n_active_cells());
4908 
4909  // Iterate over p4est and Triangulation relations
4910  // to find refined/coarsened/kept
4911  // cells. Then append cell_weight.
4912  // Note that we need to follow the p4est ordering
4913  // instead of the deal.II ordering to get the cell_weights
4914  // in the same order p4est will encounter them during repartitioning.
4915  for (const auto &quad_cell_rel : local_quadrant_cell_relations)
4916  {
4917  const auto &cell_status = std::get<1>(quad_cell_rel);
4918  const auto &cell_it = std::get<2>(quad_cell_rel);
4919 
4920  switch (cell_status)
4921  {
4923  spacedim>::CELL_PERSIST:
4924  weights.push_back(1000);
4925  weights.back() += this->signals.cell_weight(
4926  cell_it,
4928  spacedim>::CELL_PERSIST);
4929  break;
4930 
4932  spacedim>::CELL_REFINE:
4934  spacedim>::CELL_INVALID:
4935  {
4936  // calculate weight of parent cell
4937  unsigned int parent_weight = 1000;
4938  parent_weight += this->signals.cell_weight(
4939  cell_it,
4941  CELL_REFINE);
4942  // assign the weight of the parent cell equally to all
4943  // children
4944  weights.push_back(parent_weight);
4945  break;
4946  }
4947 
4949  spacedim>::CELL_COARSEN:
4950  weights.push_back(1000);
4951  weights.back() += this->signals.cell_weight(
4952  cell_it,
4954  spacedim>::CELL_COARSEN);
4955  break;
4956 
4957  default:
4958  Assert(false, ExcInternalError());
4959  break;
4960  }
4961  }
4962 
4963  return weights;
4964  }
4965 
4966 
4967 
4968  template <int spacedim>
4970  MPI_Comm mpi_communicator,
4971  const typename ::Triangulation<1, spacedim>::MeshSmoothing
4972  smooth_grid,
4973  const Settings /*settings*/)
4974  : ::parallel::DistributedTriangulationBase<1, spacedim>(
4975  mpi_communicator,
4976  smooth_grid,
4977  false)
4978  {
4979  Assert(false, ExcNotImplemented());
4980  }
4981 
4982 
4983  template <int spacedim>
4985  {
4986  AssertNothrow(false, ExcNotImplemented());
4987  }
4988 
4989 
4990 
4991  template <int spacedim>
4992  void
4994  const std::vector<bool> & /*vertex_locally_moved*/)
4995  {
4996  Assert(false, ExcNotImplemented());
4997  }
4998 
4999 
5000 
5001  template <int spacedim>
5002  unsigned int
5004  const std::function<std::vector<char>(
5005  const typename ::Triangulation<1, spacedim>::cell_iterator &,
5006  const typename ::Triangulation<1, spacedim>::CellStatus)>
5007  & /*pack_callback*/,
5008  const bool /*returns_variable_size_data*/)
5009  {
5010  Assert(false, ExcNotImplemented());
5011  return 0;
5012  }
5013 
5014 
5015 
5016  template <int spacedim>
5017  void
5019  const unsigned int /*handle*/,
5020  const std::function<
5021  void(const typename ::Triangulation<1, spacedim>::cell_iterator &,
5022  const typename ::Triangulation<1, spacedim>::CellStatus,
5023  const boost::iterator_range<std::vector<char>::const_iterator> &)>
5024  & /*unpack_callback*/)
5025  {
5026  Assert(false, ExcNotImplemented());
5027  }
5028 
5029 
5030 
5031  template <int spacedim>
5032  const std::vector<types::global_dof_index> &
5034  const
5035  {
5036  static std::vector<types::global_dof_index> a;
5037  return a;
5038  }
5039 
5040 
5041 
5042  template <int spacedim>
5043  std::map<unsigned int, std::set<::types::subdomain_id>>
5045  const unsigned int /*level*/) const
5046  {
5047  Assert(false, ExcNotImplemented());
5048 
5049  return std::map<unsigned int, std::set<::types::subdomain_id>>();
5050  }
5051 
5052 
5053 
5054  template <int spacedim>
5055  std::vector<bool>
5057  const unsigned int) const
5058  {
5059  Assert(false, ExcNotImplemented());
5060  return std::vector<bool>();
5061  }
5062 
5063 
5064 
5065  template <int spacedim>
5066  unsigned int
5068  const types::coarse_cell_id) const
5069  {
5070  Assert(false, ExcNotImplemented());
5071  return 0;
5072  }
5073 
5074 
5075 
5076  template <int spacedim>
5079  const unsigned int) const
5080  {
5081  Assert(false, ExcNotImplemented());
5082  return 0;
5083  }
5084 
5085 
5086  template <int spacedim>
5087  void
5088  Triangulation<1, spacedim>::load(const std::string &, const bool)
5089  {
5090  Assert(false, ExcNotImplemented());
5091  }
5092 
5093 
5094 
5095  template <int spacedim>
5096  void
5097  Triangulation<1, spacedim>::save(const std::string &) const
5098  {
5099  Assert(false, ExcNotImplemented());
5100  }
5101 
5102 
5103 
5104  template <int spacedim>
5105  bool
5107  {
5108  Assert(false, ExcNotImplemented());
5109  return false;
5110  }
5111 
5112  } // namespace distributed
5113 } // namespace parallel
5114 
5115 
5116 #endif // DEAL_II_WITH_P4EST
5117 
5118 
5119 
5120 /*-------------- Explicit Instantiations -------------------------------*/
5121 #include "tria.inst"
5122 
5123 
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const typename identity< BaseIterator >::type &end)
bool is_multilevel_hierarchy_constructed() const override
Definition: tria.cc:2604
unsigned int n_active_cells() const
Definition: tria.cc:11776
virtual unsigned int coarse_cell_id_to_coarse_cell_index(const types::coarse_cell_id coarse_cell_id) const
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
Definition: tria.cc:9534
unsigned int n_vertices() const
unsigned int register_data_attach(const std::function< std::vector< char >(const cell_iterator &, const CellStatus)> &pack_callback, const bool returns_variable_size_data)
Definition: tria.cc:4360
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1474
bool tree_exists_locally(const typename types< dim >::forest *parallel_forest, const typename types< dim >::topidx coarse_grid_cell)
boost::signals2::signal< void()> post_distributed_repartition
Definition: tria.h:2187
virtual bool has_hanging_nodes() const
Definition: tria.cc:11908
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
void load(const std::string &filename, const bool autopartition=true)
Definition: tria.cc:2751
void save(const typename ::internal::p4est::types< dim >::forest *parallel_forest, const std::string &filename) const
Definition: tria.cc:1765
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition: tria_base.cc:63
unsigned int n_cells() const
Definition: tria.cc:11768
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void unpack_cell_status(std::vector< quadrant_cell_relation_t > &quad_cell_relations) const
Definition: tria.cc:1562
DataTransfer(MPI_Comm mpi_communicator)
Definition: tria.cc:1125
static ::ExceptionBase & ExcIO()
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
std::vector< Point< spacedim > > vertices
Definition: tria.h:3849
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:11286
unsigned int compute_n_point_to_point_communications(const MPI_Comm &mpi_comm, const std::vector< unsigned int > &destinations)
Definition: mpi.cc:479
boost::signals2::signal< void()> pre_distributed_refinement
Definition: tria.h:2158
void pack_data(const std::vector< quadrant_cell_relation_t > &quad_cell_relations, const std::vector< typename CellAttachedData::pack_callback_t > &pack_callbacks_fixed, const std::vector< typename CellAttachedData::pack_callback_t > &pack_callbacks_variable)
Definition: tria.cc:1135
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11114
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
std::vector< unsigned int * > first_vertex_indices
Definition: tria.cc:2251
void save(const std::string &filename) const
Definition: tria.cc:2671
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11094
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
Definition: tria.cc:4472
unsigned int n_levels() const
bool is_locally_owned() const
std::vector<::Point< spacedim > * > first_vertices
Definition: tria.cc:2252
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2505
cell_iterator end() const
Definition: tria.cc:11180
virtual unsigned int coarse_cell_id_to_coarse_cell_index(const types::coarse_cell_id coarse_cell_id) const override
Definition: tria.cc:4567
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator >> &)
Definition: tria.cc:12414
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:12443
virtual std::size_t memory_consumption() const override
Definition: tria_base.cc:88
unsigned int n_active_lines() const
Definition: tria.cc:11968
virtual bool prepare_coarsening_and_refinement()
Definition: tria.cc:13139
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
Triangulation(const MeshSmoothing smooth_grid=none, const bool check_for_distorted_cells=false)
Definition: tria.cc:9213
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const override
Definition: tria.cc:4577
T sum(const T &t, const MPI_Comm &mpi_communicator)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
typename ::Triangulation< dim, spacedim >::CellStatus CellStatus
Definition: tria.h:287
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:9617
#define Assert(cond, exc)
Definition: exceptions.h:1411
Signals signals
Definition: tria.h:2223
void load(const typename ::internal::p4est::types< dim >::forest *parallel_forest, const std::string &filename, const unsigned int n_attached_deserialize_fixed, const unsigned int n_attached_deserialize_variable)
Definition: tria.cc:1928
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void load(Archive &ar, const unsigned int version)
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const
void save_coarsen_flags(std::ostream &out) const
Definition: tria.cc:10057
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
void notify_ready_to_unpack(const unsigned int handle, const std::function< void(const cell_iterator &, const CellStatus, const boost::iterator_range< std::vector< char >::const_iterator > &)> &unpack_callback)
Definition: tria.cc:4390
unsigned int level
Definition: grid_out.cc:4341
void save_refine_flags(std::ostream &out) const
Definition: tria.cc:9988
boost::signals2::signal< void()> pre_distributed_save
Definition: tria.h:2195
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1419
virtual std::size_t memory_consumption() const
Definition: tria.cc:14227
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
void save(Archive &ar, const unsigned int version) const
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1182
std::vector< typename ::internal::p4est::types< dim >::quadrant > quadrants
Definition: tria.cc:2241
std::vector< unsigned int > tree_index
Definition: tria.cc:2237
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
#define DEAL_II_MPI_CONST_CAST(expr)
Definition: mpi.h:76
void reorder_hierarchical(const DynamicSparsityPattern &sparsity, std::vector< DynamicSparsityPattern::size_type > &new_indices)
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
boost::signals2::signal< void()> post_distributed_save
Definition: tria.h:2202
void communicate_locally_moved_vertices(const std::vector< bool > &vertex_locally_moved)
Definition: tria.cc:4176
void unpack_data(const std::vector< quadrant_cell_relation_t > &quad_cell_relations, const unsigned int handle, const std::function< void(const typename ::Triangulation< dim, spacedim >::cell_iterator &, const typename ::Triangulation< dim, spacedim >::CellStatus &, const boost::iterator_range< std::vector< char >::const_iterator > &)> &unpack_callback) const
Definition: tria.cc:1598
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1692
active_cell_iterator end_active(const unsigned int level) const
Definition: tria.cc:11249
void update_periodic_face_map()
Definition: tria.cc:12502
const std::vector< bool > & get_used_vertices() const
Definition: tria.cc:12320
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
virtual unsigned int n_global_levels() const
T min(const T &t, const MPI_Comm &mpi_communicator)
virtual ~Triangulation() override
Definition: tria.cc:2142
MeshSmoothing smooth_grid
Definition: tria.h:3420
virtual ~Triangulation() override
Definition: tria.cc:9287
Triangulation<dim, spacedim>::communicate_locally_moved_vertices()
Definition: mpi_tags.h:63
boost::signals2::signal< unsigned int(const cell_iterator &, const CellStatus), CellWeightSum< unsigned int > > cell_weight
Definition: tria.h:2145
boost::signals2::signal< void()> pre_distributed_load
Definition: tria.h:2210
void execute_transfer(const typename ::internal::p4est::types< dim >::forest *parallel_forest, const typename ::internal::p4est::types< dim >::gloidx *previous_global_first_quadrant)
Definition: tria.cc:1469
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1326
std::vector< unsigned int > sizes_fixed_cumulative
Definition: tria.h:1113
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
static std::mutex mutex
Definition: subscriptor.h:261
static ::ExceptionBase & ExcNotImplemented()
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2959
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
active_cell_iterator last_active() const
Definition: tria.cc:11159
Triangulation(MPI_Comm mpi_communicator, const typename ::Triangulation< dim, spacedim >::MeshSmoothing smooth_grid=(::Triangulation< dim, spacedim >::none), const Settings settings=default_setting)
Definition: tria.cc:2112
std::map< unsigned int, std::set<::types::subdomain_id > > * vertices_with_ghost_neighbors
boost::signals2::signal< void()> post_distributed_load
Definition: tria.h:2217
global_cell_index coarse_cell_id
Definition: types.h:114
boost::signals2::signal< void()> pre_distributed_repartition
Definition: tria.h:2180
CellAttachedData cell_attached_data
Definition: tria.h:923
std::vector< bool > mark_locally_active_vertices_on_level(const int level) const
Definition: tria.cc:4492
T max(const T &t, const MPI_Comm &mpi_communicator)
virtual types::subdomain_id locally_owned_subdomain() const
Definition: tria.cc:12387
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
virtual void clear()
Definition: tria.cc:9316
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > & get_periodic_face_map() const
Definition: tria.cc:12434
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:5540
boost::signals2::signal< void()> post_distributed_refinement
Definition: tria.h:2168