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