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