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