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