45 namespace TriangulationImplementation
54 template <
int dim,
int spacedim>
61 &cell) -> std::uint8_t {
62 if (cell->refine_flag_set())
64 if (cell->coarsen_flag_set())
72 const std::uint8_t &flag) ->
void {
73 cell->clear_coarsen_flag();
74 cell->clear_refine_flag();
76 cell->set_refine_flag();
78 cell->set_coarsen_flag();
81 GridTools::exchange_cell_data_to_ghosts<std::uint8_t>(tria,
92#ifdef DEAL_II_WITH_P4EST
96 template <
int dim,
int spacedim>
98 get_vertex_to_cell_mappings(
100 std::vector<unsigned int> &vertex_touch_count,
101 std::vector<std::list<
103 unsigned int>>> &vertex_to_cell)
108 for (
const auto &cell :
triangulation.active_cell_iterators())
111 ++vertex_touch_count[cell->vertex_index(v)];
112 vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
118 template <
int dim,
int spacedim>
120 get_edge_to_cell_mappings(
122 std::vector<unsigned int> &edge_touch_count,
123 std::vector<std::list<
125 unsigned int>>> &edge_to_cell)
132 for (
const auto &cell :
triangulation.active_cell_iterators())
133 for (unsigned
int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
135 ++edge_touch_count[cell->line(l)->index()];
136 edge_to_cell[cell->line(l)->index()].emplace_back(cell, l);
146 template <
int dim,
int spacedim>
148 set_vertex_and_cell_info(
150 const std::vector<unsigned int> &vertex_touch_count,
151 const std::vector<std::list<
153 unsigned int>>> &vertex_to_cell,
154 const std::vector<types::global_dof_index>
155 &coarse_cell_to_p4est_tree_permutation,
156 const bool set_vertex_info,
172 if (set_vertex_info ==
true)
173 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
175 connectivity->vertices[3 * v] =
triangulation.get_vertices()[v][0];
176 connectivity->vertices[3 * v + 1] =
178 connectivity->vertices[3 * v + 2] =
190 for (; cell != endc; ++cell)
192 const unsigned int index =
193 coarse_cell_to_p4est_tree_permutation[cell->index()];
197 if (set_vertex_info ==
true)
200 v] = cell->vertex_index(v);
203 v] = cell->vertex_index(v);
209 if (cell->face(f)->at_boundary() == false)
212 coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->
index()];
216 coarse_cell_to_p4est_tree_permutation[cell->index()];
221 if (cell->face(f)->at_boundary() == false)
227 connectivity->tree_to_face
229 cell->neighbor_of_neighbor(f);
248 connectivity->tree_to_face[
index * 6 + f] =
249 cell->neighbor_of_neighbor(f);
251 unsigned int face_idx_list[2] = {
252 f, cell->neighbor_of_neighbor(f)};
254 cell_list[2] = {cell, cell->neighbor(f)};
255 unsigned int smaller_idx = 0;
257 if (f > cell->neighbor_of_neighbor(f))
260 unsigned int larger_idx = (smaller_idx + 1) % 2;
268 unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
270 face_idx_list[smaller_idx],
272 cell_list[smaller_idx]->face_orientation(
273 face_idx_list[smaller_idx]),
274 cell_list[smaller_idx]->face_flip(
275 face_idx_list[smaller_idx]),
276 cell_list[smaller_idx]->face_rotation(
277 face_idx_list[smaller_idx])));
281 for (
unsigned int i = 0;
282 i < GeometryInfo<dim>::vertices_per_face;
286 cell_list[larger_idx]->vertex_index(
288 face_idx_list[larger_idx], i));
297 connectivity->tree_to_face[
index * 6 + f] += 6 * v;
311 connectivity->ctt_offset[0] = 0;
312 std::partial_sum(vertex_touch_count.begin(),
313 vertex_touch_count.end(),
314 &connectivity->ctt_offset[1]);
317 std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
322 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
324 Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
328 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
329 unsigned int>>::const_iterator p =
330 vertex_to_cell[v].begin();
331 for (
unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
333 connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
334 coarse_cell_to_p4est_tree_permutation[p->first->index()];
335 connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
343 template <
int dim,
int spacedim>
349 Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
351 return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
352 (coarse_grid_cell <= parallel_forest->last_local_tree));
356 template <
int dim,
int spacedim>
358 delete_all_children_and_self(
361 if (cell->has_children())
362 for (
unsigned int c = 0; c < cell->n_children(); ++c)
363 delete_all_children_and_self<dim, spacedim>(cell->child(c));
365 cell->set_coarsen_flag();
370 template <
int dim,
int spacedim>
375 if (cell->has_children())
376 for (
unsigned int c = 0; c < cell->n_children(); ++c)
377 delete_all_children_and_self<dim, spacedim>(cell->child(c));
381 template <
int dim,
int spacedim>
383 determine_level_subdomain_id_recursively(
390 const std::vector<std::vector<bool>> &marked_vertices)
399 if (marked_vertices[dealii_cell->level()]
400 [dealii_cell->vertex_index(v)])
419 if (!used && dealii_cell->is_active() &&
420 dealii_cell->is_artificial() ==
false &&
421 dealii_cell->level() + 1 <
static_cast<int>(marked_vertices.size()))
425 if (marked_vertices[dealii_cell->level() + 1]
426 [dealii_cell->vertex_index(v)])
435 if (!used && dealii_cell->is_active() &&
436 dealii_cell->is_artificial() ==
false && dealii_cell->level() > 0)
440 if (marked_vertices[dealii_cell->level() - 1]
441 [dealii_cell->vertex_index(v)])
452 &forest, tree_index, &p4est_cell, my_subdomain);
453 Assert((owner != -2) && (owner != -1),
455 dealii_cell->set_level_subdomain_id(owner);
459 if (dealii_cell->has_children())
463 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
468 P4EST_QUADRANT_INIT(&p4est_child[c]);
471 P8EST_QUADRANT_INIT(&p4est_child[c]);
481 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
484 determine_level_subdomain_id_recursively<dim, spacedim>(
487 dealii_cell->child(c),
497 template <
int dim,
int spacedim>
499 match_tree_recursively(
507 if (sc_array_bsearch(
const_cast<sc_array_t *
>(&tree.quadrants),
513 delete_all_children<dim, spacedim>(dealii_cell);
514 if (dealii_cell->is_active())
515 dealii_cell->set_subdomain_id(my_subdomain);
524 if (dealii_cell->is_active())
525 dealii_cell->set_refine_flag();
530 for (
unsigned int c = 0;
531 c < GeometryInfo<dim>::max_children_per_cell;
536 P4EST_QUADRANT_INIT(&p4est_child[c]);
539 P8EST_QUADRANT_INIT(&p4est_child[c]);
549 for (
unsigned int c = 0;
550 c < GeometryInfo<dim>::max_children_per_cell;
555 &p4est_child[c]) ==
false)
561 delete_all_children<dim, spacedim>(dealii_cell->child(c));
562 dealii_cell->child(c)->recursively_set_subdomain_id(
569 match_tree_recursively<dim, spacedim>(tree,
570 dealii_cell->child(c),
580 template <
int dim,
int spacedim>
583 const ::Triangulation<dim, spacedim> *tria,
584 unsigned int dealii_index,
588 const int l = ghost_quadrant.level;
590 for (
int i = 0; i <
l; ++i)
595 if (cell->is_active())
597 cell->clear_coarsen_flag();
598 cell->set_refine_flag();
605 dealii_index = cell->child_index(child_id);
611 if (cell->has_children())
612 delete_all_children<dim, spacedim>(cell);
615 cell->clear_coarsen_flag();
616 cell->set_subdomain_id(ghost_owner);
620# ifdef P4EST_SEARCH_LOCAL
622 class PartitionSearch
630 PartitionSearch(
const PartitionSearch<dim> &other) =
delete;
632 PartitionSearch<dim> &
633 operator=(
const PartitionSearch<dim> &other) =
delete;
689 quad_length_on_level);
692 initialize_mapping();
695 map_real_to_unit_cell(
const Point<dim> &p)
const;
698 is_in_this_quadrant(
const Point<dim> &p)
const;
701 std::vector<Point<dim>> cell_vertices;
709 bool are_vertices_initialized;
711 bool is_reference_mapping_initialized;
717 QuadrantData quadrant_data;
724 PartitionSearch<dim>::local_quadrant_fn(
738 PartitionSearch<dim> *this_object =
739 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
743 quad_length_on_level =
745 (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
749 this_object->quadrant_data.set_cell_vertices(forest,
752 quad_length_on_level);
755 this_object->quadrant_data.initialize_mapping();
765 PartitionSearch<dim>::local_point_fn(
778 PartitionSearch<dim> *this_object =
779 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
782 double *this_point_dptr =
static_cast<double *
>(
point);
785 (dim == 2 ?
Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
786 Point<dim>(this_point_dptr[0],
788 this_point_dptr[2]));
791 const bool is_in_this_quadrant =
792 this_object->quadrant_data.is_in_this_quadrant(this_point);
796 if (!is_in_this_quadrant)
805 if (rank_begin < rank_end)
813 this_point_dptr[dim] =
static_cast<double>(rank_begin);
823 PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
826 const Point<dim> p_ref = map_real_to_unit_cell(p);
835 PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
838 Assert(is_reference_mapping_initialized,
840 "Cell vertices and mapping coefficients must be fully "
841 "initialized before transforming a point to the unit cell."));
847 for (
unsigned int alpha = 0;
848 alpha < GeometryInfo<dim>::vertices_per_cell;
854 p_out += (quadrant_mapping_matrix(alpha, 0) +
855 quadrant_mapping_matrix(alpha, 1) * p(0) +
856 quadrant_mapping_matrix(alpha, 2) * p(1) +
857 quadrant_mapping_matrix(alpha, 3) * p(0) * p(1)) *
863 for (
unsigned int alpha = 0;
864 alpha < GeometryInfo<dim>::vertices_per_cell;
870 p_out += (quadrant_mapping_matrix(alpha, 0) +
871 quadrant_mapping_matrix(alpha, 1) * p(0) +
872 quadrant_mapping_matrix(alpha, 2) * p(1) +
873 quadrant_mapping_matrix(alpha, 3) * p(2) +
874 quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
875 quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
876 quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
877 quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
887 PartitionSearch<dim>::QuadrantData::QuadrantData()
889 , quadrant_mapping_matrix(
GeometryInfo<dim>::vertices_per_cell,
891 , are_vertices_initialized(false)
892 , is_reference_mapping_initialized(false)
899 PartitionSearch<dim>::QuadrantData::initialize_mapping()
902 are_vertices_initialized,
904 "Cell vertices must be initialized before the cell mapping can be filled."));
911 for (
unsigned int alpha = 0;
912 alpha < GeometryInfo<dim>::vertices_per_cell;
916 point_matrix(0, alpha) = 1;
917 point_matrix(1, alpha) = cell_vertices[alpha](0);
918 point_matrix(2, alpha) = cell_vertices[alpha](1);
919 point_matrix(3, alpha) =
920 cell_vertices[alpha](0) * cell_vertices[alpha](1);
927 quadrant_mapping_matrix.invert(point_matrix);
931 for (
unsigned int alpha = 0;
932 alpha < GeometryInfo<dim>::vertices_per_cell;
936 point_matrix(0, alpha) = 1;
937 point_matrix(1, alpha) = cell_vertices[alpha](0);
938 point_matrix(2, alpha) = cell_vertices[alpha](1);
939 point_matrix(3, alpha) = cell_vertices[alpha](2);
940 point_matrix(4, alpha) =
941 cell_vertices[alpha](0) * cell_vertices[alpha](1);
942 point_matrix(5, alpha) =
943 cell_vertices[alpha](1) * cell_vertices[alpha](2);
944 point_matrix(6, alpha) =
945 cell_vertices[alpha](0) * cell_vertices[alpha](2);
946 point_matrix(7, alpha) = cell_vertices[alpha](0) *
947 cell_vertices[alpha](1) *
948 cell_vertices[alpha](2);
955 quadrant_mapping_matrix.invert(point_matrix);
958 is_reference_mapping_initialized =
true;
965 PartitionSearch<2>::QuadrantData::set_cell_vertices(
970 quad_length_on_level)
972 constexpr unsigned int dim = 2;
976 double corner_point[dim + 1] = {0};
979 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
981 for (
unsigned int d = 0;
d < dim; ++
d)
983 cell_vertices[vertex_index](
d) = corner_point[d];
993 unsigned int vertex_index = 0;
995 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
998 copy_vertex(vertex_index);
1005 forest->connectivity,
1007 quadrant->x + quad_length_on_level,
1012 copy_vertex(vertex_index);
1019 forest->connectivity,
1022 quadrant->y + quad_length_on_level,
1026 copy_vertex(vertex_index);
1033 forest->connectivity,
1035 quadrant->x + quad_length_on_level,
1036 quadrant->y + quad_length_on_level,
1040 copy_vertex(vertex_index);
1042 are_vertices_initialized =
true;
1049 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1054 quad_length_on_level)
1056 constexpr unsigned int dim = 3;
1058 double corner_point[dim] = {0};
1061 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1063 for (
unsigned int d = 0;
d < dim; ++
d)
1065 cell_vertices[vertex_index](
d) = corner_point[d];
1067 corner_point[
d] = 0;
1075 unsigned int vertex_index = 0;
1077 forest->connectivity,
1085 copy_vertex(vertex_index);
1093 forest->connectivity,
1095 quadrant->x + quad_length_on_level,
1101 copy_vertex(vertex_index);
1108 forest->connectivity,
1111 quadrant->y + quad_length_on_level,
1116 copy_vertex(vertex_index);
1123 forest->connectivity,
1125 quadrant->x + quad_length_on_level,
1126 quadrant->y + quad_length_on_level,
1131 copy_vertex(vertex_index);
1138 forest->connectivity,
1142 quadrant->z + quad_length_on_level,
1146 copy_vertex(vertex_index);
1153 forest->connectivity,
1155 quadrant->x + quad_length_on_level,
1157 quadrant->z + quad_length_on_level,
1161 copy_vertex(vertex_index);
1168 forest->connectivity,
1171 quadrant->y + quad_length_on_level,
1172 quadrant->z + quad_length_on_level,
1176 copy_vertex(vertex_index);
1183 forest->connectivity,
1185 quadrant->x + quad_length_on_level,
1186 quadrant->y + quad_length_on_level,
1187 quadrant->z + quad_length_on_level,
1191 copy_vertex(vertex_index);
1194 are_vertices_initialized =
true;
1204 template <
int dim,
int spacedim>
1205 class RefineAndCoarsenList
1209 const std::vector<types::global_dof_index>
1210 &p4est_tree_to_coarse_cell_permutation,
1238 pointers_are_at_end()
const;
1241 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1242 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1243 const_iterator current_refine_pointer;
1245 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1246 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1247 const_iterator current_coarsen_pointer;
1258 template <
int dim,
int spacedim>
1260 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1262 return ((current_refine_pointer == refine_list.end()) &&
1263 (current_coarsen_pointer == coarsen_list.end()));
1268 template <
int dim,
int spacedim>
1269 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1271 const std::vector<types::global_dof_index>
1272 &p4est_tree_to_coarse_cell_permutation,
1276 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1277 for (
const auto &cell :
triangulation.active_cell_iterators())
1280 if (cell->subdomain_id() != my_subdomain)
1283 if (cell->refine_flag_set())
1285 else if (cell->coarsen_flag_set())
1289 refine_list.reserve(n_refine_flags);
1290 coarsen_list.reserve(n_coarsen_flags);
1302 unsigned int coarse_cell_index =
1303 p4est_tree_to_coarse_cell_permutation[c];
1312 p4est_cell.p.which_tree = c;
1313 build_lists(cell, p4est_cell, my_subdomain);
1321 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1322 Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1324 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1325 Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1328 current_refine_pointer = refine_list.begin();
1329 current_coarsen_pointer = coarsen_list.begin();
1334 template <
int dim,
int spacedim>
1336 RefineAndCoarsenList<dim, spacedim>::build_lists(
1341 if (cell->is_active())
1343 if (cell->subdomain_id() == my_subdomain)
1345 if (cell->refine_flag_set())
1346 refine_list.push_back(p4est_cell);
1347 else if (cell->coarsen_flag_set())
1348 coarsen_list.push_back(p4est_cell);
1355 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1360 P4EST_QUADRANT_INIT(&p4est_child[c]);
1363 P8EST_QUADRANT_INIT(&p4est_child[c]);
1370 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1373 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1374 build_lists(cell->child(c), p4est_child[c], my_subdomain);
1380 template <
int dim,
int spacedim>
1382 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1387 RefineAndCoarsenList<dim, spacedim> *this_object =
1388 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1389 forest->user_pointer);
1393 if (this_object->current_refine_pointer == this_object->refine_list.end())
1396 Assert(coarse_cell_index <=
1397 this_object->current_refine_pointer->p.which_tree,
1402 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1406 Assert(coarse_cell_index <=
1407 this_object->current_refine_pointer->p.which_tree,
1413 quadrant, &*this_object->current_refine_pointer) <= 0,
1418 quadrant, &*this_object->current_refine_pointer))
1420 ++this_object->current_refine_pointer;
1430 template <
int dim,
int spacedim>
1432 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1437 RefineAndCoarsenList<dim, spacedim> *this_object =
1438 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1439 forest->user_pointer);
1443 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1446 Assert(coarse_cell_index <=
1447 this_object->current_coarsen_pointer->p.which_tree,
1452 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1456 Assert(coarse_cell_index <=
1457 this_object->current_coarsen_pointer->p.which_tree,
1463 children[0], &*this_object->current_coarsen_pointer) <= 0,
1469 children[0], &*this_object->current_coarsen_pointer))
1472 ++this_object->current_coarsen_pointer;
1476 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1480 children[c], &*this_object->current_coarsen_pointer),
1482 ++this_object->current_coarsen_pointer;
1500 template <
int dim,
int spacedim>
1501 class PartitionWeights
1509 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1524 std::vector<unsigned int> cell_weights_list;
1525 std::vector<unsigned int>::const_iterator current_pointer;
1529 template <
int dim,
int spacedim>
1530 PartitionWeights<dim, spacedim>::PartitionWeights(
1531 const std::vector<unsigned int> &cell_weights)
1532 : cell_weights_list(cell_weights)
1536 current_pointer = cell_weights_list.begin();
1540 template <
int dim,
int spacedim>
1542 PartitionWeights<dim, spacedim>::cell_weight(
1551 PartitionWeights<dim, spacedim> *this_object =
1552 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1554 Assert(this_object->current_pointer >=
1555 this_object->cell_weights_list.begin(),
1557 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1563 const unsigned int weight = *this_object->current_pointer;
1564 ++this_object->current_pointer;
1566 Assert(weight <
static_cast<unsigned int>(std::numeric_limits<int>::max()),
1567 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1568 "weights for cells. The weight provided here exceeds "
1569 "the maximum value represented as a 'signed int'."));
1570 return static_cast<int>(weight);
1573 template <
int dim,
int spacedim>
1574 using cell_relation_t =
typename std::pair<
1575 typename ::Triangulation<dim, spacedim>::cell_iterator,
1587 template <
int dim,
int spacedim>
1589 add_single_cell_relation(
1590 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1591 const typename ::internal::p4est::types<dim>::tree &tree,
1592 const unsigned int idx,
1596 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1602 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1616 template <
int dim,
int spacedim>
1618 update_cell_relations_recursively(
1619 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1620 const typename ::internal::p4est::types<dim>::tree &tree,
1622 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1625 const int idx = sc_array_bsearch(
1626 const_cast<sc_array_t *
>(&tree.quadrants),
1631 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1633 &p4est_cell) ==
false))
1638 const bool p4est_has_children = (idx == -1);
1639 if (p4est_has_children && dealii_cell->has_children())
1642 typename ::internal::p4est::types<dim>::quadrant
1645 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1650 P4EST_QUADRANT_INIT(&p4est_child[c]);
1653 P8EST_QUADRANT_INIT(&p4est_child[c]);
1660 &p4est_cell, p4est_child);
1662 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1665 update_cell_relations_recursively<dim, spacedim>(
1666 cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1669 else if (!p4est_has_children && !dealii_cell->has_children())
1673 add_single_cell_relation<dim, spacedim>(
1676 else if (p4est_has_children)
1684 typename ::internal::p4est::types<dim>::quadrant
1686 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1691 P4EST_QUADRANT_INIT(&p4est_child[c]);
1694 P8EST_QUADRANT_INIT(&p4est_child[c]);
1701 &p4est_cell, p4est_child);
1709 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1712 child_idx = sc_array_bsearch(
1713 const_cast<sc_array_t *
>(&tree.quadrants),
1720 add_single_cell_relation<dim, spacedim>(
1721 cell_rel, tree, child_idx, dealii_cell, cell_status);
1729 add_single_cell_relation<dim, spacedim>(
1743 namespace distributed
1746 template <
int dim,
int spacedim>
1758 (settings & construct_multigrid_hierarchy) ?
1762 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1765 , settings(settings)
1766 , triangulation_has_content(false)
1767 , connectivity(
nullptr)
1768 , parallel_forest(
nullptr)
1770 parallel_ghost =
nullptr;
1775 template <
int dim,
int spacedim>
1796 template <
int dim,
int spacedim>
1809 const typename ::Triangulation<dim, spacedim>::DistortedCellList
1818 this->all_reference_cells_are_hyper_cube(),
1820 "The class parallel::distributed::Triangulation only supports meshes "
1821 "consisting only of hypercube-like cells."));
1826 triangulation_has_content =
true;
1828 setup_coarse_cell_to_p4est_tree_permutation();
1830 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1834 copy_local_forest_to_triangulation();
1843 this->update_periodic_face_map();
1844 this->update_number_cache();
1849 template <
int dim,
int spacedim>
1855 (void)construction_data;
1862 template <
int dim,
int spacedim>
1866 triangulation_has_content =
false;
1868 if (parallel_ghost !=
nullptr)
1872 parallel_ghost =
nullptr;
1875 if (parallel_forest !=
nullptr)
1878 parallel_forest =
nullptr;
1881 if (connectivity !=
nullptr)
1885 connectivity =
nullptr;
1888 coarse_cell_to_p4est_tree_permutation.resize(0);
1889 p4est_tree_to_coarse_cell_permutation.resize(0);
1893 this->update_number_cache();
1898 template <
int dim,
int spacedim>
1909 template <
int dim,
int spacedim>
1920 template <
int dim,
int spacedim>
1926 *previous_global_first_quadrant)
1928 Assert(this->data_serializer.sizes_fixed_cumulative.size() > 0,
1932 this->data_serializer.dest_data_fixed.resize(
1933 parallel_forest->local_num_quadrants *
1934 this->data_serializer.sizes_fixed_cumulative.back());
1937 typename ::internal::p4est::types<dim>::transfer_context
1941 parallel_forest->global_first_quadrant,
1942 previous_global_first_quadrant,
1943 parallel_forest->mpicomm,
1945 this->data_serializer.dest_data_fixed.data(),
1946 this->data_serializer.src_data_fixed.data(),
1947 this->data_serializer.sizes_fixed_cumulative.back());
1949 if (this->data_serializer.variable_size_data_stored)
1952 this->data_serializer.dest_sizes_variable.resize(
1953 parallel_forest->local_num_quadrants);
1958 parallel_forest->global_first_quadrant,
1959 previous_global_first_quadrant,
1960 parallel_forest->mpicomm,
1962 this->data_serializer.dest_sizes_variable.data(),
1963 this->data_serializer.src_sizes_variable.data(),
1964 sizeof(
unsigned int));
1970 this->data_serializer.src_data_fixed.clear();
1971 this->data_serializer.src_data_fixed.shrink_to_fit();
1973 if (this->data_serializer.variable_size_data_stored)
1976 this->data_serializer.dest_data_variable.resize(
1977 std::accumulate(this->data_serializer.dest_sizes_variable.begin(),
1978 this->data_serializer.dest_sizes_variable.end(),
1979 std::vector<int>::size_type(0)));
1981# if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1988 if (this->data_serializer.src_sizes_variable.empty())
1989 this->data_serializer.src_sizes_variable.resize(1);
1990 if (this->data_serializer.dest_sizes_variable.empty())
1991 this->data_serializer.dest_sizes_variable.resize(1);
1996 parallel_forest->global_first_quadrant,
1997 previous_global_first_quadrant,
1998 parallel_forest->mpicomm,
2000 this->data_serializer.dest_data_variable.data(),
2001 this->data_serializer.dest_sizes_variable.data(),
2002 this->data_serializer.src_data_variable.data(),
2003 this->data_serializer.src_sizes_variable.data());
2006 this->data_serializer.src_sizes_variable.clear();
2007 this->data_serializer.src_sizes_variable.shrink_to_fit();
2008 this->data_serializer.src_data_variable.clear();
2009 this->data_serializer.src_data_variable.shrink_to_fit();
2015 template <
int dim,
int spacedim>
2018 spacedim>::setup_coarse_cell_to_p4est_tree_permutation()
2023 coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
2025 cell_connectivity, coarse_cell_to_p4est_tree_permutation);
2027 p4est_tree_to_coarse_cell_permutation =
2033 template <
int dim,
int spacedim>
2036 const
std::
string &file_basename)
const
2038 Assert(parallel_forest !=
nullptr,
2039 ExcMessage(
"Can't produce output when no forest is created yet."));
2043 "To use this function the triangulation's flag "
2044 "Settings::communicate_vertices_to_p4est must be set."));
2047 parallel_forest,
nullptr, file_basename.c_str());
2052 template <
int dim,
int spacedim>
2057 this->cell_attached_data.n_attached_deserialize == 0,
2059 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2060 Assert(this->n_cells() > 0,
2061 ExcMessage(
"Can not save() an empty Triangulation."));
2067 this->signals.pre_distributed_save();
2069 if (this->my_subdomain == 0)
2071 std::string fname = std::string(filename) +
".info";
2072 std::ofstream f(fname);
2073 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2077 << this->cell_attached_data.pack_callbacks_fixed.size() <<
" "
2078 << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
2079 << this->n_cells(0) << std::endl;
2083 for (
const auto &cell_rel : this->local_cell_relations)
2086 Assert((cell_rel.second ==
2092 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2093 parallel_forest->global_num_quadrants,
2101 this->signals.post_distributed_save();
2106 template <
int dim,
int spacedim>
2111 this->n_cells() > 0,
2113 "load() only works if the Triangulation already contains a coarse mesh!"));
2115 this->n_levels() == 1,
2117 "Triangulation may only contain coarse cells when calling load()."));
2123 this->signals.pre_distributed_load();
2125 if (parallel_ghost !=
nullptr)
2129 parallel_ghost =
nullptr;
2132 parallel_forest =
nullptr;
2135 connectivity =
nullptr;
2137 unsigned int version, numcpus, attached_count_fixed,
2138 attached_count_variable, n_coarse_cells;
2140 std::string fname = std::string(filename) +
".info";
2141 std::ifstream f(fname);
2143 std::string firstline;
2144 getline(f, firstline);
2145 f >> version >> numcpus >> attached_count_fixed >>
2146 attached_count_variable >> n_coarse_cells;
2150 ExcMessage(
"Incompatible version found in .info file."));
2151 Assert(this->n_cells(0) == n_coarse_cells,
2152 ExcMessage(
"Number of coarse cells differ!"));
2156 this->cell_attached_data.n_attached_data_sets = 0;
2157 this->cell_attached_data.n_attached_deserialize =
2158 attached_count_fixed + attached_count_variable;
2162 this->mpi_communicator,
2180 copy_local_forest_to_triangulation();
2190 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2191 parallel_forest->global_num_quadrants,
2192 parallel_forest->local_num_quadrants,
2194 attached_count_fixed,
2195 attached_count_variable);
2198 this->signals.post_distributed_load();
2200 this->update_periodic_face_map();
2201 this->update_number_cache();
2206 template <
int dim,
int spacedim>
2209 const
bool autopartition)
2211 (void)autopartition;
2217 template <
int dim,
int spacedim>
2220 const typename ::
internal::p4est::
types<dim>::forest *forest)
2222 Assert(this->n_cells() > 0,
2224 "load() only works if the Triangulation already contains "
2226 Assert(this->n_cells() == forest->trees->elem_count,
2228 "Coarse mesh of the Triangulation does not match the one "
2229 "of the provided forest!"));
2232 if (parallel_ghost !=
nullptr)
2236 parallel_ghost =
nullptr;
2239 parallel_forest =
nullptr;
2245 typename ::internal::p4est::types<dim>::forest *temp =
2246 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2250 parallel_forest->connectivity = connectivity;
2251 parallel_forest->user_pointer =
this;
2255 copy_local_forest_to_triangulation();
2264 this->update_periodic_face_map();
2265 this->update_number_cache();
2270 template <
int dim,
int spacedim>
2274 Assert(parallel_forest !=
nullptr,
2276 "Can't produce a check sum when no forest is created yet."));
2281# if !DEAL_II_P4EST_VERSION_GTE(2, 8, 6, 0)
2297 template <
int dim,
int spacedim>
2299 const typename ::internal::p4est::types<dim>::forest
2302 Assert(parallel_forest !=
nullptr,
2303 ExcMessage(
"The forest has not been allocated yet."));
2304 return parallel_forest;
2309 template <
int dim,
int spacedim>
2311 typename ::internal::p4est::types<dim>::tree
2313 const int dealii_coarse_cell_index)
const
2315 const unsigned int tree_index =
2316 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2317 typename ::internal::p4est::types<dim>::tree *tree =
2318 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2319 sc_array_index(parallel_forest->trees, tree_index));
2333 std::integral_constant<int, 2>)
2335 const unsigned int dim = 2, spacedim = 2;
2343 std::vector<unsigned int> vertex_touch_count;
2345 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2348 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2349 const ::internal::p4est::types<2>::locidx num_vtt =
2350 std::accumulate(vertex_touch_count.begin(),
2351 vertex_touch_count.end(),
2357 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2360 (set_vertex_info ==
true ? this->n_vertices() : 0),
2365 set_vertex_and_cell_info(*
this,
2368 coarse_cell_to_p4est_tree_permutation,
2372 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2377 this->mpi_communicator,
2394 std::integral_constant<int, 2>)
2396 const unsigned int dim = 2, spacedim = 3;
2404 std::vector<unsigned int> vertex_touch_count;
2406 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2409 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2410 const ::internal::p4est::types<2>::locidx num_vtt =
2411 std::accumulate(vertex_touch_count.begin(),
2412 vertex_touch_count.end(),
2418 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2421 (set_vertex_info ==
true ? this->n_vertices() : 0),
2426 set_vertex_and_cell_info(*
this,
2429 coarse_cell_to_p4est_tree_permutation,
2433 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2438 this->mpi_communicator,
2453 std::integral_constant<int, 3>)
2455 const int dim = 3, spacedim = 3;
2463 std::vector<unsigned int> vertex_touch_count;
2464 std::vector<std::list<
2465 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2467 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2468 const ::internal::p4est::types<2>::locidx num_vtt =
2469 std::accumulate(vertex_touch_count.begin(),
2470 vertex_touch_count.end(),
2473 std::vector<unsigned int> edge_touch_count;
2474 std::vector<std::list<
2475 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2477 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2478 const ::internal::p4est::types<2>::locidx num_ett =
2479 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2482 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2485 (set_vertex_info ==
true ? this->n_vertices() : 0),
2487 this->n_active_lines(),
2492 set_vertex_and_cell_info(*
this,
2495 coarse_cell_to_p4est_tree_permutation,
2524 const unsigned int deal_to_p4est_line_index[12] = {
2525 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2528 this->begin_active();
2529 cell != this->
end();
2532 const unsigned int index =
2533 coarse_cell_to_p4est_tree_permutation[cell->index()];
2534 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2536 deal_to_p4est_line_index[e]] =
2537 cell->line(e)->index();
2542 connectivity->ett_offset[0] = 0;
2543 std::partial_sum(edge_touch_count.begin(),
2544 edge_touch_count.end(),
2545 &connectivity->ett_offset[1]);
2547 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2550 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2552 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2556 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2557 unsigned int>>::const_iterator p =
2558 edge_to_cell[v].begin();
2559 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2561 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2562 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2563 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2564 deal_to_p4est_line_index[p->second];
2568 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2573 this->mpi_communicator,
2590 template <
int dim,
int spacedim>
2592 enforce_mesh_balance_over_periodic_boundaries(
2598 std::vector<bool> flags_before[2];
2602 std::vector<unsigned int> topological_vertex_numbering(
2604 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2605 topological_vertex_numbering[i] = i;
2623 using cell_iterator =
2625 typename std::map<std::pair<cell_iterator, unsigned int>,
2626 std::pair<std::pair<cell_iterator, unsigned int>,
2627 unsigned char>>::const_iterator it;
2632 const cell_iterator &cell_1 = it->first.first;
2633 const unsigned int face_no_1 = it->first.second;
2634 const cell_iterator &cell_2 = it->second.first.first;
2635 const unsigned int face_no_2 = it->second.first.second;
2636 const unsigned char combined_orientation = it->second.second;
2637 const auto [orientation, rotation, flip] =
2640 if (cell_1->level() == cell_2->level())
2642 for (
unsigned int v = 0;
2648 const unsigned int vface0 =
2650 v, orientation, flip, rotation);
2651 const unsigned int vi0 =
2652 topological_vertex_numbering[cell_1->face(face_no_1)
2653 ->vertex_index(vface0)];
2654 const unsigned int vi1 =
2655 topological_vertex_numbering[cell_2->face(face_no_2)
2657 const unsigned int min_index =
std::min(vi0, vi1);
2658 topological_vertex_numbering[cell_1->face(face_no_1)
2659 ->vertex_index(vface0)] =
2660 topological_vertex_numbering[cell_2->face(face_no_2)
2661 ->vertex_index(v)] =
2668 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2670 const unsigned int j = topological_vertex_numbering[i];
2678 bool continue_iterating =
true;
2679 std::vector<int> vertex_level(tria.
n_vertices());
2680 while (continue_iterating)
2684 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2688 for (; cell != endc; ++cell)
2690 if (cell->refine_flag_set())
2691 for (
const unsigned int vertex :
2693 vertex_level[topological_vertex_numbering
2694 [cell->vertex_index(vertex)]] =
2695 std::
max(vertex_level[topological_vertex_numbering
2696 [cell->vertex_index(vertex)]],
2698 else if (!cell->coarsen_flag_set())
2699 for (
const unsigned int vertex :
2701 vertex_level[topological_vertex_numbering
2702 [cell->vertex_index(vertex)]] =
2703 std::
max(vertex_level[topological_vertex_numbering
2704 [cell->vertex_index(vertex)]],
2715 for (
const unsigned int vertex :
2717 vertex_level[topological_vertex_numbering
2718 [cell->vertex_index(vertex)]] =
2719 std::
max(vertex_level[topological_vertex_numbering
2720 [cell->vertex_index(vertex)]],
2725 continue_iterating =
false;
2736 for (cell = tria.
last_active(); cell != endc; --cell)
2737 if (cell->refine_flag_set() ==
false)
2739 for (
const unsigned int vertex :
2741 if (vertex_level[topological_vertex_numbering
2742 [cell->vertex_index(vertex)]] >=
2746 cell->clear_coarsen_flag();
2751 if (vertex_level[topological_vertex_numbering
2752 [cell->vertex_index(vertex)]] >
2755 cell->set_refine_flag();
2756 continue_iterating =
true;
2758 for (
const unsigned int v :
2760 vertex_level[topological_vertex_numbering
2761 [cell->vertex_index(v)]] =
2763 vertex_level[topological_vertex_numbering
2764 [cell->vertex_index(v)]],
2775 for (
const auto &cell : tria.cell_iterators())
2778 if (cell->is_active())
2781 const unsigned int n_children = cell->n_children();
2782 unsigned int flagged_children = 0;
2783 for (
unsigned int child = 0; child < n_children; ++child)
2784 if (cell->child(child)->is_active() &&
2785 cell->child(child)->coarsen_flag_set())
2790 if (flagged_children < n_children)
2791 for (
unsigned int child = 0; child < n_children; ++child)
2792 if (cell->child(child)->is_active())
2793 cell->child(child)->clear_coarsen_flag();
2796 std::vector<bool> flags_after[2];
2799 return ((flags_before[0] != flags_after[0]) ||
2800 (flags_before[1] != flags_after[1]));
2806 template <
int dim,
int spacedim>
2825 template <
int dim,
int spacedim>
2846 if (settings & construct_multigrid_hierarchy)
2849 spacedim>::limit_level_difference_at_vertices;
2853 bool mesh_changed =
false;
2861 if (settings & mesh_reconstruction_after_repartitioning)
2862 while (this->n_levels() > 1)
2869 for (
const auto &cell :
2870 this->active_cell_iterators_on_level(this->n_levels() - 1))
2872 cell->set_coarsen_flag();
2890 if (parallel_ghost !=
nullptr)
2894 parallel_ghost =
nullptr;
2898 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2899 P4EST_CONNECT_CORNER) :
2900 typename ::internal::p4est::types<dim>::balance_type(
2901 P8EST_CONNECT_CORNER)));
2908 for (
const auto &cell : this->cell_iterators_on_level(0))
2913 for (
const auto &cell : this->cell_iterators_on_level(0))
2918 if (tree_exists_locally<dim, spacedim>(
2920 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2923 delete_all_children<dim, spacedim>(cell);
2924 if (cell->is_active())
2933 typename ::internal::p4est::types<dim>::quadrant
2935 typename ::internal::p4est::types<dim>::tree *tree =
2936 init_tree(cell->index());
2938 ::internal::p4est::init_coarse_quadrant<dim>(
2941 match_tree_recursively<dim, spacedim>(*tree,
2945 this->my_subdomain);
2952 typename ::internal::p4est::types<dim>::quadrant *quadr;
2954 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2956 for (
unsigned int g_idx = 0;
2957 g_idx < parallel_ghost->ghosts.elem_count;
2960 while (g_idx >=
static_cast<unsigned int>(
2961 parallel_ghost->proc_offsets[ghost_owner + 1]))
2963 while (g_idx >=
static_cast<unsigned int>(
2964 parallel_ghost->tree_offsets[ghost_tree + 1]))
2967 quadr =
static_cast<
2968 typename ::internal::p4est::types<dim>::quadrant *
>(
2969 sc_array_index(¶llel_ghost->ghosts, g_idx));
2971 unsigned int coarse_cell_index =
2972 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2974 match_quadrant<dim, spacedim>(
this,
2989 bool mesh_changed =
true;
2990 unsigned int loop_counter = 0;
2997 this->update_periodic_face_map();
3000 enforce_mesh_balance_over_periodic_boundaries(*
this);
3012 "parallel::distributed::Triangulation::copy_local_forest_to_triangulation() "
3013 "for periodic boundaries detected. Aborting."));
3015 while (mesh_changed);
3020 std::any_of(this->begin_active(),
3023 return cell.refine_flag_set() ||
3024 cell.coarsen_flag_set();
3042 while (mesh_changed);
3046 unsigned int num_ghosts = 0;
3048 for (
const auto &cell : this->active_cell_iterators())
3050 if (cell->subdomain_id() != this->my_subdomain &&
3055 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
3066 if (settings & construct_multigrid_hierarchy)
3072 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3075 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3077 if ((cell->is_active() &&
3078 cell->subdomain_id() ==
3079 this->locally_owned_subdomain()) ||
3080 (cell->has_children() &&
3081 cell->child(0)->level_subdomain_id() ==
3082 this->locally_owned_subdomain()))
3083 cell->set_level_subdomain_id(
3084 this->locally_owned_subdomain());
3088 cell->set_level_subdomain_id(
3096 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3097 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3098 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3100 for (
const auto &cell : this->cell_iterators_on_level(0))
3102 typename ::internal::p4est::types<dim>::quadrant
3104 const unsigned int tree_index =
3105 coarse_cell_to_p4est_tree_permutation[cell->index()];
3106 typename ::internal::p4est::types<dim>::tree *tree =
3107 init_tree(cell->index());
3109 ::internal::p4est::init_coarse_quadrant<dim>(
3112 determine_level_subdomain_id_recursively<dim, spacedim>(
3123 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3126 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3128 if (cell->has_children())
3129 for (
unsigned int c = 0;
3130 c < GeometryInfo<dim>::max_children_per_cell;
3133 if (cell->child(c)->level_subdomain_id() ==
3134 this->locally_owned_subdomain())
3139 cell->child(0)->level_subdomain_id();
3143 cell->set_level_subdomain_id(mark);
3159 const unsigned int total_local_cells = this->n_active_cells();
3160 (void)total_local_cells;
3164 Assert(
static_cast<unsigned int>(
3165 parallel_forest->local_num_quadrants) == total_local_cells,
3170 Assert(
static_cast<unsigned int>(
3171 parallel_forest->local_num_quadrants) <= total_local_cells,
3177 unsigned int n_owned = 0;
3178 for (
const auto &cell : this->active_cell_iterators())
3180 if (cell->subdomain_id() == this->my_subdomain)
3184 Assert(
static_cast<unsigned int>(
3185 parallel_forest->local_num_quadrants) == n_owned,
3190 this->smooth_grid = save_smooth;
3196 update_cell_relations();
3201 template <
int dim,
int spacedim>
3207 std::vector<Point<dim>> point{p};
3208 std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
3215 template <
int dim,
int spacedim>
3218 find_point_owner_rank(const
std::vector<
Point<dim>> &points)
3220# ifndef P4EST_SEARCH_LOCAL
3225 "This function is only available if p4est is version 2.2 and higher."));
3227 return std::vector<unsigned int>(1,
3231 AssertThrow(this->are_vertices_communicated_to_p4est(),
3233 "Vertices need to be communicated to p4est to use this "
3234 "function. This must explicitly be turned on in the "
3235 "settings of the triangulation's constructor."));
3238 for (
const auto &manifold_id : this->get_manifold_ids())
3243 "This function can only be used if the triangulation "
3244 "has no other manifold than a Cartesian (flat) manifold attached."));
3248 PartitionSearch<dim> partition_search;
3254 parallel_forest->user_pointer = &partition_search;
3260 sc_array_t *point_sc_array;
3264 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3267 for (
size_t i = 0; i < points.size(); ++i)
3272 double *this_sc_point =
3273 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3275 for (
unsigned int d = 0; d < dim; ++d)
3277 this_sc_point[d] = p(d);
3279 this_sc_point[dim] = -1.0;
3285 static_cast<int>(
false),
3286 &PartitionSearch<dim>::local_quadrant_fn,
3287 &PartitionSearch<dim>::local_point_fn,
3291 std::vector<types::subdomain_id> owner_rank(
3295 for (
size_t i = 0; i < points.size(); ++i)
3298 double *this_sc_point =
3299 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3300 Assert(this_sc_point[dim] >= 0. || this_sc_point[dim] == -1.,
3302 if (this_sc_point[dim] < 0.)
3310 parallel_forest->user_pointer =
this;
3313 sc_array_destroy_null(&point_sc_array);
3321 template <
int dim,
int spacedim>
3327 for (
const auto &cell : this->active_cell_iterators())
3328 if (cell->is_locally_owned() && cell->refine_flag_set())
3329 Assert(cell->refine_flag_set() ==
3332 "This class does not support anisotropic refinement"));
3337 if (this->n_levels() ==
3341 cell = this->begin_active(
3349 !(cell->refine_flag_set()),
3351 "Fatal Error: maximum refinement level of p4est reached."));
3355 this->prepare_coarsening_and_refinement();
3358 this->signals.pre_distributed_refinement();
3363 for (
const auto &cell : this->active_cell_iterators())
3364 if (cell->is_ghost() || cell->is_artificial())
3366 cell->clear_refine_flag();
3367 cell->clear_coarsen_flag();
3373 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3374 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3381 parallel_forest->user_pointer = &refine_and_coarsen_list;
3383 if (parallel_ghost !=
nullptr)
3387 parallel_ghost =
nullptr;
3392 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3397 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3404 parallel_forest->user_pointer =
this;
3410 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3411 P4EST_CONNECT_FULL) :
3412 typename ::internal::p4est::types<dim>::balance_type(
3413 P8EST_CONNECT_FULL)),
3418 update_cell_relations();
3422 this->signals.post_p4est_refinement();
3426 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3427 previous_global_first_quadrant;
3429 if (this->cell_attached_data.n_attached_data_sets > 0)
3431 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3432 std::memcpy(previous_global_first_quadrant.data(),
3433 parallel_forest->global_first_quadrant,
3435 typename ::internal::p4est::types<dim>::gloidx) *
3436 (parallel_forest->mpisize + 1));
3439 if (!(settings & no_automatic_repartitioning))
3443 if (this->signals.weight.empty())
3451 const std::vector<unsigned int> cell_weights = get_cell_weights();
3457 this->mpi_communicator) > 0,
3459 "The global sum of weights over all active cells "
3460 "is zero. Please verify how you generate weights."));
3462 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3467 parallel_forest->user_pointer = &partition_weights;
3473 &PartitionWeights<dim, spacedim>::cell_weight);
3477 parallel_forest, 0,
nullptr,
nullptr);
3479 parallel_forest->user_pointer =
this;
3484 if (this->cell_attached_data.n_attached_data_sets > 0)
3486 this->data_serializer.pack_data(
3487 this->local_cell_relations,
3488 this->cell_attached_data.pack_callbacks_fixed,
3489 this->cell_attached_data.pack_callbacks_variable,
3490 this->get_communicator());
3496 for (
const auto &cell : this->active_cell_iterators())
3498 cell->clear_refine_flag();
3499 cell->clear_coarsen_flag();
3504 copy_local_forest_to_triangulation();
3514 if (this->cell_attached_data.n_attached_data_sets > 0)
3516 this->execute_transfer(parallel_forest,
3517 previous_global_first_quadrant.data());
3520 this->data_serializer.unpack_cell_status(this->local_cell_relations);
3539 if (settings & construct_multigrid_hierarchy)
3541 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3543 std::vector<bool> active_verts =
3544 this->mark_locally_active_vertices_on_level(lvl);
3546 const unsigned int maybe_coarser_lvl =
3547 (lvl > 0) ? (lvl - 1) : lvl;
3549 cell = this->begin(maybe_coarser_lvl),
3550 endc = this->end(lvl);
3551 for (; cell != endc; ++cell)
3552 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3554 const bool is_level_artificial =
3555 (cell->level_subdomain_id() ==
3557 bool need_to_know =
false;
3558 for (
const unsigned int vertex :
3560 if (active_verts[cell->vertex_index(vertex)])
3562 need_to_know =
true;
3567 !need_to_know || !is_level_artificial,
3569 "Internal error: the owner of cell" +
3570 cell->id().to_string() +
3571 " is unknown even though it is needed for geometric multigrid."));
3577 this->update_periodic_face_map();
3578 this->update_number_cache();
3581 this->signals.post_distributed_refinement();
3586 template <
int dim,
int spacedim>
3591 for (
const auto &cell : this->active_cell_iterators())
3592 if (cell->is_locally_owned())
3594 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3596 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3600 this->signals.pre_distributed_repartition();
3604 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3605 previous_global_first_quadrant;
3607 if (this->cell_attached_data.n_attached_data_sets > 0)
3609 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3610 std::memcpy(previous_global_first_quadrant.data(),
3611 parallel_forest->global_first_quadrant,
3613 typename ::internal::p4est::types<dim>::gloidx) *
3614 (parallel_forest->mpisize + 1));
3617 if (this->signals.weight.empty())
3629 const std::vector<unsigned int> cell_weights = get_cell_weights();
3635 this->mpi_communicator) > 0,
3637 "The global sum of weights over all active cells "
3638 "is zero. Please verify how you generate weights."));
3640 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3645 parallel_forest->user_pointer = &partition_weights;
3651 &PartitionWeights<dim, spacedim>::cell_weight);
3654 parallel_forest->user_pointer =
this;
3658 if (this->cell_attached_data.n_attached_data_sets > 0)
3660 this->data_serializer.pack_data(
3661 this->local_cell_relations,
3662 this->cell_attached_data.pack_callbacks_fixed,
3663 this->cell_attached_data.pack_callbacks_variable,
3664 this->get_communicator());
3669 copy_local_forest_to_triangulation();
3679 if (this->cell_attached_data.n_attached_data_sets > 0)
3681 this->execute_transfer(parallel_forest,
3682 previous_global_first_quadrant.data());
3685 this->update_periodic_face_map();
3688 this->update_number_cache();
3691 this->signals.post_distributed_repartition();
3696 template <
int dim,
int spacedim>
3698 const std::vector<types::global_dof_index>
3702 return p4est_tree_to_coarse_cell_permutation;
3707 template <
int dim,
int spacedim>
3709 const std::vector<types::global_dof_index>
3713 return coarse_cell_to_p4est_tree_permutation;
3718 template <
int dim,
int spacedim>
3721 mark_locally_active_vertices_on_level(const
int level)
const
3725 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3726 for (
const auto &cell : this->cell_iterators_on_level(
level))
3727 if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3729 marked_vertices[cell->vertex_index(v)] =
true;
3746 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3747 for (
const auto &it : this->get_periodic_face_map())
3750 const unsigned int face_no_1 = it.first.second;
3752 const unsigned int face_no_2 = it.second.first.second;
3753 const unsigned char combined_orientation = it.second.second;
3754 const auto [orientation, rotation, flip] =
3757 if (cell_1->level() ==
level && cell_2->level() ==
level)
3759 for (
unsigned int v = 0;
3765 const unsigned int vface0 =
3767 v, orientation, flip, rotation);
3768 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3770 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3772 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3774 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3780 return marked_vertices;
3785 template <
int dim,
int spacedim>
3788 coarse_cell_id_to_coarse_cell_index(
3789 const
types::coarse_cell_id coarse_cell_id)
const
3791 return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3796 template <
int dim,
int spacedim>
3800 const unsigned int coarse_cell_index)
const
3802 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3807 template <
int dim,
int spacedim>
3811 &periodicity_vector)
3813 Assert(triangulation_has_content ==
true,
3815 Assert(this->n_levels() == 1,
3816 ExcMessage(
"The triangulation is refined!"));
3823 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
3825 for (
const auto &face_pair : periodicity_vector)
3829 const unsigned int face_left = face_pair.face_idx[0];
3830 const unsigned int face_right = face_pair.face_idx[1];
3833 const unsigned int tree_left =
3834 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3835 const unsigned int tree_right =
3836 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3844 unsigned int p4est_orientation = 0;
3849 face_pair.orientation ==
3856 const unsigned int face_idx_list[] = {face_left, face_right};
3857 const cell_iterator cell_list[] = {first_cell, second_cell};
3858 unsigned int lower_idx, higher_idx;
3859 unsigned char orientation;
3860 if (face_left <= face_right)
3865 face_reference_cell.get_inverse_combined_orientation(
3866 face_pair.orientation);
3872 orientation = face_pair.orientation;
3877 unsigned int first_p4est_idx_on_cell =
3878 p8est_face_corners[face_idx_list[lower_idx]][0];
3879 unsigned int first_dealii_idx_on_face =
3881 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3884 const unsigned int first_dealii_idx_on_cell =
3886 face_idx_list[lower_idx],
3888 cell_list[lower_idx]->face_orientation(
3889 face_idx_list[lower_idx]),
3890 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3891 cell_list[lower_idx]->face_rotation(
3892 face_idx_list[lower_idx]));
3893 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3895 first_dealii_idx_on_face = i;
3903 const unsigned int second_dealii_idx_on_face =
3904 reference_cell.standard_to_real_face_vertex(
3905 first_dealii_idx_on_face,
3906 face_idx_list[lower_idx],
3908 const unsigned int second_dealii_idx_on_cell =
3909 reference_cell.face_to_cell_vertices(
3910 face_idx_list[higher_idx],
3911 second_dealii_idx_on_face,
3912 cell_list[higher_idx]->combined_face_orientation(
3913 face_idx_list[higher_idx]));
3915 const unsigned int second_p4est_idx_on_face =
3916 p8est_corner_face_corners[second_dealii_idx_on_cell]
3917 [face_idx_list[higher_idx]];
3918 p4est_orientation = second_p4est_idx_on_face;
3938 this->mpi_communicator,
3949 copy_local_forest_to_triangulation();
3960 this->update_number_cache();
3965 template <
int dim,
int spacedim>
3976 this->cell_attached_data.n_attached_data_sets) +
3983 coarse_cell_to_p4est_tree_permutation) +
3985 p4est_tree_to_coarse_cell_permutation) +
3986 memory_consumption_p4est();
3993 template <
int dim,
int spacedim>
3997 return ::internal::p4est::functions<dim>::forest_memory_used(
4005 template <
int dim,
int spacedim>
4012 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
4014 (other_tria.n_global_levels() == 1),
4025 const typename ::Triangulation<dim, spacedim>::DistortedCellList
4033 if (const ::parallel::distributed::Triangulation<dim, spacedim>
4034 *other_distributed =
4035 dynamic_cast<const ::parallel::distributed::
4036 Triangulation<dim, spacedim> *
>(&other_tria))
4039 settings = other_distributed->settings;
4040 triangulation_has_content =
4041 other_distributed->triangulation_has_content;
4042 coarse_cell_to_p4est_tree_permutation =
4043 other_distributed->coarse_cell_to_p4est_tree_permutation;
4044 p4est_tree_to_coarse_cell_permutation =
4045 other_distributed->p4est_tree_to_coarse_cell_permutation;
4048 typename ::internal::p4est::types<dim>::connectivity
4049 *temp_connectivity =
const_cast<
4050 typename ::internal::p4est::types<dim>::connectivity *
>(
4051 other_distributed->connectivity);
4053 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
4056 typename ::internal::p4est::types<dim>::forest *temp_forest =
4057 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
4058 other_distributed->parallel_forest);
4062 parallel_forest->connectivity = connectivity;
4063 parallel_forest->user_pointer =
this;
4067 triangulation_has_content =
true;
4068 setup_coarse_cell_to_p4est_tree_permutation();
4069 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
4074 copy_local_forest_to_triangulation();
4083 this->update_periodic_face_map();
4084 this->update_number_cache();
4089 template <
int dim,
int spacedim>
4094 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
4095 this->local_cell_relations.shrink_to_fit();
4098 for (
const auto &cell : this->cell_iterators_on_level(0))
4101 if (tree_exists_locally<dim, spacedim>(
4103 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4107 typename ::internal::p4est::types<dim>::quadrant
4109 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4112 typename ::internal::p4est::types<dim>::tree *tree =
4113 init_tree(cell->index());
4115 update_cell_relations_recursively<dim, spacedim>(
4116 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4122 template <
int dim,
int spacedim>
4129 Assert(this->local_cell_relations.size() ==
4130 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4135 std::vector<unsigned int> weights;
4136 weights.reserve(this->local_cell_relations.size());
4144 for (
const auto &cell_rel : this->local_cell_relations)
4146 const auto &cell_it = cell_rel.first;
4147 const auto &cell_status = cell_rel.second;
4149 weights.push_back(this->signals.weight(cell_it, cell_status));
4157 template <
int spacedim>
4173 template <
int spacedim>
4182 template <
int spacedim>
4184 const std::vector<types::global_dof_index>
4188 static std::vector<types::global_dof_index> a;
4194 template <
int spacedim>
4196 std::map<
unsigned int,
4199 const unsigned int )
const
4203 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4208 template <
int spacedim>
4211 mark_locally_active_vertices_on_level(const
unsigned int)
const
4214 return std::vector<bool>();
4219 template <
int spacedim>
4222 coarse_cell_id_to_coarse_cell_index(const
types::coarse_cell_id)
const
4230 template <
int spacedim>
4234 const unsigned int)
const
4242 template <
int spacedim>
4251 template <
int spacedim>
4260 template <
int spacedim>
4269 template <
int spacedim>
4279 template <
int spacedim>
4289 template <
int spacedim>
4306 namespace distributed
4308 template <
int dim,
int spacedim>
4316#ifdef DEAL_II_WITH_P4EST
4326 const auto &cell = pair.first;
4327 const auto &status = pair.second;
4333 cell->clear_refine_flag();
4334 cell->clear_coarsen_flag();
4339 cell->clear_coarsen_flag();
4340 cell->set_refine_flag();
4345 for (
const auto &child : cell->child_iterators())
4347 child->clear_refine_flag();
4348 child->set_coarsen_flag();
4367 template <
int dim,
int spacedim>
4370#ifdef DEAL_II_WITH_P4EST
4371 if (distributed_tria)
4374 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4375 distributed_tria->load_refine_flags(saved_refine_flags);
4379 (void)distributed_tria;
@ children_will_be_coarsened
static constexpr unsigned char default_combined_face_orientation()
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator > > &)
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, unsigned char > > & get_periodic_face_map() const
active_cell_iterator last_active() const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
cell_iterator end() const
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const
virtual void execute_coarsening_and_refinement()
virtual bool prepare_coarsening_and_refinement()
void save_refine_flags(std::ostream &out) const
unsigned int n_vertices() const
void save_coarsen_flags(std::ostream &out) const
active_cell_iterator begin_active(const unsigned int level=0) const
virtual std::size_t memory_consumption() const override
virtual void clear() override
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
~TemporarilyMatchRefineFlags()
const SmartPointer< ::parallel::distributed::Triangulation< dim, spacedim > > distributed_tria
std::vector< bool > saved_refine_flags
std::vector< bool > saved_coarsen_flags
virtual void clear() override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
typename ::Triangulation< dim, spacedim >::active_cell_iterator active_cell_iterator
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell & get_hypercube()
VectorType::value_type * end(VectorType &V)
std::enable_if_t< is_mpi_type< T >==false, T > broadcast(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
bool tree_exists_locally(const typename types< dim >::forest *parallel_forest, const typename types< dim >::topidx coarse_grid_cell)
void exchange_refinement_flags(::parallel::distributed::Triangulation< dim, spacedim > &tria)
std::tuple< bool, bool, bool > split_face_orientation(const unsigned char combined_face_orientation)
const types::subdomain_id artificial_subdomain_id
const types::subdomain_id invalid_subdomain_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static 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)