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