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