Reference documentation for deal.II version GIT d77e5ebb0a 2023-01-27 22:35:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tria.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/logstream.h>
19 #include <deal.II/base/point.h>
20 #include <deal.II/base/utilities.h>
21 
24 
26 #include <deal.II/grid/tria.h>
29 
30 #include <algorithm>
31 #include <fstream>
32 #include <iostream>
33 #include <limits>
34 #include <numeric>
35 
36 
38 
39 
40 #ifdef DEAL_II_WITH_P4EST
41 
42 namespace
43 {
44  template <int dim, int spacedim>
45  void
46  get_vertex_to_cell_mappings(
48  std::vector<unsigned int> & vertex_touch_count,
49  std::vector<std::list<
51  unsigned int>>> & vertex_to_cell)
52  {
53  vertex_touch_count.resize(triangulation.n_vertices());
54  vertex_to_cell.resize(triangulation.n_vertices());
55 
56  for (const auto &cell : triangulation.active_cell_iterators())
57  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
58  {
59  ++vertex_touch_count[cell->vertex_index(v)];
60  vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
61  }
62  }
63 
64 
65 
66  template <int dim, int spacedim>
67  void
68  get_edge_to_cell_mappings(
70  std::vector<unsigned int> & edge_touch_count,
71  std::vector<std::list<
73  unsigned int>>> & edge_to_cell)
74  {
75  Assert(triangulation.n_levels() == 1, ExcInternalError());
76 
77  edge_touch_count.resize(triangulation.n_active_lines());
78  edge_to_cell.resize(triangulation.n_active_lines());
79 
80  for (const auto &cell : triangulation.active_cell_iterators())
81  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
82  {
83  ++edge_touch_count[cell->line(l)->index()];
84  edge_to_cell[cell->line(l)->index()].emplace_back(cell, l);
85  }
86  }
87 
88 
89 
94  template <int dim, int spacedim>
95  void
96  set_vertex_and_cell_info(
98  const std::vector<unsigned int> & vertex_touch_count,
99  const std::vector<std::list<
101  unsigned int>>> & vertex_to_cell,
102  const std::vector<types::global_dof_index>
103  & coarse_cell_to_p4est_tree_permutation,
104  const bool set_vertex_info,
105  typename internal::p4est::types<dim>::connectivity *connectivity)
106  {
107  // copy the vertices into the connectivity structure. the triangulation
108  // exports the array of vertices, but some of the entries are sometimes
109  // unused; this shouldn't be the case for a newly created triangulation,
110  // but make sure
111  //
112  // note that p4est stores coordinates as a triplet of values even in 2d
113  Assert(triangulation.get_used_vertices().size() ==
114  triangulation.get_vertices().size(),
115  ExcInternalError());
116  Assert(std::find(triangulation.get_used_vertices().begin(),
117  triangulation.get_used_vertices().end(),
118  false) == triangulation.get_used_vertices().end(),
119  ExcInternalError());
120  if (set_vertex_info == true)
121  for (unsigned int v = 0; v < triangulation.n_vertices(); ++v)
122  {
123  connectivity->vertices[3 * v] = triangulation.get_vertices()[v][0];
124  connectivity->vertices[3 * v + 1] =
125  triangulation.get_vertices()[v][1];
126  connectivity->vertices[3 * v + 2] =
127  (spacedim == 2 ? 0 : triangulation.get_vertices()[v][2]);
128  }
129 
130  // next store the tree_to_vertex indices (each tree is here only a single
131  // cell in the coarse mesh). p4est requires vertex numbering in clockwise
132  // orientation
133  //
134  // while we're at it, also copy the neighborship information between cells
136  cell = triangulation.begin_active(),
137  endc = triangulation.end();
138  for (; cell != endc; ++cell)
139  {
140  const unsigned int index =
141  coarse_cell_to_p4est_tree_permutation[cell->index()];
142 
143  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
144  {
145  if (set_vertex_info == true)
146  connectivity
147  ->tree_to_vertex[index * GeometryInfo<dim>::vertices_per_cell +
148  v] = cell->vertex_index(v);
149  connectivity
150  ->tree_to_corner[index * GeometryInfo<dim>::vertices_per_cell +
151  v] = cell->vertex_index(v);
152  }
153 
154  // neighborship information. if a cell is at a boundary, then enter
155  // the index of the cell itself here
156  for (auto f : GeometryInfo<dim>::face_indices())
157  if (cell->face(f)->at_boundary() == false)
158  connectivity
159  ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
160  coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->index()];
161  else
162  connectivity
163  ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
164  coarse_cell_to_p4est_tree_permutation[cell->index()];
165 
166  // fill tree_to_face, which is essentially neighbor_to_neighbor;
167  // however, we have to remap the resulting face number as well
168  for (auto f : GeometryInfo<dim>::face_indices())
169  if (cell->face(f)->at_boundary() == false)
170  {
171  switch (dim)
172  {
173  case 2:
174  {
175  connectivity->tree_to_face
177  cell->neighbor_of_neighbor(f);
178  break;
179  }
180 
181  case 3:
182  {
183  /*
184  * The values for tree_to_face are in 0..23 where ttf % 6
185  * gives the face number and ttf / 4 the face orientation
186  * code. The orientation is determined as follows. Let
187  * my_face and other_face be the two face numbers of the
188  * connecting trees in 0..5. Then the first face vertex
189  * of the lower of my_face and other_face connects to a
190  * face vertex numbered 0..3 in the higher of my_face and
191  * other_face. The face orientation is defined as this
192  * number. If my_face == other_face, treating either of
193  * both faces as the lower one leads to the same result.
194  */
195 
196  connectivity->tree_to_face[index * 6 + f] =
197  cell->neighbor_of_neighbor(f);
198 
199  unsigned int face_idx_list[2] = {
200  f, cell->neighbor_of_neighbor(f)};
202  cell_list[2] = {cell, cell->neighbor(f)};
203  unsigned int smaller_idx = 0;
204 
205  if (f > cell->neighbor_of_neighbor(f))
206  smaller_idx = 1;
207 
208  unsigned int larger_idx = (smaller_idx + 1) % 2;
209  // smaller = *_list[smaller_idx]
210  // larger = *_list[larger_idx]
211 
212  unsigned int v = 0;
213 
214  // global vertex index of vertex 0 on face of cell with
215  // smaller local face index
216  unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
218  face_idx_list[smaller_idx],
219  0,
220  cell_list[smaller_idx]->face_orientation(
221  face_idx_list[smaller_idx]),
222  cell_list[smaller_idx]->face_flip(
223  face_idx_list[smaller_idx]),
224  cell_list[smaller_idx]->face_rotation(
225  face_idx_list[smaller_idx])));
226 
227  // loop over vertices on face from other cell and compare
228  // global vertex numbers
229  for (unsigned int i = 0;
230  i < GeometryInfo<dim>::vertices_per_face;
231  ++i)
232  {
233  unsigned int idx =
234  cell_list[larger_idx]->vertex_index(
236  face_idx_list[larger_idx], i));
237 
238  if (idx == g_idx)
239  {
240  v = i;
241  break;
242  }
243  }
244 
245  connectivity->tree_to_face[index * 6 + f] += 6 * v;
246  break;
247  }
248 
249  default:
250  Assert(false, ExcNotImplemented());
251  }
252  }
253  else
254  connectivity
255  ->tree_to_face[index * GeometryInfo<dim>::faces_per_cell + f] = f;
256  }
257 
258  // now fill the vertex information
259  connectivity->ctt_offset[0] = 0;
260  std::partial_sum(vertex_touch_count.begin(),
261  vertex_touch_count.end(),
262  &connectivity->ctt_offset[1]);
263 
264  const typename internal::p4est::types<dim>::locidx num_vtt =
265  std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
266  (void)num_vtt;
267  Assert(connectivity->ctt_offset[triangulation.n_vertices()] == num_vtt,
268  ExcInternalError());
269 
270  for (unsigned int v = 0; v < triangulation.n_vertices(); ++v)
271  {
272  Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
273  ExcInternalError());
274 
275  typename std::list<
276  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
277  unsigned int>>::const_iterator p =
278  vertex_to_cell[v].begin();
279  for (unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
280  {
281  connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
282  coarse_cell_to_p4est_tree_permutation[p->first->index()];
283  connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
284  p->second;
285  }
286  }
287  }
288 
289 
290 
291  template <int dim, int spacedim>
292  bool
294  const typename internal::p4est::types<dim>::forest *parallel_forest,
295  const typename internal::p4est::types<dim>::topidx coarse_grid_cell)
296  {
297  Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
298  ExcInternalError());
299  return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
300  (coarse_grid_cell <= parallel_forest->last_local_tree));
301  }
302 
303 
304  template <int dim, int spacedim>
305  void
306  delete_all_children_and_self(
307  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
308  {
309  if (cell->has_children())
310  for (unsigned int c = 0; c < cell->n_children(); ++c)
311  delete_all_children_and_self<dim, spacedim>(cell->child(c));
312  else
313  cell->set_coarsen_flag();
314  }
315 
316 
317 
318  template <int dim, int spacedim>
319  void
320  delete_all_children(
321  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
322  {
323  if (cell->has_children())
324  for (unsigned int c = 0; c < cell->n_children(); ++c)
325  delete_all_children_and_self<dim, spacedim>(cell->child(c));
326  }
327 
328 
329  template <int dim, int spacedim>
330  void
331  determine_level_subdomain_id_recursively(
332  const typename internal::p4est::types<dim>::tree & tree,
333  const typename internal::p4est::types<dim>::locidx & tree_index,
334  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
335  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
336  typename internal::p4est::types<dim>::forest & forest,
337  const types::subdomain_id my_subdomain,
338  const std::vector<std::vector<bool>> & marked_vertices)
339  {
340  if (dealii_cell->level_subdomain_id() == numbers::artificial_subdomain_id)
341  {
342  // important: only assign the level_subdomain_id if it is a ghost cell
343  // even though we could fill in all.
344  bool used = false;
345  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
346  {
347  if (marked_vertices[dealii_cell->level()]
348  [dealii_cell->vertex_index(v)])
349  {
350  used = true;
351  break;
352  }
353  }
354 
355  // Special case: if this cell is active we might be a ghost neighbor
356  // to a locally owned cell across a vertex that is finer.
357  // Example (M= my, O=dealii_cell, owned by somebody else):
358  // *------*
359  // | |
360  // | O |
361  // | |
362  // *---*---*------*
363  // | M | M |
364  // *---*---*
365  // | | M |
366  // *---*---*
367  if (!used && dealii_cell->is_active() &&
368  dealii_cell->is_artificial() == false &&
369  dealii_cell->level() + 1 < static_cast<int>(marked_vertices.size()))
370  {
371  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
372  {
373  if (marked_vertices[dealii_cell->level() + 1]
374  [dealii_cell->vertex_index(v)])
375  {
376  used = true;
377  break;
378  }
379  }
380  }
381 
382  // Like above, but now the other way around
383  if (!used && dealii_cell->is_active() &&
384  dealii_cell->is_artificial() == false && dealii_cell->level() > 0)
385  {
386  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
387  {
388  if (marked_vertices[dealii_cell->level() - 1]
389  [dealii_cell->vertex_index(v)])
390  {
391  used = true;
392  break;
393  }
394  }
395  }
396 
397  if (used)
398  {
400  &forest, tree_index, &p4est_cell, my_subdomain);
401  Assert((owner != -2) && (owner != -1),
402  ExcMessage("p4est should know the owner."));
403  dealii_cell->set_level_subdomain_id(owner);
404  }
405  }
406 
407  if (dealii_cell->has_children())
408  {
411  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
412  ++c)
413  switch (dim)
414  {
415  case 2:
416  P4EST_QUADRANT_INIT(&p4est_child[c]);
417  break;
418  case 3:
419  P8EST_QUADRANT_INIT(&p4est_child[c]);
420  break;
421  default:
422  Assert(false, ExcNotImplemented());
423  }
424 
425 
427  p4est_child);
428 
429  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
430  ++c)
431  {
432  determine_level_subdomain_id_recursively<dim, spacedim>(
433  tree,
434  tree_index,
435  dealii_cell->child(c),
436  p4est_child[c],
437  forest,
438  my_subdomain,
439  marked_vertices);
440  }
441  }
442  }
443 
444 
445  template <int dim, int spacedim>
446  void
447  match_tree_recursively(
448  const typename internal::p4est::types<dim>::tree & tree,
449  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
450  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
451  const typename internal::p4est::types<dim>::forest & forest,
452  const types::subdomain_id my_subdomain)
453  {
454  // check if this cell exists in the local p4est cell
455  if (sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
456  &p4est_cell,
458  -1)
459  {
460  // yes, cell found in local part of p4est
461  delete_all_children<dim, spacedim>(dealii_cell);
462  if (dealii_cell->is_active())
463  dealii_cell->set_subdomain_id(my_subdomain);
464  }
465  else
466  {
467  // no, cell not found in local part of p4est. this means that the
468  // local part is more refined than the current cell. if this cell has
469  // no children of its own, we need to refine it, and if it does
470  // already have children then loop over all children and see if they
471  // are locally available as well
472  if (dealii_cell->is_active())
473  dealii_cell->set_refine_flag();
474  else
475  {
478  for (unsigned int c = 0;
479  c < GeometryInfo<dim>::max_children_per_cell;
480  ++c)
481  switch (dim)
482  {
483  case 2:
484  P4EST_QUADRANT_INIT(&p4est_child[c]);
485  break;
486  case 3:
487  P8EST_QUADRANT_INIT(&p4est_child[c]);
488  break;
489  default:
490  Assert(false, ExcNotImplemented());
491  }
492 
493 
495  p4est_child);
496 
497  for (unsigned int c = 0;
498  c < GeometryInfo<dim>::max_children_per_cell;
499  ++c)
501  const_cast<typename internal::p4est::types<dim>::tree *>(
502  &tree),
503  &p4est_child[c]) == false)
504  {
505  // no, this child is locally not available in the p4est.
506  // delete all its children but, because this may not be
507  // successful, make sure to mark all children recursively
508  // as not local.
509  delete_all_children<dim, spacedim>(dealii_cell->child(c));
510  dealii_cell->child(c)->recursively_set_subdomain_id(
512  }
513  else
514  {
515  // at least some part of the tree rooted in this child is
516  // locally available
517  match_tree_recursively<dim, spacedim>(tree,
518  dealii_cell->child(c),
519  p4est_child[c],
520  forest,
521  my_subdomain);
522  }
523  }
524  }
525  }
526 
527 
528  template <int dim, int spacedim>
529  void
530  match_quadrant(
531  const ::Triangulation<dim, spacedim> * tria,
532  unsigned int dealii_index,
533  const typename internal::p4est::types<dim>::quadrant &ghost_quadrant,
534  types::subdomain_id ghost_owner)
535  {
536  const int l = ghost_quadrant.level;
537 
538  for (int i = 0; i < l; ++i)
539  {
541  i,
542  dealii_index);
543  if (cell->is_active())
544  {
545  cell->clear_coarsen_flag();
546  cell->set_refine_flag();
547  return;
548  }
549 
550  const int child_id =
552  i + 1);
553  dealii_index = cell->child_index(child_id);
554  }
555 
557  l,
558  dealii_index);
559  if (cell->has_children())
560  delete_all_children<dim, spacedim>(cell);
561  else
562  {
563  cell->clear_coarsen_flag();
564  cell->set_subdomain_id(ghost_owner);
565  }
566  }
567 
568 
569 # ifdef P4EST_SEARCH_LOCAL
570  template <int dim>
571  class PartitionSearch
572  {
573  public:
574  PartitionSearch()
575  {
576  Assert(dim > 1, ExcNotImplemented());
577  }
578 
579  PartitionSearch(const PartitionSearch<dim> &other) = delete;
580 
581  PartitionSearch<dim> &
582  operator=(const PartitionSearch<dim> &other) = delete;
583 
584  public:
594  static int
595  local_quadrant_fn(typename internal::p4est::types<dim>::forest *forest,
596  typename internal::p4est::types<dim>::topidx which_tree,
597  typename internal::p4est::types<dim>::quadrant *quadrant,
598  int rank_begin,
599  int rank_end,
600  void *point);
601 
614  static int
615  local_point_fn(typename internal::p4est::types<dim>::forest * forest,
616  typename internal::p4est::types<dim>::topidx which_tree,
617  typename internal::p4est::types<dim>::quadrant *quadrant,
618  int rank_begin,
619  int rank_end,
620  void * point);
621 
622  private:
627  class QuadrantData
628  {
629  public:
630  QuadrantData();
631 
632  void
633  set_cell_vertices(
634  typename internal::p4est::types<dim>::forest * forest,
635  typename internal::p4est::types<dim>::topidx which_tree,
636  typename internal::p4est::types<dim>::quadrant *quadrant,
638  quad_length_on_level);
639 
640  void
641  initialize_mapping();
642 
643  Point<dim>
644  map_real_to_unit_cell(const Point<dim> &p) const;
645 
646  bool
647  is_in_this_quadrant(const Point<dim> &p) const;
648 
649  private:
650  std::vector<Point<dim>> cell_vertices;
651 
656  FullMatrix<double> quadrant_mapping_matrix;
657 
658  bool are_vertices_initialized;
659 
660  bool is_reference_mapping_initialized;
661  };
662 
666  QuadrantData quadrant_data;
667  }; // class PartitionSearch
668 
669 
670 
671  template <int dim>
672  int
673  PartitionSearch<dim>::local_quadrant_fn(
674  typename internal::p4est::types<dim>::forest * forest,
675  typename internal::p4est::types<dim>::topidx which_tree,
676  typename internal::p4est::types<dim>::quadrant *quadrant,
677  int /* rank_begin */,
678  int /* rank_end */,
679  void * /* this is always nullptr */ point)
680  {
681  // point must be nullptr here
682  (void)point;
683  Assert(point == nullptr, ::ExcInternalError());
684 
685  // we need the user pointer
686  // note that this is not available since function is static
687  PartitionSearch<dim> *this_object =
688  reinterpret_cast<PartitionSearch<dim> *>(forest->user_pointer);
689 
690  // Avoid p4est macros, instead do bitshifts manually with fixed size types
692  quad_length_on_level =
693  1 << (static_cast<typename internal::p4est::types<dim>::quadrant_coord>(
694  (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
695  static_cast<typename internal::p4est::types<dim>::quadrant_coord>(
696  quadrant->level));
697 
698  this_object->quadrant_data.set_cell_vertices(forest,
699  which_tree,
700  quadrant,
701  quad_length_on_level);
702 
703  // from cell vertices we can initialize the mapping
704  this_object->quadrant_data.initialize_mapping();
705 
706  // always return true since we must decide by point
707  return /* true */ 1;
708  }
709 
710 
711 
712  template <int dim>
713  int
714  PartitionSearch<dim>::local_point_fn(
715  typename internal::p4est::types<dim>::forest *forest,
716  typename internal::p4est::types<dim>::topidx /* which_tree */,
717  typename internal::p4est::types<dim>::quadrant * /* quadrant */,
718  int rank_begin,
719  int rank_end,
720  void *point)
721  {
722  // point must NOT be be nullptr here
723  Assert(point != nullptr, ::ExcInternalError());
724 
725  // we need the user pointer
726  // note that this is not available since function is static
727  PartitionSearch<dim> *this_object =
728  reinterpret_cast<PartitionSearch<dim> *>(forest->user_pointer);
729 
730  // point with rank as double pointer
731  double *this_point_dptr = static_cast<double *>(point);
732 
733  Point<dim> this_point =
734  (dim == 2 ? Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
735  Point<dim>(this_point_dptr[0],
736  this_point_dptr[1],
737  this_point_dptr[2]));
738 
739  // use reference mapping to decide whether this point is in this quadrant
740  const bool is_in_this_quadrant =
741  this_object->quadrant_data.is_in_this_quadrant(this_point);
742 
743 
744 
745  if (!is_in_this_quadrant)
746  {
747  // no need to search further, stop recursion
748  return /* false */ 0;
749  }
750 
751 
752 
753  // From here we have a candidate
754  if (rank_begin < rank_end)
755  {
756  // continue recursion
757  return /* true */ 1;
758  }
759 
760  // Now, we know that the point is found (rank_begin==rank_end) and we have
761  // the MPI rank, so no need to search further.
762  this_point_dptr[dim] = static_cast<double>(rank_begin);
763 
764  // stop recursion.
765  return /* false */ 0;
766  }
767 
768 
769 
770  template <int dim>
771  bool
772  PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
773  const Point<dim> &p) const
774  {
775  const Point<dim> p_ref = map_real_to_unit_cell(p);
776 
778  }
779 
780 
781 
782  template <int dim>
783  Point<dim>
784  PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
785  const Point<dim> &p) const
786  {
787  Assert(is_reference_mapping_initialized,
788  ::ExcMessage(
789  "Cell vertices and mapping coefficients must be fully "
790  "initialized before transforming a point to the unit cell."));
791 
792  Point<dim> p_out;
793 
794  if (dim == 2)
795  {
796  for (unsigned int alpha = 0;
797  alpha < GeometryInfo<dim>::vertices_per_cell;
798  ++alpha)
799  {
800  const Point<dim> &p_ref =
802 
803  p_out += (quadrant_mapping_matrix(alpha, 0) +
804  quadrant_mapping_matrix(alpha, 1) * p(0) +
805  quadrant_mapping_matrix(alpha, 2) * p(1) +
806  quadrant_mapping_matrix(alpha, 3) * p(0) * p(1)) *
807  p_ref;
808  }
809  }
810  else
811  {
812  for (unsigned int alpha = 0;
813  alpha < GeometryInfo<dim>::vertices_per_cell;
814  ++alpha)
815  {
816  const Point<dim> &p_ref =
818 
819  p_out += (quadrant_mapping_matrix(alpha, 0) +
820  quadrant_mapping_matrix(alpha, 1) * p(0) +
821  quadrant_mapping_matrix(alpha, 2) * p(1) +
822  quadrant_mapping_matrix(alpha, 3) * p(2) +
823  quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
824  quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
825  quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
826  quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
827  p_ref;
828  }
829  }
830 
831  return p_out;
832  }
833 
834 
835  template <int dim>
836  PartitionSearch<dim>::QuadrantData::QuadrantData()
837  : cell_vertices(GeometryInfo<dim>::vertices_per_cell)
838  , quadrant_mapping_matrix(GeometryInfo<dim>::vertices_per_cell,
839  GeometryInfo<dim>::vertices_per_cell)
840  , are_vertices_initialized(false)
841  , is_reference_mapping_initialized(false)
842  {}
843 
844 
845 
846  template <int dim>
847  void
848  PartitionSearch<dim>::QuadrantData::initialize_mapping()
849  {
850  Assert(
851  are_vertices_initialized,
852  ::ExcMessage(
853  "Cell vertices must be initialized before the cell mapping can be filled."));
854 
857 
858  if (dim == 2)
859  {
860  for (unsigned int alpha = 0;
861  alpha < GeometryInfo<dim>::vertices_per_cell;
862  ++alpha)
863  {
864  // point matrix to be inverted
865  point_matrix(0, alpha) = 1;
866  point_matrix(1, alpha) = cell_vertices[alpha](0);
867  point_matrix(2, alpha) = cell_vertices[alpha](1);
868  point_matrix(3, alpha) =
869  cell_vertices[alpha](0) * cell_vertices[alpha](1);
870  }
871 
872  /*
873  * Rows of quadrant_mapping_matrix are the coefficients of the basis
874  * on the physical cell
875  */
876  quadrant_mapping_matrix.invert(point_matrix);
877  }
878  else
879  {
880  for (unsigned int alpha = 0;
881  alpha < GeometryInfo<dim>::vertices_per_cell;
882  ++alpha)
883  {
884  // point matrix to be inverted
885  point_matrix(0, alpha) = 1;
886  point_matrix(1, alpha) = cell_vertices[alpha](0);
887  point_matrix(2, alpha) = cell_vertices[alpha](1);
888  point_matrix(3, alpha) = cell_vertices[alpha](2);
889  point_matrix(4, alpha) =
890  cell_vertices[alpha](0) * cell_vertices[alpha](1);
891  point_matrix(5, alpha) =
892  cell_vertices[alpha](1) * cell_vertices[alpha](2);
893  point_matrix(6, alpha) =
894  cell_vertices[alpha](0) * cell_vertices[alpha](2);
895  point_matrix(7, alpha) = cell_vertices[alpha](0) *
896  cell_vertices[alpha](1) *
897  cell_vertices[alpha](2);
898  }
899 
900  /*
901  * Rows of quadrant_mapping_matrix are the coefficients of the basis
902  * on the physical cell
903  */
904  quadrant_mapping_matrix.invert(point_matrix);
905  }
906 
907  is_reference_mapping_initialized = true;
908  }
909 
910 
911 
912  template <>
913  void
914  PartitionSearch<2>::QuadrantData::set_cell_vertices(
915  typename internal::p4est::types<2>::forest * forest,
916  typename internal::p4est::types<2>::topidx which_tree,
917  typename internal::p4est::types<2>::quadrant *quadrant,
919  quad_length_on_level)
920  {
921  constexpr unsigned int dim = 2;
922 
923  // p4est for some reason always needs double vxyz[3] as last argument to
924  // quadrant_coord_to_vertex
925  double corner_point[dim + 1] = {0};
926 
927  // A lambda to avoid code duplication.
928  const auto copy_vertex = [&](unsigned int vertex_index) -> void {
929  // copy into local struct
930  for (unsigned int d = 0; d < dim; ++d)
931  {
932  cell_vertices[vertex_index](d) = corner_point[d];
933  // reset
934  corner_point[d] = 0;
935  }
936  };
937 
938  // Fill points of QuadrantData in lexicographic order
939  /*
940  * Corner #0
941  */
942  unsigned int vertex_index = 0;
944  forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
945 
946  // copy into local struct
947  copy_vertex(vertex_index);
948 
949  /*
950  * Corner #1
951  */
952  vertex_index = 1;
954  forest->connectivity,
955  which_tree,
956  quadrant->x + quad_length_on_level,
957  quadrant->y,
958  corner_point);
959 
960  // copy into local struct
961  copy_vertex(vertex_index);
962 
963  /*
964  * Corner #2
965  */
966  vertex_index = 2;
968  forest->connectivity,
969  which_tree,
970  quadrant->x,
971  quadrant->y + quad_length_on_level,
972  corner_point);
973 
974  // copy into local struct
975  copy_vertex(vertex_index);
976 
977  /*
978  * Corner #3
979  */
980  vertex_index = 3;
982  forest->connectivity,
983  which_tree,
984  quadrant->x + quad_length_on_level,
985  quadrant->y + quad_length_on_level,
986  corner_point);
987 
988  // copy into local struct
989  copy_vertex(vertex_index);
990 
991  are_vertices_initialized = true;
992  }
993 
994 
995 
996  template <>
997  void
998  PartitionSearch<3>::QuadrantData::set_cell_vertices(
999  typename internal::p4est::types<3>::forest * forest,
1000  typename internal::p4est::types<3>::topidx which_tree,
1001  typename internal::p4est::types<3>::quadrant *quadrant,
1003  quad_length_on_level)
1004  {
1005  constexpr unsigned int dim = 3;
1006 
1007  double corner_point[dim] = {0};
1008 
1009  // A lambda to avoid code duplication.
1010  auto copy_vertex = [&](unsigned int vertex_index) -> void {
1011  // copy into local struct
1012  for (unsigned int d = 0; d < dim; ++d)
1013  {
1014  cell_vertices[vertex_index](d) = corner_point[d];
1015  // reset
1016  corner_point[d] = 0;
1017  }
1018  };
1019 
1020  // Fill points of QuadrantData in lexicographic order
1021  /*
1022  * Corner #0
1023  */
1024  unsigned int vertex_index = 0;
1026  forest->connectivity,
1027  which_tree,
1028  quadrant->x,
1029  quadrant->y,
1030  quadrant->z,
1031  corner_point);
1032 
1033  // copy into local struct
1034  copy_vertex(vertex_index);
1035 
1036 
1037  /*
1038  * Corner #1
1039  */
1040  vertex_index = 1;
1042  forest->connectivity,
1043  which_tree,
1044  quadrant->x + quad_length_on_level,
1045  quadrant->y,
1046  quadrant->z,
1047  corner_point);
1048 
1049  // copy into local struct
1050  copy_vertex(vertex_index);
1051 
1052  /*
1053  * Corner #2
1054  */
1055  vertex_index = 2;
1057  forest->connectivity,
1058  which_tree,
1059  quadrant->x,
1060  quadrant->y + quad_length_on_level,
1061  quadrant->z,
1062  corner_point);
1063 
1064  // copy into local struct
1065  copy_vertex(vertex_index);
1066 
1067  /*
1068  * Corner #3
1069  */
1070  vertex_index = 3;
1072  forest->connectivity,
1073  which_tree,
1074  quadrant->x + quad_length_on_level,
1075  quadrant->y + quad_length_on_level,
1076  quadrant->z,
1077  corner_point);
1078 
1079  // copy into local struct
1080  copy_vertex(vertex_index);
1081 
1082  /*
1083  * Corner #4
1084  */
1085  vertex_index = 4;
1087  forest->connectivity,
1088  which_tree,
1089  quadrant->x,
1090  quadrant->y,
1091  quadrant->z + quad_length_on_level,
1092  corner_point);
1093 
1094  // copy into local struct
1095  copy_vertex(vertex_index);
1096 
1097  /*
1098  * Corner #5
1099  */
1100  vertex_index = 5;
1102  forest->connectivity,
1103  which_tree,
1104  quadrant->x + quad_length_on_level,
1105  quadrant->y,
1106  quadrant->z + quad_length_on_level,
1107  corner_point);
1108 
1109  // copy into local struct
1110  copy_vertex(vertex_index);
1111 
1112  /*
1113  * Corner #6
1114  */
1115  vertex_index = 6;
1117  forest->connectivity,
1118  which_tree,
1119  quadrant->x,
1120  quadrant->y + quad_length_on_level,
1121  quadrant->z + quad_length_on_level,
1122  corner_point);
1123 
1124  // copy into local struct
1125  copy_vertex(vertex_index);
1126 
1127  /*
1128  * Corner #7
1129  */
1130  vertex_index = 7;
1132  forest->connectivity,
1133  which_tree,
1134  quadrant->x + quad_length_on_level,
1135  quadrant->y + quad_length_on_level,
1136  quadrant->z + quad_length_on_level,
1137  corner_point);
1138 
1139  // copy into local struct
1140  copy_vertex(vertex_index);
1141 
1142 
1143  are_vertices_initialized = true;
1144  }
1145 # endif // P4EST_SEARCH_LOCAL defined
1146 
1147 
1153  template <int dim, int spacedim>
1154  class RefineAndCoarsenList
1155  {
1156  public:
1157  RefineAndCoarsenList(const Triangulation<dim, spacedim> &triangulation,
1158  const std::vector<types::global_dof_index>
1159  &p4est_tree_to_coarse_cell_permutation,
1160  const types::subdomain_id my_subdomain);
1161 
1170  static int
1171  refine_callback(
1172  typename internal::p4est::types<dim>::forest * forest,
1173  typename internal::p4est::types<dim>::topidx coarse_cell_index,
1174  typename internal::p4est::types<dim>::quadrant *quadrant);
1175 
1180  static int
1181  coarsen_callback(
1182  typename internal::p4est::types<dim>::forest * forest,
1183  typename internal::p4est::types<dim>::topidx coarse_cell_index,
1184  typename internal::p4est::types<dim>::quadrant *children[]);
1185 
1186  bool
1187  pointers_are_at_end() const;
1188 
1189  private:
1190  std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1191  typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1192  const_iterator current_refine_pointer;
1193 
1194  std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1195  typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1196  const_iterator current_coarsen_pointer;
1197 
1198  void
1199  build_lists(
1200  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1201  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
1202  const types::subdomain_id myid);
1203  };
1204 
1205 
1206 
1207  template <int dim, int spacedim>
1208  bool
1209  RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end() const
1210  {
1211  return ((current_refine_pointer == refine_list.end()) &&
1212  (current_coarsen_pointer == coarsen_list.end()));
1213  }
1214 
1215 
1216 
1217  template <int dim, int spacedim>
1218  RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1220  const std::vector<types::global_dof_index>
1221  & p4est_tree_to_coarse_cell_permutation,
1222  const types::subdomain_id my_subdomain)
1223  {
1224  // count how many flags are set and allocate that much memory
1225  unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1226  for (const auto &cell : triangulation.active_cell_iterators())
1227  {
1228  // skip cells that are not local
1229  if (cell->subdomain_id() != my_subdomain)
1230  continue;
1231 
1232  if (cell->refine_flag_set())
1233  ++n_refine_flags;
1234  else if (cell->coarsen_flag_set())
1235  ++n_coarsen_flags;
1236  }
1237 
1238  refine_list.reserve(n_refine_flags);
1239  coarsen_list.reserve(n_coarsen_flags);
1240 
1241 
1242  // now build the lists of cells that are flagged. note that p4est will
1243  // traverse its cells in the order in which trees appear in the
1244  // forest. this order is not the same as the order of coarse cells in the
1245  // deal.II Triangulation because we have translated everything by the
1246  // coarse_cell_to_p4est_tree_permutation permutation. in order to make
1247  // sure that the output array is already in the correct order, traverse
1248  // our coarse cells in the same order in which p4est will:
1249  for (unsigned int c = 0; c < triangulation.n_cells(0); ++c)
1250  {
1251  unsigned int coarse_cell_index =
1252  p4est_tree_to_coarse_cell_permutation[c];
1253 
1254  const typename Triangulation<dim, spacedim>::cell_iterator cell(
1255  &triangulation, 0, coarse_cell_index);
1256 
1257  typename internal::p4est::types<dim>::quadrant p4est_cell;
1259  /*level=*/0,
1260  /*index=*/0);
1261  p4est_cell.p.which_tree = c;
1262  build_lists(cell, p4est_cell, my_subdomain);
1263  }
1264 
1265 
1266  Assert(refine_list.size() == n_refine_flags, ExcInternalError());
1267  Assert(coarsen_list.size() == n_coarsen_flags, ExcInternalError());
1268 
1269  // make sure that our ordering in fact worked
1270  for (unsigned int i = 1; i < refine_list.size(); ++i)
1271  Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1272  ExcInternalError());
1273  for (unsigned int i = 1; i < coarsen_list.size(); ++i)
1274  Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1275  ExcInternalError());
1276 
1277  current_refine_pointer = refine_list.begin();
1278  current_coarsen_pointer = coarsen_list.begin();
1279  }
1280 
1281 
1282 
1283  template <int dim, int spacedim>
1284  void
1285  RefineAndCoarsenList<dim, spacedim>::build_lists(
1286  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1287  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
1288  const types::subdomain_id my_subdomain)
1289  {
1290  if (cell->is_active())
1291  {
1292  if (cell->subdomain_id() == my_subdomain)
1293  {
1294  if (cell->refine_flag_set())
1295  refine_list.push_back(p4est_cell);
1296  else if (cell->coarsen_flag_set())
1297  coarsen_list.push_back(p4est_cell);
1298  }
1299  }
1300  else
1301  {
1304  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1305  ++c)
1306  switch (dim)
1307  {
1308  case 2:
1309  P4EST_QUADRANT_INIT(&p4est_child[c]);
1310  break;
1311  case 3:
1312  P8EST_QUADRANT_INIT(&p4est_child[c]);
1313  break;
1314  default:
1315  Assert(false, ExcNotImplemented());
1316  }
1318  p4est_child);
1319  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1320  ++c)
1321  {
1322  p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1323  build_lists(cell->child(c), p4est_child[c], my_subdomain);
1324  }
1325  }
1326  }
1327 
1328 
1329  template <int dim, int spacedim>
1330  int
1331  RefineAndCoarsenList<dim, spacedim>::refine_callback(
1332  typename internal::p4est::types<dim>::forest * forest,
1333  typename internal::p4est::types<dim>::topidx coarse_cell_index,
1334  typename internal::p4est::types<dim>::quadrant *quadrant)
1335  {
1336  RefineAndCoarsenList<dim, spacedim> *this_object =
1337  reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *>(
1338  forest->user_pointer);
1339 
1340  // if there are no more cells in our list the current cell can't be
1341  // flagged for refinement
1342  if (this_object->current_refine_pointer == this_object->refine_list.end())
1343  return 0;
1344 
1345  Assert(coarse_cell_index <=
1346  this_object->current_refine_pointer->p.which_tree,
1347  ExcInternalError());
1348 
1349  // if p4est hasn't yet reached the tree of the next flagged cell the
1350  // current cell can't be flagged for refinement
1351  if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1352  return 0;
1353 
1354  // now we're in the right tree in the forest
1355  Assert(coarse_cell_index <=
1356  this_object->current_refine_pointer->p.which_tree,
1357  ExcInternalError());
1358 
1359  // make sure that the p4est loop over cells hasn't gotten ahead of our own
1360  // pointer
1362  quadrant, &*this_object->current_refine_pointer) <= 0,
1363  ExcInternalError());
1364 
1365  // now, if the p4est cell is one in the list, it is supposed to be refined
1367  quadrant, &*this_object->current_refine_pointer))
1368  {
1369  ++this_object->current_refine_pointer;
1370  return 1;
1371  }
1372 
1373  // p4est cell is not in list
1374  return 0;
1375  }
1376 
1377 
1378 
1379  template <int dim, int spacedim>
1380  int
1381  RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1382  typename internal::p4est::types<dim>::forest * forest,
1383  typename internal::p4est::types<dim>::topidx coarse_cell_index,
1384  typename internal::p4est::types<dim>::quadrant *children[])
1385  {
1386  RefineAndCoarsenList<dim, spacedim> *this_object =
1387  reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *>(
1388  forest->user_pointer);
1389 
1390  // if there are no more cells in our list the current cell can't be
1391  // flagged for coarsening
1392  if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1393  return 0;
1394 
1395  Assert(coarse_cell_index <=
1396  this_object->current_coarsen_pointer->p.which_tree,
1397  ExcInternalError());
1398 
1399  // if p4est hasn't yet reached the tree of the next flagged cell the
1400  // current cell can't be flagged for coarsening
1401  if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1402  return 0;
1403 
1404  // now we're in the right tree in the forest
1405  Assert(coarse_cell_index <=
1406  this_object->current_coarsen_pointer->p.which_tree,
1407  ExcInternalError());
1408 
1409  // make sure that the p4est loop over cells hasn't gotten ahead of our own
1410  // pointer
1412  children[0], &*this_object->current_coarsen_pointer) <= 0,
1413  ExcInternalError());
1414 
1415  // now, if the p4est cell is one in the list, it is supposed to be
1416  // coarsened
1418  children[0], &*this_object->current_coarsen_pointer))
1419  {
1420  // move current pointer one up
1421  ++this_object->current_coarsen_pointer;
1422 
1423  // note that the next 3 cells in our list need to correspond to the
1424  // other siblings of the cell we have just found
1425  for (unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1426  ++c)
1427  {
1429  children[c], &*this_object->current_coarsen_pointer),
1430  ExcInternalError());
1431  ++this_object->current_coarsen_pointer;
1432  }
1433 
1434  return 1;
1435  }
1436 
1437  // p4est cell is not in list
1438  return 0;
1439  }
1440 
1441 
1442 
1449  template <int dim, int spacedim>
1450  class PartitionWeights
1451  {
1452  public:
1458  explicit PartitionWeights(const std::vector<unsigned int> &cell_weights);
1459 
1467  static int
1468  cell_weight(typename internal::p4est::types<dim>::forest *forest,
1469  typename internal::p4est::types<dim>::topidx coarse_cell_index,
1470  typename internal::p4est::types<dim>::quadrant *quadrant);
1471 
1472  private:
1473  std::vector<unsigned int> cell_weights_list;
1474  std::vector<unsigned int>::const_iterator current_pointer;
1475  };
1476 
1477 
1478  template <int dim, int spacedim>
1479  PartitionWeights<dim, spacedim>::PartitionWeights(
1480  const std::vector<unsigned int> &cell_weights)
1481  : cell_weights_list(cell_weights)
1482  {
1483  // set the current pointer to the first element of the list, given that
1484  // we will walk through it sequentially
1485  current_pointer = cell_weights_list.begin();
1486  }
1487 
1488 
1489  template <int dim, int spacedim>
1490  int
1491  PartitionWeights<dim, spacedim>::cell_weight(
1492  typename internal::p4est::types<dim>::forest *forest,
1495  {
1496  // the function gets two additional arguments, but we don't need them
1497  // since we know in which order p4est will walk through the cells
1498  // and have already built our weight lists in this order
1499 
1500  PartitionWeights<dim, spacedim> *this_object =
1501  reinterpret_cast<PartitionWeights<dim, spacedim> *>(forest->user_pointer);
1502 
1503  Assert(this_object->current_pointer >=
1504  this_object->cell_weights_list.begin(),
1505  ExcInternalError());
1506  Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1507  ExcInternalError());
1508 
1509  // Get the weight, increment the pointer, and return the weight. Also
1510  // make sure that we don't exceed the 'int' data type that p4est uses
1511  // to represent weights
1512  const unsigned int weight = *this_object->current_pointer;
1513  ++this_object->current_pointer;
1514 
1515  Assert(weight < static_cast<unsigned int>(std::numeric_limits<int>::max()),
1516  ExcMessage("p4est uses 'signed int' to represent the partition "
1517  "weights for cells. The weight provided here exceeds "
1518  "the maximum value represented as a 'signed int'."));
1519  return static_cast<int>(weight);
1520  }
1521 
1522  template <int dim, int spacedim>
1523  using cell_relation_t = typename std::pair<
1524  typename ::Triangulation<dim, spacedim>::cell_iterator,
1525  typename ::Triangulation<dim, spacedim>::CellStatus>;
1526 
1536  template <int dim, int spacedim>
1537  inline void
1538  add_single_cell_relation(
1539  std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
1540  const typename ::internal::p4est::types<dim>::tree & tree,
1541  const unsigned int idx,
1542  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
1543  const typename Triangulation<dim, spacedim>::CellStatus status)
1544  {
1545  const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1546 
1547  // check if we will be writing into valid memory
1548  Assert(local_quadrant_index < cell_rel.size(), ExcInternalError());
1549 
1550  // store relation
1551  cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1552  }
1553 
1554 
1555 
1565  template <int dim, int spacedim>
1566  void
1567  update_cell_relations_recursively(
1568  std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
1569  const typename ::internal::p4est::types<dim>::tree & tree,
1570  const typename Triangulation<dim, spacedim>::cell_iterator & dealii_cell,
1571  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1572  {
1573  // find index of p4est_cell in the quadrants array of the corresponding tree
1574  const int idx = sc_array_bsearch(
1575  const_cast<sc_array_t *>(&tree.quadrants),
1576  &p4est_cell,
1578  if (idx == -1 &&
1580  const_cast<typename ::internal::p4est::types<dim>::tree *>(
1581  &tree),
1582  &p4est_cell) == false))
1583  // this quadrant and none of its children belong to us.
1584  return;
1585 
1586  // recurse further if both p4est and dealii still have children
1587  const bool p4est_has_children = (idx == -1);
1588  if (p4est_has_children && dealii_cell->has_children())
1589  {
1590  // recurse further
1591  typename ::internal::p4est::types<dim>::quadrant
1593 
1594  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1595  ++c)
1596  switch (dim)
1597  {
1598  case 2:
1599  P4EST_QUADRANT_INIT(&p4est_child[c]);
1600  break;
1601  case 3:
1602  P8EST_QUADRANT_INIT(&p4est_child[c]);
1603  break;
1604  default:
1605  Assert(false, ExcNotImplemented());
1606  }
1607 
1609  &p4est_cell, p4est_child);
1610 
1611  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1612  ++c)
1613  {
1614  update_cell_relations_recursively<dim, spacedim>(
1615  cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1616  }
1617  }
1618  else if (!p4est_has_children && !dealii_cell->has_children())
1619  {
1620  // this active cell didn't change
1621  // save pair into corresponding position
1622  add_single_cell_relation<dim, spacedim>(
1623  cell_rel,
1624  tree,
1625  idx,
1626  dealii_cell,
1628  }
1629  else if (p4est_has_children) // based on the conditions above, we know that
1630  // dealii_cell has no children
1631  {
1632  // this cell got refined in p4est, but the dealii_cell has not yet been
1633  // refined
1634 
1635  // this quadrant is not active
1636  // generate its children, and store information in those
1637  typename ::internal::p4est::types<dim>::quadrant
1639  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1640  ++c)
1641  switch (dim)
1642  {
1643  case 2:
1644  P4EST_QUADRANT_INIT(&p4est_child[c]);
1645  break;
1646  case 3:
1647  P8EST_QUADRANT_INIT(&p4est_child[c]);
1648  break;
1649  default:
1650  Assert(false, ExcNotImplemented());
1651  }
1652 
1654  &p4est_cell, p4est_child);
1655 
1656  // mark first child with CELL_REFINE and the remaining children with
1657  // CELL_INVALID, but associate them all with the parent cell unpack
1658  // algorithm will be called only on CELL_REFINE flagged quadrant
1659  int child_idx;
1660  typename Triangulation<dim, spacedim>::CellStatus cell_status;
1661  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1662  ++i)
1663  {
1664  child_idx = sc_array_bsearch(
1665  const_cast<sc_array_t *>(&tree.quadrants),
1666  &p4est_child[i],
1668 
1669  cell_status = (i == 0) ? Triangulation<dim, spacedim>::CELL_REFINE :
1671 
1672  add_single_cell_relation<dim, spacedim>(
1673  cell_rel, tree, child_idx, dealii_cell, cell_status);
1674  }
1675  }
1676  else // based on the conditions above, we know that p4est_cell has no
1677  // children, and the dealii_cell does
1678  {
1679  // its children got coarsened into this cell in p4est,
1680  // but the dealii_cell still has its children
1681  add_single_cell_relation<dim, spacedim>(
1682  cell_rel,
1683  tree,
1684  idx,
1685  dealii_cell,
1687  }
1688  }
1689 } // namespace
1690 
1691 
1692 
1693 namespace parallel
1694 {
1695  namespace distributed
1696  {
1697  /*----------------- class Triangulation<dim,spacedim> ---------------\*/
1698  template <int dim, int spacedim>
1700  const MPI_Comm &mpi_communicator,
1701  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
1702  smooth_grid,
1703  const Settings settings)
1704  : // Do not check for distorted cells.
1705  // For multigrid, we need limit_level_difference_at_vertices
1706  // to make sure the transfer operators only need to consider two levels.
1707  ::parallel::DistributedTriangulationBase<dim, spacedim>(
1708  mpi_communicator,
1710  static_cast<
1711  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
1712  smooth_grid |
1713  Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1714  smooth_grid,
1715  false)
1716  , settings(settings)
1717  , triangulation_has_content(false)
1718  , connectivity(nullptr)
1719  , parallel_forest(nullptr)
1720  {
1721  parallel_ghost = nullptr;
1722  }
1723 
1724 
1725 
1726  template <int dim, int spacedim>
1728  {
1729  // virtual functions called in constructors and destructors never use the
1730  // override in a derived class
1731  // for clarity be explicit on which function is called
1732  try
1733  {
1735  }
1736  catch (...)
1737  {}
1738 
1739  AssertNothrow(triangulation_has_content == false, ExcInternalError());
1740  AssertNothrow(connectivity == nullptr, ExcInternalError());
1741  AssertNothrow(parallel_forest == nullptr, ExcInternalError());
1742  }
1743 
1744 
1745 
1746  template <int dim, int spacedim>
1747  void
1749  const std::vector<Point<spacedim>> &vertices,
1750  const std::vector<CellData<dim>> & cells,
1751  const SubCellData & subcelldata)
1752  {
1753  try
1754  {
1756  vertices, cells, subcelldata);
1757  }
1758  catch (
1759  const typename ::Triangulation<dim, spacedim>::DistortedCellList
1760  &)
1761  {
1762  // the underlying triangulation should not be checking for distorted
1763  // cells
1764  Assert(false, ExcInternalError());
1765  }
1766 
1767  Assert(
1768  this->all_reference_cells_are_hyper_cube(),
1769  ExcMessage(
1770  "The class parallel::distributed::Triangulation only supports meshes "
1771  "consisting only of hypercube-like cells."));
1772 
1773  // note that now we have some content in the p4est objects and call the
1774  // functions that do the actual work (which are dimension dependent, so
1775  // separate)
1776  triangulation_has_content = true;
1777 
1778  setup_coarse_cell_to_p4est_tree_permutation();
1779 
1780  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1781 
1782  try
1783  {
1784  copy_local_forest_to_triangulation();
1785  }
1786  catch (const typename Triangulation<dim>::DistortedCellList &)
1787  {
1788  // the underlying triangulation should not be checking for distorted
1789  // cells
1790  Assert(false, ExcInternalError());
1791  }
1792 
1793  this->update_periodic_face_map();
1794  this->update_number_cache();
1795  }
1796 
1797 
1798 
1799  template <int dim, int spacedim>
1800  void
1803  &construction_data)
1804  {
1805  (void)construction_data;
1806 
1807  Assert(false, ExcInternalError());
1808  }
1809 
1810 
1811 
1812  template <int dim, int spacedim>
1813  void
1815  {
1816  triangulation_has_content = false;
1817 
1818  if (parallel_ghost != nullptr)
1819  {
1821  parallel_ghost);
1822  parallel_ghost = nullptr;
1823  }
1824 
1825  if (parallel_forest != nullptr)
1826  {
1828  parallel_forest = nullptr;
1829  }
1830 
1831  if (connectivity != nullptr)
1832  {
1834  connectivity);
1835  connectivity = nullptr;
1836  }
1837 
1838  coarse_cell_to_p4est_tree_permutation.resize(0);
1839  p4est_tree_to_coarse_cell_permutation.resize(0);
1840 
1842 
1843  this->update_number_cache();
1844  }
1845 
1846 
1847 
1848  template <int dim, int spacedim>
1849  bool
1851  {
1852  return settings &
1854  }
1855 
1856 
1857 
1858  template <int dim, int spacedim>
1859  bool
1861  {
1862  return settings &
1864  }
1865 
1866 
1867 
1868  template <int dim, int spacedim>
1869  void
1871  const typename ::internal::p4est::types<dim>::forest
1872  *parallel_forest,
1873  const typename ::internal::p4est::types<dim>::gloidx
1874  *previous_global_first_quadrant)
1875  {
1876  Assert(this->data_transfer.sizes_fixed_cumulative.size() > 0,
1877  ExcMessage("No data has been packed!"));
1878 
1879  // Resize memory according to the data that we will receive.
1880  this->data_transfer.dest_data_fixed.resize(
1881  parallel_forest->local_num_quadrants *
1882  this->data_transfer.sizes_fixed_cumulative.back());
1883 
1884  // Execute non-blocking fixed size transfer.
1885  typename ::internal::p4est::types<dim>::transfer_context
1886  *tf_context;
1887  tf_context =
1889  parallel_forest->global_first_quadrant,
1890  previous_global_first_quadrant,
1891  parallel_forest->mpicomm,
1892  0,
1893  this->data_transfer.dest_data_fixed.data(),
1894  this->data_transfer.src_data_fixed.data(),
1895  this->data_transfer.sizes_fixed_cumulative.back());
1896 
1897  if (this->data_transfer.variable_size_data_stored)
1898  {
1899  // Resize memory according to the data that we will receive.
1900  this->data_transfer.dest_sizes_variable.resize(
1901  parallel_forest->local_num_quadrants);
1902 
1903  // Execute fixed size transfer of data sizes for variable size
1904  // transfer.
1906  parallel_forest->global_first_quadrant,
1907  previous_global_first_quadrant,
1908  parallel_forest->mpicomm,
1909  1,
1910  this->data_transfer.dest_sizes_variable.data(),
1911  this->data_transfer.src_sizes_variable.data(),
1912  sizeof(unsigned int));
1913  }
1914 
1916 
1917  // Release memory of previously packed data.
1918  this->data_transfer.src_data_fixed.clear();
1919  this->data_transfer.src_data_fixed.shrink_to_fit();
1920 
1921  if (this->data_transfer.variable_size_data_stored)
1922  {
1923  // Resize memory according to the data that we will receive.
1924  this->data_transfer.dest_data_variable.resize(
1925  std::accumulate(this->data_transfer.dest_sizes_variable.begin(),
1926  this->data_transfer.dest_sizes_variable.end(),
1928 
1929 # if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1930 # else
1931  // ----- WORKAROUND -----
1932  // An assertion in p4est prevents us from sending/receiving no data
1933  // at all, which is mandatory if one of our processes does not own
1934  // any quadrant. This bypasses the assertion from being triggered.
1935  // - see: https://github.com/cburstedde/p4est/issues/48
1936  if (this->data_transfer.src_sizes_variable.size() == 0)
1937  this->data_transfer.src_sizes_variable.resize(1);
1938  if (this->data_transfer.dest_sizes_variable.size() == 0)
1939  this->data_transfer.dest_sizes_variable.resize(1);
1940 # endif
1941 
1942  // Execute variable size transfer.
1944  parallel_forest->global_first_quadrant,
1945  previous_global_first_quadrant,
1946  parallel_forest->mpicomm,
1947  1,
1948  this->data_transfer.dest_data_variable.data(),
1949  this->data_transfer.dest_sizes_variable.data(),
1950  this->data_transfer.src_data_variable.data(),
1951  this->data_transfer.src_sizes_variable.data());
1952 
1953  // Release memory of previously packed data.
1954  this->data_transfer.src_sizes_variable.clear();
1955  this->data_transfer.src_sizes_variable.shrink_to_fit();
1956  this->data_transfer.src_data_variable.clear();
1957  this->data_transfer.src_data_variable.shrink_to_fit();
1958  }
1959  }
1960 
1961 
1962 
1963  template <int dim, int spacedim>
1964  void
1966  {
1967  DynamicSparsityPattern cell_connectivity;
1969  cell_connectivity);
1970  coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
1972  cell_connectivity, coarse_cell_to_p4est_tree_permutation);
1973 
1974  p4est_tree_to_coarse_cell_permutation =
1975  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
1976  }
1977 
1978 
1979 
1980  template <int dim, int spacedim>
1981  void
1983  const std::string &file_basename) const
1984  {
1985  Assert(parallel_forest != nullptr,
1986  ExcMessage("Can't produce output when no forest is created yet."));
1987 
1988  AssertThrow(are_vertices_communicated_to_p4est(),
1989  ExcMessage(
1990  "To use this function the triangulation's flag "
1991  "Settings::communicate_vertices_to_p4est must be set."));
1992 
1994  parallel_forest, nullptr, file_basename.c_str());
1995  }
1996 
1997 
1998 
1999  template <int dim, int spacedim>
2000  void
2001  Triangulation<dim, spacedim>::save(const std::string &filename) const
2002  {
2003  Assert(
2004  this->cell_attached_data.n_attached_deserialize == 0,
2005  ExcMessage(
2006  "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2007  Assert(this->n_cells() > 0,
2008  ExcMessage("Can not save() an empty Triangulation."));
2009 
2010  const int myrank =
2011  Utilities::MPI::this_mpi_process(this->mpi_communicator);
2012 
2013  // signal that serialization is going to happen
2014  this->signals.pre_distributed_save();
2015 
2016  if (this->my_subdomain == 0)
2017  {
2018  std::string fname = std::string(filename) + ".info";
2019  std::ofstream f(fname.c_str());
2020  f << "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2021  << std::endl
2022  << 5 << " "
2023  << Utilities::MPI::n_mpi_processes(this->mpi_communicator) << " "
2024  << this->cell_attached_data.pack_callbacks_fixed.size() << " "
2025  << this->cell_attached_data.pack_callbacks_variable.size() << " "
2026  << this->n_cells(0) << std::endl;
2027  }
2028 
2029  // each cell should have been flagged `CELL_PERSIST`
2030  for (const auto &cell_rel : this->local_cell_relations)
2031  {
2032  (void)cell_rel;
2033  Assert(
2034  (cell_rel.second == // cell_status
2036  ExcInternalError());
2037  }
2038 
2039  // Save cell attached data.
2040  this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2041  parallel_forest->global_num_quadrants,
2042  filename);
2043 
2044  ::internal::p4est::functions<dim>::save(filename.c_str(),
2045  parallel_forest,
2046  false);
2047 
2048  // signal that serialization has finished
2049  this->signals.post_distributed_save();
2050  }
2051 
2052 
2053 
2054  template <int dim, int spacedim>
2055  void
2056  Triangulation<dim, spacedim>::load(const std::string &filename)
2057  {
2058  Assert(
2059  this->n_cells() > 0,
2060  ExcMessage(
2061  "load() only works if the Triangulation already contains a coarse mesh!"));
2062  Assert(
2063  this->n_levels() == 1,
2064  ExcMessage(
2065  "Triangulation may only contain coarse cells when calling load()."));
2066 
2067  const int myrank =
2068  Utilities::MPI::this_mpi_process(this->mpi_communicator);
2069 
2070  // signal that de-serialization is going to happen
2071  this->signals.pre_distributed_load();
2072 
2073  if (parallel_ghost != nullptr)
2074  {
2076  parallel_ghost);
2077  parallel_ghost = nullptr;
2078  }
2080  parallel_forest = nullptr;
2082  connectivity);
2083  connectivity = nullptr;
2084 
2085  unsigned int version, numcpus, attached_count_fixed,
2086  attached_count_variable, n_coarse_cells;
2087  {
2088  std::string fname = std::string(filename) + ".info";
2089  std::ifstream f(fname.c_str());
2090  AssertThrow(f.fail() == false, ExcIO());
2091  std::string firstline;
2092  getline(f, firstline); // skip first line
2093  f >> version >> numcpus >> attached_count_fixed >>
2094  attached_count_variable >> n_coarse_cells;
2095  }
2096 
2097  AssertThrow(version == 5,
2098  ExcMessage("Incompatible version found in .info file."));
2099  Assert(this->n_cells(0) == n_coarse_cells,
2100  ExcMessage("Number of coarse cells differ!"));
2101 
2102  // clear all of the callback data, as explained in the documentation of
2103  // register_data_attach()
2104  this->cell_attached_data.n_attached_data_sets = 0;
2105  this->cell_attached_data.n_attached_deserialize =
2106  attached_count_fixed + attached_count_variable;
2107 
2109  filename.c_str(),
2110  this->mpi_communicator,
2111  0,
2112  0,
2113  1,
2114  0,
2115  this,
2116  &connectivity);
2117 
2118  // We partition the p4est mesh that it conforms to the requirements of the
2119  // deal.II mesh, i.e., partition for coarsening.
2120  // This function call is optional.
2122  parallel_forest,
2123  /* prepare coarsening */ 1,
2124  /* weight_callback */ nullptr);
2125 
2126  try
2127  {
2128  copy_local_forest_to_triangulation();
2129  }
2130  catch (const typename Triangulation<dim>::DistortedCellList &)
2131  {
2132  // the underlying triangulation should not be checking for distorted
2133  // cells
2134  Assert(false, ExcInternalError());
2135  }
2136 
2137  // Load attached cell data, if any was stored.
2138  this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2139  parallel_forest->global_num_quadrants,
2140  parallel_forest->local_num_quadrants,
2141  filename,
2142  attached_count_fixed,
2143  attached_count_variable);
2144 
2145  // signal that de-serialization is finished
2146  this->signals.post_distributed_load();
2147 
2148  this->update_periodic_face_map();
2149  this->update_number_cache();
2150  }
2151 
2152 
2153 
2154  template <int dim, int spacedim>
2155  void
2156  Triangulation<dim, spacedim>::load(const std::string &filename,
2157  const bool autopartition)
2158  {
2159  (void)autopartition;
2160  load(filename);
2161  }
2162 
2163 
2164 
2165  template <int dim, int spacedim>
2166  void
2168  const typename ::internal::p4est::types<dim>::forest *forest)
2169  {
2170  Assert(this->n_cells() > 0,
2171  ExcMessage(
2172  "load() only works if the Triangulation already contains "
2173  "a coarse mesh!"));
2174  Assert(this->n_cells() == forest->trees->elem_count,
2175  ExcMessage(
2176  "Coarse mesh of the Triangulation does not match the one "
2177  "of the provided forest!"));
2178 
2179  // clear the old forest
2180  if (parallel_ghost != nullptr)
2181  {
2183  parallel_ghost);
2184  parallel_ghost = nullptr;
2185  }
2187  parallel_forest = nullptr;
2188 
2189  // note: we can keep the connectivity, since the coarse grid does not
2190  // change
2191 
2192  // create deep copy of the new forest
2193  typename ::internal::p4est::types<dim>::forest *temp =
2194  const_cast<typename ::internal::p4est::types<dim>::forest *>(
2195  forest);
2196  parallel_forest =
2198  parallel_forest->connectivity = connectivity;
2199  parallel_forest->user_pointer = this;
2200 
2201  try
2202  {
2203  copy_local_forest_to_triangulation();
2204  }
2205  catch (const typename Triangulation<dim>::DistortedCellList &)
2206  {
2207  // the underlying triangulation should not be checking for distorted
2208  // cells
2209  Assert(false, ExcInternalError());
2210  }
2211 
2212  this->update_periodic_face_map();
2213  this->update_number_cache();
2214  }
2215 
2216 
2217 
2218  template <int dim, int spacedim>
2219  unsigned int
2221  {
2222  Assert(parallel_forest != nullptr,
2223  ExcMessage(
2224  "Can't produce a check sum when no forest is created yet."));
2225  return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2226  }
2227 
2228 
2229 
2230  template <int dim, int spacedim>
2231  const typename ::internal::p4est::types<dim>::forest *
2233  {
2234  Assert(parallel_forest != nullptr,
2235  ExcMessage("The forest has not been allocated yet."));
2236  return parallel_forest;
2237  }
2238 
2239 
2240 
2241  template <int dim, int spacedim>
2242  typename ::internal::p4est::types<dim>::tree *
2244  const int dealii_coarse_cell_index) const
2245  {
2246  const unsigned int tree_index =
2247  coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2248  typename ::internal::p4est::types<dim>::tree *tree =
2249  static_cast<typename ::internal::p4est::types<dim>::tree *>(
2250  sc_array_index(parallel_forest->trees, tree_index));
2251 
2252  return tree;
2253  }
2254 
2255 
2256 
2257  // Note: this has been added here to prevent that these functions
2258  // appear in the Doxygen documentation of ::Triangulation
2259 # ifndef DOXYGEN
2260 
2261  template <>
2263  std::integral_constant<int, 2>)
2264  {
2265  const unsigned int dim = 2, spacedim = 2;
2266  Assert(this->n_cells(0) > 0, ExcInternalError());
2267  Assert(this->n_levels() == 1, ExcInternalError());
2268 
2269  // data structures that counts how many cells touch each vertex
2270  // (vertex_touch_count), and which cells touch a given vertex (together
2271  // with the local numbering of that vertex within the cells that touch
2272  // it)
2273  std::vector<unsigned int> vertex_touch_count;
2274  std::vector<
2275  std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2276  unsigned int>>>
2277  vertex_to_cell;
2278  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
2279  const ::internal::p4est::types<2>::locidx num_vtt =
2280  std::accumulate(vertex_touch_count.begin(),
2281  vertex_touch_count.end(),
2282  0u);
2283 
2284  // now create a connectivity object with the right sizes for all
2285  // arrays. set vertex information only in debug mode (saves a few bytes
2286  // in optimized mode)
2287  const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2288 
2290  (set_vertex_info == true ? this->n_vertices() : 0),
2291  this->n_cells(0),
2292  this->n_vertices(),
2293  num_vtt);
2294 
2295  set_vertex_and_cell_info(*this,
2296  vertex_touch_count,
2297  vertex_to_cell,
2298  coarse_cell_to_p4est_tree_permutation,
2299  set_vertex_info,
2300  connectivity);
2301 
2302  Assert(p4est_connectivity_is_valid(connectivity) == 1,
2303  ExcInternalError());
2304 
2305  // now create a forest out of the connectivity data structure
2307  this->mpi_communicator,
2308  connectivity,
2309  /* minimum initial number of quadrants per tree */ 0,
2310  /* minimum level of upfront refinement */ 0,
2311  /* use uniform upfront refinement */ 1,
2312  /* user_data_size = */ 0,
2313  /* user_data_constructor = */ nullptr,
2314  /* user_pointer */ this);
2315  }
2316 
2317 
2318 
2319  // TODO: This is a verbatim copy of the 2,2 case. However, we can't just
2320  // specialize the dim template argument, but let spacedim open
2321  template <>
2323  std::integral_constant<int, 2>)
2324  {
2325  const unsigned int dim = 2, spacedim = 3;
2326  Assert(this->n_cells(0) > 0, ExcInternalError());
2327  Assert(this->n_levels() == 1, ExcInternalError());
2328 
2329  // data structures that counts how many cells touch each vertex
2330  // (vertex_touch_count), and which cells touch a given vertex (together
2331  // with the local numbering of that vertex within the cells that touch
2332  // it)
2333  std::vector<unsigned int> vertex_touch_count;
2334  std::vector<
2335  std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2336  unsigned int>>>
2337  vertex_to_cell;
2338  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
2339  const ::internal::p4est::types<2>::locidx num_vtt =
2340  std::accumulate(vertex_touch_count.begin(),
2341  vertex_touch_count.end(),
2342  0u);
2343 
2344  // now create a connectivity object with the right sizes for all
2345  // arrays. set vertex information only in debug mode (saves a few bytes
2346  // in optimized mode)
2347  const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2348 
2350  (set_vertex_info == true ? this->n_vertices() : 0),
2351  this->n_cells(0),
2352  this->n_vertices(),
2353  num_vtt);
2354 
2355  set_vertex_and_cell_info(*this,
2356  vertex_touch_count,
2357  vertex_to_cell,
2358  coarse_cell_to_p4est_tree_permutation,
2359  set_vertex_info,
2360  connectivity);
2361 
2362  Assert(p4est_connectivity_is_valid(connectivity) == 1,
2363  ExcInternalError());
2364 
2365  // now create a forest out of the connectivity data structure
2367  this->mpi_communicator,
2368  connectivity,
2369  /* minimum initial number of quadrants per tree */ 0,
2370  /* minimum level of upfront refinement */ 0,
2371  /* use uniform upfront refinement */ 1,
2372  /* user_data_size = */ 0,
2373  /* user_data_constructor = */ nullptr,
2374  /* user_pointer */ this);
2375  }
2376 
2377 
2378 
2379  template <>
2381  std::integral_constant<int, 3>)
2382  {
2383  const int dim = 3, spacedim = 3;
2384  Assert(this->n_cells(0) > 0, ExcInternalError());
2385  Assert(this->n_levels() == 1, ExcInternalError());
2386 
2387  // data structures that counts how many cells touch each vertex
2388  // (vertex_touch_count), and which cells touch a given vertex (together
2389  // with the local numbering of that vertex within the cells that touch
2390  // it)
2391  std::vector<unsigned int> vertex_touch_count;
2392  std::vector<std::list<
2393  std::pair<Triangulation<3>::active_cell_iterator, unsigned int>>>
2394  vertex_to_cell;
2395  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
2396  const ::internal::p4est::types<2>::locidx num_vtt =
2397  std::accumulate(vertex_touch_count.begin(),
2398  vertex_touch_count.end(),
2399  0u);
2400 
2401  std::vector<unsigned int> edge_touch_count;
2402  std::vector<std::list<
2403  std::pair<Triangulation<3>::active_cell_iterator, unsigned int>>>
2404  edge_to_cell;
2405  get_edge_to_cell_mappings(*this, edge_touch_count, edge_to_cell);
2406  const ::internal::p4est::types<2>::locidx num_ett =
2407  std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2408 
2409  // now create a connectivity object with the right sizes for all arrays
2410  const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2411 
2413  (set_vertex_info == true ? this->n_vertices() : 0),
2414  this->n_cells(0),
2415  this->n_active_lines(),
2416  num_ett,
2417  this->n_vertices(),
2418  num_vtt);
2419 
2420  set_vertex_and_cell_info(*this,
2421  vertex_touch_count,
2422  vertex_to_cell,
2423  coarse_cell_to_p4est_tree_permutation,
2424  set_vertex_info,
2425  connectivity);
2426 
2427  // next to tree-to-edge
2428  // data. note that in p4est lines
2429  // are ordered as follows
2430  // *---3---* *---3---*
2431  // /| | / /|
2432  // 6 | 11 6 7 11
2433  // / 10 | / / |
2434  // * | | *---2---* |
2435  // | *---1---* | | *
2436  // | / / | 9 /
2437  // 8 4 5 8 | 5
2438  // |/ / | |/
2439  // *---0---* *---0---*
2440  // whereas in deal.II they are like this:
2441  // *---7---* *---7---*
2442  // /| | / /|
2443  // 4 | 11 4 5 11
2444  // / 10 | / / |
2445  // * | | *---6---* |
2446  // | *---3---* | | *
2447  // | / / | 9 /
2448  // 8 0 1 8 | 1
2449  // |/ / | |/
2450  // *---2---* *---2---*
2451 
2452  const unsigned int deal_to_p4est_line_index[12] = {
2453  4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2454 
2456  this->begin_active();
2457  cell != this->end();
2458  ++cell)
2459  {
2460  const unsigned int index =
2461  coarse_cell_to_p4est_tree_permutation[cell->index()];
2462  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++e)
2463  connectivity->tree_to_edge[index * GeometryInfo<3>::lines_per_cell +
2464  deal_to_p4est_line_index[e]] =
2465  cell->line(e)->index();
2466  }
2467 
2468  // now also set edge-to-tree
2469  // information
2470  connectivity->ett_offset[0] = 0;
2471  std::partial_sum(edge_touch_count.begin(),
2472  edge_touch_count.end(),
2473  &connectivity->ett_offset[1]);
2474 
2475  Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2476  ExcInternalError());
2477 
2478  for (unsigned int v = 0; v < this->n_active_lines(); ++v)
2479  {
2480  Assert(edge_to_cell[v].size() == edge_touch_count[v],
2481  ExcInternalError());
2482 
2483  std::list<
2484  std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2485  unsigned int>>::const_iterator p =
2486  edge_to_cell[v].begin();
2487  for (unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2488  {
2489  connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2490  coarse_cell_to_p4est_tree_permutation[p->first->index()];
2491  connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2492  deal_to_p4est_line_index[p->second];
2493  }
2494  }
2495 
2496  Assert(p8est_connectivity_is_valid(connectivity) == 1,
2497  ExcInternalError());
2498 
2499  // now create a forest out of the connectivity data structure
2501  this->mpi_communicator,
2502  connectivity,
2503  /* minimum initial number of quadrants per tree */ 0,
2504  /* minimum level of upfront refinement */ 0,
2505  /* use uniform upfront refinement */ 1,
2506  /* user_data_size = */ 0,
2507  /* user_data_constructor = */ nullptr,
2508  /* user_pointer */ this);
2509  }
2510 # endif
2511 
2512 
2513 
2514  namespace
2515  {
2516  // ensures the 2:1 mesh balance for periodic boundary conditions in the
2517  // artificial cell layer (the active cells are taken care of by p4est)
2518  template <int dim, int spacedim>
2519  bool
2520  enforce_mesh_balance_over_periodic_boundaries(
2522  {
2523  if (tria.get_periodic_face_map().size() == 0)
2524  return false;
2525 
2526  std::vector<bool> flags_before[2];
2527  tria.save_coarsen_flags(flags_before[0]);
2528  tria.save_refine_flags(flags_before[1]);
2529 
2530  std::vector<unsigned int> topological_vertex_numbering(
2531  tria.n_vertices());
2532  for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2533  topological_vertex_numbering[i] = i;
2534  // combine vertices that have different locations (and thus, different
2535  // vertex_index) but represent the same topological entity over
2536  // periodic boundaries. The vector topological_vertex_numbering
2537  // contains a linear map from 0 to n_vertices at input and at output
2538  // relates periodic vertices with only one vertex index. The output is
2539  // used to always identify the same vertex according to the
2540  // periodicity, e.g. when finding the maximum cell level around a
2541  // vertex.
2542  //
2543  // Example: On a 3D cell with vertices numbered from 0 to 7 and
2544  // periodic boundary conditions in x direction, the vector
2545  // topological_vertex_numbering will contain the numbers
2546  // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5},
2547  // {6,7} belong together, respectively). If periodicity is set in x
2548  // and z direction, the output is {0,0,2,2,0,0,2,2}, and if
2549  // periodicity is in all directions, the output is simply
2550  // {0,0,0,0,0,0,0,0}.
2551  using cell_iterator =
2553  typename std::map<std::pair<cell_iterator, unsigned int>,
2554  std::pair<std::pair<cell_iterator, unsigned int>,
2555  std::bitset<3>>>::const_iterator it;
2556  for (it = tria.get_periodic_face_map().begin();
2557  it != tria.get_periodic_face_map().end();
2558  ++it)
2559  {
2560  const cell_iterator &cell_1 = it->first.first;
2561  const unsigned int face_no_1 = it->first.second;
2562  const cell_iterator &cell_2 = it->second.first.first;
2563  const unsigned int face_no_2 = it->second.first.second;
2564  const std::bitset<3> face_orientation = it->second.second;
2565 
2566  if (cell_1->level() == cell_2->level())
2567  {
2568  for (unsigned int v = 0;
2569  v < GeometryInfo<dim - 1>::vertices_per_cell;
2570  ++v)
2571  {
2572  // take possible non-standard orientation of face on
2573  // cell[0] into account
2574  const unsigned int vface0 =
2576  v,
2577  face_orientation[0],
2578  face_orientation[1],
2579  face_orientation[2]);
2580  const unsigned int vi0 =
2581  topological_vertex_numbering[cell_1->face(face_no_1)
2582  ->vertex_index(vface0)];
2583  const unsigned int vi1 =
2584  topological_vertex_numbering[cell_2->face(face_no_2)
2585  ->vertex_index(v)];
2586  const unsigned int min_index = std::min(vi0, vi1);
2587  topological_vertex_numbering[cell_1->face(face_no_1)
2588  ->vertex_index(vface0)] =
2589  topological_vertex_numbering[cell_2->face(face_no_2)
2590  ->vertex_index(v)] =
2591  min_index;
2592  }
2593  }
2594  }
2595 
2596  // There must not be any chains!
2597  for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2598  {
2599  const unsigned int j = topological_vertex_numbering[i];
2600  if (j != i)
2601  Assert(topological_vertex_numbering[j] == j, ExcInternalError());
2602  }
2603 
2604 
2605  // this code is replicated from grid/tria.cc but using an indirection
2606  // for periodic boundary conditions
2607  bool continue_iterating = true;
2608  std::vector<int> vertex_level(tria.n_vertices());
2609  while (continue_iterating)
2610  {
2611  // store highest level one of the cells adjacent to a vertex
2612  // belongs to
2613  std::fill(vertex_level.begin(), vertex_level.end(), 0);
2615  cell = tria.begin_active(),
2616  endc = tria.end();
2617  for (; cell != endc; ++cell)
2618  {
2619  if (cell->refine_flag_set())
2620  for (const unsigned int vertex :
2622  vertex_level[topological_vertex_numbering
2623  [cell->vertex_index(vertex)]] =
2624  std::max(vertex_level[topological_vertex_numbering
2625  [cell->vertex_index(vertex)]],
2626  cell->level() + 1);
2627  else if (!cell->coarsen_flag_set())
2628  for (const unsigned int vertex :
2630  vertex_level[topological_vertex_numbering
2631  [cell->vertex_index(vertex)]] =
2632  std::max(vertex_level[topological_vertex_numbering
2633  [cell->vertex_index(vertex)]],
2634  cell->level());
2635  else
2636  {
2637  // if coarsen flag is set then tentatively assume
2638  // that the cell will be coarsened. this isn't
2639  // always true (the coarsen flag could be removed
2640  // again) and so we may make an error here. we try
2641  // to correct this by iterating over the entire
2642  // process until we are converged
2643  Assert(cell->coarsen_flag_set(), ExcInternalError());
2644  for (const unsigned int vertex :
2646  vertex_level[topological_vertex_numbering
2647  [cell->vertex_index(vertex)]] =
2648  std::max(vertex_level[topological_vertex_numbering
2649  [cell->vertex_index(vertex)]],
2650  cell->level() - 1);
2651  }
2652  }
2653 
2654  continue_iterating = false;
2655 
2656  // loop over all cells in reverse order. do so because we
2657  // can then update the vertex levels on the adjacent
2658  // vertices and maybe already flag additional cells in this
2659  // loop
2660  //
2661  // note that not only may we have to add additional
2662  // refinement flags, but we will also have to remove
2663  // coarsening flags on cells adjacent to vertices that will
2664  // see refinement
2665  for (cell = tria.last_active(); cell != endc; --cell)
2666  if (cell->refine_flag_set() == false)
2667  {
2668  for (const unsigned int vertex :
2670  if (vertex_level[topological_vertex_numbering
2671  [cell->vertex_index(vertex)]] >=
2672  cell->level() + 1)
2673  {
2674  // remove coarsen flag...
2675  cell->clear_coarsen_flag();
2676 
2677  // ...and if necessary also refine the current
2678  // cell, at the same time updating the level
2679  // information about vertices
2680  if (vertex_level[topological_vertex_numbering
2681  [cell->vertex_index(vertex)]] >
2682  cell->level() + 1)
2683  {
2684  cell->set_refine_flag();
2685  continue_iterating = true;
2686 
2687  for (const unsigned int v :
2689  vertex_level[topological_vertex_numbering
2690  [cell->vertex_index(v)]] =
2691  std::max(
2692  vertex_level[topological_vertex_numbering
2693  [cell->vertex_index(v)]],
2694  cell->level() + 1);
2695  }
2696 
2697  // continue and see whether we may, for example,
2698  // go into the inner 'if' above based on a
2699  // different vertex
2700  }
2701  }
2702 
2703  // clear coarsen flag if not all children were marked
2704  for (const auto &cell : tria.cell_iterators())
2705  {
2706  // nothing to do if we are already on the finest level
2707  if (cell->is_active())
2708  continue;
2709 
2710  const unsigned int n_children = cell->n_children();
2711  unsigned int flagged_children = 0;
2712  for (unsigned int child = 0; child < n_children; ++child)
2713  if (cell->child(child)->is_active() &&
2714  cell->child(child)->coarsen_flag_set())
2715  ++flagged_children;
2716 
2717  // if not all children were flagged for coarsening, remove
2718  // coarsen flags
2719  if (flagged_children < n_children)
2720  for (unsigned int child = 0; child < n_children; ++child)
2721  if (cell->child(child)->is_active())
2722  cell->child(child)->clear_coarsen_flag();
2723  }
2724  }
2725  std::vector<bool> flags_after[2];
2726  tria.save_coarsen_flags(flags_after[0]);
2727  tria.save_refine_flags(flags_after[1]);
2728  return ((flags_before[0] != flags_after[0]) ||
2729  (flags_before[1] != flags_after[1]));
2730  }
2731  } // namespace
2732 
2733 
2734 
2735  template <int dim, int spacedim>
2736  bool
2738  {
2739  bool mesh_changed = false;
2740  unsigned int loop_counter = 0;
2741  unsigned int n_changes = 0;
2742  do
2743  {
2744  n_changes += this->::Triangulation<dim, spacedim>::
2746  this->update_periodic_face_map();
2747  // enforce 2:1 mesh balance over periodic boundaries
2748  mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this);
2749  n_changes += mesh_changed;
2750 
2751  // We can't be sure that we won't run into a situation where we can
2752  // not reconcile mesh smoothing and balancing of periodic faces. As
2753  // we don't know what else to do, at least abort with an error
2754  // message.
2755  ++loop_counter;
2756  AssertThrow(
2757  loop_counter < 32,
2758  ExcMessage(
2759  "Infinite loop in "
2760  "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2761  "for periodic boundaries detected. Aborting."));
2762  }
2763  while (mesh_changed);
2764 
2765  // report if we observed changes in any of the sub-functions
2766  return n_changes > 0;
2767  }
2768 
2769 
2770 
2771  template <int dim, int spacedim>
2772  void
2774  {
2775  // Disable mesh smoothing for recreating the deal.II triangulation,
2776  // otherwise we might not be able to reproduce the p4est mesh
2777  // exactly. We restore the original smoothing at the end of this
2778  // function. Note that the smoothing flag is used in the normal
2779  // refinement process.
2780  typename Triangulation<dim, spacedim>::MeshSmoothing save_smooth =
2781  this->smooth_grid;
2782 
2783  // We will refine manually to match the p4est further down, which
2784  // obeys a level difference of 2 at each vertex (see the balance call
2785  // to p4est). We can disable this here so we store fewer artificial
2786  // cells (in some cases).
2787  // For geometric multigrid it turns out that
2788  // we will miss level cells at shared vertices if we ignore this.
2789  // See tests/mpi/mg_06. In particular, the flag is still necessary
2790  // even though we force it for the original smooth_grid in the
2791  // constructor.
2793  this->smooth_grid =
2794  ::Triangulation<dim,
2795  spacedim>::limit_level_difference_at_vertices;
2796  else
2797  this->smooth_grid = ::Triangulation<dim, spacedim>::none;
2798 
2799  bool mesh_changed = false;
2800 
2801  // Remove all deal.II refinements. Note that we could skip this and
2802  // start from our current state, because the algorithm later coarsens as
2803  // necessary. This has the advantage of being faster when large parts
2804  // of the local partition changes (likely) and gives a deterministic
2805  // ordering of the cells (useful for snapshot/resume).
2806  // TODO: is there a more efficient way to do this?
2807  if (settings & mesh_reconstruction_after_repartitioning)
2808  while (this->n_levels() > 1)
2809  {
2810  // Instead of marking all active cells, we slice off the finest
2811  // level, one level at a time. This takes the same number of
2812  // iterations but solves an issue where not all cells on a
2813  // periodic boundary are indeed coarsened and we run into an
2814  // irrelevant Assert() in update_periodic_face_map().
2815  for (const auto &cell :
2816  this->active_cell_iterators_on_level(this->n_levels() - 1))
2817  {
2818  cell->set_coarsen_flag();
2819  }
2820  try
2821  {
2824  }
2825  catch (
2827  {
2828  // the underlying triangulation should not be checking for
2829  // distorted cells
2830  Assert(false, ExcInternalError());
2831  }
2832  }
2833 
2834 
2835  // query p4est for the ghost cells
2836  if (parallel_ghost != nullptr)
2837  {
2839  parallel_ghost);
2840  parallel_ghost = nullptr;
2841  }
2843  parallel_forest,
2844  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2845  P4EST_CONNECT_CORNER) :
2846  typename ::internal::p4est::types<dim>::balance_type(
2847  P8EST_CONNECT_CORNER)));
2848 
2849  Assert(parallel_ghost, ExcInternalError());
2850 
2851 
2852  // set all cells to artificial. we will later set it to the correct
2853  // subdomain in match_tree_recursively
2854  for (const auto &cell : this->cell_iterators_on_level(0))
2855  cell->recursively_set_subdomain_id(numbers::artificial_subdomain_id);
2856 
2857  do
2858  {
2859  for (const auto &cell : this->cell_iterators_on_level(0))
2860  {
2861  // if this processor stores no part of the forest that comes out
2862  // of this coarse grid cell, then we need to delete all children
2863  // of this cell (the coarse grid cell remains)
2864  if (tree_exists_locally<dim, spacedim>(
2865  parallel_forest,
2866  coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2867  false)
2868  {
2869  delete_all_children<dim, spacedim>(cell);
2870  if (cell->is_active())
2871  cell->set_subdomain_id(numbers::artificial_subdomain_id);
2872  }
2873 
2874  else
2875  {
2876  // this processor stores at least a part of the tree that
2877  // comes out of this cell.
2878 
2879  typename ::internal::p4est::types<dim>::quadrant
2880  p4est_coarse_cell;
2881  typename ::internal::p4est::types<dim>::tree *tree =
2882  init_tree(cell->index());
2883 
2884  ::internal::p4est::init_coarse_quadrant<dim>(
2885  p4est_coarse_cell);
2886 
2887  match_tree_recursively<dim, spacedim>(*tree,
2888  cell,
2889  p4est_coarse_cell,
2890  *parallel_forest,
2891  this->my_subdomain);
2892  }
2893  }
2894 
2895  // check mesh for ghost cells, refine as necessary. iterate over
2896  // every ghostquadrant, find corresponding deal coarsecell and
2897  // recurse.
2898  typename ::internal::p4est::types<dim>::quadrant *quadr;
2899  types::subdomain_id ghost_owner = 0;
2900  typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2901 
2902  for (unsigned int g_idx = 0;
2903  g_idx < parallel_ghost->ghosts.elem_count;
2904  ++g_idx)
2905  {
2906  while (g_idx >= static_cast<unsigned int>(
2907  parallel_ghost->proc_offsets[ghost_owner + 1]))
2908  ++ghost_owner;
2909  while (g_idx >= static_cast<unsigned int>(
2910  parallel_ghost->tree_offsets[ghost_tree + 1]))
2911  ++ghost_tree;
2912 
2913  quadr = static_cast<
2914  typename ::internal::p4est::types<dim>::quadrant *>(
2915  sc_array_index(&parallel_ghost->ghosts, g_idx));
2916 
2917  unsigned int coarse_cell_index =
2918  p4est_tree_to_coarse_cell_permutation[ghost_tree];
2919 
2920  match_quadrant<dim, spacedim>(this,
2921  coarse_cell_index,
2922  *quadr,
2923  ghost_owner);
2924  }
2925 
2926  // fix all the flags to make sure we have a consistent mesh
2927  this->prepare_coarsening_and_refinement();
2928 
2929  // see if any flags are still set
2930  mesh_changed =
2931  std::any_of(this->begin_active(),
2932  active_cell_iterator{this->end()},
2933  [](const CellAccessor<dim, spacedim> &cell) {
2934  return cell.refine_flag_set() ||
2935  cell.coarsen_flag_set();
2936  });
2937 
2938  // actually do the refinement to change the local mesh by
2939  // calling the base class refinement function directly
2940  try
2941  {
2944  }
2945  catch (
2947  {
2948  // the underlying triangulation should not be checking for
2949  // distorted cells
2950  Assert(false, ExcInternalError());
2951  }
2952  }
2953  while (mesh_changed);
2954 
2955 # ifdef DEBUG
2956  // check if correct number of ghosts is created
2957  unsigned int num_ghosts = 0;
2958 
2959  for (const auto &cell : this->active_cell_iterators())
2960  {
2961  if (cell->subdomain_id() != this->my_subdomain &&
2962  cell->subdomain_id() != numbers::artificial_subdomain_id)
2963  ++num_ghosts;
2964  }
2965 
2966  Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2967  ExcInternalError());
2968 # endif
2969 
2970 
2971 
2972  // fill level_subdomain_ids for geometric multigrid
2973  // the level ownership of a cell is defined as the owner if the cell is
2974  // active or as the owner of child(0) we need this information for all
2975  // our ancestors and the same-level neighbors of our own cells (=level
2976  // ghosts)
2978  {
2979  // step 1: We set our own ids all the way down and all the others to
2980  // -1. Note that we do not fill other cells we could figure out the
2981  // same way, because we might accidentally set an id for a cell that
2982  // is not a ghost cell.
2983  for (unsigned int lvl = this->n_levels(); lvl > 0;)
2984  {
2985  --lvl;
2986  for (const auto &cell : this->cell_iterators_on_level(lvl))
2987  {
2988  if ((cell->is_active() &&
2989  cell->subdomain_id() ==
2990  this->locally_owned_subdomain()) ||
2991  (cell->has_children() &&
2992  cell->child(0)->level_subdomain_id() ==
2993  this->locally_owned_subdomain()))
2994  cell->set_level_subdomain_id(
2995  this->locally_owned_subdomain());
2996  else
2997  {
2998  // not our cell
2999  cell->set_level_subdomain_id(
3001  }
3002  }
3003  }
3004 
3005  // step 2: make sure all the neighbors to our level_cells exist.
3006  // Need to look up in p4est...
3007  std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3008  for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3009  marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3010 
3011  for (const auto &cell : this->cell_iterators_on_level(0))
3012  {
3013  typename ::internal::p4est::types<dim>::quadrant
3014  p4est_coarse_cell;
3015  const unsigned int tree_index =
3016  coarse_cell_to_p4est_tree_permutation[cell->index()];
3017  typename ::internal::p4est::types<dim>::tree *tree =
3018  init_tree(cell->index());
3019 
3020  ::internal::p4est::init_coarse_quadrant<dim>(
3021  p4est_coarse_cell);
3022 
3023  determine_level_subdomain_id_recursively<dim, spacedim>(
3024  *tree,
3025  tree_index,
3026  cell,
3027  p4est_coarse_cell,
3028  *parallel_forest,
3029  this->my_subdomain,
3030  marked_vertices);
3031  }
3032 
3033  // step 3: make sure we have the parent of our level cells
3034  for (unsigned int lvl = this->n_levels(); lvl > 0;)
3035  {
3036  --lvl;
3037  for (const auto &cell : this->cell_iterators_on_level(lvl))
3038  {
3039  if (cell->has_children())
3040  for (unsigned int c = 0;
3041  c < GeometryInfo<dim>::max_children_per_cell;
3042  ++c)
3043  {
3044  if (cell->child(c)->level_subdomain_id() ==
3045  this->locally_owned_subdomain())
3046  {
3047  // at least one of the children belongs to us, so
3048  // make sure we set the level subdomain id
3049  const types::subdomain_id mark =
3050  cell->child(0)->level_subdomain_id();
3052  ExcInternalError()); // we should know the
3053  // child(0)
3054  cell->set_level_subdomain_id(mark);
3055  break;
3056  }
3057  }
3058  }
3059  }
3060  }
3061 
3062 
3063 
3064  // check that our local copy has exactly as many cells as the p4est
3065  // original (at least if we are on only one processor); for parallel
3066  // computations, we want to check that we have at least as many as p4est
3067  // stores locally (in the future we should check that we have exactly as
3068  // many non-artificial cells as parallel_forest->local_num_quadrants)
3069  {
3070  const unsigned int total_local_cells = this->n_active_cells();
3071  (void)total_local_cells;
3072 
3073  if (Utilities::MPI::n_mpi_processes(this->mpi_communicator) == 1)
3074  {
3075  Assert(static_cast<unsigned int>(
3076  parallel_forest->local_num_quadrants) == total_local_cells,
3077  ExcInternalError());
3078  }
3079  else
3080  {
3081  Assert(static_cast<unsigned int>(
3082  parallel_forest->local_num_quadrants) <= total_local_cells,
3083  ExcInternalError());
3084  }
3085 
3086  // count the number of owned, active cells and compare with p4est.
3087  unsigned int n_owned = 0;
3088  for (const auto &cell : this->active_cell_iterators())
3089  {
3090  if (cell->subdomain_id() == this->my_subdomain)
3091  ++n_owned;
3092  }
3093 
3094  Assert(static_cast<unsigned int>(
3095  parallel_forest->local_num_quadrants) == n_owned,
3096  ExcInternalError());
3097  }
3098 
3099  this->smooth_grid = save_smooth;
3100 
3101  // finally, after syncing the parallel_forest with the triangulation,
3102  // also update the cell_relations, which will be used for
3103  // repartitioning, further refinement/coarsening, and unpacking
3104  // of stored or transferred data.
3105  update_cell_relations();
3106  }
3107 
3108 
3109 
3110  template <int dim, int spacedim>
3113  {
3114  // Call the other function
3115  std::vector<Point<dim>> point{p};
3116  std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
3117 
3118  return owner[0];
3119  }
3120 
3121 
3122 
3123  template <int dim, int spacedim>
3124  std::vector<types::subdomain_id>
3126  const std::vector<Point<dim>> &points)
3127  {
3128 # ifndef P4EST_SEARCH_LOCAL
3129  (void)points;
3130  AssertThrow(
3131  false,
3132  ExcMessage(
3133  "This function is only available if p4est is version 2.2 and higher."));
3134  // Just return to satisfy compiler
3135  return std::vector<unsigned int>(1,
3137 # else
3138  // We can only use this function if vertices are communicated to p4est
3139  AssertThrow(this->are_vertices_communicated_to_p4est(),
3140  ExcMessage(
3141  "Vertices need to be communicated to p4est to use this "
3142  "function. This must explicitly be turned on in the "
3143  "settings of the triangulation's constructor."));
3144 
3145  // We can only use this function if all manifolds are flat
3146  for (const auto &manifold_id : this->get_manifold_ids())
3147  {
3148  AssertThrow(
3150  ExcMessage(
3151  "This function can only be used if the triangulation "
3152  "has no other manifold than a Cartesian (flat) manifold attached."));
3153  }
3154 
3155  // Create object for callback
3156  PartitionSearch<dim> partition_search;
3157 
3158  // Pointer should be this triangulation before we set it to something else
3159  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3160 
3161  // re-assign p4est's user pointer
3162  parallel_forest->user_pointer = &partition_search;
3163 
3164  //
3165  // Copy points into p4est internal array data struct
3166  //
3167  // pointer to an array of points.
3168  sc_array_t *point_sc_array;
3169  // allocate memory for a number of dim-dimensional points including their
3170  // MPI rank, i.e., dim + 1 fields
3171  point_sc_array =
3172  sc_array_new_count(sizeof(double[dim + 1]), points.size());
3173 
3174  // now assign the actual value
3175  for (size_t i = 0; i < points.size(); ++i)
3176  {
3177  // alias
3178  const Point<dim> &p = points[i];
3179  // get a non-const view of the array
3180  double *this_sc_point =
3181  static_cast<double *>(sc_array_index_ssize_t(point_sc_array, i));
3182  // fill this with the point data
3183  for (unsigned int d = 0; d < dim; ++d)
3184  {
3185  this_sc_point[d] = p(d);
3186  }
3187  this_sc_point[dim] = -1.0; // owner rank
3188  }
3189 
3191  parallel_forest,
3192  /* execute quadrant function when leaving quadrant */
3193  static_cast<int>(false),
3194  &PartitionSearch<dim>::local_quadrant_fn,
3195  &PartitionSearch<dim>::local_point_fn,
3196  point_sc_array);
3197 
3198  // copy the points found to an std::array
3199  std::vector<types::subdomain_id> owner_rank(
3200  points.size(), numbers::invalid_subdomain_id);
3201 
3202  // fill the array
3203  for (size_t i = 0; i < points.size(); ++i)
3204  {
3205  // get a non-const view of the array
3206  double *this_sc_point =
3207  static_cast<double *>(sc_array_index_ssize_t(point_sc_array, i));
3208  owner_rank[i] = static_cast<types::subdomain_id>(this_sc_point[dim]);
3209  }
3210 
3211  // reset the internal pointer to this triangulation
3212  parallel_forest->user_pointer = this;
3213 
3214  // release the memory (otherwise p4est will complain)
3215  sc_array_destroy_null(&point_sc_array);
3216 
3217  return owner_rank;
3218 # endif // P4EST_SEARCH_LOCAL defined
3219  }
3220 
3221 
3222 
3223  template <int dim, int spacedim>
3224  void
3226  {
3227  // do not allow anisotropic refinement
3228 # ifdef DEBUG
3229  for (const auto &cell : this->active_cell_iterators())
3230  if (cell->is_locally_owned() && cell->refine_flag_set())
3231  Assert(cell->refine_flag_set() ==
3233  ExcMessage(
3234  "This class does not support anisotropic refinement"));
3235 # endif
3236 
3237 
3238  // safety check: p4est has an upper limit on the level of a cell
3239  if (this->n_levels() ==
3241  {
3243  cell = this->begin_active(
3245  cell !=
3247  1);
3248  ++cell)
3249  {
3250  AssertThrow(
3251  !(cell->refine_flag_set()),
3252  ExcMessage(
3253  "Fatal Error: maximum refinement level of p4est reached."));
3254  }
3255  }
3256 
3257  this->prepare_coarsening_and_refinement();
3258 
3259  // signal that refinement is going to happen
3260  this->signals.pre_distributed_refinement();
3261 
3262  // now do the work we're supposed to do when we are in charge
3263  // make sure all flags are cleared on cells we don't own, since nothing
3264  // good can come of that if they are still around
3265  for (const auto &cell : this->active_cell_iterators())
3266  if (cell->is_ghost() || cell->is_artificial())
3267  {
3268  cell->clear_refine_flag();
3269  cell->clear_coarsen_flag();
3270  }
3271 
3272 
3273  // count how many cells will be refined and coarsened, and allocate that
3274  // much memory
3275  RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3276  *this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3277 
3278  // copy refine and coarsen flags into p4est and execute the refinement
3279  // and coarsening. this uses the refine_and_coarsen_list just built,
3280  // which is communicated to the callback functions through
3281  // p4est's user_pointer object
3282  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3283  parallel_forest->user_pointer = &refine_and_coarsen_list;
3284 
3285  if (parallel_ghost != nullptr)
3286  {
3288  parallel_ghost);
3289  parallel_ghost = nullptr;
3290  }
3292  parallel_forest,
3293  /* refine_recursive */ false,
3294  &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3295  /*init_callback=*/nullptr);
3297  parallel_forest,
3298  /* coarsen_recursive */ false,
3299  &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3300  /*init_callback=*/nullptr);
3301 
3302  // make sure all cells in the lists have been consumed
3303  Assert(refine_and_coarsen_list.pointers_are_at_end(), ExcInternalError());
3304 
3305  // reset the pointer
3306  parallel_forest->user_pointer = this;
3307 
3308  // enforce 2:1 hanging node condition
3310  parallel_forest,
3311  /* face and corner balance */
3312  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3313  P4EST_CONNECT_FULL) :
3314  typename ::internal::p4est::types<dim>::balance_type(
3315  P8EST_CONNECT_FULL)),
3316  /*init_callback=*/nullptr);
3317 
3318  // since refinement and/or coarsening on the parallel forest
3319  // has happened, we need to update the quadrant cell relations
3320  update_cell_relations();
3321 
3322  // signals that parallel_forest has been refined and cell relations have
3323  // been updated
3324  this->signals.post_p4est_refinement();
3325 
3326  // before repartitioning the mesh, save a copy of the current positions
3327  // of quadrants only if data needs to be transferred later
3328  std::vector<typename ::internal::p4est::types<dim>::gloidx>
3329  previous_global_first_quadrant;
3330 
3331  if (this->cell_attached_data.n_attached_data_sets > 0)
3332  {
3333  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3334  std::memcpy(previous_global_first_quadrant.data(),
3335  parallel_forest->global_first_quadrant,
3336  sizeof(
3337  typename ::internal::p4est::types<dim>::gloidx) *
3338  (parallel_forest->mpisize + 1));
3339  }
3340 
3341  if (!(settings & no_automatic_repartitioning))
3342  {
3343  // partition the new mesh between all processors. If cell weights
3344  // have not been given balance the number of cells.
3345  if (this->signals.weight.empty())
3347  parallel_forest,
3348  /* prepare coarsening */ 1,
3349  /* weight_callback */ nullptr);
3350  else
3351  {
3352  // get cell weights for a weighted repartitioning.
3353  const std::vector<unsigned int> cell_weights = get_cell_weights();
3354 
3355  // verify that the global sum of weights is larger than 0
3356  Assert(Utilities::MPI::sum(std::accumulate(cell_weights.begin(),
3357  cell_weights.end(),
3358  std::uint64_t(0)),
3359  this->mpi_communicator) > 0,
3360  ExcMessage(
3361  "The global sum of weights over all active cells "
3362  "is zero. Please verify how you generate weights."));
3363 
3364  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3365 
3366  // attach (temporarily) a pointer to the cell weights through
3367  // p4est's user_pointer object
3368  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3369  parallel_forest->user_pointer = &partition_weights;
3370 
3372  parallel_forest,
3373  /* prepare coarsening */ 1,
3374  /* weight_callback */
3375  &PartitionWeights<dim, spacedim>::cell_weight);
3376 
3377  // release data
3379  parallel_forest, 0, nullptr, nullptr);
3380  // reset the user pointer to its previous state
3381  parallel_forest->user_pointer = this;
3382  }
3383  }
3384 
3385  // pack data before triangulation gets updated
3386  if (this->cell_attached_data.n_attached_data_sets > 0)
3387  {
3388  this->data_transfer.pack_data(
3389  this->local_cell_relations,
3390  this->cell_attached_data.pack_callbacks_fixed,
3391  this->cell_attached_data.pack_callbacks_variable);
3392  }
3393 
3394  // finally copy back from local part of tree to deal.II
3395  // triangulation. before doing so, make sure there are no refine or
3396  // coarsen flags pending
3397  for (const auto &cell : this->active_cell_iterators())
3398  {
3399  cell->clear_refine_flag();
3400  cell->clear_coarsen_flag();
3401  }
3402 
3403  try
3404  {
3405  copy_local_forest_to_triangulation();
3406  }
3407  catch (const typename Triangulation<dim>::DistortedCellList &)
3408  {
3409  // the underlying triangulation should not be checking for distorted
3410  // cells
3411  Assert(false, ExcInternalError());
3412  }
3413 
3414  // transfer data after triangulation got updated
3415  if (this->cell_attached_data.n_attached_data_sets > 0)
3416  {
3417  this->execute_transfer(parallel_forest,
3418  previous_global_first_quadrant.data());
3419 
3420  // also update the CellStatus information on the new mesh
3421  this->data_transfer.unpack_cell_status(this->local_cell_relations);
3422  }
3423 
3424 # ifdef DEBUG
3425  // Check that we know the level subdomain ids of all our neighbors. This
3426  // also involves coarser cells that share a vertex if they are active.
3427  //
3428  // Example (M= my, O=other):
3429  // *------*
3430  // | |
3431  // | O |
3432  // | |
3433  // *---*---*------*
3434  // | M | M |
3435  // *---*---*
3436  // | | M |
3437  // *---*---*
3438  // ^- the parent can be owned by somebody else, so O is not a neighbor
3439  // one level coarser
3441  {
3442  for (unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3443  {
3444  std::vector<bool> active_verts =
3445  this->mark_locally_active_vertices_on_level(lvl);
3446 
3447  const unsigned int maybe_coarser_lvl =
3448  (lvl > 0) ? (lvl - 1) : lvl;
3450  cell = this->begin(maybe_coarser_lvl),
3451  endc = this->end(lvl);
3452  for (; cell != endc; ++cell)
3453  if (cell->level() == static_cast<int>(lvl) || cell->is_active())
3454  {
3455  const bool is_level_artificial =
3456  (cell->level_subdomain_id() ==
3458  bool need_to_know = false;
3459  for (const unsigned int vertex :
3461  if (active_verts[cell->vertex_index(vertex)])
3462  {
3463  need_to_know = true;
3464  break;
3465  }
3466 
3467  Assert(
3468  !need_to_know || !is_level_artificial,
3469  ExcMessage(
3470  "Internal error: the owner of cell" +
3471  cell->id().to_string() +
3472  " is unknown even though it is needed for geometric multigrid."));
3473  }
3474  }
3475  }
3476 # endif
3477 
3478  this->update_periodic_face_map();
3479  this->update_number_cache();
3480 
3481  // signal that refinement is finished
3482  this->signals.post_distributed_refinement();
3483  }
3484 
3485 
3486 
3487  template <int dim, int spacedim>
3488  void
3490  {
3491 # ifdef DEBUG
3492  for (const auto &cell : this->active_cell_iterators())
3493  if (cell->is_locally_owned())
3494  Assert(
3495  !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3496  ExcMessage(
3497  "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3498 # endif
3499 
3500  // signal that repartitioning is going to happen
3501  this->signals.pre_distributed_repartition();
3502 
3503  // before repartitioning the mesh, save a copy of the current positions
3504  // of quadrants only if data needs to be transferred later
3505  std::vector<typename ::internal::p4est::types<dim>::gloidx>
3506  previous_global_first_quadrant;
3507 
3508  if (this->cell_attached_data.n_attached_data_sets > 0)
3509  {
3510  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3511  std::memcpy(previous_global_first_quadrant.data(),
3512  parallel_forest->global_first_quadrant,
3513  sizeof(
3514  typename ::internal::p4est::types<dim>::gloidx) *
3515  (parallel_forest->mpisize + 1));
3516  }
3517 
3518  if (this->signals.weight.empty())
3519  {
3520  // no cell weights given -- call p4est's 'partition' without a
3521  // callback for cell weights
3523  parallel_forest,
3524  /* prepare coarsening */ 1,
3525  /* weight_callback */ nullptr);
3526  }
3527  else
3528  {
3529  // get cell weights for a weighted repartitioning.
3530  const std::vector<unsigned int> cell_weights = get_cell_weights();
3531 
3532  // verify that the global sum of weights is larger than 0
3533  Assert(Utilities::MPI::sum(std::accumulate(cell_weights.begin(),
3534  cell_weights.end(),
3535  std::uint64_t(0)),
3536  this->mpi_communicator) > 0,
3537  ExcMessage(
3538  "The global sum of weights over all active cells "
3539  "is zero. Please verify how you generate weights."));
3540 
3541  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3542 
3543  // attach (temporarily) a pointer to the cell weights through
3544  // p4est's user_pointer object
3545  Assert(parallel_forest->user_pointer == this, ExcInternalError());
3546  parallel_forest->user_pointer = &partition_weights;
3547 
3549  parallel_forest,
3550  /* prepare coarsening */ 1,
3551  /* weight_callback */
3552  &PartitionWeights<dim, spacedim>::cell_weight);
3553 
3554  // reset the user pointer to its previous state
3555  parallel_forest->user_pointer = this;
3556  }
3557 
3558  // pack data before triangulation gets updated
3559  if (this->cell_attached_data.n_attached_data_sets > 0)
3560  {
3561  this->data_transfer.pack_data(
3562  this->local_cell_relations,
3563  this->cell_attached_data.pack_callbacks_fixed,
3564  this->cell_attached_data.pack_callbacks_variable);
3565  }
3566 
3567  try
3568  {
3569  copy_local_forest_to_triangulation();
3570  }
3571  catch (const typename Triangulation<dim>::DistortedCellList &)
3572  {
3573  // the underlying triangulation should not be checking for distorted
3574  // cells
3575  Assert(false, ExcInternalError());
3576  }
3577 
3578  // transfer data after triangulation got updated
3579  if (this->cell_attached_data.n_attached_data_sets > 0)
3580  {
3581  this->execute_transfer(parallel_forest,
3582  previous_global_first_quadrant.data());
3583  }
3584 
3585  this->update_periodic_face_map();
3586 
3587  // update how many cells, edges, etc, we store locally
3588  this->update_number_cache();
3589 
3590  // signal that repartitioning is finished
3591  this->signals.post_distributed_repartition();
3592  }
3593 
3594 
3595 
3596  template <int dim, int spacedim>
3597  const std::vector<types::global_dof_index> &
3599  const
3600  {
3601  return p4est_tree_to_coarse_cell_permutation;
3602  }
3603 
3604 
3605 
3606  template <int dim, int spacedim>
3607  const std::vector<types::global_dof_index> &
3609  const
3610  {
3611  return coarse_cell_to_p4est_tree_permutation;
3612  }
3613 
3614 
3615 
3616  template <int dim, int spacedim>
3617  std::vector<bool>
3619  const int level) const
3620  {
3621  Assert(dim > 1, ExcNotImplemented());
3622 
3623  std::vector<bool> marked_vertices(this->n_vertices(), false);
3624  for (const auto &cell : this->cell_iterators_on_level(level))
3625  if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3626  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
3627  marked_vertices[cell->vertex_index(v)] = true;
3628 
3634  // When a connectivity in the code below is detected, the assignment
3635  // 'marked_vertices[v1] = marked_vertices[v2] = true' makes sure that
3636  // the information about the periodicity propagates back to vertices on
3637  // cells that are not owned locally. However, in the worst case we want
3638  // to connect to a vertex that is 'dim' hops away from the locally owned
3639  // cell. Depending on the order of the periodic face map, we might
3640  // connect to that point by chance or miss it. However, after looping
3641  // through all the periodic directions (which are at most as many as
3642  // the number of space dimensions) we can be sure that all connections
3643  // to vertices have been created.
3644  for (unsigned int repetition = 0; repetition < dim; ++repetition)
3645  for (const auto &it : this->get_periodic_face_map())
3646  {
3647  const cell_iterator & cell_1 = it.first.first;
3648  const unsigned int face_no_1 = it.first.second;
3649  const cell_iterator & cell_2 = it.second.first.first;
3650  const unsigned int face_no_2 = it.second.first.second;
3651  const std::bitset<3> &face_orientation = it.second.second;
3652 
3653  if (cell_1->level() == level && cell_2->level() == level)
3654  {
3655  for (unsigned int v = 0;
3656  v < GeometryInfo<dim - 1>::vertices_per_cell;
3657  ++v)
3658  {
3659  // take possible non-standard orientation of faces into
3660  // account
3661  const unsigned int vface0 =
3663  v,
3664  face_orientation[0],
3665  face_orientation[1],
3666  face_orientation[2]);
3667  if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3668  vface0)] ||
3669  marked_vertices[cell_2->face(face_no_2)->vertex_index(
3670  v)])
3671  marked_vertices[cell_1->face(face_no_1)->vertex_index(
3672  vface0)] =
3673  marked_vertices[cell_2->face(face_no_2)->vertex_index(
3674  v)] = true;
3675  }
3676  }
3677  }
3678 
3679  return marked_vertices;
3680  }
3681 
3682 
3683 
3684  template <int dim, int spacedim>
3685  unsigned int
3688  {
3689  return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3690  }
3691 
3692 
3693 
3694  template <int dim, int spacedim>
3697  const unsigned int coarse_cell_index) const
3698  {
3699  return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3700  }
3701 
3702 
3703 
3704  template <int dim, int spacedim>
3705  void
3707  const std::vector<::GridTools::PeriodicFacePair<cell_iterator>>
3708  &periodicity_vector)
3709  {
3710  Assert(triangulation_has_content == true,
3711  ExcMessage("The triangulation is empty!"));
3712  Assert(this->n_levels() == 1,
3713  ExcMessage("The triangulation is refined!"));
3714 
3715  // call the base class for storing the periodicity information; we must
3716  // do this before going to p4est and rebuilding the triangulation to get
3717  // the level subdomain ids correct in the multigrid case
3719 
3720  for (const auto &face_pair : periodicity_vector)
3721  {
3722  const cell_iterator first_cell = face_pair.cell[0];
3723  const cell_iterator second_cell = face_pair.cell[1];
3724  const unsigned int face_left = face_pair.face_idx[0];
3725  const unsigned int face_right = face_pair.face_idx[1];
3726 
3727  // respective cells of the matching faces in p4est
3728  const unsigned int tree_left =
3729  coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3730  const unsigned int tree_right =
3731  coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3732 
3733  // p4est wants to know which corner the first corner on
3734  // the face with the lower id is mapped to on the face with
3735  // with the higher id. For d==2 there are only two possibilities
3736  // that are determined by it->orientation[1].
3737  // For d==3 we have to use GridTools::OrientationLookupTable.
3738  // The result is given below.
3739 
3740  unsigned int p4est_orientation = 0;
3741  if (dim == 2)
3742  p4est_orientation = face_pair.orientation[1];
3743  else
3744  {
3745  const unsigned int face_idx_list[] = {face_left, face_right};
3746  const cell_iterator cell_list[] = {first_cell, second_cell};
3747  unsigned int lower_idx, higher_idx;
3748  if (face_left <= face_right)
3749  {
3750  higher_idx = 1;
3751  lower_idx = 0;
3752  }
3753  else
3754  {
3755  higher_idx = 0;
3756  lower_idx = 1;
3757  }
3758 
3759  // get the cell index of the first index on the face with the
3760  // lower id
3761  unsigned int first_p4est_idx_on_cell =
3762  p8est_face_corners[face_idx_list[lower_idx]][0];
3763  unsigned int first_dealii_idx_on_face =
3765  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3766  ++i)
3767  {
3768  const unsigned int first_dealii_idx_on_cell =
3770  face_idx_list[lower_idx],
3771  i,
3772  cell_list[lower_idx]->face_orientation(
3773  face_idx_list[lower_idx]),
3774  cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3775  cell_list[lower_idx]->face_rotation(
3776  face_idx_list[lower_idx]));
3777  if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3778  {
3779  first_dealii_idx_on_face = i;
3780  break;
3781  }
3782  }
3783  Assert(first_dealii_idx_on_face != numbers::invalid_unsigned_int,
3784  ExcInternalError());
3785  // Now map dealii_idx_on_face according to the orientation
3786  constexpr unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
3787  {0, 1, 2, 3},
3788  {3, 1, 2, 0},
3789  {3, 2, 1, 0},
3790  {2, 3, 0, 1},
3791  {1, 3, 0, 2},
3792  {1, 0, 3, 2},
3793  {2, 0, 3, 1}};
3794  constexpr unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
3795  {0, 1, 2, 3},
3796  {3, 1, 2, 0},
3797  {3, 2, 1, 0},
3798  {2, 3, 0, 1},
3799  {2, 0, 3, 1},
3800  {1, 0, 3, 2},
3801  {1, 3, 0, 2}};
3802  const unsigned int second_dealii_idx_on_face =
3803  lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
3804  [first_dealii_idx_on_face] :
3805  right_to_left[face_pair.orientation.to_ulong()]
3806  [first_dealii_idx_on_face];
3807  const unsigned int second_dealii_idx_on_cell =
3809  face_idx_list[higher_idx],
3810  second_dealii_idx_on_face,
3811  cell_list[higher_idx]->face_orientation(
3812  face_idx_list[higher_idx]),
3813  cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
3814  cell_list[higher_idx]->face_rotation(
3815  face_idx_list[higher_idx]));
3816  // map back to p4est
3817  const unsigned int second_p4est_idx_on_face =
3818  p8est_corner_face_corners[second_dealii_idx_on_cell]
3819  [face_idx_list[higher_idx]];
3820  p4est_orientation = second_p4est_idx_on_face;
3821  }
3822 
3824  connectivity,
3825  tree_left,
3826  tree_right,
3827  face_left,
3828  face_right,
3829  p4est_orientation);
3830  }
3831 
3832 
3834  connectivity) == 1,
3835  ExcInternalError());
3836 
3837  // now create a forest out of the connectivity data structure
3840  this->mpi_communicator,
3841  connectivity,
3842  /* minimum initial number of quadrants per tree */ 0,
3843  /* minimum level of upfront refinement */ 0,
3844  /* use uniform upfront refinement */ 1,
3845  /* user_data_size = */ 0,
3846  /* user_data_constructor = */ nullptr,
3847  /* user_pointer */ this);
3848 
3849  try
3850  {
3851  copy_local_forest_to_triangulation();
3852  }
3853  catch (const typename Triangulation<dim>::DistortedCellList &)
3854  {
3855  // the underlying triangulation should not be checking for distorted
3856  // cells
3857  Assert(false, ExcInternalError());
3858  }
3859 
3860  // The range of ghost_owners might have changed so update that
3861  // information
3862  this->update_number_cache();
3863  }
3864 
3865 
3866 
3867  template <int dim, int spacedim>
3868  std::size_t
3870  {
3871  std::size_t mem =
3874  MemoryConsumption::memory_consumption(triangulation_has_content) +
3876  MemoryConsumption::memory_consumption(parallel_forest) +
3878  this->cell_attached_data.n_attached_data_sets) +
3879  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_fixed)
3880  // +
3881  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_variable)
3882  // +
3883  // TODO[TH]: how?
3885  coarse_cell_to_p4est_tree_permutation) +
3887  p4est_tree_to_coarse_cell_permutation) +
3888  memory_consumption_p4est();
3889 
3890  return mem;
3891  }
3892 
3893 
3894 
3895  template <int dim, int spacedim>
3896  std::size_t
3898  {
3899  return ::internal::p4est::functions<dim>::forest_memory_used(
3900  parallel_forest) +
3902  connectivity);
3903  }
3904 
3905 
3906 
3907  template <int dim, int spacedim>
3908  void
3910  const ::Triangulation<dim, spacedim> &other_tria)
3911  {
3912  Assert(
3913  (dynamic_cast<
3914  const ::parallel::distributed::Triangulation<dim, spacedim> *>(
3915  &other_tria)) ||
3916  (other_tria.n_global_levels() == 1),
3917  ExcNotImplemented());
3918 
3920 
3921  try
3922  {
3924  copy_triangulation(other_tria);
3925  }
3926  catch (
3927  const typename ::Triangulation<dim, spacedim>::DistortedCellList
3928  &)
3929  {
3930  // the underlying triangulation should not be checking for distorted
3931  // cells
3932  Assert(false, ExcInternalError());
3933  }
3934 
3935  if (const ::parallel::distributed::Triangulation<dim, spacedim>
3936  *other_distributed =
3937  dynamic_cast<const ::parallel::distributed::
3938  Triangulation<dim, spacedim> *>(&other_tria))
3939  {
3940  // copy parallel distributed specifics
3941  settings = other_distributed->settings;
3942  triangulation_has_content =
3943  other_distributed->triangulation_has_content;
3944  coarse_cell_to_p4est_tree_permutation =
3945  other_distributed->coarse_cell_to_p4est_tree_permutation;
3946  p4est_tree_to_coarse_cell_permutation =
3947  other_distributed->p4est_tree_to_coarse_cell_permutation;
3948 
3949  // create deep copy of connectivity graph
3950  typename ::internal::p4est::types<dim>::connectivity
3951  *temp_connectivity = const_cast<
3952  typename ::internal::p4est::types<dim>::connectivity *>(
3953  other_distributed->connectivity);
3954  connectivity =
3955  ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
3956 
3957  // create deep copy of parallel forest
3958  typename ::internal::p4est::types<dim>::forest *temp_forest =
3959  const_cast<typename ::internal::p4est::types<dim>::forest *>(
3960  other_distributed->parallel_forest);
3961  parallel_forest =
3963  false);
3964  parallel_forest->connectivity = connectivity;
3965  parallel_forest->user_pointer = this;
3966  }
3967  else
3968  {
3969  triangulation_has_content = true;
3970  setup_coarse_cell_to_p4est_tree_permutation();
3971  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
3972  }
3973 
3974  try
3975  {
3976  copy_local_forest_to_triangulation();
3977  }
3978  catch (const typename Triangulation<dim>::DistortedCellList &)
3979  {
3980  // the underlying triangulation should not be checking for distorted
3981  // cells
3982  Assert(false, ExcInternalError());
3983  }
3984 
3985  this->update_periodic_face_map();
3986  this->update_number_cache();
3987  }
3988 
3989 
3990 
3991  template <int dim, int spacedim>
3992  void
3994  {
3995  // reorganize memory for local_cell_relations
3996  this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
3997  this->local_cell_relations.shrink_to_fit();
3998 
3999  // recurse over p4est
4000  for (const auto &cell : this->cell_iterators_on_level(0))
4001  {
4002  // skip coarse cells that are not ours
4003  if (tree_exists_locally<dim, spacedim>(
4004  parallel_forest,
4005  coarse_cell_to_p4est_tree_permutation[cell->index()]) == false)
4006  continue;
4007 
4008  // initialize auxiliary top level p4est quadrant
4009  typename ::internal::p4est::types<dim>::quadrant
4010  p4est_coarse_cell;
4011  ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4012 
4013  // determine tree to start recursion on
4014  typename ::internal::p4est::types<dim>::tree *tree =
4015  init_tree(cell->index());
4016 
4017  update_cell_relations_recursively<dim, spacedim>(
4018  this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4019  }
4020  }
4021 
4022 
4023 
4024  template <int dim, int spacedim>
4025  std::vector<unsigned int>
4027  {
4028  // check if local_cell_relations have been previously gathered
4029  // correctly
4030  Assert(this->local_cell_relations.size() ==
4031  static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4032  ExcInternalError());
4033 
4034  // Allocate the space for the weights. We reserve an integer for each
4035  // locally owned quadrant on the already refined p4est object.
4036  std::vector<unsigned int> weights;
4037  weights.reserve(this->local_cell_relations.size());
4038 
4039  // Iterate over p4est and Triangulation relations
4040  // to find refined/coarsened/kept
4041  // cells. Then append weight.
4042  // Note that we need to follow the p4est ordering
4043  // instead of the deal.II ordering to get the weights
4044  // in the same order p4est will encounter them during repartitioning.
4045  for (const auto &cell_rel : this->local_cell_relations)
4046  {
4047  const auto &cell_it = cell_rel.first;
4048  const auto &cell_status = cell_rel.second;
4049 
4050  weights.push_back(this->signals.weight(cell_it, cell_status));
4051  }
4052 
4053  return weights;
4054  }
4055 
4056 
4057 
4058  template <int spacedim>
4060  const MPI_Comm &mpi_communicator,
4061  const typename ::Triangulation<1, spacedim>::MeshSmoothing
4062  smooth_grid,
4063  const Settings /*settings*/)
4064  : ::parallel::DistributedTriangulationBase<1, spacedim>(
4065  mpi_communicator,
4066  smooth_grid,
4067  false)
4068  {
4069  Assert(false, ExcNotImplemented());
4070  }
4071 
4072 
4073  template <int spacedim>
4075  {
4076  AssertNothrow(false, ExcNotImplemented());
4077  }
4078 
4079 
4080 
4081  template <int spacedim>
4082  const std::vector<types::global_dof_index> &
4084  const
4085  {
4086  static std::vector<types::global_dof_index> a;
4087  return a;
4088  }
4089 
4090 
4091 
4092  template <int spacedim>
4093  std::map<unsigned int, std::set<::types::subdomain_id>>
4095  const unsigned int /*level*/) const
4096  {
4097  Assert(false, ExcNotImplemented());
4098 
4099  return std::map<unsigned int, std::set<::types::subdomain_id>>();
4100  }
4101 
4102 
4103 
4104  template <int spacedim>
4105  std::vector<bool>
4107  const unsigned int) const
4108  {
4109  Assert(false, ExcNotImplemented());
4110  return std::vector<bool>();
4111  }
4112 
4113 
4114 
4115  template <int spacedim>
4116  unsigned int
4118  const types::coarse_cell_id) const
4119  {
4120  Assert(false, ExcNotImplemented());
4121  return 0;
4122  }
4123 
4124 
4125 
4126  template <int spacedim>
4129  const unsigned int) const
4130  {
4131  Assert(false, ExcNotImplemented());
4132  return 0;
4133  }
4134 
4135 
4136 
4137  template <int spacedim>
4138  void
4140  {
4141  Assert(false, ExcNotImplemented());
4142  }
4143 
4144 
4145 
4146  template <int spacedim>
4147  void
4148  Triangulation<1, spacedim>::load(const std::string &, const bool)
4149  {
4150  Assert(false, ExcNotImplemented());
4151  }
4152 
4153 
4154 
4155  template <int spacedim>
4156  void
4157  Triangulation<1, spacedim>::save(const std::string &) const
4158  {
4159  Assert(false, ExcNotImplemented());
4160  }
4161 
4162 
4163 
4164  template <int spacedim>
4165  bool
4167  {
4168  Assert(false, ExcNotImplemented());
4169  return false;
4170  }
4171 
4172 
4173 
4174  template <int spacedim>
4175  bool
4177  {
4178  Assert(false, ExcNotImplemented());
4179  return false;
4180  }
4181 
4182 
4183 
4184  template <int spacedim>
4185  void
4187  {
4188  Assert(false, ExcNotImplemented());
4189  }
4190 
4191  } // namespace distributed
4192 } // namespace parallel
4193 
4194 
4195 #endif // DEAL_II_WITH_P4EST
4196 
4197 
4198 
4199 namespace parallel
4200 {
4201  namespace distributed
4202  {
4203  template <int dim, int spacedim>
4206  : distributed_tria(
4207  dynamic_cast<
4208  ::parallel::distributed::Triangulation<dim, spacedim> *>(
4209  &tria))
4210  {
4211 #ifdef DEAL_II_WITH_P4EST
4212  if (distributed_tria != nullptr)
4213  {
4214  // Save the current set of refinement flags, and adjust the
4215  // refinement flags to be consistent with the p4est oracle.
4216  distributed_tria->save_coarsen_flags(saved_coarsen_flags);
4217  distributed_tria->save_refine_flags(saved_refine_flags);
4218 
4219  for (const auto &pair : distributed_tria->local_cell_relations)
4220  {
4221  const auto &cell = pair.first;
4222  const auto &status = pair.second;
4223 
4224  switch (status)
4225  {
4226  case ::Triangulation<dim, spacedim>::CELL_PERSIST:
4227  // cell remains unchanged
4228  cell->clear_refine_flag();
4229  cell->clear_coarsen_flag();
4230  break;
4231 
4232  case ::Triangulation<dim, spacedim>::CELL_REFINE:
4233  // cell will be refined
4234  cell->clear_coarsen_flag();
4235  cell->set_refine_flag();
4236  break;
4237 
4238  case ::Triangulation<dim, spacedim>::CELL_COARSEN:
4239  // children of this cell will be coarsened
4240  for (const auto &child : cell->child_iterators())
4241  {
4242  child->clear_refine_flag();
4243  child->set_coarsen_flag();
4244  }
4245  break;
4246 
4247  case ::Triangulation<dim, spacedim>::CELL_INVALID:
4248  // do nothing as cell does not exist yet
4249  break;
4250 
4251  default:
4252  Assert(false, ExcInternalError());
4253  break;
4254  }
4255  }
4256  }
4257 #endif
4258  }
4259 
4260 
4261 
4262  template <int dim, int spacedim>
4264  {
4265 #ifdef DEAL_II_WITH_P4EST
4266  if (distributed_tria)
4267  {
4268  // Undo the refinement flags modification.
4269  distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4270  distributed_tria->load_refine_flags(saved_refine_flags);
4271  }
4272 #else
4273  // pretend that this destructor does something to silence clang-tidy
4274  (void)distributed_tria;
4275 #endif
4276  }
4277  } // namespace distributed
4278 } // namespace parallel
4279 
4280 
4281 
4282 /*-------------- Explicit Instantiations -------------------------------*/
4283 #include "tria.inst"
4284 
4285 
Definition: point.h:111
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
active_cell_iterator last_active() const
virtual types::subdomain_id locally_owned_subdomain() const
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator >> &)
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
virtual bool prepare_coarsening_and_refinement()
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
void save_refine_flags(std::ostream &out) const
unsigned int n_vertices() const
void save_coarsen_flags(std::ostream &out) const
active_cell_iterator begin_active(const unsigned int level=0) const
virtual void clear() override
Definition: tria_base.cc:655
virtual std::size_t memory_consumption() const override
Definition: tria_base.cc:92
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
Definition: tria_base.cc:67
TemporarilyMatchRefineFlags(::Triangulation< dim, spacedim > &tria)
Definition: tria.cc:4204
const SmartPointer< ::parallel::distributed::Triangulation< dim, spacedim > > distributed_tria
Definition: tria.h:1133
virtual void execute_coarsening_and_refinement() override
Definition: tria.cc:3225
void setup_coarse_cell_to_p4est_tree_permutation()
Definition: tria.cc:1965
bool are_vertices_communicated_to_p4est() const
Definition: tria.cc:1860
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const override
Definition: tria.cc:3696
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &other_tria) override
Definition: tria.cc:3909
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
Definition: tria.cc:3598
typename ::internal::p4est::types< dim >::ghost * parallel_ghost
Definition: tria.h:731
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:1699
std::vector< unsigned int > get_cell_weights() const
Definition: tria.cc:4026
unsigned int get_checksum() const
Definition: tria.cc:2220
virtual ~Triangulation() override
Definition: tria.cc:1727
virtual void update_cell_relations() override
Definition: tria.cc:3993
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:1870
typename ::internal::p4est::types< dim >::tree * init_tree(const int dealii_coarse_cell_index) const
Definition: tria.cc:2243
types::subdomain_id find_point_owner_rank(const Point< dim > &p)
Definition: tria.cc:3112
virtual bool prepare_coarsening_and_refinement() override
Definition: tria.cc:2737
std::vector< bool > mark_locally_active_vertices_on_level(const int level) const
Definition: tria.cc:3618
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator >> &) override
Definition: tria.cc:3706
const ::internal::p4est::types< dim >::forest * get_p4est() const
Definition: tria.cc:2232
const std::vector< types::global_dof_index > & get_coarse_cell_to_p4est_tree_permutation() const
Definition: tria.cc:3608
void write_mesh_vtk(const std::string &file_basename) const
Definition: tria.cc:1982
virtual unsigned int coarse_cell_id_to_coarse_cell_index(const types::coarse_cell_id coarse_cell_id) const override
Definition: tria.cc:3686
void copy_new_triangulation_to_p4est(std::integral_constant< int, 2 >)
virtual void save(const std::string &filename) const override
Definition: tria.cc:2001
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata) override
Definition: tria.cc:1748
virtual std::size_t memory_consumption_p4est() const
Definition: tria.cc:3897
virtual void clear() override
Definition: tria.cc:1814
bool is_multilevel_hierarchy_constructed() const override
Definition: tria.cc:1850
virtual void load(const std::string &filename) override
Definition: tria.cc:2056
virtual std::size_t memory_consumption() const override
Definition: tria.cc:3869
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4608
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1509
static ::ExceptionBase & ExcNotImplemented()
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1572
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1619
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
Definition: tria.h:267
typename ::Triangulation< dim, spacedim >::active_cell_iterator active_cell_iterator
Definition: tria.h:288
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3921
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:189
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void partition(const SparsityPattern &sparsity_pattern, const unsigned int n_partitions, std::vector< unsigned int > &partition_indices, const Partitioner partitioner=Partitioner::metis)
void reorder_hierarchical(const DynamicSparsityPattern &sparsity, std::vector< DynamicSparsityPattern::size_type > &new_indices)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:157
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:146
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1655
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 3 > &c)
Definition: tria.cc:13786
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 3 > &c)
Definition: tria.cc:13793
bool tree_exists_locally(const typename types< dim >::forest *parallel_forest, const typename types< dim >::topidx coarse_grid_cell)
const types::subdomain_id artificial_subdomain_id
Definition: types.h:308
const types::subdomain_id invalid_subdomain_id
Definition: types.h:291
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::manifold_id flat_manifold_id
Definition: types.h:279
unsigned int manifold_id
Definition: types.h:146
global_cell_index coarse_cell_id
Definition: types.h:119
unsigned int subdomain_id
Definition: types.h:43
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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)
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)
static Point< dim > unit_cell_vertex(const unsigned int vertex)
static bool is_inside_unit_cell(const Point< dim > &p)
const TriangulationDescription::Settings settings
const ::Triangulation< dim, spacedim > & tria