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