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