Reference documentation for deal.II version 9.2.0
\(\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  do
3443  {
3446  this->update_periodic_face_map();
3447  // enforce 2:1 mesh balance over periodic boundaries
3448  mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this);
3449  }
3450  while (mesh_changed);
3451 
3452  // check if any of the refinement flags were changed during this
3453  // function and return that value
3454  std::vector<bool> flags_after[2];
3455  this->save_coarsen_flags(flags_after[0]);
3456  this->save_refine_flags(flags_after[1]);
3457  return ((flags_before[0] != flags_after[0]) ||
3458  (flags_before[1] != flags_after[1]));
3459  }
3460 
3461 
3462 
3463  template <int dim, int spacedim>
3464  void
3466  {
3467  // disable mesh smoothing for recreating the deal.II triangulation,
3468  // otherwise we might not be able to reproduce the p4est mesh
3469  // exactly. We restore the original smoothing at the end of this
3470  // function. Note that the smoothing flag is used in the normal
3471  // refinement process.
3472  typename Triangulation<dim, spacedim>::MeshSmoothing save_smooth =
3473  this->smooth_grid;
3474 
3475  // We will refine manually to match the p4est further down, which
3476  // obeys a level difference of 2 at each vertex (see the balance call
3477  // to p4est). We can disable this here so we store fewer artificial
3478  // cells (in some cases).
3479  // For geometric multigrid it turns out that
3480  // we will miss level cells at shared vertices if we ignore this.
3481  // See tests/mpi/mg_06. In particular, the flag is still necessary
3482  // even though we force it for the original smooth_grid in the
3483  // constructor.
3484  if (settings & construct_multigrid_hierarchy)
3485  this->smooth_grid =
3486  ::Triangulation<dim,
3488  else
3490 
3491  bool mesh_changed = false;
3492 
3493  // remove all deal.II refinements. Note that we could skip this and
3494  // start from our current state, because the algorithm later coarsens as
3495  // necessary. This has the advantage of being faster when large parts
3496  // of the local partition changes (likely) and gives a deterministic
3497  // ordering of the cells (useful for snapshot/resume).
3498  // TODO: is there a more efficient way to do this?
3499  if (settings & mesh_reconstruction_after_repartitioning)
3500  while (this->begin_active()->level() > 0)
3501  {
3502  for (const auto &cell : this->active_cell_iterators())
3503  {
3504  cell->set_coarsen_flag();
3505  }
3506 
3508  try
3509  {
3512  }
3513  catch (
3515  {
3516  // the underlying triangulation should not be checking for
3517  // distorted cells
3518  Assert(false, ExcInternalError());
3519  }
3520  }
3521 
3522 
3523  // query p4est for the ghost cells
3524  if (parallel_ghost != nullptr)
3525  {
3527  parallel_ghost);
3528  parallel_ghost = nullptr;
3529  }
3531  parallel_forest,
3532  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3533  P4EST_CONNECT_CORNER) :
3534  typename ::internal::p4est::types<dim>::balance_type(
3535  P8EST_CONNECT_CORNER)));
3536 
3537  Assert(parallel_ghost, ExcInternalError());
3538 
3539 
3540  // set all cells to artificial. we will later set it to the correct
3541  // subdomain in match_tree_recursively
3542  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3543  this->begin(0);
3544  cell != this->end(0);
3545  ++cell)
3546  cell->recursively_set_subdomain_id(numbers::artificial_subdomain_id);
3547 
3548  do
3549  {
3550  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3551  this->begin(0);
3552  cell != this->end(0);
3553  ++cell)
3554  {
3555  // if this processor stores no part of the forest that comes out
3556  // of this coarse grid cell, then we need to delete all children
3557  // of this cell (the coarse grid cell remains)
3558  if (tree_exists_locally<dim, spacedim>(
3559  parallel_forest,
3560  coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
3561  false)
3562  {
3563  delete_all_children<dim, spacedim>(cell);
3564  if (cell->is_active())
3565  cell->set_subdomain_id(numbers::artificial_subdomain_id);
3566  }
3567 
3568  else
3569  {
3570  // this processor stores at least a part of the tree that
3571  // comes out of this cell.
3572 
3573  typename ::internal::p4est::types<dim>::quadrant
3574  p4est_coarse_cell;
3575  typename ::internal::p4est::types<dim>::tree *tree =
3576  init_tree(cell->index());
3577 
3578  ::internal::p4est::init_coarse_quadrant<dim>(
3579  p4est_coarse_cell);
3580 
3581  match_tree_recursively<dim, spacedim>(*tree,
3582  cell,
3583  p4est_coarse_cell,
3584  *parallel_forest,
3585  this->my_subdomain);
3586  }
3587  }
3588 
3589  // check mesh for ghost cells, refine as necessary. iterate over
3590  // every ghostquadrant, find corresponding deal coarsecell and
3591  // recurse.
3592  typename ::internal::p4est::types<dim>::quadrant *quadr;
3593  types::subdomain_id ghost_owner = 0;
3594  typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
3595 
3596  for (unsigned int g_idx = 0;
3597  g_idx < parallel_ghost->ghosts.elem_count;
3598  ++g_idx)
3599  {
3600  while (g_idx >= static_cast<unsigned int>(
3601  parallel_ghost->proc_offsets[ghost_owner + 1]))
3602  ++ghost_owner;
3603  while (g_idx >= static_cast<unsigned int>(
3604  parallel_ghost->tree_offsets[ghost_tree + 1]))
3605  ++ghost_tree;
3606 
3607  quadr = static_cast<
3608  typename ::internal::p4est::types<dim>::quadrant *>(
3609  sc_array_index(&parallel_ghost->ghosts, g_idx));
3610 
3611  unsigned int coarse_cell_index =
3612  p4est_tree_to_coarse_cell_permutation[ghost_tree];
3613 
3614  match_quadrant<dim, spacedim>(this,
3615  coarse_cell_index,
3616  *quadr,
3617  ghost_owner);
3618  }
3619 
3620  // fix all the flags to make sure we have a consistent mesh
3622 
3623  // see if any flags are still set
3624  mesh_changed =
3625  std::any_of(this->begin_active(),
3626  active_cell_iterator{this->end()},
3627  [](const CellAccessor<dim, spacedim> &cell) {
3628  return cell.refine_flag_set() ||
3629  cell.coarsen_flag_set();
3630  });
3631 
3632  // actually do the refinement to change the local mesh by
3633  // calling the base class refinement function directly
3634  try
3635  {
3638  }
3639  catch (
3641  {
3642  // the underlying triangulation should not be checking for
3643  // distorted cells
3644  Assert(false, ExcInternalError());
3645  }
3646  }
3647  while (mesh_changed);
3648 
3649 # ifdef DEBUG
3650  // check if correct number of ghosts is created
3651  unsigned int num_ghosts = 0;
3652 
3653  for (const auto &cell : this->active_cell_iterators())
3654  {
3655  if (cell->subdomain_id() != this->my_subdomain &&
3656  cell->subdomain_id() != numbers::artificial_subdomain_id)
3657  ++num_ghosts;
3658  }
3659 
3660  Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
3661  ExcInternalError());
3662 # endif
3663 
3664 
3665 
3666  // fill level_subdomain_ids for geometric multigrid
3667  // the level ownership of a cell is defined as the owner if the cell is
3668  // active or as the owner of child(0) we need this information for all our
3669  // ancestors and the same-level neighbors of our own cells (=level ghosts)
3670  if (settings & construct_multigrid_hierarchy)
3671  {
3672  // step 1: We set our own ids all the way down and all the others to
3673  // -1. Note that we do not fill other cells we could figure out the
3674  // same way, because we might accidentally set an id for a cell that
3675  // is not a ghost cell.
3676  for (unsigned int lvl = this->n_levels(); lvl > 0;)
3677  {
3678  --lvl;
3680  endc = this->end(lvl);
3681  for (cell = this->begin(lvl); cell != endc; ++cell)
3682  {
3683  if ((cell->is_active() &&
3684  cell->subdomain_id() ==
3685  this->locally_owned_subdomain()) ||
3686  (cell->has_children() &&
3687  cell->child(0)->level_subdomain_id() ==
3688  this->locally_owned_subdomain()))
3689  cell->set_level_subdomain_id(
3690  this->locally_owned_subdomain());
3691  else
3692  {
3693  // not our cell
3694  cell->set_level_subdomain_id(
3696  }
3697  }
3698  }
3699 
3700  // step 2: make sure all the neighbors to our level_cells exist. Need
3701  // to look up in p4est...
3702  std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3703  for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3704  marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3705 
3706  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
3707  this->begin(0);
3708  cell != this->end(0);
3709  ++cell)
3710  {
3711  typename ::internal::p4est::types<dim>::quadrant
3712  p4est_coarse_cell;
3713  const unsigned int tree_index =
3714  coarse_cell_to_p4est_tree_permutation[cell->index()];
3715  typename ::internal::p4est::types<dim>::tree *tree =
3716  init_tree(cell->index());
3717 
3718  ::internal::p4est::init_coarse_quadrant<dim>(
3719  p4est_coarse_cell);
3720 
3721  determine_level_subdomain_id_recursively<dim, spacedim>(
3722  *tree,
3723  tree_index,
3724  cell,
3725  p4est_coarse_cell,
3726  *parallel_forest,
3727  this->my_subdomain,
3728  marked_vertices);
3729  }
3730 
3731  // step 3: make sure we have the parent of our level cells
3732  for (unsigned int lvl = this->n_levels(); lvl > 0;)
3733  {
3734  --lvl;
3736  endc = this->end(lvl);
3737  for (cell = this->begin(lvl); cell != endc; ++cell)
3738  {
3739  if (cell->has_children())
3740  for (unsigned int c = 0;
3741  c < GeometryInfo<dim>::max_children_per_cell;
3742  ++c)
3743  {
3744  if (cell->child(c)->level_subdomain_id() ==
3745  this->locally_owned_subdomain())
3746  {
3747  // at least one of the children belongs to us, so
3748  // make sure we set the level subdomain id
3749  const types::subdomain_id mark =
3750  cell->child(0)->level_subdomain_id();
3752  ExcInternalError()); // we should know the
3753  // child(0)
3754  cell->set_level_subdomain_id(mark);
3755  break;
3756  }
3757  }
3758  }
3759  }
3760  }
3761 
3762 
3763 
3764  // check that our local copy has exactly as many cells as the p4est
3765  // original (at least if we are on only one processor); for parallel
3766  // computations, we want to check that we have at least as many as p4est
3767  // stores locally (in the future we should check that we have exactly as
3768  // many non-artificial cells as parallel_forest->local_num_quadrants)
3769  {
3770  const unsigned int total_local_cells = this->n_active_cells();
3771  (void)total_local_cells;
3772 
3773  if (Utilities::MPI::n_mpi_processes(this->mpi_communicator) == 1)
3774  {
3775  Assert(static_cast<unsigned int>(
3776  parallel_forest->local_num_quadrants) == total_local_cells,
3777  ExcInternalError());
3778  }
3779  else
3780  {
3781  Assert(static_cast<unsigned int>(
3782  parallel_forest->local_num_quadrants) <= total_local_cells,
3783  ExcInternalError());
3784  }
3785 
3786  // count the number of owned, active cells and compare with p4est.
3787  unsigned int n_owned = 0;
3788  for (const auto &cell : this->active_cell_iterators())
3789  {
3790  if (cell->subdomain_id() == this->my_subdomain)
3791  ++n_owned;
3792  }
3793 
3794  Assert(static_cast<unsigned int>(
3795  parallel_forest->local_num_quadrants) == n_owned,
3796  ExcInternalError());
3797  }
3798 
3799  this->smooth_grid = save_smooth;
3800 
3801  // finally, after syncing the parallel_forest with the triangulation,
3802  // also update the quadrant_cell_relations, which will be used for
3803  // repartitioning, further refinement/coarsening, and unpacking
3804  // of stored or transferred data.
3805  update_quadrant_cell_relations();
3806  }
3807 
3808 
3809 
3810  template <int dim, int spacedim>
3811  void
3813  {
3814  // do not allow anisotropic refinement
3815 # ifdef DEBUG
3816  for (const auto &cell : this->active_cell_iterators())
3817  if (cell->is_locally_owned() && cell->refine_flag_set())
3818  Assert(cell->refine_flag_set() ==
3820  ExcMessage(
3821  "This class does not support anisotropic refinement"));
3822 # endif
3823 
3824 
3825  // safety check: p4est has an upper limit on the level of a cell
3826  if (this->n_levels() ==
3828  {
3830  cell = this->begin_active(
3832  cell !=
3834  1);
3835  ++cell)
3836  {
3837  AssertThrow(
3838  !(cell->refine_flag_set()),
3839  ExcMessage(
3840  "Fatal Error: maximum refinement level of p4est reached."));
3841  }
3842  }
3843 
3845 
3846  // signal that refinement is going to happen
3848 
3849  // now do the work we're supposed to do when we are in charge
3850  // make sure all flags are cleared on cells we don't own, since nothing
3851  // good can come of that if they are still around
3852  for (const auto &cell : this->active_cell_iterators())
3853  if (cell->is_ghost() || cell->is_artificial())
3854  {
3855  cell->clear_refine_flag();
3856  cell->clear_coarsen_flag();
3857  }
3858 
3859 
3860  // count how many cells will be refined and coarsened, and allocate that
3861  // much memory
3862  RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3863  *this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3864 
3865  // copy refine and coarsen flags into p4est and execute the refinement
3866  // and coarsening. this uses the refine_and_coarsen_list just built,
3867  // which is communicated to the callback functions through
3868  // p4est's user_pointer object
3869  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3870  parallel_forest->user_pointer = &refine_and_coarsen_list;
3871 
3872  if (parallel_ghost != nullptr)
3873  {
3875  parallel_ghost);
3876  parallel_ghost = nullptr;
3877  }
3879  parallel_forest,
3880  /* refine_recursive */ false,
3881  &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3882  /*init_callback=*/nullptr);
3884  parallel_forest,
3885  /* coarsen_recursive */ false,
3886  &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3887  /*init_callback=*/nullptr);
3888 
3889  // make sure all cells in the lists have been consumed
3890  Assert(refine_and_coarsen_list.pointers_are_at_end(), ExcInternalError());
3891 
3892  // reset the pointer
3893  parallel_forest->user_pointer = this;
3894 
3895  // enforce 2:1 hanging node condition
3897  parallel_forest,
3898  /* face and corner balance */
3899  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3900  P4EST_CONNECT_FULL) :
3901  typename ::internal::p4est::types<dim>::balance_type(
3902  P8EST_CONNECT_FULL)),
3903  /*init_callback=*/nullptr);
3904 
3905  // since refinement and/or coarsening on the parallel forest
3906  // has happened, we need to update the quadrant cell relations
3907  update_quadrant_cell_relations();
3908 
3909  // before repartitioning the mesh, store the current distribution
3910  // of the p4est quadrants and let others attach mesh related info
3911  // (such as SolutionTransfer data)
3912  std::vector<typename ::internal::p4est::types<dim>::gloidx>
3913  previous_global_first_quadrant;
3914 
3915  // pack data only if anything has been attached
3916  if (cell_attached_data.n_attached_data_sets > 0)
3917  {
3918  data_transfer.pack_data(local_quadrant_cell_relations,
3919  cell_attached_data.pack_callbacks_fixed,
3920  cell_attached_data.pack_callbacks_variable);
3921 
3922  // before repartitioning the p4est object, save a copy of the
3923  // positions of the global first quadrants for data transfer later
3924  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3925  std::memcpy(previous_global_first_quadrant.data(),
3926  parallel_forest->global_first_quadrant,
3927  sizeof(
3928  typename ::internal::p4est::types<dim>::gloidx) *
3929  (parallel_forest->mpisize + 1));
3930  }
3931 
3932  if (!(settings & no_automatic_repartitioning))
3933  {
3934  // partition the new mesh between all processors. If cell weights have
3935  // not been given balance the number of cells.
3936  if (this->signals.cell_weight.num_slots() == 0)
3938  parallel_forest,
3939  /* prepare coarsening */ 1,
3940  /* weight_callback */ nullptr);
3941  else
3942  {
3943  // get cell weights for a weighted repartitioning.
3944  const std::vector<unsigned int> cell_weights = get_cell_weights();
3945 
3946  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3947 
3948  // attach (temporarily) a pointer to the cell weights through
3949  // p4est's user_pointer object
3950  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3951  parallel_forest->user_pointer = &partition_weights;
3952 
3954  parallel_forest,
3955  /* prepare coarsening */ 1,
3956  /* weight_callback */
3957  &PartitionWeights<dim, spacedim>::cell_weight);
3958 
3959  // release data
3961  parallel_forest, 0, nullptr, nullptr);
3962  // reset the user pointer to its previous state
3963  parallel_forest->user_pointer = this;
3964  }
3965  }
3966 
3967  // finally copy back from local part of tree to deal.II
3968  // triangulation. before doing so, make sure there are no refine or
3969  // coarsen flags pending
3970  for (const auto &cell : this->active_cell_iterators())
3971  {
3972  cell->clear_refine_flag();
3973  cell->clear_coarsen_flag();
3974  }
3975 
3976  try
3977  {
3978  copy_local_forest_to_triangulation();
3979  }
3980  catch (const typename Triangulation<dim>::DistortedCellList &)
3981  {
3982  // the underlying triangulation should not be checking for distorted
3983  // cells
3984  Assert(false, ExcInternalError());
3985  }
3986 
3987  // transfer data
3988  // only if anything has been attached
3989  if (cell_attached_data.n_attached_data_sets > 0)
3990  {
3991  // execute transfer after triangulation got updated
3992  data_transfer.execute_transfer(parallel_forest,
3993  previous_global_first_quadrant.data());
3994 
3995  // also update the CellStatus information on the new mesh
3996  data_transfer.unpack_cell_status(local_quadrant_cell_relations);
3997  }
3998 
3999 # ifdef DEBUG
4000  // Check that we know the level subdomain ids of all our neighbors. This
4001  // also involves coarser cells that share a vertex if they are active.
4002  //
4003  // Example (M= my, O=other):
4004  // *------*
4005  // | |
4006  // | O |
4007  // | |
4008  // *---*---*------*
4009  // | M | M |
4010  // *---*---*
4011  // | | M |
4012  // *---*---*
4013  // ^- the parent can be owned by somebody else, so O is not a neighbor
4014  // one level coarser
4015  if (settings & construct_multigrid_hierarchy)
4016  {
4017  for (unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
4018  {
4019  std::vector<bool> active_verts =
4020  this->mark_locally_active_vertices_on_level(lvl);
4021 
4022  const unsigned int maybe_coarser_lvl =
4023  (lvl > 0) ? (lvl - 1) : lvl;
4025  cell = this->begin(maybe_coarser_lvl),
4026  endc = this->end(lvl);
4027  for (; cell != endc; ++cell)
4028  if (cell->level() == static_cast<int>(lvl) || cell->is_active())
4029  {
4030  const bool is_level_artificial =
4031  (cell->level_subdomain_id() ==
4033  bool need_to_know = false;
4034  for (const unsigned int vertex :
4036  if (active_verts[cell->vertex_index(vertex)])
4037  {
4038  need_to_know = true;
4039  break;
4040  }
4041 
4042  Assert(
4043  !need_to_know || !is_level_artificial,
4044  ExcMessage(
4045  "Internal error: the owner of cell" +
4046  cell->id().to_string() +
4047  " is unknown even though it is needed for geometric multigrid."));
4048  }
4049  }
4050  }
4051 # endif
4052 
4053  this->update_periodic_face_map();
4054  this->update_number_cache();
4055 
4056  // signal that refinement is finished
4058  }
4059 
4060 
4061 
4062  template <int dim, int spacedim>
4063  void
4065  {
4066 # ifdef DEBUG
4067  for (const auto &cell : this->active_cell_iterators())
4068  if (cell->is_locally_owned())
4069  Assert(
4070  !cell->refine_flag_set() && !cell->coarsen_flag_set(),
4071  ExcMessage(
4072  "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
4073 # endif
4074 
4075  // signal that repartitioning is going to happen
4077 
4078  // before repartitioning the mesh let others attach mesh related info
4079  // (such as SolutionTransfer data) to the p4est
4080  std::vector<typename ::internal::p4est::types<dim>::gloidx>
4081  previous_global_first_quadrant;
4082 
4083  // pack data only if anything has been attached
4084  if (cell_attached_data.n_attached_data_sets > 0)
4085  {
4086  data_transfer.pack_data(local_quadrant_cell_relations,
4087  cell_attached_data.pack_callbacks_fixed,
4088  cell_attached_data.pack_callbacks_variable);
4089 
4090  // before repartitioning the p4est object, save a copy of the
4091  // positions of quadrant for data transfer later
4092  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
4093  std::memcpy(previous_global_first_quadrant.data(),
4094  parallel_forest->global_first_quadrant,
4095  sizeof(
4096  typename ::internal::p4est::types<dim>::gloidx) *
4097  (parallel_forest->mpisize + 1));
4098  }
4099 
4100  if (this->signals.cell_weight.num_slots() == 0)
4101  {
4102  // no cell weights given -- call p4est's 'partition' without a
4103  // callback for cell weights
4105  parallel_forest,
4106  /* prepare coarsening */ 1,
4107  /* weight_callback */ nullptr);
4108  }
4109  else
4110  {
4111  // get cell weights for a weighted repartitioning.
4112  const std::vector<unsigned int> cell_weights = get_cell_weights();
4113 
4114  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
4115 
4116  // attach (temporarily) a pointer to the cell weights through p4est's
4117  // user_pointer object
4118  Assert(parallel_forest->user_pointer == this, ExcInternalError());
4119  parallel_forest->user_pointer = &partition_weights;
4120 
4122  parallel_forest,
4123  /* prepare coarsening */ 1,
4124  /* weight_callback */
4125  &PartitionWeights<dim, spacedim>::cell_weight);
4126 
4127  // reset the user pointer to its previous state
4128  parallel_forest->user_pointer = this;
4129  }
4130 
4131  try
4132  {
4133  copy_local_forest_to_triangulation();
4134  }
4135  catch (const typename Triangulation<dim>::DistortedCellList &)
4136  {
4137  // the underlying triangulation should not be checking for distorted
4138  // cells
4139  Assert(false, ExcInternalError());
4140  }
4141 
4142  // transfer data
4143  // only if anything has been attached
4144  if (cell_attached_data.n_attached_data_sets > 0)
4145  {
4146  // execute transfer after triangulation got updated
4147  data_transfer.execute_transfer(parallel_forest,
4148  previous_global_first_quadrant.data());
4149  }
4150 
4151  this->update_periodic_face_map();
4152 
4153  // update how many cells, edges, etc, we store locally
4154  this->update_number_cache();
4155 
4156  // signal that repartitioning is finished
4158  }
4159 
4160 
4161 
4162  template <int dim, int spacedim>
4163  void
4165  const std::vector<bool> &vertex_locally_moved)
4166  {
4167  Assert(vertex_locally_moved.size() == this->n_vertices(),
4168  ExcDimensionMismatch(vertex_locally_moved.size(),
4169  this->n_vertices()));
4170 # ifdef DEBUG
4171  {
4172  const std::vector<bool> locally_owned_vertices =
4174  for (unsigned int i = 0; i < locally_owned_vertices.size(); ++i)
4175  Assert((vertex_locally_moved[i] == false) ||
4176  (locally_owned_vertices[i] == true),
4177  ExcMessage("The vertex_locally_moved argument must not "
4178  "contain vertices that are not locally owned"));
4179  }
4180 # endif
4181 
4182  // First find out which process should receive which vertices.
4183  // These are specifically the ones that are located on cells at the
4184  // boundary of the subdomain this process owns and the receiving
4185  // process taking periodic faces into account.
4186  // Here, it is sufficient to collect all vertices that are located
4187  // at that boundary.
4188  const std::map<unsigned int, std::set<::types::subdomain_id>>
4191 
4192  // now collect cells and their vertices
4193  // for the interested neighbors
4194  using cellmap_t =
4195  std::map<::types::subdomain_id,
4196  CommunicateLocallyMovedVertices::CellInfo<dim, spacedim>>;
4197  cellmap_t needs_to_get_cells;
4198 
4199  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
4200  this->begin(0);
4201  cell != this->end(0);
4202  ++cell)
4203  {
4204  typename ::internal::p4est::types<dim>::quadrant
4205  p4est_coarse_cell;
4206  ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4207 
4208  CommunicateLocallyMovedVertices::fill_vertices_recursively<dim,
4209  spacedim>(
4210  *this,
4211  this->get_coarse_cell_to_p4est_tree_permutation()[cell->index()],
4212  cell,
4213  p4est_coarse_cell,
4215  vertex_locally_moved,
4216  needs_to_get_cells);
4217  }
4218 
4219  // Send information
4220 
4221  // We need to protect this communication below using a mutex:
4224  mutex, this->get_communicator());
4225 
4226  const int mpi_tag = Utilities::MPI::internal::Tags::
4228 
4229  std::vector<std::vector<char>> sendbuffers(needs_to_get_cells.size());
4230  std::vector<std::vector<char>>::iterator buffer = sendbuffers.begin();
4231  std::vector<MPI_Request> requests(needs_to_get_cells.size());
4232  std::vector<unsigned int> destinations;
4233 
4234  unsigned int idx = 0;
4235 
4236  for (typename cellmap_t::iterator it = needs_to_get_cells.begin();
4237  it != needs_to_get_cells.end();
4238  ++it, ++buffer, ++idx)
4239  {
4240  const unsigned int num_cells = it->second.tree_index.size();
4241  (void)num_cells;
4242  destinations.push_back(it->first);
4243 
4244  Assert(num_cells == it->second.quadrants.size(), ExcInternalError());
4245  Assert(num_cells > 0, ExcInternalError());
4246 
4247  // pack all the data into
4248  // the buffer for this
4249  // recipient and send
4250  // it. keep data around
4251  // till we can make sure
4252  // that the packet has been
4253  // received
4254  it->second.pack_data(*buffer);
4255  const int ierr = MPI_Isend(buffer->data(),
4256  buffer->size(),
4257  MPI_BYTE,
4258  it->first,
4259  mpi_tag,
4260  this->get_communicator(),
4261  &requests[idx]);
4262  AssertThrowMPI(ierr);
4263  }
4264 
4265  Assert(destinations.size() == needs_to_get_cells.size(),
4266  ExcInternalError());
4267 
4268  // collect the neighbors
4269  // that are going to send stuff to us
4270  const unsigned int n_senders =
4272  this->get_communicator(), destinations);
4273 
4274  // receive ghostcelldata
4275  std::vector<char> receive;
4276  CommunicateLocallyMovedVertices::CellInfo<dim, spacedim> cellinfo;
4277  for (unsigned int i = 0; i < n_senders; ++i)
4278  {
4279  MPI_Status status;
4280  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4281  mpi_tag,
4282  this->get_communicator(),
4283  &status);
4284  AssertThrowMPI(ierr);
4285 
4286  int len;
4287  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4288  AssertThrowMPI(ierr);
4289  receive.resize(len);
4290 
4291  char *ptr = receive.data();
4292  ierr = MPI_Recv(ptr,
4293  len,
4294  MPI_BYTE,
4295  status.MPI_SOURCE,
4296  status.MPI_TAG,
4297  this->get_communicator(),
4298  &status);
4299  AssertThrowMPI(ierr);
4300 
4301  cellinfo.unpack_data(receive);
4302  const unsigned int cells = cellinfo.tree_index.size();
4303  for (unsigned int c = 0; c < cells; ++c)
4304  {
4305  typename ::parallel::distributed::
4306  Triangulation<dim, spacedim>::cell_iterator cell(
4307  this,
4308  0,
4309  this->get_p4est_tree_to_coarse_cell_permutation()
4310  [cellinfo.tree_index[c]]);
4311 
4312  typename ::internal::p4est::types<dim>::quadrant
4313  p4est_coarse_cell;
4314  ::internal::p4est::init_coarse_quadrant<dim>(
4315  p4est_coarse_cell);
4316 
4317  CommunicateLocallyMovedVertices::set_vertices_recursively<
4318  dim,
4319  spacedim>(*this,
4320  p4est_coarse_cell,
4321  cell,
4322  cellinfo.quadrants[c],
4323  cellinfo.first_vertices[c],
4324  cellinfo.first_vertex_indices[c]);
4325  }
4326  }
4327 
4328  // complete all sends, so that we can
4329  // safely destroy the buffers.
4330  if (requests.size() > 0)
4331  {
4332  const int ierr =
4333  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4334  AssertThrowMPI(ierr);
4335  }
4336 
4337  // check all msgs got sent and received
4338  Assert(Utilities::MPI::sum(needs_to_get_cells.size(),
4339  this->get_communicator()) ==
4340  Utilities::MPI::sum(n_senders, this->get_communicator()),
4341  ExcInternalError());
4342  }
4343 
4344 
4345 
4346  template <int dim, int spacedim>
4347  unsigned int
4349  const std::function<std::vector<char>(const cell_iterator &,
4350  const CellStatus)> &pack_callback,
4351  const bool returns_variable_size_data)
4352  {
4353  unsigned int handle = numbers::invalid_unsigned_int;
4354 
4355  // Add new callback function to the corresponding register.
4356  // Encode handles according to returns_variable_size_data.
4357  if (returns_variable_size_data)
4358  {
4359  handle = 2 * cell_attached_data.pack_callbacks_variable.size();
4360  cell_attached_data.pack_callbacks_variable.push_back(pack_callback);
4361  }
4362  else
4363  {
4364  handle = 2 * cell_attached_data.pack_callbacks_fixed.size() + 1;
4365  cell_attached_data.pack_callbacks_fixed.push_back(pack_callback);
4366  }
4367 
4368  // Increase overall counter.
4369  ++cell_attached_data.n_attached_data_sets;
4370 
4371  return handle;
4372  }
4373 
4374 
4375 
4376  template <int dim, int spacedim>
4377  void
4379  const unsigned int handle,
4380  const std::function<
4381  void(const cell_iterator &,
4382  const CellStatus,
4383  const boost::iterator_range<std::vector<char>::const_iterator> &)>
4384  &unpack_callback)
4385  {
4386  Assert(cell_attached_data.n_attached_data_sets > 0,
4387  ExcMessage("The notify_ready_to_unpack() has already been called "
4388  "once for each registered callback."));
4389 
4390  // check if local_quadrant_cell_relations have been previously gathered
4391  // correctly
4392  Assert(local_quadrant_cell_relations.size() ==
4393  static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4394  ExcInternalError());
4395 
4396 # ifdef DEBUG
4397  // check validity of handle and deregister pack_callback function.
4398  // first reset with invalid entries to preserve ambiguity of
4399  // handles, then free memory when all were unpacked (see below).
4400  const unsigned int callback_index = handle / 2;
4401  if (handle % 2 == 0)
4402  {
4403  Assert(callback_index <
4404  cell_attached_data.pack_callbacks_variable.size(),
4405  ExcMessage("Invalid handle."));
4406 
4407  Assert(cell_attached_data.pack_callbacks_variable[callback_index] !=
4408  nullptr,
4409  ExcInternalError());
4410  cell_attached_data.pack_callbacks_variable[callback_index] = nullptr;
4411  }
4412  else
4413  {
4414  Assert(callback_index <
4415  cell_attached_data.pack_callbacks_fixed.size(),
4416  ExcMessage("Invalid handle."));
4417 
4418  Assert(cell_attached_data.pack_callbacks_fixed[callback_index] !=
4419  nullptr,
4420  ExcInternalError());
4421  cell_attached_data.pack_callbacks_fixed[callback_index] = nullptr;
4422  }
4423 # endif
4424 
4425  // perform unpacking
4426  data_transfer.unpack_data(local_quadrant_cell_relations,
4427  handle,
4428  unpack_callback);
4429 
4430  // decrease counters
4431  --cell_attached_data.n_attached_data_sets;
4432  if (cell_attached_data.n_attached_deserialize > 0)
4433  --cell_attached_data.n_attached_deserialize;
4434 
4435  // important: only remove data if we are not in the deserialization
4436  // process. There, each SolutionTransfer registers and unpacks before
4437  // the next one does this, so n_attached_data_sets is only 1 here. This
4438  // would destroy the saved data before the second SolutionTransfer can
4439  // get it. This created a bug that is documented in
4440  // tests/mpi/p4est_save_03 with more than one SolutionTransfer.
4441  if (cell_attached_data.n_attached_data_sets == 0 &&
4442  cell_attached_data.n_attached_deserialize == 0)
4443  {
4444  // everybody got their data, time for cleanup!
4445  cell_attached_data.pack_callbacks_fixed.clear();
4446  cell_attached_data.pack_callbacks_variable.clear();
4447  data_transfer.clear();
4448 
4449  // reset all cell_status entries after coarsening/refinement
4450  for (auto &quad_cell_rel : local_quadrant_cell_relations)
4451  std::get<1>(quad_cell_rel) =
4453  }
4454  }
4455 
4456 
4457 
4458  template <int dim, int spacedim>
4459  const std::vector<types::global_dof_index> &
4461  const
4462  {
4463  return p4est_tree_to_coarse_cell_permutation;
4464  }
4465 
4466 
4467 
4468  template <int dim, int spacedim>
4469  const std::vector<types::global_dof_index> &
4471  const
4472  {
4473  return coarse_cell_to_p4est_tree_permutation;
4474  }
4475 
4476 
4477 
4478  template <int dim, int spacedim>
4479  std::vector<bool>
4481  const int level) const
4482  {
4483  Assert(dim > 1, ExcNotImplemented());
4484 
4485  std::vector<bool> marked_vertices(this->n_vertices(), false);
4486  cell_iterator cell = this->begin(level), endc = this->end(level);
4487  for (; cell != endc; ++cell)
4488  if (cell->level_subdomain_id() == this->locally_owned_subdomain())
4489  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
4490  marked_vertices[cell->vertex_index(v)] = true;
4491 
4497  typename std::map<std::pair<cell_iterator, unsigned int>,
4498  std::pair<std::pair<cell_iterator, unsigned int>,
4499  std::bitset<3>>>::const_iterator it;
4500 
4501  // When a connectivity in the code below is detected, the assignment
4502  // 'marked_vertices[v1] = marked_vertices[v2] = true' makes sure that
4503  // the information about the periodicity propagates back to vertices on
4504  // cells that are not owned locally. However, in the worst case we want
4505  // to connect to a vertex that is 'dim' hops away from the locally owned
4506  // cell. Depending on the order of the periodic face map, we might
4507  // connect to that point by chance or miss it. However, after looping
4508  // through all the periodic directions (which are at most as many as
4509  // the number of space dimensions) we can be sure that all connections
4510  // to vertices have been created.
4511  for (unsigned int repetition = 0; repetition < dim; ++repetition)
4512  for (it = this->get_periodic_face_map().begin();
4513  it != this->get_periodic_face_map().end();
4514  ++it)
4515  {
4516  const cell_iterator & cell_1 = it->first.first;
4517  const unsigned int face_no_1 = it->first.second;
4518  const cell_iterator & cell_2 = it->second.first.first;
4519  const unsigned int face_no_2 = it->second.first.second;
4520  const std::bitset<3> &face_orientation = it->second.second;
4521 
4522  if (cell_1->level() == level && cell_2->level() == level)
4523  {
4524  for (unsigned int v = 0;
4525  v < GeometryInfo<dim - 1>::vertices_per_cell;
4526  ++v)
4527  {
4528  // take possible non-standard orientation of faces into
4529  // account
4530  const unsigned int vface0 =
4532  v,
4533  face_orientation[0],
4534  face_orientation[1],
4535  face_orientation[2]);
4536  if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
4537  vface0)] ||
4538  marked_vertices[cell_2->face(face_no_2)->vertex_index(
4539  v)])
4540  marked_vertices[cell_1->face(face_no_1)->vertex_index(
4541  vface0)] =
4542  marked_vertices[cell_2->face(face_no_2)->vertex_index(
4543  v)] = true;
4544  }
4545  }
4546  }
4547 
4548  return marked_vertices;
4549  }
4550 
4551 
4552 
4553  template <int dim, int spacedim>
4554  unsigned int
4557  {
4558  return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
4559  }
4560 
4561 
4562 
4563  template <int dim, int spacedim>
4566  const unsigned int coarse_cell_index) const
4567  {
4568  return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
4569  }
4570 
4571 
4572 
4573  template <int dim, int spacedim>
4574  void
4576  const std::vector<::GridTools::PeriodicFacePair<cell_iterator>>
4577  &periodicity_vector)
4578  {
4579  Assert(triangulation_has_content == true,
4580  ExcMessage("The triangulation is empty!"));
4581  Assert(this->n_levels() == 1,
4582  ExcMessage("The triangulation is refined!"));
4583 
4584  // call the base class for storing the periodicity information; we must
4585  // do this before going to p4est and rebuilding the triangulation to get
4586  // the level subdomain ids correct in the multigrid case
4588 
4589  for (const auto &face_pair : periodicity_vector)
4590  {
4591  const cell_iterator first_cell = face_pair.cell[0];
4592  const cell_iterator second_cell = face_pair.cell[1];
4593  const unsigned int face_left = face_pair.face_idx[0];
4594  const unsigned int face_right = face_pair.face_idx[1];
4595 
4596  // respective cells of the matching faces in p4est
4597  const unsigned int tree_left =
4598  coarse_cell_to_p4est_tree_permutation[first_cell->index()];
4599  const unsigned int tree_right =
4600  coarse_cell_to_p4est_tree_permutation[second_cell->index()];
4601 
4602  // p4est wants to know which corner the first corner on
4603  // the face with the lower id is mapped to on the face with
4604  // with the higher id. For d==2 there are only two possibilities
4605  // that are determined by it->orientation[1].
4606  // For d==3 we have to use GridTools::OrientationLookupTable.
4607  // The result is given below.
4608 
4609  unsigned int p4est_orientation = 0;
4610  if (dim == 2)
4611  p4est_orientation = face_pair.orientation[1];
4612  else
4613  {
4614  const unsigned int face_idx_list[] = {face_left, face_right};
4615  const cell_iterator cell_list[] = {first_cell, second_cell};
4616  unsigned int lower_idx, higher_idx;
4617  if (face_left <= face_right)
4618  {
4619  higher_idx = 1;
4620  lower_idx = 0;
4621  }
4622  else
4623  {
4624  higher_idx = 0;
4625  lower_idx = 1;
4626  }
4627 
4628  // get the cell index of the first index on the face with the
4629  // lower id
4630  unsigned int first_p4est_idx_on_cell =
4631  p8est_face_corners[face_idx_list[lower_idx]][0];
4632  unsigned int first_dealii_idx_on_face =
4634  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
4635  ++i)
4636  {
4637  const unsigned int first_dealii_idx_on_cell =
4639  face_idx_list[lower_idx],
4640  i,
4641  cell_list[lower_idx]->face_orientation(
4642  face_idx_list[lower_idx]),
4643  cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
4644  cell_list[lower_idx]->face_rotation(
4645  face_idx_list[lower_idx]));
4646  if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
4647  {
4648  first_dealii_idx_on_face = i;
4649  break;
4650  }
4651  }
4652  Assert(first_dealii_idx_on_face != numbers::invalid_unsigned_int,
4653  ExcInternalError());
4654  // Now map dealii_idx_on_face according to the orientation
4655  constexpr unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
4656  {0, 1, 2, 3},
4657  {3, 1, 2, 0},
4658  {3, 2, 1, 0},
4659  {2, 3, 0, 1},
4660  {1, 3, 0, 2},
4661  {1, 0, 3, 2},
4662  {2, 0, 3, 1}};
4663  constexpr unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
4664  {0, 1, 2, 3},
4665  {3, 1, 2, 0},
4666  {3, 2, 1, 0},
4667  {2, 3, 0, 1},
4668  {2, 0, 3, 1},
4669  {1, 0, 3, 2},
4670  {1, 3, 0, 2}};
4671  const unsigned int second_dealii_idx_on_face =
4672  lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
4673  [first_dealii_idx_on_face] :
4674  right_to_left[face_pair.orientation.to_ulong()]
4675  [first_dealii_idx_on_face];
4676  const unsigned int second_dealii_idx_on_cell =
4678  face_idx_list[higher_idx],
4679  second_dealii_idx_on_face,
4680  cell_list[higher_idx]->face_orientation(
4681  face_idx_list[higher_idx]),
4682  cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
4683  cell_list[higher_idx]->face_rotation(
4684  face_idx_list[higher_idx]));
4685  // map back to p4est
4686  const unsigned int second_p4est_idx_on_face =
4687  p8est_corner_face_corners[second_dealii_idx_on_cell]
4688  [face_idx_list[higher_idx]];
4689  p4est_orientation = second_p4est_idx_on_face;
4690  }
4691 
4693  connectivity,
4694  tree_left,
4695  tree_right,
4696  face_left,
4697  face_right,
4698  p4est_orientation);
4699  }
4700 
4701 
4703  connectivity) == 1,
4704  ExcInternalError());
4705 
4706  // now create a forest out of the connectivity data structure
4709  this->mpi_communicator,
4710  connectivity,
4711  /* minimum initial number of quadrants per tree */ 0,
4712  /* minimum level of upfront refinement */ 0,
4713  /* use uniform upfront refinement */ 1,
4714  /* user_data_size = */ 0,
4715  /* user_data_constructor = */ nullptr,
4716  /* user_pointer */ this);
4717 
4718  try
4719  {
4720  copy_local_forest_to_triangulation();
4721  }
4722  catch (const typename Triangulation<dim>::DistortedCellList &)
4723  {
4724  // the underlying triangulation should not be checking for distorted
4725  // cells
4726  Assert(false, ExcInternalError());
4727  }
4728 
4729  // The range of ghost_owners might have changed so update that information
4730  this->update_number_cache();
4731  }
4732 
4733 
4734 
4735  template <int dim, int spacedim>
4736  std::size_t
4738  {
4739  std::size_t mem =
4742  MemoryConsumption::memory_consumption(triangulation_has_content) +
4744  MemoryConsumption::memory_consumption(parallel_forest) +
4746  cell_attached_data.n_attached_data_sets) +
4747  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_fixed)
4748  // +
4749  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_variable)
4750  // +
4751  // TODO[TH]: how?
4753  coarse_cell_to_p4est_tree_permutation) +
4755  p4est_tree_to_coarse_cell_permutation) +
4756  memory_consumption_p4est();
4757 
4758  return mem;
4759  }
4760 
4761 
4762 
4763  template <int dim, int spacedim>
4764  std::size_t
4766  {
4767  return ::internal::p4est::functions<dim>::forest_memory_used(
4768  parallel_forest) +
4770  connectivity);
4771  }
4772 
4773 
4774 
4775  template <int dim, int spacedim>
4776  void
4778  const ::Triangulation<dim, spacedim> &other_tria)
4779  {
4780  try
4781  {
4783  copy_triangulation(other_tria);
4784  }
4785  catch (
4786  const typename ::Triangulation<dim, spacedim>::DistortedCellList
4787  &)
4788  {
4789  // the underlying triangulation should not be checking for distorted
4790  // cells
4791  Assert(false, ExcInternalError());
4792  }
4793 
4794  // note that now we have some content in the p4est objects and call the
4795  // functions that do the actual work (which are dimension dependent, so
4796  // separate)
4797  triangulation_has_content = true;
4798 
4799  Assert(other_tria.n_levels() == 1,
4800  ExcMessage(
4801  "Parallel distributed triangulations can only be copied, "
4802  "if they are not refined!"));
4803 
4804  if (const ::parallel::distributed::Triangulation<dim, spacedim>
4805  *other_tria_x =
4806  dynamic_cast<const ::parallel::distributed::
4807  Triangulation<dim, spacedim> *>(&other_tria))
4808  {
4809  coarse_cell_to_p4est_tree_permutation =
4810  other_tria_x->coarse_cell_to_p4est_tree_permutation;
4811  p4est_tree_to_coarse_cell_permutation =
4812  other_tria_x->p4est_tree_to_coarse_cell_permutation;
4813  cell_attached_data = other_tria_x->cell_attached_data;
4814  data_transfer = other_tria_x->data_transfer;
4815 
4816  settings = other_tria_x->settings;
4817  }
4818  else
4819  {
4820  setup_coarse_cell_to_p4est_tree_permutation();
4821  }
4822 
4823  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
4824 
4825  try
4826  {
4827  copy_local_forest_to_triangulation();
4828  }
4829  catch (const typename Triangulation<dim>::DistortedCellList &)
4830  {
4831  // the underlying triangulation should not be checking for distorted
4832  // cells
4833  Assert(false, ExcInternalError());
4834  }
4835 
4836  this->update_periodic_face_map();
4837  this->update_number_cache();
4838  }
4839 
4840 
4841 
4842  template <int dim, int spacedim>
4843  void
4845  {
4846  // reorganize memory for local_quadrant_cell_relations
4847  local_quadrant_cell_relations.resize(
4848  parallel_forest->local_num_quadrants);
4849  local_quadrant_cell_relations.shrink_to_fit();
4850 
4851  // recurse over p4est
4852  for (typename Triangulation<dim, spacedim>::cell_iterator cell =
4853  this->begin(0);
4854  cell != this->end(0);
4855  ++cell)
4856  {
4857  // skip coarse cells that are not ours
4858  if (tree_exists_locally<dim, spacedim>(
4859  parallel_forest,
4860  coarse_cell_to_p4est_tree_permutation[cell->index()]) == false)
4861  continue;
4862 
4863  // initialize auxiliary top level p4est quadrant
4864  typename ::internal::p4est::types<dim>::quadrant
4865  p4est_coarse_cell;
4866  ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4867 
4868  // determine tree to start recursion on
4869  typename ::internal::p4est::types<dim>::tree *tree =
4870  init_tree(cell->index());
4871 
4872  update_quadrant_cell_relations_recursively<dim, spacedim>(
4873  local_quadrant_cell_relations, *tree, cell, p4est_coarse_cell);
4874  }
4875  }
4876 
4877 
4878 
4879  template <int dim, int spacedim>
4880  std::vector<unsigned int>
4882  {
4883  // check if local_quadrant_cell_relations have been previously gathered
4884  // correctly
4885  Assert(local_quadrant_cell_relations.size() ==
4886  static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4887  ExcInternalError());
4888 
4889  // Allocate the space for the weights. In fact we do not know yet, how
4890  // many cells we own after the refinement (only p4est knows that
4891  // at this point). We simply reserve n_active_cells space and if many
4892  // more cells are refined than coarsened than additional reallocation
4893  // will be done inside get_cell_weights_recursively.
4894  std::vector<unsigned int> weights;
4895  weights.reserve(this->n_active_cells());
4896 
4897  // Iterate over p4est and Triangulation relations
4898  // to find refined/coarsened/kept
4899  // cells. Then append cell_weight.
4900  // Note that we need to follow the p4est ordering
4901  // instead of the deal.II ordering to get the cell_weights
4902  // in the same order p4est will encounter them during repartitioning.
4903  for (const auto &quad_cell_rel : local_quadrant_cell_relations)
4904  {
4905  const auto &cell_status = std::get<1>(quad_cell_rel);
4906  const auto &cell_it = std::get<2>(quad_cell_rel);
4907 
4908  switch (cell_status)
4909  {
4911  spacedim>::CELL_PERSIST:
4912  weights.push_back(1000);
4913  weights.back() += this->signals.cell_weight(
4914  cell_it,
4916  spacedim>::CELL_PERSIST);
4917  break;
4918 
4920  spacedim>::CELL_REFINE:
4922  spacedim>::CELL_INVALID:
4923  {
4924  // calculate weight of parent cell
4925  unsigned int parent_weight = 1000;
4926  parent_weight += this->signals.cell_weight(
4927  cell_it,
4929  CELL_REFINE);
4930  // assign the weight of the parent cell equally to all
4931  // children
4932  weights.push_back(parent_weight);
4933  break;
4934  }
4935 
4937  spacedim>::CELL_COARSEN:
4938  weights.push_back(1000);
4939  weights.back() += this->signals.cell_weight(
4940  cell_it,
4942  spacedim>::CELL_COARSEN);
4943  break;
4944 
4945  default:
4946  Assert(false, ExcInternalError());
4947  break;
4948  }
4949  }
4950 
4951  return weights;
4952  }
4953 
4954 
4955 
4956  template <int spacedim>
4958  MPI_Comm mpi_communicator,
4959  const typename ::Triangulation<1, spacedim>::MeshSmoothing
4960  smooth_grid,
4961  const Settings /*settings*/)
4962  : ::parallel::DistributedTriangulationBase<1, spacedim>(
4963  mpi_communicator,
4964  smooth_grid,
4965  false)
4966  {
4967  Assert(false, ExcNotImplemented());
4968  }
4969 
4970 
4971  template <int spacedim>
4973  {
4974  AssertNothrow(false, ExcNotImplemented());
4975  }
4976 
4977 
4978 
4979  template <int spacedim>
4980  void
4982  const std::vector<bool> & /*vertex_locally_moved*/)
4983  {
4984  Assert(false, ExcNotImplemented());
4985  }
4986 
4987 
4988 
4989  template <int spacedim>
4990  unsigned int
4992  const std::function<std::vector<char>(
4993  const typename ::Triangulation<1, spacedim>::cell_iterator &,
4994  const typename ::Triangulation<1, spacedim>::CellStatus)>
4995  & /*pack_callback*/,
4996  const bool /*returns_variable_size_data*/)
4997  {
4998  Assert(false, ExcNotImplemented());
4999  return 0;
5000  }
5001 
5002 
5003 
5004  template <int spacedim>
5005  void
5007  const unsigned int /*handle*/,
5008  const std::function<
5009  void(const typename ::Triangulation<1, spacedim>::cell_iterator &,
5010  const typename ::Triangulation<1, spacedim>::CellStatus,
5011  const boost::iterator_range<std::vector<char>::const_iterator> &)>
5012  & /*unpack_callback*/)
5013  {
5014  Assert(false, ExcNotImplemented());
5015  }
5016 
5017 
5018 
5019  template <int spacedim>
5020  const std::vector<types::global_dof_index> &
5022  const
5023  {
5024  static std::vector<types::global_dof_index> a;
5025  return a;
5026  }
5027 
5028 
5029 
5030  template <int spacedim>
5031  std::map<unsigned int, std::set<::types::subdomain_id>>
5033  const unsigned int /*level*/) const
5034  {
5035  Assert(false, ExcNotImplemented());
5036 
5037  return std::map<unsigned int, std::set<::types::subdomain_id>>();
5038  }
5039 
5040 
5041 
5042  template <int spacedim>
5043  std::vector<bool>
5045  const unsigned int) const
5046  {
5047  Assert(false, ExcNotImplemented());
5048  return std::vector<bool>();
5049  }
5050 
5051 
5052 
5053  template <int spacedim>
5054  unsigned int
5056  const types::coarse_cell_id) const
5057  {
5058  Assert(false, ExcNotImplemented());
5059  return 0;
5060  }
5061 
5062 
5063 
5064  template <int spacedim>
5067  const unsigned int) const
5068  {
5069  Assert(false, ExcNotImplemented());
5070  return 0;
5071  }
5072 
5073 
5074  template <int spacedim>
5075  void
5076  Triangulation<1, spacedim>::load(const std::string &, const bool)
5077  {
5078  Assert(false, ExcNotImplemented());
5079  }
5080 
5081 
5082 
5083  template <int spacedim>
5084  void
5085  Triangulation<1, spacedim>::save(const std::string &) const
5086  {
5087  Assert(false, ExcNotImplemented());
5088  }
5089 
5090 
5091 
5092  template <int spacedim>
5093  bool
5095  {
5096  Assert(false, ExcNotImplemented());
5097  return false;
5098  }
5099 
5100  } // namespace distributed
5101 } // namespace parallel
5102 
5103 
5104 #endif // DEAL_II_WITH_P4EST
5105 
5106 
5107 
5108 /*-------------- Explicit Instantiations -------------------------------*/
5109 #include "tria.inst"
5110 
5111 
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:12675
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:10438
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:4348
static const unsigned int invalid_unsigned_int
Definition: types.h:191
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1483
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:2186
virtual bool has_hanging_nodes() const
Definition: tria.cc:12807
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:12667
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:3848
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12185
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:2157
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:12013
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
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:11993
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
Definition: tria.cc:4460
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:2524
cell_iterator end() const
Definition: tria.cc:12079
virtual unsigned int coarse_cell_id_to_coarse_cell_index(const types::coarse_cell_id coarse_cell_id) const override
Definition: tria.cc:4555
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator >> &)
Definition: tria.cc:13356
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:13385
virtual std::size_t memory_consumption() const override
Definition: tria_base.cc:88
unsigned int n_active_lines() const
Definition: tria.cc:12910
virtual bool prepare_coarsening_and_refinement()
Definition: tria.cc:14083
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:10106
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const override
Definition: tria.cc:4565
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:288
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
Definition: tria.cc:10521
#define Assert(cond, exc)
Definition: exceptions.h:1419
Signals signals
Definition: tria.h:2222
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:10953
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
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:4378
unsigned int level
Definition: grid_out.cc:4355
void save_refine_flags(std::ostream &out) const
Definition: tria.cc:10885
boost::signals2::signal< void()> pre_distributed_save
Definition: tria.h:2194
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1446
virtual std::size_t memory_consumption() const
Definition: tria.cc:15157
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:1209
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:302
boost::signals2::signal< void()> post_distributed_save
Definition: tria.h:2201
void communicate_locally_moved_vertices(const std::vector< bool > &vertex_locally_moved)
Definition: tria.cc:4164
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:1707
active_cell_iterator end_active(const unsigned int level) const
Definition: tria.cc:12148
void update_periodic_face_map()
Definition: tria.cc:13444
const std::vector< bool > & get_used_vertices() const
Definition: tria.cc:13262
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
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:3419
virtual ~Triangulation() override
Definition: tria.cc:10180
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:2144
boost::signals2::signal< void()> pre_distributed_load
Definition: tria.h:2209
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:1353
std::vector< unsigned int > sizes_fixed_cumulative
Definition: tria.h:1114
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
static std::mutex mutex
Definition: subscriptor.h:263
static ::ExceptionBase & ExcNotImplemented()
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2978
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
active_cell_iterator last_active() const
Definition: tria.cc:12058
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:2216
global_cell_index coarse_cell_id
Definition: types.h:114
boost::signals2::signal< void()> pre_distributed_repartition
Definition: tria.h:2179
CellAttachedData cell_attached_data
Definition: tria.h:924
std::vector< bool > mark_locally_active_vertices_on_level(const int level) const
Definition: tria.cc:4480
T max(const T &t, const MPI_Comm &mpi_communicator)
virtual types::subdomain_id locally_owned_subdomain() const
Definition: tria.cc:13329
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:10209
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:13376
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:5517
boost::signals2::signal< void()> post_distributed_refinement
Definition: tria.h:2167