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