Reference documentation for deal.II version Git 713825e468 2021-05-17 16:05:53 -0400
\(\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/utilities.h>
20 
23 
25 #include <deal.II/grid/tria.h>
28 
31 
32 #include <algorithm>
33 #include <fstream>
34 #include <iostream>
35 #include <numeric>
36 
37 
39 
40 
41 #ifdef DEAL_II_WITH_P4EST
42 
43 namespace
44 {
45  template <int dim, int spacedim>
46  void
47  get_vertex_to_cell_mappings(
49  std::vector<unsigned int> & vertex_touch_count,
50  std::vector<std::list<
52  unsigned int>>> & vertex_to_cell)
53  {
54  vertex_touch_count.resize(triangulation.n_vertices());
55  vertex_to_cell.resize(triangulation.n_vertices());
56 
57  for (const auto &cell : triangulation.active_cell_iterators())
58  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
59  {
60  ++vertex_touch_count[cell->vertex_index(v)];
61  vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
62  }
63  }
64 
65 
66 
67  template <int dim, int spacedim>
68  void
69  get_edge_to_cell_mappings(
70  const Triangulation<dim, spacedim> &triangulation,
71  std::vector<unsigned int> & edge_touch_count,
72  std::vector<std::list<
74  unsigned int>>> & edge_to_cell)
75  {
76  Assert(triangulation.n_levels() == 1, ExcInternalError());
77 
78  edge_touch_count.resize(triangulation.n_active_lines());
79  edge_to_cell.resize(triangulation.n_active_lines());
80 
81  for (const auto &cell : triangulation.active_cell_iterators())
82  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
83  {
84  ++edge_touch_count[cell->line(l)->index()];
85  edge_to_cell[cell->line(l)->index()].emplace_back(cell, l);
86  }
87  }
88 
89 
90 
95  template <int dim, int spacedim>
96  void
97  set_vertex_and_cell_info(
98  const Triangulation<dim, spacedim> &triangulation,
99  const std::vector<unsigned int> & vertex_touch_count,
100  const std::vector<std::list<
102  unsigned int>>> & vertex_to_cell,
103  const std::vector<types::global_dof_index>
104  & coarse_cell_to_p4est_tree_permutation,
105  const bool set_vertex_info,
106  typename internal::p4est::types<dim>::connectivity *connectivity)
107  {
108  // copy the vertices into the connectivity structure. the triangulation
109  // exports the array of vertices, but some of the entries are sometimes
110  // unused; this shouldn't be the case for a newly created triangulation,
111  // but make sure
112  //
113  // note that p4est stores coordinates as a triplet of values even in 2d
114  Assert(triangulation.get_used_vertices().size() ==
115  triangulation.get_vertices().size(),
116  ExcInternalError());
117  Assert(std::find(triangulation.get_used_vertices().begin(),
118  triangulation.get_used_vertices().end(),
119  false) == triangulation.get_used_vertices().end(),
120  ExcInternalError());
121  if (set_vertex_info == true)
122  for (unsigned int v = 0; v < triangulation.n_vertices(); ++v)
123  {
124  connectivity->vertices[3 * v] = triangulation.get_vertices()[v][0];
125  connectivity->vertices[3 * v + 1] =
126  triangulation.get_vertices()[v][1];
127  connectivity->vertices[3 * v + 2] =
128  (spacedim == 2 ? 0 : triangulation.get_vertices()[v][2]);
129  }
130 
131  // next store the tree_to_vertex indices (each tree is here only a single
132  // cell in the coarse mesh). p4est requires vertex numbering in clockwise
133  // orientation
134  //
135  // while we're at it, also copy the neighborship information between cells
137  cell = triangulation.begin_active(),
138  endc = triangulation.end();
139  for (; cell != endc; ++cell)
140  {
141  const unsigned int index =
142  coarse_cell_to_p4est_tree_permutation[cell->index()];
143 
144  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
145  {
146  if (set_vertex_info == true)
147  connectivity
148  ->tree_to_vertex[index * GeometryInfo<dim>::vertices_per_cell +
149  v] = cell->vertex_index(v);
150  connectivity
151  ->tree_to_corner[index * GeometryInfo<dim>::vertices_per_cell +
152  v] = cell->vertex_index(v);
153  }
154 
155  // neighborship information. if a cell is at a boundary, then enter
156  // the index of the cell itself here
157  for (auto f : GeometryInfo<dim>::face_indices())
158  if (cell->face(f)->at_boundary() == false)
159  connectivity
160  ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
161  coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->index()];
162  else
163  connectivity
164  ->tree_to_tree[index * GeometryInfo<dim>::faces_per_cell + f] =
165  coarse_cell_to_p4est_tree_permutation[cell->index()];
166 
167  // fill tree_to_face, which is essentially neighbor_to_neighbor;
168  // however, we have to remap the resulting face number as well
169  for (auto f : GeometryInfo<dim>::face_indices())
170  if (cell->face(f)->at_boundary() == false)
171  {
172  switch (dim)
173  {
174  case 2:
175  {
176  connectivity->tree_to_face
177  [index * GeometryInfo<dim>::faces_per_cell + f] =
178  cell->neighbor_of_neighbor(f);
179  break;
180  }
181 
182  case 3:
183  {
184  /*
185  * The values for tree_to_face are in 0..23 where ttf % 6
186  * gives the face number and ttf / 4 the face orientation
187  * code. The orientation is determined as follows. Let
188  * my_face and other_face be the two face numbers of the
189  * connecting trees in 0..5. Then the first face vertex
190  * of the lower of my_face and other_face connects to a
191  * face vertex numbered 0..3 in the higher of my_face and
192  * other_face. The face orientation is defined as this
193  * number. If my_face == other_face, treating either of
194  * both faces as the lower one leads to the same result.
195  */
196 
197  connectivity->tree_to_face[index * 6 + f] =
198  cell->neighbor_of_neighbor(f);
199 
200  unsigned int face_idx_list[2] = {
201  f, cell->neighbor_of_neighbor(f)};
203  cell_list[2] = {cell, cell->neighbor(f)};
204  unsigned int smaller_idx = 0;
205 
206  if (f > cell->neighbor_of_neighbor(f))
207  smaller_idx = 1;
208 
209  unsigned int larger_idx = (smaller_idx + 1) % 2;
210  // smaller = *_list[smaller_idx]
211  // larger = *_list[larger_idx]
212 
213  unsigned int v = 0;
214 
215  // global vertex index of vertex 0 on face of cell with
216  // smaller local face index
217  unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
219  face_idx_list[smaller_idx],
220  0,
221  cell_list[smaller_idx]->face_orientation(
222  face_idx_list[smaller_idx]),
223  cell_list[smaller_idx]->face_flip(
224  face_idx_list[smaller_idx]),
225  cell_list[smaller_idx]->face_rotation(
226  face_idx_list[smaller_idx])));
227 
228  // loop over vertices on face from other cell and compare
229  // global vertex numbers
230  for (unsigned int i = 0;
231  i < GeometryInfo<dim>::vertices_per_face;
232  ++i)
233  {
234  unsigned int idx =
235  cell_list[larger_idx]->vertex_index(
237  face_idx_list[larger_idx], i));
238 
239  if (idx == g_idx)
240  {
241  v = i;
242  break;
243  }
244  }
245 
246  connectivity->tree_to_face[index * 6 + f] += 6 * v;
247  break;
248  }
249 
250  default:
251  Assert(false, ExcNotImplemented());
252  }
253  }
254  else
255  connectivity
256  ->tree_to_face[index * GeometryInfo<dim>::faces_per_cell + f] = f;
257  }
258 
259  // now fill the vertex information
260  connectivity->ctt_offset[0] = 0;
261  std::partial_sum(vertex_touch_count.begin(),
262  vertex_touch_count.end(),
263  &connectivity->ctt_offset[1]);
264 
265  const typename internal::p4est::types<dim>::locidx num_vtt =
266  std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
267  (void)num_vtt;
268  Assert(connectivity->ctt_offset[triangulation.n_vertices()] == num_vtt,
269  ExcInternalError());
270 
271  for (unsigned int v = 0; v < triangulation.n_vertices(); ++v)
272  {
273  Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
274  ExcInternalError());
275 
276  typename std::list<
277  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
278  unsigned int>>::const_iterator p =
279  vertex_to_cell[v].begin();
280  for (unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
281  {
282  connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
283  coarse_cell_to_p4est_tree_permutation[p->first->index()];
284  connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
285  p->second;
286  }
287  }
288  }
289 
290 
291 
292  template <int dim, int spacedim>
293  bool
295  const typename internal::p4est::types<dim>::forest *parallel_forest,
296  const typename internal::p4est::types<dim>::topidx coarse_grid_cell)
297  {
298  Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
299  ExcInternalError());
300  return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
301  (coarse_grid_cell <= parallel_forest->last_local_tree));
302  }
303 
304 
305  template <int dim, int spacedim>
306  void
307  delete_all_children_and_self(
308  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
309  {
310  if (cell->has_children())
311  for (unsigned int c = 0; c < cell->n_children(); ++c)
312  delete_all_children_and_self<dim, spacedim>(cell->child(c));
313  else
314  cell->set_coarsen_flag();
315  }
316 
317 
318 
319  template <int dim, int spacedim>
320  void
321  delete_all_children(
322  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
323  {
324  if (cell->has_children())
325  for (unsigned int c = 0; c < cell->n_children(); ++c)
326  delete_all_children_and_self<dim, spacedim>(cell->child(c));
327  }
328 
329 
330  template <int dim, int spacedim>
331  void
332  determine_level_subdomain_id_recursively(
333  const typename internal::p4est::types<dim>::tree & tree,
334  const typename internal::p4est::types<dim>::locidx & tree_index,
335  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
336  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
337  typename internal::p4est::types<dim>::forest & forest,
338  const types::subdomain_id my_subdomain,
339  const std::vector<std::vector<bool>> & marked_vertices)
340  {
341  if (dealii_cell->level_subdomain_id() == numbers::artificial_subdomain_id)
342  {
343  // important: only assign the level_subdomain_id if it is a ghost cell
344  // even though we could fill in all.
345  bool used = false;
346  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
347  {
348  if (marked_vertices[dealii_cell->level()]
349  [dealii_cell->vertex_index(v)])
350  {
351  used = true;
352  break;
353  }
354  }
355 
356  // Special case: if this cell is active we might be a ghost neighbor
357  // to a locally owned cell across a vertex that is finer.
358  // Example (M= my, O=dealii_cell, owned by somebody else):
359  // *------*
360  // | |
361  // | O |
362  // | |
363  // *---*---*------*
364  // | M | M |
365  // *---*---*
366  // | | M |
367  // *---*---*
368  if (!used && dealii_cell->is_active() &&
369  dealii_cell->is_artificial() == false &&
370  dealii_cell->level() + 1 < static_cast<int>(marked_vertices.size()))
371  {
372  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
373  {
374  if (marked_vertices[dealii_cell->level() + 1]
375  [dealii_cell->vertex_index(v)])
376  {
377  used = true;
378  break;
379  }
380  }
381  }
382 
383  // Like above, but now the other way around
384  if (!used && dealii_cell->is_active() &&
385  dealii_cell->is_artificial() == false && dealii_cell->level() > 0)
386  {
387  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
388  {
389  if (marked_vertices[dealii_cell->level() - 1]
390  [dealii_cell->vertex_index(v)])
391  {
392  used = true;
393  break;
394  }
395  }
396  }
397 
398  if (used)
399  {
401  &forest, tree_index, &p4est_cell, my_subdomain);
402  Assert((owner != -2) && (owner != -1),
403  ExcMessage("p4est should know the owner."));
404  dealii_cell->set_level_subdomain_id(owner);
405  }
406  }
407 
408  if (dealii_cell->has_children())
409  {
412  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
413  ++c)
414  switch (dim)
415  {
416  case 2:
417  P4EST_QUADRANT_INIT(&p4est_child[c]);
418  break;
419  case 3:
420  P8EST_QUADRANT_INIT(&p4est_child[c]);
421  break;
422  default:
423  Assert(false, ExcNotImplemented());
424  }
425 
426 
428  p4est_child);
429 
430  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
431  ++c)
432  {
433  determine_level_subdomain_id_recursively<dim, spacedim>(
434  tree,
435  tree_index,
436  dealii_cell->child(c),
437  p4est_child[c],
438  forest,
439  my_subdomain,
440  marked_vertices);
441  }
442  }
443  }
444 
445 
446  template <int dim, int spacedim>
447  void
448  match_tree_recursively(
449  const typename internal::p4est::types<dim>::tree & tree,
450  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
451  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
452  const typename internal::p4est::types<dim>::forest & forest,
453  const types::subdomain_id my_subdomain)
454  {
455  // check if this cell exists in the local p4est cell
456  if (sc_array_bsearch(const_cast<sc_array_t *>(&tree.quadrants),
457  &p4est_cell,
459  -1)
460  {
461  // yes, cell found in local part of p4est
462  delete_all_children<dim, spacedim>(dealii_cell);
463  if (dealii_cell->is_active())
464  dealii_cell->set_subdomain_id(my_subdomain);
465  }
466  else
467  {
468  // no, cell not found in local part of p4est. this means that the
469  // local part is more refined than the current cell. if this cell has
470  // no children of its own, we need to refine it, and if it does
471  // already have children then loop over all children and see if they
472  // are locally available as well
473  if (dealii_cell->is_active())
474  dealii_cell->set_refine_flag();
475  else
476  {
479  for (unsigned int c = 0;
480  c < GeometryInfo<dim>::max_children_per_cell;
481  ++c)
482  switch (dim)
483  {
484  case 2:
485  P4EST_QUADRANT_INIT(&p4est_child[c]);
486  break;
487  case 3:
488  P8EST_QUADRANT_INIT(&p4est_child[c]);
489  break;
490  default:
491  Assert(false, ExcNotImplemented());
492  }
493 
494 
496  p4est_child);
497 
498  for (unsigned int c = 0;
499  c < GeometryInfo<dim>::max_children_per_cell;
500  ++c)
502  const_cast<typename internal::p4est::types<dim>::tree *>(
503  &tree),
504  &p4est_child[c]) == false)
505  {
506  // no, this child is locally not available in the p4est.
507  // delete all its children but, because this may not be
508  // successful, make sure to mark all children recursively
509  // as not local.
510  delete_all_children<dim, spacedim>(dealii_cell->child(c));
511  dealii_cell->child(c)->recursively_set_subdomain_id(
513  }
514  else
515  {
516  // at least some part of the tree rooted in this child is
517  // locally available
518  match_tree_recursively<dim, spacedim>(tree,
519  dealii_cell->child(c),
520  p4est_child[c],
521  forest,
522  my_subdomain);
523  }
524  }
525  }
526  }
527 
528 
529  template <int dim, int spacedim>
530  void
531  match_quadrant(
532  const ::Triangulation<dim, spacedim> * tria,
533  unsigned int dealii_index,
534  const typename internal::p4est::types<dim>::quadrant &ghost_quadrant,
535  types::subdomain_id ghost_owner)
536  {
537  const int l = ghost_quadrant.level;
538 
539  for (int i = 0; i < l; ++i)
540  {
542  i,
543  dealii_index);
544  if (cell->is_active())
545  {
546  cell->clear_coarsen_flag();
547  cell->set_refine_flag();
548  return;
549  }
550 
551  const int child_id =
553  i + 1);
554  dealii_index = cell->child_index(child_id);
555  }
556 
558  l,
559  dealii_index);
560  if (cell->has_children())
561  delete_all_children<dim, spacedim>(cell);
562  else
563  {
564  cell->clear_coarsen_flag();
565  cell->set_subdomain_id(ghost_owner);
566  }
567  }
568 
569 
570 
576  template <int dim, int spacedim>
577  class RefineAndCoarsenList
578  {
579  public:
580  RefineAndCoarsenList(const Triangulation<dim, spacedim> &triangulation,
581  const std::vector<types::global_dof_index>
582  &p4est_tree_to_coarse_cell_permutation,
583  const types::subdomain_id my_subdomain);
584 
593  static int
594  refine_callback(
595  typename internal::p4est::types<dim>::forest * forest,
596  typename internal::p4est::types<dim>::topidx coarse_cell_index,
597  typename internal::p4est::types<dim>::quadrant *quadrant);
598 
603  static int
604  coarsen_callback(
605  typename internal::p4est::types<dim>::forest * forest,
606  typename internal::p4est::types<dim>::topidx coarse_cell_index,
607  typename internal::p4est::types<dim>::quadrant *children[]);
608 
609  bool
610  pointers_are_at_end() const;
611 
612  private:
613  std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
614  typename std::vector<typename internal::p4est::types<dim>::quadrant>::
615  const_iterator current_refine_pointer;
616 
617  std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
618  typename std::vector<typename internal::p4est::types<dim>::quadrant>::
619  const_iterator current_coarsen_pointer;
620 
621  void
622  build_lists(
623  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
624  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
625  const types::subdomain_id myid);
626  };
627 
628 
629 
630  template <int dim, int spacedim>
631  bool
632  RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end() const
633  {
634  return ((current_refine_pointer == refine_list.end()) &&
635  (current_coarsen_pointer == coarsen_list.end()));
636  }
637 
638 
639 
640  template <int dim, int spacedim>
641  RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
642  const Triangulation<dim, spacedim> &triangulation,
643  const std::vector<types::global_dof_index>
644  & p4est_tree_to_coarse_cell_permutation,
645  const types::subdomain_id my_subdomain)
646  {
647  // count how many flags are set and allocate that much memory
648  unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
649  for (const auto &cell : triangulation.active_cell_iterators())
650  {
651  // skip cells that are not local
652  if (cell->subdomain_id() != my_subdomain)
653  continue;
654 
655  if (cell->refine_flag_set())
656  ++n_refine_flags;
657  else if (cell->coarsen_flag_set())
658  ++n_coarsen_flags;
659  }
660 
661  refine_list.reserve(n_refine_flags);
662  coarsen_list.reserve(n_coarsen_flags);
663 
664 
665  // now build the lists of cells that are flagged. note that p4est will
666  // traverse its cells in the order in which trees appear in the
667  // forest. this order is not the same as the order of coarse cells in the
668  // deal.II Triangulation because we have translated everything by the
669  // coarse_cell_to_p4est_tree_permutation permutation. in order to make
670  // sure that the output array is already in the correct order, traverse
671  // our coarse cells in the same order in which p4est will:
672  for (unsigned int c = 0; c < triangulation.n_cells(0); ++c)
673  {
674  unsigned int coarse_cell_index =
675  p4est_tree_to_coarse_cell_permutation[c];
676 
678  &triangulation, 0, coarse_cell_index);
679 
680  typename internal::p4est::types<dim>::quadrant p4est_cell;
682  /*level=*/0,
683  /*index=*/0);
684  p4est_cell.p.which_tree = c;
685  build_lists(cell, p4est_cell, my_subdomain);
686  }
687 
688 
689  Assert(refine_list.size() == n_refine_flags, ExcInternalError());
690  Assert(coarsen_list.size() == n_coarsen_flags, ExcInternalError());
691 
692  // make sure that our ordering in fact worked
693  for (unsigned int i = 1; i < refine_list.size(); ++i)
694  Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
695  ExcInternalError());
696  for (unsigned int i = 1; i < coarsen_list.size(); ++i)
697  Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
698  ExcInternalError());
699 
700  current_refine_pointer = refine_list.begin();
701  current_coarsen_pointer = coarsen_list.begin();
702  }
703 
704 
705 
706  template <int dim, int spacedim>
707  void
708  RefineAndCoarsenList<dim, spacedim>::build_lists(
709  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
710  const typename internal::p4est::types<dim>::quadrant & p4est_cell,
711  const types::subdomain_id my_subdomain)
712  {
713  if (cell->is_active())
714  {
715  if (cell->subdomain_id() == my_subdomain)
716  {
717  if (cell->refine_flag_set())
718  refine_list.push_back(p4est_cell);
719  else if (cell->coarsen_flag_set())
720  coarsen_list.push_back(p4est_cell);
721  }
722  }
723  else
724  {
727  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
728  ++c)
729  switch (dim)
730  {
731  case 2:
732  P4EST_QUADRANT_INIT(&p4est_child[c]);
733  break;
734  case 3:
735  P8EST_QUADRANT_INIT(&p4est_child[c]);
736  break;
737  default:
738  Assert(false, ExcNotImplemented());
739  }
741  p4est_child);
742  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
743  ++c)
744  {
745  p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
746  build_lists(cell->child(c), p4est_child[c], my_subdomain);
747  }
748  }
749  }
750 
751 
752  template <int dim, int spacedim>
753  int
754  RefineAndCoarsenList<dim, spacedim>::refine_callback(
755  typename internal::p4est::types<dim>::forest * forest,
756  typename internal::p4est::types<dim>::topidx coarse_cell_index,
757  typename internal::p4est::types<dim>::quadrant *quadrant)
758  {
759  RefineAndCoarsenList<dim, spacedim> *this_object =
760  reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *>(
761  forest->user_pointer);
762 
763  // if there are no more cells in our list the current cell can't be
764  // flagged for refinement
765  if (this_object->current_refine_pointer == this_object->refine_list.end())
766  return false;
767 
768  Assert(coarse_cell_index <=
769  this_object->current_refine_pointer->p.which_tree,
770  ExcInternalError());
771 
772  // if p4est hasn't yet reached the tree of the next flagged cell the
773  // current cell can't be flagged for refinement
774  if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
775  return false;
776 
777  // now we're in the right tree in the forest
778  Assert(coarse_cell_index <=
779  this_object->current_refine_pointer->p.which_tree,
780  ExcInternalError());
781 
782  // make sure that the p4est loop over cells hasn't gotten ahead of our own
783  // pointer
785  quadrant, &*this_object->current_refine_pointer) <= 0,
786  ExcInternalError());
787 
788  // now, if the p4est cell is one in the list, it is supposed to be refined
790  quadrant, &*this_object->current_refine_pointer))
791  {
792  ++this_object->current_refine_pointer;
793  return true;
794  }
795 
796  // p4est cell is not in list
797  return false;
798  }
799 
800 
801 
802  template <int dim, int spacedim>
803  int
804  RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
805  typename internal::p4est::types<dim>::forest * forest,
806  typename internal::p4est::types<dim>::topidx coarse_cell_index,
807  typename internal::p4est::types<dim>::quadrant *children[])
808  {
809  RefineAndCoarsenList<dim, spacedim> *this_object =
810  reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *>(
811  forest->user_pointer);
812 
813  // if there are no more cells in our list the current cell can't be
814  // flagged for coarsening
815  if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
816  return false;
817 
818  Assert(coarse_cell_index <=
819  this_object->current_coarsen_pointer->p.which_tree,
820  ExcInternalError());
821 
822  // if p4est hasn't yet reached the tree of the next flagged cell the
823  // current cell can't be flagged for coarsening
824  if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
825  return false;
826 
827  // now we're in the right tree in the forest
828  Assert(coarse_cell_index <=
829  this_object->current_coarsen_pointer->p.which_tree,
830  ExcInternalError());
831 
832  // make sure that the p4est loop over cells hasn't gotten ahead of our own
833  // pointer
835  children[0], &*this_object->current_coarsen_pointer) <= 0,
836  ExcInternalError());
837 
838  // now, if the p4est cell is one in the list, it is supposed to be
839  // coarsened
841  children[0], &*this_object->current_coarsen_pointer))
842  {
843  // move current pointer one up
844  ++this_object->current_coarsen_pointer;
845 
846  // note that the next 3 cells in our list need to correspond to the
847  // other siblings of the cell we have just found
848  for (unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
849  ++c)
850  {
852  children[c], &*this_object->current_coarsen_pointer),
853  ExcInternalError());
854  ++this_object->current_coarsen_pointer;
855  }
856 
857  return true;
858  }
859 
860  // p4est cell is not in list
861  return false;
862  }
863 
864 
865 
872  template <int dim, int spacedim>
873  class PartitionWeights
874  {
875  public:
881  explicit PartitionWeights(const std::vector<unsigned int> &cell_weights);
882 
890  static int
891  cell_weight(typename internal::p4est::types<dim>::forest *forest,
892  typename internal::p4est::types<dim>::topidx coarse_cell_index,
893  typename internal::p4est::types<dim>::quadrant *quadrant);
894 
895  private:
896  std::vector<unsigned int> cell_weights_list;
897  std::vector<unsigned int>::const_iterator current_pointer;
898  };
899 
900 
901  template <int dim, int spacedim>
902  PartitionWeights<dim, spacedim>::PartitionWeights(
903  const std::vector<unsigned int> &cell_weights)
904  : cell_weights_list(cell_weights)
905  {
906  // set the current pointer to the first element of the list, given that
907  // we will walk through it sequentially
908  current_pointer = cell_weights_list.begin();
909  }
910 
911 
912  template <int dim, int spacedim>
913  int
914  PartitionWeights<dim, spacedim>::cell_weight(
915  typename internal::p4est::types<dim>::forest *forest,
918  {
919  // the function gets two additional arguments, but we don't need them
920  // since we know in which order p4est will walk through the cells
921  // and have already built our weight lists in this order
922 
923  PartitionWeights<dim, spacedim> *this_object =
924  reinterpret_cast<PartitionWeights<dim, spacedim> *>(forest->user_pointer);
925 
926  Assert(this_object->current_pointer >=
927  this_object->cell_weights_list.begin(),
928  ExcInternalError());
929  Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
930  ExcInternalError());
931 
932  // get the weight, increment the pointer, and return the weight
933  return *this_object->current_pointer++;
934  }
935 
936  template <int dim, int spacedim>
937  using cell_relation_t = typename std::pair<
938  typename ::Triangulation<dim, spacedim>::cell_iterator,
939  typename ::Triangulation<dim, spacedim>::CellStatus>;
940 
950  template <int dim, int spacedim>
951  inline void
952  add_single_cell_relation(
953  std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
954  const typename ::internal::p4est::types<dim>::tree & tree,
955  const unsigned int idx,
956  const typename Triangulation<dim, spacedim>::cell_iterator &dealii_cell,
957  const typename Triangulation<dim, spacedim>::CellStatus status)
958  {
959  const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
960 
961  // check if we will be writing into valid memory
962  Assert(local_quadrant_index < cell_rel.size(), ExcInternalError());
963 
964  // store relation
965  cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
966  }
967 
968 
969 
979  template <int dim, int spacedim>
980  void
981  update_cell_relations_recursively(
982  std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
983  const typename ::internal::p4est::types<dim>::tree & tree,
984  const typename Triangulation<dim, spacedim>::cell_iterator & dealii_cell,
985  const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
986  {
987  // find index of p4est_cell in the quadrants array of the corresponding tree
988  const int idx = sc_array_bsearch(
989  const_cast<sc_array_t *>(&tree.quadrants),
990  &p4est_cell,
992  if (idx == -1 &&
994  const_cast<typename ::internal::p4est::types<dim>::tree *>(
995  &tree),
996  &p4est_cell) == false))
997  // this quadrant and none of its children belong to us.
998  return;
999 
1000  // recurse further if both p4est and dealii still have children
1001  const bool p4est_has_children = (idx == -1);
1002  if (p4est_has_children && dealii_cell->has_children())
1003  {
1004  // recurse further
1005  typename ::internal::p4est::types<dim>::quadrant
1007 
1008  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1009  ++c)
1010  switch (dim)
1011  {
1012  case 2:
1013  P4EST_QUADRANT_INIT(&p4est_child[c]);
1014  break;
1015  case 3:
1016  P8EST_QUADRANT_INIT(&p4est_child[c]);
1017  break;
1018  default:
1019  Assert(false, ExcNotImplemented());
1020  }
1021 
1023  &p4est_cell, p4est_child);
1024 
1025  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1026  ++c)
1027  {
1028  update_cell_relations_recursively<dim, spacedim>(
1029  cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1030  }
1031  }
1032  else if (!p4est_has_children && !dealii_cell->has_children())
1033  {
1034  // this active cell didn't change
1035  // save pair into corresponding position
1036  add_single_cell_relation<dim, spacedim>(
1037  cell_rel,
1038  tree,
1039  idx,
1040  dealii_cell,
1042  }
1043  else if (p4est_has_children) // based on the conditions above, we know that
1044  // dealii_cell has no children
1045  {
1046  // this cell got refined in p4est, but the dealii_cell has not yet been
1047  // refined
1048 
1049  // this quadrant is not active
1050  // generate its children, and store information in those
1051  typename ::internal::p4est::types<dim>::quadrant
1053  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1054  ++c)
1055  switch (dim)
1056  {
1057  case 2:
1058  P4EST_QUADRANT_INIT(&p4est_child[c]);
1059  break;
1060  case 3:
1061  P8EST_QUADRANT_INIT(&p4est_child[c]);
1062  break;
1063  default:
1064  Assert(false, ExcNotImplemented());
1065  }
1066 
1068  &p4est_cell, p4est_child);
1069 
1070  // mark first child with CELL_REFINE and the remaining children with
1071  // CELL_INVALID, but associate them all with the parent cell unpack
1072  // algorithm will be called only on CELL_REFINE flagged quadrant
1073  int child_idx;
1074  typename Triangulation<dim, spacedim>::CellStatus cell_status;
1075  for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1076  ++i)
1077  {
1078  child_idx = sc_array_bsearch(
1079  const_cast<sc_array_t *>(&tree.quadrants),
1080  &p4est_child[i],
1082 
1083  cell_status = (i == 0) ? Triangulation<dim, spacedim>::CELL_REFINE :
1085 
1086  add_single_cell_relation<dim, spacedim>(
1087  cell_rel, tree, child_idx, dealii_cell, cell_status);
1088  }
1089  }
1090  else // based on the conditions above, we know that p4est_cell has no
1091  // children, and the dealii_cell does
1092  {
1093  // its children got coarsened into this cell in p4est,
1094  // but the dealii_cell still has its children
1095  add_single_cell_relation<dim, spacedim>(
1096  cell_rel,
1097  tree,
1098  idx,
1099  dealii_cell,
1101  }
1102  }
1103 } // namespace
1104 
1105 
1106 
1107 namespace parallel
1108 {
1109  namespace distributed
1110  {
1111  /*----------------- class Triangulation<dim,spacedim> ---------------\*/
1112  template <int dim, int spacedim>
1114  const MPI_Comm &mpi_communicator,
1115  const typename ::Triangulation<dim, spacedim>::MeshSmoothing
1116  smooth_grid,
1117  const Settings settings)
1118  : // Do not check for distorted cells.
1119  // For multigrid, we need limit_level_difference_at_vertices
1120  // to make sure the transfer operators only need to consider two levels.
1121  ::parallel::DistributedTriangulationBase<dim, spacedim>(
1122  mpi_communicator,
1123  (settings & construct_multigrid_hierarchy) ?
1124  static_cast<
1125  typename ::Triangulation<dim, spacedim>::MeshSmoothing>(
1126  smooth_grid |
1127  Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1128  smooth_grid,
1129  false)
1130  , settings(settings)
1131  , triangulation_has_content(false)
1132  , connectivity(nullptr)
1133  , parallel_forest(nullptr)
1134  {
1135  parallel_ghost = nullptr;
1136  }
1137 
1138 
1139 
1140  template <int dim, int spacedim>
1142  {
1143  // virtual functions called in constructors and destructors never use the
1144  // override in a derived class
1145  // for clarity be explicit on which function is called
1146  try
1147  {
1149  }
1150  catch (...)
1151  {}
1152 
1153  AssertNothrow(triangulation_has_content == false, ExcInternalError());
1154  AssertNothrow(connectivity == nullptr, ExcInternalError());
1155  AssertNothrow(parallel_forest == nullptr, ExcInternalError());
1156  }
1157 
1158 
1159 
1160  template <int dim, int spacedim>
1161  void
1163  const std::vector<Point<spacedim>> &vertices,
1164  const std::vector<CellData<dim>> & cells,
1165  const SubCellData & subcelldata)
1166  {
1167  try
1168  {
1170  vertices, cells, subcelldata);
1171  }
1172  catch (
1173  const typename ::Triangulation<dim, spacedim>::DistortedCellList
1174  &)
1175  {
1176  // the underlying triangulation should not be checking for distorted
1177  // cells
1178  Assert(false, ExcInternalError());
1179  }
1180 
1181  Assert(
1183  ExcMessage(
1184  "The class parallel::distributed::Triangulation only supports meshes "
1185  "consisting only of hypercube-like cells."));
1186 
1187  // note that now we have some content in the p4est objects and call the
1188  // functions that do the actual work (which are dimension dependent, so
1189  // separate)
1190  triangulation_has_content = true;
1191 
1192  setup_coarse_cell_to_p4est_tree_permutation();
1193 
1194  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1195 
1196  try
1197  {
1198  copy_local_forest_to_triangulation();
1199  }
1200  catch (const typename Triangulation<dim>::DistortedCellList &)
1201  {
1202  // the underlying triangulation should not be checking for distorted
1203  // cells
1204  Assert(false, ExcInternalError());
1205  }
1206 
1207  this->update_periodic_face_map();
1208  this->update_number_cache();
1209  }
1210 
1211 
1212 
1213  template <int dim, int spacedim>
1214  void
1217  &construction_data)
1218  {
1219  (void)construction_data;
1220 
1221  Assert(false, ExcInternalError());
1222  }
1223 
1224 
1225 
1226  template <int dim, int spacedim>
1227  void
1229  {
1230  triangulation_has_content = false;
1231 
1232  if (parallel_ghost != nullptr)
1233  {
1235  parallel_ghost);
1236  parallel_ghost = nullptr;
1237  }
1238 
1239  if (parallel_forest != nullptr)
1240  {
1242  parallel_forest = nullptr;
1243  }
1244 
1245  if (connectivity != nullptr)
1246  {
1248  connectivity);
1249  connectivity = nullptr;
1250  }
1251 
1252  coarse_cell_to_p4est_tree_permutation.resize(0);
1253  p4est_tree_to_coarse_cell_permutation.resize(0);
1254 
1256 
1257  this->update_number_cache();
1258  }
1259 
1260 
1261 
1262  template <int dim, int spacedim>
1263  bool
1265  {
1266  return settings &
1268  }
1269 
1270 
1271 
1272  template <int dim, int spacedim>
1273  void
1275  const typename ::internal::p4est::types<dim>::forest
1276  *parallel_forest,
1277  const typename ::internal::p4est::types<dim>::gloidx
1278  *previous_global_first_quadrant)
1279  {
1280  Assert(this->data_transfer.sizes_fixed_cumulative.size() > 0,
1281  ExcMessage("No data has been packed!"));
1282 
1283  // Resize memory according to the data that we will receive.
1284  this->data_transfer.dest_data_fixed.resize(
1285  parallel_forest->local_num_quadrants *
1286  this->data_transfer.sizes_fixed_cumulative.back());
1287 
1288  // Execute non-blocking fixed size transfer.
1289  typename ::internal::p4est::types<dim>::transfer_context
1290  *tf_context;
1291  tf_context =
1293  parallel_forest->global_first_quadrant,
1294  previous_global_first_quadrant,
1295  parallel_forest->mpicomm,
1296  0,
1297  this->data_transfer.dest_data_fixed.data(),
1298  this->data_transfer.src_data_fixed.data(),
1299  this->data_transfer.sizes_fixed_cumulative.back());
1300 
1301  if (this->data_transfer.variable_size_data_stored)
1302  {
1303  // Resize memory according to the data that we will receive.
1304  this->data_transfer.dest_sizes_variable.resize(
1305  parallel_forest->local_num_quadrants);
1306 
1307  // Execute fixed size transfer of data sizes for variable size
1308  // transfer.
1310  parallel_forest->global_first_quadrant,
1311  previous_global_first_quadrant,
1312  parallel_forest->mpicomm,
1313  1,
1314  this->data_transfer.dest_sizes_variable.data(),
1315  this->data_transfer.src_sizes_variable.data(),
1316  sizeof(unsigned int));
1317  }
1318 
1320 
1321  // Release memory of previously packed data.
1322  this->data_transfer.src_data_fixed.clear();
1323  this->data_transfer.src_data_fixed.shrink_to_fit();
1324 
1325  if (this->data_transfer.variable_size_data_stored)
1326  {
1327  // Resize memory according to the data that we will receive.
1328  this->data_transfer.dest_data_variable.resize(
1329  std::accumulate(this->data_transfer.dest_sizes_variable.begin(),
1330  this->data_transfer.dest_sizes_variable.end(),
1332 
1333 # if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1334 # else
1335  // ----- WORKAROUND -----
1336  // An assertion in p4est prevents us from sending/receiving no data
1337  // at all, which is mandatory if one of our processes does not own
1338  // any quadrant. This bypasses the assertion from being triggered.
1339  // - see: https://github.com/cburstedde/p4est/issues/48
1340  if (this->data_transfer.src_sizes_variable.size() == 0)
1341  this->data_transfer.src_sizes_variable.resize(1);
1342  if (this->data_transfer.dest_sizes_variable.size() == 0)
1343  this->data_transfer.dest_sizes_variable.resize(1);
1344 # endif
1345 
1346  // Execute variable size transfer.
1348  parallel_forest->global_first_quadrant,
1349  previous_global_first_quadrant,
1350  parallel_forest->mpicomm,
1351  1,
1352  this->data_transfer.dest_data_variable.data(),
1353  this->data_transfer.dest_sizes_variable.data(),
1354  this->data_transfer.src_data_variable.data(),
1355  this->data_transfer.src_sizes_variable.data());
1356 
1357  // Release memory of previously packed data.
1358  this->data_transfer.src_sizes_variable.clear();
1359  this->data_transfer.src_sizes_variable.shrink_to_fit();
1360  this->data_transfer.src_data_variable.clear();
1361  this->data_transfer.src_data_variable.shrink_to_fit();
1362  }
1363  }
1364 
1365 
1366  template <int dim, int spacedim>
1367  bool
1369  {
1370  if (this->n_global_levels() <= 1)
1371  return false; // can not have hanging nodes without refined cells
1372 
1373  // if there are any active cells with level less than n_global_levels()-1,
1374  // then there is obviously also one with level n_global_levels()-1, and
1375  // consequently there must be a hanging node somewhere.
1376  //
1377  // The problem is that we cannot just ask for the first active cell, but
1378  // instead need to filter over locally owned cells.
1379  const bool have_coarser_cell =
1380  std::any_of(this->begin_active(this->n_global_levels() - 2),
1381  this->end_active(this->n_global_levels() - 2),
1382  [](const CellAccessor<dim, spacedim> &cell) {
1383  return cell.is_locally_owned();
1384  });
1385 
1386  // return true if at least one process has a coarser cell
1387  return 0 < Utilities::MPI::max(have_coarser_cell ? 1 : 0,
1388  this->mpi_communicator);
1389  }
1390 
1391 
1392 
1393  template <int dim, int spacedim>
1394  void
1396  {
1397  DynamicSparsityPattern cell_connectivity;
1399  cell_connectivity);
1400  coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
1402  cell_connectivity, coarse_cell_to_p4est_tree_permutation);
1403 
1404  p4est_tree_to_coarse_cell_permutation =
1405  Utilities::invert_permutation(coarse_cell_to_p4est_tree_permutation);
1406  }
1407 
1408 
1409 
1410  template <int dim, int spacedim>
1411  void
1413  const std::string &file_basename) const
1414  {
1415  Assert(parallel_forest != nullptr,
1416  ExcMessage("Can't produce output when no forest is created yet."));
1418  parallel_forest, nullptr, file_basename.c_str());
1419  }
1420 
1421 
1422 
1423  template <int dim, int spacedim>
1424  void
1425  Triangulation<dim, spacedim>::save(const std::string &filename) const
1426  {
1427  Assert(
1428  this->cell_attached_data.n_attached_deserialize == 0,
1429  ExcMessage(
1430  "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
1431  Assert(this->n_cells() > 0,
1432  ExcMessage("Can not save() an empty Triangulation."));
1433 
1434  const int myrank =
1435  Utilities::MPI::this_mpi_process(this->mpi_communicator);
1436 
1437  // signal that serialization is going to happen
1438  this->signals.pre_distributed_save();
1439 
1440  if (this->my_subdomain == 0)
1441  {
1442  std::string fname = std::string(filename) + ".info";
1443  std::ofstream f(fname.c_str());
1444  f << "version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
1445  << std::endl
1446  << 4 << " "
1447  << Utilities::MPI::n_mpi_processes(this->mpi_communicator) << " "
1448  << this->cell_attached_data.pack_callbacks_fixed.size() << " "
1449  << this->cell_attached_data.pack_callbacks_variable.size() << " "
1450  << this->n_cells(0) << std::endl;
1451  }
1452 
1453  // each cell should have been flagged `CELL_PERSIST`
1454  for (const auto &cell_rel : this->local_cell_relations)
1455  {
1456  (void)cell_rel;
1457  Assert(
1458  (cell_rel.second == // cell_status
1460  ExcInternalError());
1461  }
1462 
1463  // Save cell attached data.
1464  this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
1465  parallel_forest->global_num_quadrants,
1466  filename);
1467 
1468  ::internal::p4est::functions<dim>::save(filename.c_str(),
1469  parallel_forest,
1470  false);
1471 
1472  // signal that serialization has finished
1474  }
1475 
1476 
1477 
1478  template <int dim, int spacedim>
1479  void
1480  Triangulation<dim, spacedim>::load(const std::string &filename,
1481  const bool autopartition)
1482  {
1483  Assert(
1484  this->n_cells() > 0,
1485  ExcMessage(
1486  "load() only works if the Triangulation already contains a coarse mesh!"));
1487  Assert(
1488  this->n_levels() == 1,
1489  ExcMessage(
1490  "Triangulation may only contain coarse cells when calling load()."));
1491 
1492  const int myrank =
1493  Utilities::MPI::this_mpi_process(this->mpi_communicator);
1494 
1495  // signal that de-serialization is going to happen
1496  this->signals.pre_distributed_load();
1497 
1498  if (parallel_ghost != nullptr)
1499  {
1501  parallel_ghost);
1502  parallel_ghost = nullptr;
1503  }
1505  parallel_forest = nullptr;
1507  connectivity);
1508  connectivity = nullptr;
1509 
1510  unsigned int version, numcpus, attached_count_fixed,
1511  attached_count_variable, n_coarse_cells;
1512  {
1513  std::string fname = std::string(filename) + ".info";
1514  std::ifstream f(fname.c_str());
1515  AssertThrow(f, ExcIO());
1516  std::string firstline;
1517  getline(f, firstline); // skip first line
1518  f >> version >> numcpus >> attached_count_fixed >>
1519  attached_count_variable >> n_coarse_cells;
1520  }
1521 
1522  AssertThrow(version == 4,
1523  ExcMessage("Incompatible version found in .info file."));
1524  Assert(this->n_cells(0) == n_coarse_cells,
1525  ExcMessage("Number of coarse cells differ!"));
1526 
1527  // clear all of the callback data, as explained in the documentation of
1528  // register_data_attach()
1529  this->cell_attached_data.n_attached_data_sets = 0;
1530  this->cell_attached_data.n_attached_deserialize =
1531  attached_count_fixed + attached_count_variable;
1532 
1534  filename.c_str(),
1535  this->mpi_communicator,
1536  0,
1537  false,
1538  autopartition,
1539  0,
1540  this,
1541  &connectivity);
1542 
1543  if (numcpus != Utilities::MPI::n_mpi_processes(this->mpi_communicator))
1544  {
1545  // We are changing the number of CPUs so we need to repartition.
1546  // Note that p4est actually distributes the cells between the changed
1547  // number of CPUs and so everything works without this call, but
1548  // this command changes the distribution for some reason, so we
1549  // will leave it in here.
1550  if (this->signals.cell_weight.num_slots() == 0)
1551  {
1552  // no cell weights given -- call p4est's 'partition' without a
1553  // callback for cell weights
1555  parallel_forest,
1556  /* prepare coarsening */ 1,
1557  /* weight_callback */ nullptr);
1558  }
1559  else
1560  {
1561  // get cell weights for a weighted repartitioning.
1562  const std::vector<unsigned int> cell_weights = get_cell_weights();
1563 
1564  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
1565 
1566  // attach (temporarily) a pointer to the cell weights through
1567  // p4est's user_pointer object
1568  Assert(parallel_forest->user_pointer == this, ExcInternalError());
1569  parallel_forest->user_pointer = &partition_weights;
1570 
1572  parallel_forest,
1573  /* prepare coarsening */ 1,
1574  /* weight_callback */
1575  &PartitionWeights<dim, spacedim>::cell_weight);
1576 
1577  // reset the user pointer to its previous state
1578  parallel_forest->user_pointer = this;
1579  }
1580  }
1581 
1582  try
1583  {
1584  copy_local_forest_to_triangulation();
1585  }
1586  catch (const typename Triangulation<dim>::DistortedCellList &)
1587  {
1588  // the underlying
1589  // triangulation should not
1590  // be checking for
1591  // distorted cells
1592  Assert(false, ExcInternalError());
1593  }
1594 
1595  // Load attached cell data, if any was stored.
1596  this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
1597  parallel_forest->global_num_quadrants,
1598  parallel_forest->local_num_quadrants,
1599  filename,
1600  attached_count_fixed,
1601  attached_count_variable);
1602 
1603  // signal that de-serialization is finished
1605 
1606  this->update_periodic_face_map();
1607  this->update_number_cache();
1608  }
1609 
1610 
1611 
1612  template <int dim, int spacedim>
1613  void
1615  const typename ::internal::p4est::types<dim>::forest *forest)
1616  {
1617  Assert(this->n_cells() > 0,
1618  ExcMessage(
1619  "load() only works if the Triangulation already contains "
1620  "a coarse mesh!"));
1621  Assert(this->n_cells() == forest->trees->elem_count,
1622  ExcMessage(
1623  "Coarse mesh of the Triangulation does not match the one "
1624  "of the provided forest!"));
1625 
1626  // clear the old forest
1627  if (parallel_ghost != nullptr)
1628  {
1630  parallel_ghost);
1631  parallel_ghost = nullptr;
1632  }
1634  parallel_forest = nullptr;
1635 
1636  // note: we can keep the connectivity, since the coarse grid does not
1637  // change
1638 
1639  // create deep copy of the new forest
1640  typename ::internal::p4est::types<dim>::forest *temp =
1641  const_cast<typename ::internal::p4est::types<dim>::forest *>(
1642  forest);
1643  parallel_forest =
1645  parallel_forest->connectivity = connectivity;
1646  parallel_forest->user_pointer = this;
1647 
1648  try
1649  {
1650  copy_local_forest_to_triangulation();
1651  }
1652  catch (const typename Triangulation<dim>::DistortedCellList &)
1653  {
1654  // the underlying triangulation should not be checking for distorted
1655  // cells
1656  Assert(false, ExcInternalError());
1657  }
1658 
1659  this->update_periodic_face_map();
1660  this->update_number_cache();
1661  }
1662 
1663 
1664 
1665  template <int dim, int spacedim>
1666  unsigned int
1668  {
1669  Assert(parallel_forest != nullptr,
1670  ExcMessage(
1671  "Can't produce a check sum when no forest is created yet."));
1672  return ::internal::p4est::functions<dim>::checksum(parallel_forest);
1673  }
1674 
1675 
1676 
1677  template <int dim, int spacedim>
1678  const typename ::internal::p4est::types<dim>::forest *
1680  {
1681  Assert(parallel_forest != nullptr,
1682  ExcMessage("The forest has not been allocated yet."));
1683  return parallel_forest;
1684  }
1685 
1686 
1687 
1688  template <int dim, int spacedim>
1689  typename ::internal::p4est::types<dim>::tree *
1691  const int dealii_coarse_cell_index) const
1692  {
1693  const unsigned int tree_index =
1694  coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
1695  typename ::internal::p4est::types<dim>::tree *tree =
1696  static_cast<typename ::internal::p4est::types<dim>::tree *>(
1697  sc_array_index(parallel_forest->trees, tree_index));
1698 
1699  return tree;
1700  }
1701 
1702 
1703 
1704  // Note: this has been added here to prevent that these functions
1705  // appear in the Doxygen documentation of ::Triangulation
1706 # ifndef DOXYGEN
1707 
1708  template <>
1709  void
1711  std::integral_constant<int, 2>)
1712  {
1713  const unsigned int dim = 2, spacedim = 2;
1714  Assert(this->n_cells(0) > 0, ExcInternalError());
1715  Assert(this->n_levels() == 1, ExcInternalError());
1716 
1717  // data structures that counts how many cells touch each vertex
1718  // (vertex_touch_count), and which cells touch a given vertex (together
1719  // with the local numbering of that vertex within the cells that touch
1720  // it)
1721  std::vector<unsigned int> vertex_touch_count;
1722  std::vector<
1723  std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
1724  unsigned int>>>
1725  vertex_to_cell;
1726  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
1727  const ::internal::p4est::types<2>::locidx num_vtt =
1728  std::accumulate(vertex_touch_count.begin(),
1729  vertex_touch_count.end(),
1730  0u);
1731 
1732  // now create a connectivity object with the right sizes for all
1733  // arrays. set vertex information only in debug mode (saves a few bytes
1734  // in optimized mode)
1735  const bool set_vertex_info
1736 # ifdef DEBUG
1737  = true
1738 # else
1739  = false
1740 # endif
1741  ;
1742 
1744  (set_vertex_info == true ? this->n_vertices() : 0),
1745  this->n_cells(0),
1746  this->n_vertices(),
1747  num_vtt);
1748 
1749  set_vertex_and_cell_info(*this,
1750  vertex_touch_count,
1751  vertex_to_cell,
1752  coarse_cell_to_p4est_tree_permutation,
1753  set_vertex_info,
1754  connectivity);
1755 
1756  Assert(p4est_connectivity_is_valid(connectivity) == 1,
1757  ExcInternalError());
1758 
1759  // now create a forest out of the connectivity data structure
1761  this->mpi_communicator,
1762  connectivity,
1763  /* minimum initial number of quadrants per tree */ 0,
1764  /* minimum level of upfront refinement */ 0,
1765  /* use uniform upfront refinement */ 1,
1766  /* user_data_size = */ 0,
1767  /* user_data_constructor = */ nullptr,
1768  /* user_pointer */ this);
1769  }
1770 
1771 
1772 
1773  // TODO: This is a verbatim copy of the 2,2 case. However, we can't just
1774  // specialize the dim template argument, but let spacedim open
1775  template <>
1776  void
1778  std::integral_constant<int, 2>)
1779  {
1780  const unsigned int dim = 2, spacedim = 3;
1781  Assert(this->n_cells(0) > 0, ExcInternalError());
1782  Assert(this->n_levels() == 1, ExcInternalError());
1783 
1784  // data structures that counts how many cells touch each vertex
1785  // (vertex_touch_count), and which cells touch a given vertex (together
1786  // with the local numbering of that vertex within the cells that touch
1787  // it)
1788  std::vector<unsigned int> vertex_touch_count;
1789  std::vector<
1790  std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
1791  unsigned int>>>
1792  vertex_to_cell;
1793  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
1794  const ::internal::p4est::types<2>::locidx num_vtt =
1795  std::accumulate(vertex_touch_count.begin(),
1796  vertex_touch_count.end(),
1797  0u);
1798 
1799  // now create a connectivity object with the right sizes for all
1800  // arrays. set vertex information only in debug mode (saves a few bytes
1801  // in optimized mode)
1802  const bool set_vertex_info
1803 # ifdef DEBUG
1804  = true
1805 # else
1806  = false
1807 # endif
1808  ;
1809 
1811  (set_vertex_info == true ? this->n_vertices() : 0),
1812  this->n_cells(0),
1813  this->n_vertices(),
1814  num_vtt);
1815 
1816  set_vertex_and_cell_info(*this,
1817  vertex_touch_count,
1818  vertex_to_cell,
1819  coarse_cell_to_p4est_tree_permutation,
1820  set_vertex_info,
1821  connectivity);
1822 
1823  Assert(p4est_connectivity_is_valid(connectivity) == 1,
1824  ExcInternalError());
1825 
1826  // now create a forest out of the connectivity data structure
1828  this->mpi_communicator,
1829  connectivity,
1830  /* minimum initial number of quadrants per tree */ 0,
1831  /* minimum level of upfront refinement */ 0,
1832  /* use uniform upfront refinement */ 1,
1833  /* user_data_size = */ 0,
1834  /* user_data_constructor = */ nullptr,
1835  /* user_pointer */ this);
1836  }
1837 
1838 
1839 
1840  template <>
1841  void
1843  std::integral_constant<int, 3>)
1844  {
1845  const int dim = 3, spacedim = 3;
1846  Assert(this->n_cells(0) > 0, ExcInternalError());
1847  Assert(this->n_levels() == 1, ExcInternalError());
1848 
1849  // data structures that counts how many cells touch each vertex
1850  // (vertex_touch_count), and which cells touch a given vertex (together
1851  // with the local numbering of that vertex within the cells that touch
1852  // it)
1853  std::vector<unsigned int> vertex_touch_count;
1854  std::vector<std::list<
1855  std::pair<Triangulation<3>::active_cell_iterator, unsigned int>>>
1856  vertex_to_cell;
1857  get_vertex_to_cell_mappings(*this, vertex_touch_count, vertex_to_cell);
1858  const ::internal::p4est::types<2>::locidx num_vtt =
1859  std::accumulate(vertex_touch_count.begin(),
1860  vertex_touch_count.end(),
1861  0u);
1862 
1863  std::vector<unsigned int> edge_touch_count;
1864  std::vector<std::list<
1865  std::pair<Triangulation<3>::active_cell_iterator, unsigned int>>>
1866  edge_to_cell;
1867  get_edge_to_cell_mappings(*this, edge_touch_count, edge_to_cell);
1868  const ::internal::p4est::types<2>::locidx num_ett =
1869  std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
1870 
1871  // now create a connectivity object with the right sizes for all arrays
1872  const bool set_vertex_info
1873 # ifdef DEBUG
1874  = true
1875 # else
1876  = false
1877 # endif
1878  ;
1879 
1881  (set_vertex_info == true ? this->n_vertices() : 0),
1882  this->n_cells(0),
1883  this->n_active_lines(),
1884  num_ett,
1885  this->n_vertices(),
1886  num_vtt);
1887 
1888  set_vertex_and_cell_info(*this,
1889  vertex_touch_count,
1890  vertex_to_cell,
1891  coarse_cell_to_p4est_tree_permutation,
1892  set_vertex_info,
1893  connectivity);
1894 
1895  // next to tree-to-edge
1896  // data. note that in p4est lines
1897  // are ordered as follows
1898  // *---3---* *---3---*
1899  // /| | / /|
1900  // 6 | 11 6 7 11
1901  // / 10 | / / |
1902  // * | | *---2---* |
1903  // | *---1---* | | *
1904  // | / / | 9 /
1905  // 8 4 5 8 | 5
1906  // |/ / | |/
1907  // *---0---* *---0---*
1908  // whereas in deal.II they are like this:
1909  // *---7---* *---7---*
1910  // /| | / /|
1911  // 4 | 11 4 5 11
1912  // / 10 | / / |
1913  // * | | *---6---* |
1914  // | *---3---* | | *
1915  // | / / | 9 /
1916  // 8 0 1 8 | 1
1917  // |/ / | |/
1918  // *---2---* *---2---*
1919 
1920  const unsigned int deal_to_p4est_line_index[12] = {
1921  4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
1922 
1924  this->begin_active();
1925  cell != this->end();
1926  ++cell)
1927  {
1928  const unsigned int index =
1929  coarse_cell_to_p4est_tree_permutation[cell->index()];
1930  for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++e)
1931  connectivity->tree_to_edge[index * GeometryInfo<3>::lines_per_cell +
1932  deal_to_p4est_line_index[e]] =
1933  cell->line(e)->index();
1934  }
1935 
1936  // now also set edge-to-tree
1937  // information
1938  connectivity->ett_offset[0] = 0;
1939  std::partial_sum(edge_touch_count.begin(),
1940  edge_touch_count.end(),
1941  &connectivity->ett_offset[1]);
1942 
1943  Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
1944  ExcInternalError());
1945 
1946  for (unsigned int v = 0; v < this->n_active_lines(); ++v)
1947  {
1948  Assert(edge_to_cell[v].size() == edge_touch_count[v],
1949  ExcInternalError());
1950 
1951  std::list<
1952  std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
1953  unsigned int>>::const_iterator p =
1954  edge_to_cell[v].begin();
1955  for (unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
1956  {
1957  connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
1958  coarse_cell_to_p4est_tree_permutation[p->first->index()];
1959  connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
1960  deal_to_p4est_line_index[p->second];
1961  }
1962  }
1963 
1964  Assert(p8est_connectivity_is_valid(connectivity) == 1,
1965  ExcInternalError());
1966 
1967  // now create a forest out of the connectivity data structure
1969  this->mpi_communicator,
1970  connectivity,
1971  /* minimum initial number of quadrants per tree */ 0,
1972  /* minimum level of upfront refinement */ 0,
1973  /* use uniform upfront refinement */ 1,
1974  /* user_data_size = */ 0,
1975  /* user_data_constructor = */ nullptr,
1976  /* user_pointer */ this);
1977  }
1978 # endif
1979 
1980 
1981 
1982  namespace
1983  {
1984  // ensures the 2:1 mesh balance for periodic boundary conditions in the
1985  // artificial cell layer (the active cells are taken care of by p4est)
1986  template <int dim, int spacedim>
1987  bool
1988  enforce_mesh_balance_over_periodic_boundaries(
1990  {
1991  if (tria.get_periodic_face_map().size() == 0)
1992  return false;
1993 
1994  std::vector<bool> flags_before[2];
1995  tria.save_coarsen_flags(flags_before[0]);
1996  tria.save_refine_flags(flags_before[1]);
1997 
1998  std::vector<unsigned int> topological_vertex_numbering(
1999  tria.n_vertices());
2000  for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2001  topological_vertex_numbering[i] = i;
2002  // combine vertices that have different locations (and thus, different
2003  // vertex_index) but represent the same topological entity over periodic
2004  // boundaries. The vector topological_vertex_numbering contains a linear
2005  // map from 0 to n_vertices at input and at output relates periodic
2006  // vertices with only one vertex index. The output is used to always
2007  // identify the same vertex according to the periodicity, e.g. when
2008  // finding the maximum cell level around a vertex.
2009  //
2010  // Example: On a 3D cell with vertices numbered from 0 to 7 and periodic
2011  // boundary conditions in x direction, the vector
2012  // topological_vertex_numbering will contain the numbers
2013  // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5},
2014  // {6,7} belong together, respectively). If periodicity is set in x and
2015  // z direction, the output is {0,0,2,2,0,0,2,2}, and if periodicity is
2016  // in all directions, the output is simply {0,0,0,0,0,0,0,0}.
2017  using cell_iterator =
2019  typename std::map<std::pair<cell_iterator, unsigned int>,
2020  std::pair<std::pair<cell_iterator, unsigned int>,
2021  std::bitset<3>>>::const_iterator it;
2022  for (it = tria.get_periodic_face_map().begin();
2023  it != tria.get_periodic_face_map().end();
2024  ++it)
2025  {
2026  const cell_iterator &cell_1 = it->first.first;
2027  const unsigned int face_no_1 = it->first.second;
2028  const cell_iterator &cell_2 = it->second.first.first;
2029  const unsigned int face_no_2 = it->second.first.second;
2030  const std::bitset<3> face_orientation = it->second.second;
2031 
2032  if (cell_1->level() == cell_2->level())
2033  {
2034  for (unsigned int v = 0;
2035  v < GeometryInfo<dim - 1>::vertices_per_cell;
2036  ++v)
2037  {
2038  // take possible non-standard orientation of face on cell[0]
2039  // into account
2040  const unsigned int vface0 =
2042  v,
2043  face_orientation[0],
2044  face_orientation[1],
2045  face_orientation[2]);
2046  const unsigned int vi0 =
2047  topological_vertex_numbering[cell_1->face(face_no_1)
2048  ->vertex_index(vface0)];
2049  const unsigned int vi1 =
2050  topological_vertex_numbering[cell_2->face(face_no_2)
2051  ->vertex_index(v)];
2052  const unsigned int min_index = std::min(vi0, vi1);
2053  topological_vertex_numbering[cell_1->face(face_no_1)
2054  ->vertex_index(vface0)] =
2055  topological_vertex_numbering[cell_2->face(face_no_2)
2056  ->vertex_index(v)] =
2057  min_index;
2058  }
2059  }
2060  }
2061 
2062  // There must not be any chains!
2063  for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2064  {
2065  const unsigned int j = topological_vertex_numbering[i];
2066  if (j != i)
2067  Assert(topological_vertex_numbering[j] == j, ExcInternalError());
2068  }
2069 
2070 
2071  // this code is replicated from grid/tria.cc but using an indirection
2072  // for periodic boundary conditions
2073  bool continue_iterating = true;
2074  std::vector<int> vertex_level(tria.n_vertices());
2075  while (continue_iterating)
2076  {
2077  // store highest level one of the cells adjacent to a vertex
2078  // belongs to
2079  std::fill(vertex_level.begin(), vertex_level.end(), 0);
2081  cell = tria.begin_active(),
2082  endc = tria.end();
2083  for (; cell != endc; ++cell)
2084  {
2085  if (cell->refine_flag_set())
2086  for (const unsigned int vertex :
2088  vertex_level[topological_vertex_numbering
2089  [cell->vertex_index(vertex)]] =
2090  std::max(vertex_level[topological_vertex_numbering
2091  [cell->vertex_index(vertex)]],
2092  cell->level() + 1);
2093  else if (!cell->coarsen_flag_set())
2094  for (const unsigned int vertex :
2096  vertex_level[topological_vertex_numbering
2097  [cell->vertex_index(vertex)]] =
2098  std::max(vertex_level[topological_vertex_numbering
2099  [cell->vertex_index(vertex)]],
2100  cell->level());
2101  else
2102  {
2103  // if coarsen flag is set then tentatively assume
2104  // that the cell will be coarsened. this isn't
2105  // always true (the coarsen flag could be removed
2106  // again) and so we may make an error here. we try
2107  // to correct this by iterating over the entire
2108  // process until we are converged
2109  Assert(cell->coarsen_flag_set(), ExcInternalError());
2110  for (const unsigned int vertex :
2112  vertex_level[topological_vertex_numbering
2113  [cell->vertex_index(vertex)]] =
2114  std::max(vertex_level[topological_vertex_numbering
2115  [cell->vertex_index(vertex)]],
2116  cell->level() - 1);
2117  }
2118  }
2119 
2120  continue_iterating = false;
2121 
2122  // loop over all cells in reverse order. do so because we
2123  // can then update the vertex levels on the adjacent
2124  // vertices and maybe already flag additional cells in this
2125  // loop
2126  //
2127  // note that not only may we have to add additional
2128  // refinement flags, but we will also have to remove
2129  // coarsening flags on cells adjacent to vertices that will
2130  // see refinement
2131  for (cell = tria.last_active(); cell != endc; --cell)
2132  if (cell->refine_flag_set() == false)
2133  {
2134  for (const unsigned int vertex :
2136  if (vertex_level[topological_vertex_numbering
2137  [cell->vertex_index(vertex)]] >=
2138  cell->level() + 1)
2139  {
2140  // remove coarsen flag...
2141  cell->clear_coarsen_flag();
2142 
2143  // ...and if necessary also refine the current
2144  // cell, at the same time updating the level
2145  // information about vertices
2146  if (vertex_level[topological_vertex_numbering
2147  [cell->vertex_index(vertex)]] >
2148  cell->level() + 1)
2149  {
2150  cell->set_refine_flag();
2151  continue_iterating = true;
2152 
2153  for (const unsigned int v :
2155  vertex_level[topological_vertex_numbering
2156  [cell->vertex_index(v)]] =
2157  std::max(
2158  vertex_level[topological_vertex_numbering
2159  [cell->vertex_index(v)]],
2160  cell->level() + 1);
2161  }
2162 
2163  // continue and see whether we may, for example,
2164  // go into the inner 'if' above based on a
2165  // different vertex
2166  }
2167  }
2168 
2169  // clear coarsen flag if not all children were marked
2170  for (const auto &cell : tria.cell_iterators())
2171  {
2172  // nothing to do if we are already on the finest level
2173  if (cell->is_active())
2174  continue;
2175 
2176  const unsigned int n_children = cell->n_children();
2177  unsigned int flagged_children = 0;
2178  for (unsigned int child = 0; child < n_children; ++child)
2179  if (cell->child(child)->is_active() &&
2180  cell->child(child)->coarsen_flag_set())
2181  ++flagged_children;
2182 
2183  // if not all children were flagged for coarsening, remove
2184  // coarsen flags
2185  if (flagged_children < n_children)
2186  for (unsigned int child = 0; child < n_children; ++child)
2187  if (cell->child(child)->is_active())
2188  cell->child(child)->clear_coarsen_flag();
2189  }
2190  }
2191  std::vector<bool> flags_after[2];
2192  tria.save_coarsen_flags(flags_after[0]);
2193  tria.save_refine_flags(flags_after[1]);
2194  return ((flags_before[0] != flags_after[0]) ||
2195  (flags_before[1] != flags_after[1]));
2196  }
2197  } // namespace
2198 
2199 
2200 
2201  template <int dim, int spacedim>
2202  bool
2204  {
2205  std::vector<bool> flags_before[2];
2206  this->save_coarsen_flags(flags_before[0]);
2207  this->save_refine_flags(flags_before[1]);
2208 
2209  bool mesh_changed = false;
2210  unsigned int loop_counter = 0;
2211  do
2212  {
2215  this->update_periodic_face_map();
2216  // enforce 2:1 mesh balance over periodic boundaries
2217  mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this);
2218 
2219  // We can't be sure that we won't run into a situation where we can
2220  // not reconcile mesh smoothing and balancing of periodic faces. As we
2221  // don't know what else to do, at least abort with an error message.
2222  ++loop_counter;
2223  AssertThrow(
2224  loop_counter < 32,
2225  ExcMessage(
2226  "Infinite loop in "
2227  "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2228  "for periodic boundaries detected. Aborting."));
2229  }
2230  while (mesh_changed);
2231 
2232  // check if any of the refinement flags were changed during this
2233  // function and return that value
2234  std::vector<bool> flags_after[2];
2235  this->save_coarsen_flags(flags_after[0]);
2236  this->save_refine_flags(flags_after[1]);
2237  return ((flags_before[0] != flags_after[0]) ||
2238  (flags_before[1] != flags_after[1]));
2239  }
2240 
2241 
2242 
2243  template <int dim, int spacedim>
2244  void
2246  {
2247  // Disable mesh smoothing for recreating the deal.II triangulation,
2248  // otherwise we might not be able to reproduce the p4est mesh
2249  // exactly. We restore the original smoothing at the end of this
2250  // function. Note that the smoothing flag is used in the normal
2251  // refinement process.
2252  typename Triangulation<dim, spacedim>::MeshSmoothing save_smooth =
2253  this->smooth_grid;
2254 
2255  // We will refine manually to match the p4est further down, which
2256  // obeys a level difference of 2 at each vertex (see the balance call
2257  // to p4est). We can disable this here so we store fewer artificial
2258  // cells (in some cases).
2259  // For geometric multigrid it turns out that
2260  // we will miss level cells at shared vertices if we ignore this.
2261  // See tests/mpi/mg_06. In particular, the flag is still necessary
2262  // even though we force it for the original smooth_grid in the
2263  // constructor.
2264  if (settings & construct_multigrid_hierarchy)
2265  this->smooth_grid =
2266  ::Triangulation<dim,
2268  else
2270 
2271  bool mesh_changed = false;
2272 
2273  // Remove all deal.II refinements. Note that we could skip this and
2274  // start from our current state, because the algorithm later coarsens as
2275  // necessary. This has the advantage of being faster when large parts
2276  // of the local partition changes (likely) and gives a deterministic
2277  // ordering of the cells (useful for snapshot/resume).
2278  // TODO: is there a more efficient way to do this?
2279  if (settings & mesh_reconstruction_after_repartitioning)
2280  while (this->n_levels() > 1)
2281  {
2282  // Instead of marking all active cells, we slice off the finest
2283  // level, one level at a time. This takes the same number of
2284  // iterations but solves an issue where not all cells on a
2285  // periodic boundary are indeed coarsened and we run into an
2286  // irrelevant Assert() in update_periodic_face_map().
2287  for (const auto &cell :
2288  this->active_cell_iterators_on_level(this->n_levels() - 1))
2289  {
2290  cell->set_coarsen_flag();
2291  }
2292  try
2293  {
2296  }
2297  catch (
2299  {
2300  // the underlying triangulation should not be checking for
2301  // distorted cells
2302  Assert(false, ExcInternalError());
2303  }
2304  }
2305 
2306 
2307  // query p4est for the ghost cells
2308  if (parallel_ghost != nullptr)
2309  {
2311  parallel_ghost);
2312  parallel_ghost = nullptr;
2313  }
2315  parallel_forest,
2316  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2317  P4EST_CONNECT_CORNER) :
2318  typename ::internal::p4est::types<dim>::balance_type(
2319  P8EST_CONNECT_CORNER)));
2320 
2321  Assert(parallel_ghost, ExcInternalError());
2322 
2323 
2324  // set all cells to artificial. we will later set it to the correct
2325  // subdomain in match_tree_recursively
2326  for (const auto &cell : this->cell_iterators_on_level(0))
2327  cell->recursively_set_subdomain_id(numbers::artificial_subdomain_id);
2328 
2329  do
2330  {
2331  for (const auto &cell : this->cell_iterators_on_level(0))
2332  {
2333  // if this processor stores no part of the forest that comes out
2334  // of this coarse grid cell, then we need to delete all children
2335  // of this cell (the coarse grid cell remains)
2336  if (tree_exists_locally<dim, spacedim>(
2337  parallel_forest,
2338  coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2339  false)
2340  {
2341  delete_all_children<dim, spacedim>(cell);
2342  if (cell->is_active())
2343  cell->set_subdomain_id(numbers::artificial_subdomain_id);
2344  }
2345 
2346  else
2347  {
2348  // this processor stores at least a part of the tree that
2349  // comes out of this cell.
2350 
2351  typename ::internal::p4est::types<dim>::quadrant
2352  p4est_coarse_cell;
2353  typename ::internal::p4est::types<dim>::tree *tree =
2354  init_tree(cell->index());
2355 
2356  ::internal::p4est::init_coarse_quadrant<dim>(
2357  p4est_coarse_cell);
2358 
2359  match_tree_recursively<dim, spacedim>(*tree,
2360  cell,
2361  p4est_coarse_cell,
2362  *parallel_forest,
2363  this->my_subdomain);
2364  }
2365  }
2366 
2367  // check mesh for ghost cells, refine as necessary. iterate over
2368  // every ghostquadrant, find corresponding deal coarsecell and
2369  // recurse.
2370  typename ::internal::p4est::types<dim>::quadrant *quadr;
2371  types::subdomain_id ghost_owner = 0;
2372  typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2373 
2374  for (unsigned int g_idx = 0;
2375  g_idx < parallel_ghost->ghosts.elem_count;
2376  ++g_idx)
2377  {
2378  while (g_idx >= static_cast<unsigned int>(
2379  parallel_ghost->proc_offsets[ghost_owner + 1]))
2380  ++ghost_owner;
2381  while (g_idx >= static_cast<unsigned int>(
2382  parallel_ghost->tree_offsets[ghost_tree + 1]))
2383  ++ghost_tree;
2384 
2385  quadr = static_cast<
2386  typename ::internal::p4est::types<dim>::quadrant *>(
2387  sc_array_index(&parallel_ghost->ghosts, g_idx));
2388 
2389  unsigned int coarse_cell_index =
2390  p4est_tree_to_coarse_cell_permutation[ghost_tree];
2391 
2392  match_quadrant<dim, spacedim>(this,
2393  coarse_cell_index,
2394  *quadr,
2395  ghost_owner);
2396  }
2397 
2398  // fix all the flags to make sure we have a consistent mesh
2400 
2401  // see if any flags are still set
2402  mesh_changed =
2403  std::any_of(this->begin_active(),
2404  active_cell_iterator{this->end()},
2405  [](const CellAccessor<dim, spacedim> &cell) {
2406  return cell.refine_flag_set() ||
2407  cell.coarsen_flag_set();
2408  });
2409 
2410  // actually do the refinement to change the local mesh by
2411  // calling the base class refinement function directly
2412  try
2413  {
2416  }
2417  catch (
2419  {
2420  // the underlying triangulation should not be checking for
2421  // distorted cells
2422  Assert(false, ExcInternalError());
2423  }
2424  }
2425  while (mesh_changed);
2426 
2427 # ifdef DEBUG
2428  // check if correct number of ghosts is created
2429  unsigned int num_ghosts = 0;
2430 
2431  for (const auto &cell : this->active_cell_iterators())
2432  {
2433  if (cell->subdomain_id() != this->my_subdomain &&
2434  cell->subdomain_id() != numbers::artificial_subdomain_id)
2435  ++num_ghosts;
2436  }
2437 
2438  Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2439  ExcInternalError());
2440 # endif
2441 
2442 
2443 
2444  // fill level_subdomain_ids for geometric multigrid
2445  // the level ownership of a cell is defined as the owner if the cell is
2446  // active or as the owner of child(0) we need this information for all our
2447  // ancestors and the same-level neighbors of our own cells (=level ghosts)
2448  if (settings & construct_multigrid_hierarchy)
2449  {
2450  // step 1: We set our own ids all the way down and all the others to
2451  // -1. Note that we do not fill other cells we could figure out the
2452  // same way, because we might accidentally set an id for a cell that
2453  // is not a ghost cell.
2454  for (unsigned int lvl = this->n_levels(); lvl > 0;)
2455  {
2456  --lvl;
2457  for (const auto &cell : this->cell_iterators_on_level(lvl))
2458  {
2459  if ((cell->is_active() &&
2460  cell->subdomain_id() ==
2461  this->locally_owned_subdomain()) ||
2462  (cell->has_children() &&
2463  cell->child(0)->level_subdomain_id() ==
2464  this->locally_owned_subdomain()))
2465  cell->set_level_subdomain_id(
2466  this->locally_owned_subdomain());
2467  else
2468  {
2469  // not our cell
2470  cell->set_level_subdomain_id(
2472  }
2473  }
2474  }
2475 
2476  // step 2: make sure all the neighbors to our level_cells exist. Need
2477  // to look up in p4est...
2478  std::vector<std::vector<bool>> marked_vertices(this->n_levels());
2479  for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
2480  marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
2481 
2482  for (const auto &cell : this->cell_iterators_on_level(0))
2483  {
2484  typename ::internal::p4est::types<dim>::quadrant
2485  p4est_coarse_cell;
2486  const unsigned int tree_index =
2487  coarse_cell_to_p4est_tree_permutation[cell->index()];
2488  typename ::internal::p4est::types<dim>::tree *tree =
2489  init_tree(cell->index());
2490 
2491  ::internal::p4est::init_coarse_quadrant<dim>(
2492  p4est_coarse_cell);
2493 
2494  determine_level_subdomain_id_recursively<dim, spacedim>(
2495  *tree,
2496  tree_index,
2497  cell,
2498  p4est_coarse_cell,
2499  *parallel_forest,
2500  this->my_subdomain,
2501  marked_vertices);
2502  }
2503 
2504  // step 3: make sure we have the parent of our level cells
2505  for (unsigned int lvl = this->n_levels(); lvl > 0;)
2506  {
2507  --lvl;
2508  for (const auto &cell : this->cell_iterators_on_level(lvl))
2509  {
2510  if (cell->has_children())
2511  for (unsigned int c = 0;
2512  c < GeometryInfo<dim>::max_children_per_cell;
2513  ++c)
2514  {
2515  if (cell->child(c)->level_subdomain_id() ==
2516  this->locally_owned_subdomain())
2517  {
2518  // at least one of the children belongs to us, so
2519  // make sure we set the level subdomain id
2520  const types::subdomain_id mark =
2521  cell->child(0)->level_subdomain_id();
2523  ExcInternalError()); // we should know the
2524  // child(0)
2525  cell->set_level_subdomain_id(mark);
2526  break;
2527  }
2528  }
2529  }
2530  }
2531  }
2532 
2533 
2534 
2535  // check that our local copy has exactly as many cells as the p4est
2536  // original (at least if we are on only one processor); for parallel
2537  // computations, we want to check that we have at least as many as p4est
2538  // stores locally (in the future we should check that we have exactly as
2539  // many non-artificial cells as parallel_forest->local_num_quadrants)
2540  {
2541  const unsigned int total_local_cells = this->n_active_cells();
2542  (void)total_local_cells;
2543 
2544  if (Utilities::MPI::n_mpi_processes(this->mpi_communicator) == 1)
2545  {
2546  Assert(static_cast<unsigned int>(
2547  parallel_forest->local_num_quadrants) == total_local_cells,
2548  ExcInternalError());
2549  }
2550  else
2551  {
2552  Assert(static_cast<unsigned int>(
2553  parallel_forest->local_num_quadrants) <= total_local_cells,
2554  ExcInternalError());
2555  }
2556 
2557  // count the number of owned, active cells and compare with p4est.
2558  unsigned int n_owned = 0;
2559  for (const auto &cell : this->active_cell_iterators())
2560  {
2561  if (cell->subdomain_id() == this->my_subdomain)
2562  ++n_owned;
2563  }
2564 
2565  Assert(static_cast<unsigned int>(
2566  parallel_forest->local_num_quadrants) == n_owned,
2567  ExcInternalError());
2568  }
2569 
2570  this->smooth_grid = save_smooth;
2571 
2572  // finally, after syncing the parallel_forest with the triangulation,
2573  // also update the cell_relations, which will be used for
2574  // repartitioning, further refinement/coarsening, and unpacking
2575  // of stored or transferred data.
2576  update_cell_relations();
2577  }
2578 
2579 
2580 
2581  template <int dim, int spacedim>
2582  void
2584  {
2585  // do not allow anisotropic refinement
2586 # ifdef DEBUG
2587  for (const auto &cell : this->active_cell_iterators())
2588  if (cell->is_locally_owned() && cell->refine_flag_set())
2589  Assert(cell->refine_flag_set() ==
2591  ExcMessage(
2592  "This class does not support anisotropic refinement"));
2593 # endif
2594 
2595 
2596  // safety check: p4est has an upper limit on the level of a cell
2597  if (this->n_levels() ==
2599  {
2601  cell = this->begin_active(
2603  cell !=
2605  1);
2606  ++cell)
2607  {
2608  AssertThrow(
2609  !(cell->refine_flag_set()),
2610  ExcMessage(
2611  "Fatal Error: maximum refinement level of p4est reached."));
2612  }
2613  }
2614 
2616 
2617  // signal that refinement is going to happen
2619 
2620  // now do the work we're supposed to do when we are in charge
2621  // make sure all flags are cleared on cells we don't own, since nothing
2622  // good can come of that if they are still around
2623  for (const auto &cell : this->active_cell_iterators())
2624  if (cell->is_ghost() || cell->is_artificial())
2625  {
2626  cell->clear_refine_flag();
2627  cell->clear_coarsen_flag();
2628  }
2629 
2630 
2631  // count how many cells will be refined and coarsened, and allocate that
2632  // much memory
2633  RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
2634  *this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
2635 
2636  // copy refine and coarsen flags into p4est and execute the refinement
2637  // and coarsening. this uses the refine_and_coarsen_list just built,
2638  // which is communicated to the callback functions through
2639  // p4est's user_pointer object
2640  Assert(parallel_forest->user_pointer == this, ExcInternalError());
2641  parallel_forest->user_pointer = &refine_and_coarsen_list;
2642 
2643  if (parallel_ghost != nullptr)
2644  {
2646  parallel_ghost);
2647  parallel_ghost = nullptr;
2648  }
2650  parallel_forest,
2651  /* refine_recursive */ false,
2652  &RefineAndCoarsenList<dim, spacedim>::refine_callback,
2653  /*init_callback=*/nullptr);
2655  parallel_forest,
2656  /* coarsen_recursive */ false,
2657  &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
2658  /*init_callback=*/nullptr);
2659 
2660  // make sure all cells in the lists have been consumed
2661  Assert(refine_and_coarsen_list.pointers_are_at_end(), ExcInternalError());
2662 
2663  // reset the pointer
2664  parallel_forest->user_pointer = this;
2665 
2666  // enforce 2:1 hanging node condition
2668  parallel_forest,
2669  /* face and corner balance */
2670  (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2671  P4EST_CONNECT_FULL) :
2672  typename ::internal::p4est::types<dim>::balance_type(
2673  P8EST_CONNECT_FULL)),
2674  /*init_callback=*/nullptr);
2675 
2676  // since refinement and/or coarsening on the parallel forest
2677  // has happened, we need to update the quadrant cell relations
2678  update_cell_relations();
2679 
2680  // signals that parallel_forest has been refined and cell relations have
2681  // been updated
2683 
2684  // before repartitioning the mesh, save a copy of the current positions of
2685  // quadrants
2686  // only if data needs to be transferred later
2687  std::vector<typename ::internal::p4est::types<dim>::gloidx>
2688  previous_global_first_quadrant;
2689 
2690  if (this->cell_attached_data.n_attached_data_sets > 0)
2691  {
2692  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
2693  std::memcpy(previous_global_first_quadrant.data(),
2694  parallel_forest->global_first_quadrant,
2695  sizeof(
2696  typename ::internal::p4est::types<dim>::gloidx) *
2697  (parallel_forest->mpisize + 1));
2698  }
2699 
2700  if (!(settings & no_automatic_repartitioning))
2701  {
2702  // partition the new mesh between all processors. If cell weights have
2703  // not been given balance the number of cells.
2704  if (this->signals.cell_weight.num_slots() == 0)
2706  parallel_forest,
2707  /* prepare coarsening */ 1,
2708  /* weight_callback */ nullptr);
2709  else
2710  {
2711  // get cell weights for a weighted repartitioning.
2712  const std::vector<unsigned int> cell_weights = get_cell_weights();
2713 
2714  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
2715 
2716  // attach (temporarily) a pointer to the cell weights through
2717  // p4est's user_pointer object
2718  Assert(parallel_forest->user_pointer == this, ExcInternalError());
2719  parallel_forest->user_pointer = &partition_weights;
2720 
2722  parallel_forest,
2723  /* prepare coarsening */ 1,
2724  /* weight_callback */
2725  &PartitionWeights<dim, spacedim>::cell_weight);
2726 
2727  // release data
2729  parallel_forest, 0, nullptr, nullptr);
2730  // reset the user pointer to its previous state
2731  parallel_forest->user_pointer = this;
2732  }
2733  }
2734 
2735  // pack data before triangulation gets updated
2736  if (this->cell_attached_data.n_attached_data_sets > 0)
2737  {
2738  this->data_transfer.pack_data(
2739  this->local_cell_relations,
2740  this->cell_attached_data.pack_callbacks_fixed,
2741  this->cell_attached_data.pack_callbacks_variable);
2742  }
2743 
2744  // finally copy back from local part of tree to deal.II
2745  // triangulation. before doing so, make sure there are no refine or
2746  // coarsen flags pending
2747  for (const auto &cell : this->active_cell_iterators())
2748  {
2749  cell->clear_refine_flag();
2750  cell->clear_coarsen_flag();
2751  }
2752 
2753  try
2754  {
2755  copy_local_forest_to_triangulation();
2756  }
2757  catch (const typename Triangulation<dim>::DistortedCellList &)
2758  {
2759  // the underlying triangulation should not be checking for distorted
2760  // cells
2761  Assert(false, ExcInternalError());
2762  }
2763 
2764  // transfer data after triangulation got updated
2765  if (this->cell_attached_data.n_attached_data_sets > 0)
2766  {
2767  this->execute_transfer(parallel_forest,
2768  previous_global_first_quadrant.data());
2769 
2770  // also update the CellStatus information on the new mesh
2771  this->data_transfer.unpack_cell_status(this->local_cell_relations);
2772  }
2773 
2774 # ifdef DEBUG
2775  // Check that we know the level subdomain ids of all our neighbors. This
2776  // also involves coarser cells that share a vertex if they are active.
2777  //
2778  // Example (M= my, O=other):
2779  // *------*
2780  // | |
2781  // | O |
2782  // | |
2783  // *---*---*------*
2784  // | M | M |
2785  // *---*---*
2786  // | | M |
2787  // *---*---*
2788  // ^- the parent can be owned by somebody else, so O is not a neighbor
2789  // one level coarser
2790  if (settings & construct_multigrid_hierarchy)
2791  {
2792  for (unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
2793  {
2794  std::vector<bool> active_verts =
2795  this->mark_locally_active_vertices_on_level(lvl);
2796 
2797  const unsigned int maybe_coarser_lvl =
2798  (lvl > 0) ? (lvl - 1) : lvl;
2800  cell = this->begin(maybe_coarser_lvl),
2801  endc = this->end(lvl);
2802  for (; cell != endc; ++cell)
2803  if (cell->level() == static_cast<int>(lvl) || cell->is_active())
2804  {
2805  const bool is_level_artificial =
2806  (cell->level_subdomain_id() ==
2808  bool need_to_know = false;
2809  for (const unsigned int vertex :
2811  if (active_verts[cell->vertex_index(vertex)])
2812  {
2813  need_to_know = true;
2814  break;
2815  }
2816 
2817  Assert(
2818  !need_to_know || !is_level_artificial,
2819  ExcMessage(
2820  "Internal error: the owner of cell" +
2821  cell->id().to_string() +
2822  " is unknown even though it is needed for geometric multigrid."));
2823  }
2824  }
2825  }
2826 # endif
2827 
2828  this->update_periodic_face_map();
2829  this->update_number_cache();
2830 
2831  // signal that refinement is finished
2833  }
2834 
2835 
2836 
2837  template <int dim, int spacedim>
2838  void
2840  {
2841 # ifdef DEBUG
2842  for (const auto &cell : this->active_cell_iterators())
2843  if (cell->is_locally_owned())
2844  Assert(
2845  !cell->refine_flag_set() && !cell->coarsen_flag_set(),
2846  ExcMessage(
2847  "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
2848 # endif
2849 
2850  // signal that repartitioning is going to happen
2852 
2853  // before repartitioning the mesh, save a copy of the current positions of
2854  // quadrants
2855  // only if data needs to be transferred later
2856  std::vector<typename ::internal::p4est::types<dim>::gloidx>
2857  previous_global_first_quadrant;
2858 
2859  if (this->cell_attached_data.n_attached_data_sets > 0)
2860  {
2861  previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
2862  std::memcpy(previous_global_first_quadrant.data(),
2863  parallel_forest->global_first_quadrant,
2864  sizeof(
2865  typename ::internal::p4est::types<dim>::gloidx) *
2866  (parallel_forest->mpisize + 1));
2867  }
2868 
2869  if (this->signals.cell_weight.num_slots() == 0)
2870  {
2871  // no cell weights given -- call p4est's 'partition' without a
2872  // callback for cell weights
2874  parallel_forest,
2875  /* prepare coarsening */ 1,
2876  /* weight_callback */ nullptr);
2877  }
2878  else
2879  {
2880  // get cell weights for a weighted repartitioning.
2881  const std::vector<unsigned int> cell_weights = get_cell_weights();
2882 
2883  PartitionWeights<dim, spacedim> partition_weights(cell_weights);
2884 
2885  // attach (temporarily) a pointer to the cell weights through p4est's
2886  // user_pointer object
2887  Assert(parallel_forest->user_pointer == this, ExcInternalError());
2888  parallel_forest->user_pointer = &partition_weights;
2889 
2891  parallel_forest,
2892  /* prepare coarsening */ 1,
2893  /* weight_callback */
2894  &PartitionWeights<dim, spacedim>::cell_weight);
2895 
2896  // reset the user pointer to its previous state
2897  parallel_forest->user_pointer = this;
2898  }
2899 
2900  // pack data before triangulation gets updated
2901  if (this->cell_attached_data.n_attached_data_sets > 0)
2902  {
2903  this->data_transfer.pack_data(
2904  this->local_cell_relations,
2905  this->cell_attached_data.pack_callbacks_fixed,
2906  this->cell_attached_data.pack_callbacks_variable);
2907  }
2908 
2909  try
2910  {
2911  copy_local_forest_to_triangulation();
2912  }
2913  catch (const typename Triangulation<dim>::DistortedCellList &)
2914  {
2915  // the underlying triangulation should not be checking for distorted
2916  // cells
2917  Assert(false, ExcInternalError());
2918  }
2919 
2920  // transfer data after triangulation got updated
2921  if (this->cell_attached_data.n_attached_data_sets > 0)
2922  {
2923  this->execute_transfer(parallel_forest,
2924  previous_global_first_quadrant.data());
2925  }
2926 
2927  this->update_periodic_face_map();
2928 
2929  // update how many cells, edges, etc, we store locally
2930  this->update_number_cache();
2931 
2932  // signal that repartitioning is finished
2934  }
2935 
2936 
2937 
2938  template <int dim, int spacedim>
2939  const std::vector<types::global_dof_index> &
2941  const
2942  {
2943  return p4est_tree_to_coarse_cell_permutation;
2944  }
2945 
2946 
2947 
2948  template <int dim, int spacedim>
2949  const std::vector<types::global_dof_index> &
2951  const
2952  {
2953  return coarse_cell_to_p4est_tree_permutation;
2954  }
2955 
2956 
2957 
2958  template <int dim, int spacedim>
2959  std::vector<bool>
2961  const int level) const
2962  {
2963  Assert(dim > 1, ExcNotImplemented());
2964 
2965  std::vector<bool> marked_vertices(this->n_vertices(), false);
2966  for (const auto &cell : this->cell_iterators_on_level(level))
2967  if (cell->level_subdomain_id() == this->locally_owned_subdomain())
2968  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2969  marked_vertices[cell->vertex_index(v)] = true;
2970 
2976  // When a connectivity in the code below is detected, the assignment
2977  // 'marked_vertices[v1] = marked_vertices[v2] = true' makes sure that
2978  // the information about the periodicity propagates back to vertices on
2979  // cells that are not owned locally. However, in the worst case we want
2980  // to connect to a vertex that is 'dim' hops away from the locally owned
2981  // cell. Depending on the order of the periodic face map, we might
2982  // connect to that point by chance or miss it. However, after looping
2983  // through all the periodic directions (which are at most as many as
2984  // the number of space dimensions) we can be sure that all connections
2985  // to vertices have been created.
2986  for (unsigned int repetition = 0; repetition < dim; ++repetition)
2987  for (const auto &it : this->get_periodic_face_map())
2988  {
2989  const cell_iterator & cell_1 = it.first.first;
2990  const unsigned int face_no_1 = it.first.second;
2991  const cell_iterator & cell_2 = it.second.first.first;
2992  const unsigned int face_no_2 = it.second.first.second;
2993  const std::bitset<3> &face_orientation = it.second.second;
2994 
2995  if (cell_1->level() == level && cell_2->level() == level)
2996  {
2997  for (unsigned int v = 0;
2998  v < GeometryInfo<dim - 1>::vertices_per_cell;
2999  ++v)
3000  {
3001  // take possible non-standard orientation of faces into
3002  // account
3003  const unsigned int vface0 =
3005  v,
3006  face_orientation[0],
3007  face_orientation[1],
3008  face_orientation[2]);
3009  if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3010  vface0)] ||
3011  marked_vertices[cell_2->face(face_no_2)->vertex_index(
3012  v)])
3013  marked_vertices[cell_1->face(face_no_1)->vertex_index(
3014  vface0)] =
3015  marked_vertices[cell_2->face(face_no_2)->vertex_index(
3016  v)] = true;
3017  }
3018  }
3019  }
3020 
3021  return marked_vertices;
3022  }
3023 
3024 
3025 
3026  template <int dim, int spacedim>
3027  unsigned int
3030  {
3031  return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3032  }
3033 
3034 
3035 
3036  template <int dim, int spacedim>
3039  const unsigned int coarse_cell_index) const
3040  {
3041  return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3042  }
3043 
3044 
3045 
3046  template <int dim, int spacedim>
3047  void
3049  const std::vector<::GridTools::PeriodicFacePair<cell_iterator>>
3050  &periodicity_vector)
3051  {
3052  Assert(triangulation_has_content == true,
3053  ExcMessage("The triangulation is empty!"));
3054  Assert(this->n_levels() == 1,
3055  ExcMessage("The triangulation is refined!"));
3056 
3057  // call the base class for storing the periodicity information; we must
3058  // do this before going to p4est and rebuilding the triangulation to get
3059  // the level subdomain ids correct in the multigrid case
3061 
3062  for (const auto &face_pair : periodicity_vector)
3063  {
3064  const cell_iterator first_cell = face_pair.cell[0];
3065  const cell_iterator second_cell = face_pair.cell[1];
3066  const unsigned int face_left = face_pair.face_idx[0];
3067  const unsigned int face_right = face_pair.face_idx[1];
3068 
3069  // respective cells of the matching faces in p4est
3070  const unsigned int tree_left =
3071  coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3072  const unsigned int tree_right =
3073  coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3074 
3075  // p4est wants to know which corner the first corner on
3076  // the face with the lower id is mapped to on the face with
3077  // with the higher id. For d==2 there are only two possibilities
3078  // that are determined by it->orientation[1].
3079  // For d==3 we have to use GridTools::OrientationLookupTable.
3080  // The result is given below.
3081 
3082  unsigned int p4est_orientation = 0;
3083  if (dim == 2)
3084  p4est_orientation = face_pair.orientation[1];
3085  else
3086  {
3087  const unsigned int face_idx_list[] = {face_left, face_right};
3088  const cell_iterator cell_list[] = {first_cell, second_cell};
3089  unsigned int lower_idx, higher_idx;
3090  if (face_left <= face_right)
3091  {
3092  higher_idx = 1;
3093  lower_idx = 0;
3094  }
3095  else
3096  {
3097  higher_idx = 0;
3098  lower_idx = 1;
3099  }
3100 
3101  // get the cell index of the first index on the face with the
3102  // lower id
3103  unsigned int first_p4est_idx_on_cell =
3104  p8est_face_corners[face_idx_list[lower_idx]][0];
3105  unsigned int first_dealii_idx_on_face =
3107  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3108  ++i)
3109  {
3110  const unsigned int first_dealii_idx_on_cell =
3112  face_idx_list[lower_idx],
3113  i,
3114  cell_list[lower_idx]->face_orientation(
3115  face_idx_list[lower_idx]),
3116  cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3117  cell_list[lower_idx]->face_rotation(
3118  face_idx_list[lower_idx]));
3119  if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3120  {
3121  first_dealii_idx_on_face = i;
3122  break;
3123  }
3124  }
3125  Assert(first_dealii_idx_on_face != numbers::invalid_unsigned_int,
3126  ExcInternalError());
3127  // Now map dealii_idx_on_face according to the orientation
3128  constexpr unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
3129  {0, 1, 2, 3},
3130  {3, 1, 2, 0},
3131  {3, 2, 1, 0},
3132  {2, 3, 0, 1},
3133  {1, 3, 0, 2},
3134  {1, 0, 3, 2},
3135  {2, 0, 3, 1}};
3136  constexpr unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
3137  {0, 1, 2, 3},
3138  {3, 1, 2, 0},
3139  {3, 2, 1, 0},
3140  {2, 3, 0, 1},
3141  {2, 0, 3, 1},
3142  {1, 0, 3, 2},
3143  {1, 3, 0, 2}};
3144  const unsigned int second_dealii_idx_on_face =
3145  lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
3146  [first_dealii_idx_on_face] :
3147  right_to_left[face_pair.orientation.to_ulong()]
3148  [first_dealii_idx_on_face];
3149  const unsigned int second_dealii_idx_on_cell =
3151  face_idx_list[higher_idx],
3152  second_dealii_idx_on_face,
3153  cell_list[higher_idx]->face_orientation(
3154  face_idx_list[higher_idx]),
3155  cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
3156  cell_list[higher_idx]->face_rotation(
3157  face_idx_list[higher_idx]));
3158  // map back to p4est
3159  const unsigned int second_p4est_idx_on_face =
3160  p8est_corner_face_corners[second_dealii_idx_on_cell]
3161  [face_idx_list[higher_idx]];
3162  p4est_orientation = second_p4est_idx_on_face;
3163  }
3164 
3166  connectivity,
3167  tree_left,
3168  tree_right,
3169  face_left,
3170  face_right,
3171  p4est_orientation);
3172  }
3173 
3174 
3176  connectivity) == 1,
3177  ExcInternalError());
3178 
3179  // now create a forest out of the connectivity data structure
3182  this->mpi_communicator,
3183  connectivity,
3184  /* minimum initial number of quadrants per tree */ 0,
3185  /* minimum level of upfront refinement */ 0,
3186  /* use uniform upfront refinement */ 1,
3187  /* user_data_size = */ 0,
3188  /* user_data_constructor = */ nullptr,
3189  /* user_pointer */ this);
3190 
3191  try
3192  {
3193  copy_local_forest_to_triangulation();
3194  }
3195  catch (const typename Triangulation<dim>::DistortedCellList &)
3196  {
3197  // the underlying triangulation should not be checking for distorted
3198  // cells
3199  Assert(false, ExcInternalError());
3200  }
3201 
3202  // The range of ghost_owners might have changed so update that information
3203  this->update_number_cache();
3204  }
3205 
3206 
3207 
3208  template <int dim, int spacedim>
3209  std::size_t
3211  {
3212  std::size_t mem =
3215  MemoryConsumption::memory_consumption(triangulation_has_content) +
3217  MemoryConsumption::memory_consumption(parallel_forest) +
3219  this->cell_attached_data.n_attached_data_sets) +
3220  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_fixed)
3221  // +
3222  // MemoryConsumption::memory_consumption(cell_attached_data.pack_callbacks_variable)
3223  // +
3224  // TODO[TH]: how?
3226  coarse_cell_to_p4est_tree_permutation) +
3228  p4est_tree_to_coarse_cell_permutation) +
3229  memory_consumption_p4est();
3230 
3231  return mem;
3232  }
3233 
3234 
3235 
3236  template <int dim, int spacedim>
3237  std::size_t
3239  {
3240  return ::internal::p4est::functions<dim>::forest_memory_used(
3241  parallel_forest) +
3243  connectivity);
3244  }
3245 
3246 
3247 
3248  template <int dim, int spacedim>
3249  void
3251  const ::Triangulation<dim, spacedim> &other_tria)
3252  {
3253  Assert(
3254  (dynamic_cast<
3255  const ::parallel::distributed::Triangulation<dim, spacedim> *>(
3256  &other_tria)) ||
3257  (other_tria.n_global_levels() == 1),
3258  ExcNotImplemented());
3259 
3261 
3262  try
3263  {
3265  copy_triangulation(other_tria);
3266  }
3267  catch (
3268  const typename ::Triangulation<dim, spacedim>::DistortedCellList
3269  &)
3270  {
3271  // the underlying triangulation should not be checking for distorted
3272  // cells
3273  Assert(false, ExcInternalError());
3274  }
3275 
3276  if (const ::parallel::distributed::Triangulation<dim, spacedim>
3277  *other_distributed =
3278  dynamic_cast<const ::parallel::distributed::
3279  Triangulation<dim, spacedim> *>(&other_tria))
3280  {
3281  // copy parallel distributed specifics
3282  settings = other_distributed->settings;
3283  triangulation_has_content =
3284  other_distributed->triangulation_has_content;
3285  coarse_cell_to_p4est_tree_permutation =
3286  other_distributed->coarse_cell_to_p4est_tree_permutation;
3287  p4est_tree_to_coarse_cell_permutation =
3288  other_distributed->p4est_tree_to_coarse_cell_permutation;
3289  this->cell_attached_data = other_distributed->cell_attached_data;
3290  this->data_transfer = other_distributed->data_transfer;
3291 
3292  // create deep copy of connectivity graph
3293  typename ::internal::p4est::types<dim>::connectivity
3294  *temp_connectivity = const_cast<
3295  typename ::internal::p4est::types<dim>::connectivity *>(
3296  other_distributed->connectivity);
3297  connectivity =
3298  ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
3299 
3300  // create deep copy of parallel forest
3301  typename ::internal::p4est::types<dim>::forest *temp_forest =
3302  const_cast<typename ::internal::p4est::types<dim>::forest *>(
3303  other_distributed->parallel_forest);
3304  parallel_forest =
3306  false);
3307  parallel_forest->connectivity = connectivity;
3308  parallel_forest->user_pointer = this;
3309  }
3310  else
3311  {
3312  triangulation_has_content = true;
3313  setup_coarse_cell_to_p4est_tree_permutation();
3314  copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
3315  }
3316 
3317  try
3318  {
3319  copy_local_forest_to_triangulation();
3320  }
3321  catch (const typename Triangulation<dim>::DistortedCellList &)
3322  {
3323  // the underlying triangulation should not be checking for distorted
3324  // cells
3325  Assert(false, ExcInternalError());
3326  }
3327 
3328  this->update_periodic_face_map();
3329  this->update_number_cache();
3330  }
3331 
3332 
3333 
3334  template <int dim, int spacedim>
3335  void
3337  {
3338  // reorganize memory for local_cell_relations
3339  this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
3340  this->local_cell_relations.shrink_to_fit();
3341 
3342  // recurse over p4est
3343  for (const auto &cell : this->cell_iterators_on_level(0))
3344  {
3345  // skip coarse cells that are not ours
3346  if (tree_exists_locally<dim, spacedim>(
3347  parallel_forest,
3348  coarse_cell_to_p4est_tree_permutation[cell->index()]) == false)
3349  continue;
3350 
3351  // initialize auxiliary top level p4est quadrant
3352  typename ::internal::p4est::types<dim>::quadrant
3353  p4est_coarse_cell;
3354  ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
3355 
3356  // determine tree to start recursion on
3357  typename ::internal::p4est::types<dim>::tree *tree =
3358  init_tree(cell->index());
3359 
3360  update_cell_relations_recursively<dim, spacedim>(
3361  this->local_cell_relations, *tree, cell, p4est_coarse_cell);
3362  }
3363  }
3364 
3365 
3366 
3367  template <int dim, int spacedim>
3368  std::vector<unsigned int>
3370  {
3371  // check if local_cell_relations have been previously gathered
3372  // correctly
3373  Assert(this->local_cell_relations.size() ==
3374  static_cast<unsigned int>(parallel_forest->local_num_quadrants),
3375  ExcInternalError());
3376 
3377  // Allocate the space for the weights. In fact we do not know yet, how
3378  // many cells we own after the refinement (only p4est knows that
3379  // at this point). We simply reserve n_active_cells space and if many
3380  // more cells are refined than coarsened than additional reallocation
3381  // will be done inside get_cell_weights_recursively.
3382  std::vector<unsigned int> weights;
3383  weights.reserve(this->n_active_cells());
3384 
3385  // Iterate over p4est and Triangulation relations
3386  // to find refined/coarsened/kept
3387  // cells. Then append cell_weight.
3388  // Note that we need to follow the p4est ordering
3389  // instead of the deal.II ordering to get the cell_weights
3390  // in the same order p4est will encounter them during repartitioning.
3391  for (const auto &cell_rel : this->local_cell_relations)
3392  {
3393  const auto &cell_it = cell_rel.first;
3394  const auto &cell_status = cell_rel.second;
3395 
3396  switch (cell_status)
3397  {
3399  spacedim>::CELL_PERSIST:
3400  weights.push_back(1000);
3401  weights.back() += this->signals.cell_weight(
3402  cell_it,
3404  spacedim>::CELL_PERSIST);
3405  break;
3406 
3408  spacedim>::CELL_REFINE:
3410  spacedim>::CELL_INVALID:
3411  {
3412  // calculate weight of parent cell
3413  unsigned int parent_weight = 1000;
3414  parent_weight += this->signals.cell_weight(
3415  cell_it,
3417  CELL_REFINE);
3418  // assign the weight of the parent cell equally to all
3419  // children
3420  weights.push_back(parent_weight);
3421  break;
3422  }
3423 
3425  spacedim>::CELL_COARSEN:
3426  weights.push_back(1000);
3427  weights.back() += this->signals.cell_weight(
3428  cell_it,
3430  spacedim>::CELL_COARSEN);
3431  break;
3432 
3433  default:
3434  Assert(false, ExcInternalError());
3435  break;
3436  }
3437  }
3438 
3439  return weights;
3440  }
3441 
3442 
3443 
3444  template <int spacedim>
3446  const MPI_Comm &mpi_communicator,
3447  const typename ::Triangulation<1, spacedim>::MeshSmoothing
3448  smooth_grid,
3449  const Settings /*settings*/)
3450  : ::parallel::DistributedTriangulationBase<1, spacedim>(
3451  mpi_communicator,
3452  smooth_grid,
3453  false)
3454  {
3455  Assert(false, ExcNotImplemented());
3456  }
3457 
3458 
3459  template <int spacedim>
3461  {
3462  AssertNothrow(false, ExcNotImplemented());
3463  }
3464 
3465 
3466 
3467  template <int spacedim>
3468  const std::vector<types::global_dof_index> &
3470  const
3471  {
3472  static std::vector<types::global_dof_index> a;
3473  return a;
3474  }
3475 
3476 
3477 
3478  template <int spacedim>
3479  std::map<unsigned int, std::set<::types::subdomain_id>>
3481  const unsigned int /*level*/) const
3482  {
3483  Assert(false, ExcNotImplemented());
3484 
3485  return std::map<unsigned int, std::set<::types::subdomain_id>>();
3486  }
3487 
3488 
3489 
3490  template <int spacedim>
3491  std::vector<bool>
3493  const unsigned int) const
3494  {
3495  Assert(false, ExcNotImplemented());
3496  return std::vector<bool>();
3497  }
3498 
3499 
3500 
3501  template <int spacedim>
3502  unsigned int
3504  const types::coarse_cell_id) const
3505  {
3506  Assert(false, ExcNotImplemented());
3507  return 0;
3508  }
3509 
3510 
3511 
3512  template <int spacedim>
3515  const unsigned int) const
3516  {
3517  Assert(false, ExcNotImplemented());
3518  return 0;
3519  }
3520 
3521 
3522 
3523  template <int spacedim>
3524  void
3525  Triangulation<1, spacedim>::load(const std::string &, const bool)
3526  {
3527  Assert(false, ExcNotImplemented());
3528  }
3529 
3530 
3531 
3532  template <int spacedim>
3533  void
3534  Triangulation<1, spacedim>::save(const std::string &) const
3535  {
3536  Assert(false, ExcNotImplemented());
3537  }
3538 
3539 
3540 
3541  template <int spacedim>
3542  bool
3544  {
3545  Assert(false, ExcNotImplemented());
3546  return false;
3547  }
3548 
3549 
3550 
3551  template <int spacedim>
3552  void
3554  {
3555  Assert(false, ExcNotImplemented());
3556  }
3557 
3558  } // namespace distributed
3559 } // namespace parallel
3560 
3561 
3562 #endif // DEAL_II_WITH_P4EST
3563 
3564 
3565 
3566 namespace parallel
3567 {
3568  namespace distributed
3569  {
3570  template <int dim, int spacedim>
3573  : distributed_tria(
3574  dynamic_cast<
3575  ::parallel::distributed::Triangulation<dim, spacedim> *>(
3576  &tria))
3577  {
3578 #ifdef DEAL_II_WITH_P4EST
3579  if (distributed_tria != nullptr)
3580  {
3581  // Save the current set of refinement flags, and adjust the
3582  // refinement flags to be consistent with the p4est oracle.
3583  distributed_tria->save_coarsen_flags(saved_coarsen_flags);
3584  distributed_tria->save_refine_flags(saved_refine_flags);
3585 
3586  for (const auto &pair : distributed_tria->local_cell_relations)
3587  {
3588  const auto &cell = pair.first;
3589  const auto &status = pair.second;
3590 
3591  switch (status)
3592  {
3593  case ::Triangulation<dim, spacedim>::CELL_PERSIST:
3594  // cell remains unchanged
3595  cell->clear_refine_flag();
3596  cell->clear_coarsen_flag();
3597  break;
3598 
3599  case ::Triangulation<dim, spacedim>::CELL_REFINE:
3600  // cell will be refined
3601  cell->clear_coarsen_flag();
3602  cell->set_refine_flag();
3603  break;
3604 
3605  case ::Triangulation<dim, spacedim>::CELL_COARSEN:
3606  // children of this cell will be coarsened
3607  for (const auto &child : cell->child_iterators())
3608  {
3609  child->clear_refine_flag();
3610  child->set_coarsen_flag();
3611  }
3612  break;
3613 
3614  case ::Triangulation<dim, spacedim>::CELL_INVALID:
3615  // do nothing as cell does not exist yet
3616  break;
3617 
3618  default:
3619  Assert(false, ExcInternalError());
3620  break;
3621  }
3622  }
3623  }
3624 #endif
3625  }
3626 
3627 
3628 
3629  template <int dim, int spacedim>
3631  {
3632 #ifdef DEAL_II_WITH_P4EST
3633  if (distributed_tria)
3634  {
3635  // Undo the refinement flags modification.
3636  distributed_tria->load_coarsen_flags(saved_coarsen_flags);
3637  distributed_tria->load_refine_flags(saved_refine_flags);
3638  }
3639 #else
3640  // pretend that this destructor does something to silence clang-tidy
3641  (void)distributed_tria;
3642 #endif
3643  }
3644  } // namespace distributed
3645 } // namespace parallel
3646 
3647 
3648 
3649 /*-------------- Explicit Instantiations -------------------------------*/
3650 #include "tria.inst"
3651 
3652 
bool is_multilevel_hierarchy_constructed() const override
Definition: tria.cc:1264
unsigned int n_active_cells() const
Definition: tria.cc:12638
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:10329
unsigned int n_vertices() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1528
bool tree_exists_locally(const typename types< dim >::forest *parallel_forest, const typename types< dim >::topidx coarse_grid_cell)
boost::signals2::signal< void()> post_distributed_repartition
Definition: tria.h:2259
bool all_reference_cells_are_hyper_cube() const
Definition: tria.cc:13506
virtual bool has_hanging_nodes() const
Definition: tria.cc:12770
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:65
unsigned int n_cells() const
Definition: tria.cc:12630
types::global_dof_index size_type
Definition: cuda_kernels.h:45
static ::ExceptionBase & ExcIO()
MeshSmoothing smooth_grid
Definition: tria.h:3525
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< Point< spacedim > > vertices
Definition: tria.h:3992
boost::signals2::signal< void()> post_p4est_refinement
Definition: tria.h:2230
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12148
boost::signals2::signal< void()> pre_distributed_refinement
Definition: tria.h:2221
virtual void clear() override
Definition: tria_base.cc:680
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11951
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
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:1113
virtual void save(const std::string &filename) const override
Definition: tria.cc:1425
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11931
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
Definition: tria.cc:2940
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:3697
cell_iterator end() const
Definition: tria.cc:12042
virtual unsigned int coarse_cell_id_to_coarse_cell_index(const types::coarse_cell_id coarse_cell_id) const override
Definition: tria.cc:3028
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator >> &)
Definition: tria.cc:13276
virtual void execute_coarsening_and_refinement()
Definition: tria.cc:13305
virtual std::size_t memory_consumption() const override
Definition: tria_base.cc:90
unsigned int n_active_lines() const
Definition: tria.cc:12830
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition: tria.cc:12159
virtual bool prepare_coarsening_and_refinement()
Definition: tria.cc:14095
static ::ExceptionBase & ExcMessage(std::string arg1)
Triangulation(const MeshSmoothing smooth_grid=none, const bool check_for_distorted_cells=false)
Definition: tria.cc:9998
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const override
Definition: tria.cc:3038
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:10440
#define Assert(cond, exc)
Definition: exceptions.h:1465
Signals signals
Definition: tria.h:2295
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:10894
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
unsigned int level
Definition: grid_out.cc:4590
void save_refine_flags(std::ostream &out) const
Definition: tria.cc:10825
boost::signals2::signal< void()> pre_distributed_save
Definition: tria.h:2267
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1482
virtual std::size_t memory_consumption() const
Definition: tria.cc:15181
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)
virtual void load(const std::string &filename, const bool autopartition=true) override
Definition: tria.cc:1480
virtual void update_cell_relations() override
Definition: tria.cc:3336
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:1069
boost::signals2::signal< void()> post_distributed_save
Definition: tria.h:2274
active_cell_iterator end_active(const unsigned int level) const
Definition: tria.cc:12111
const std::vector< bool > & get_used_vertices() const
Definition: tria.cc:13182
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
virtual unsigned int n_global_levels() const
virtual ~Triangulation() override
Definition: tria.cc:1141
void update_periodic_face_map()
Definition: tria.cc:13417
virtual ~Triangulation() override
Definition: tria.cc:10074
boost::signals2::signal< unsigned int(const cell_iterator &, const CellStatus), CellWeightSum< unsigned int > > cell_weight
Definition: tria.h:2208
boost::signals2::signal< void()> pre_distributed_load
Definition: tria.h:2282
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
Definition: tria.cc:12170
TemporarilyMatchRefineFlags(::Triangulation< dim, spacedim > &tria)
Definition: tria.cc:3571
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
static ::ExceptionBase & ExcNotImplemented()
IteratorRange< cell_iterator > cell_iterators() const
Definition: tria.cc:12139
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
active_cell_iterator last_active() const
Definition: tria.cc:11996
boost::signals2::signal< void()> post_distributed_load
Definition: tria.h:2289
global_cell_index coarse_cell_id
Definition: types.h:114
boost::signals2::signal< void()> pre_distributed_repartition
Definition: tria.h:2252
std::vector< bool > mark_locally_active_vertices_on_level(const int level) const
Definition: tria.cc:2960
T max(const T &t, const MPI_Comm &mpi_communicator)
virtual types::subdomain_id locally_owned_subdomain() const
Definition: tria.cc:13249
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:10103
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:13296
virtual void clear() override
Definition: tria.cc:1228
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:2240