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