39#ifdef DEAL_II_WITH_P4EST
43 template <
int dim,
int spacedim>
45 get_vertex_to_cell_mappings(
47 std::vector<unsigned int> &vertex_touch_count,
48 std::vector<std::list<
50 unsigned int>>> &vertex_to_cell)
55 for (
const auto &cell :
triangulation.active_cell_iterators())
58 ++vertex_touch_count[cell->vertex_index(v)];
59 vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
65 template <
int dim,
int spacedim>
67 get_edge_to_cell_mappings(
69 std::vector<unsigned int> &edge_touch_count,
70 std::vector<std::list<
72 unsigned int>>> &edge_to_cell)
79 for (
const auto &cell :
triangulation.active_cell_iterators())
80 for (unsigned
int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
82 ++edge_touch_count[cell->line(l)->index()];
83 edge_to_cell[cell->line(l)->index()].emplace_back(cell, l);
93 template <
int dim,
int spacedim>
95 set_vertex_and_cell_info(
97 const std::vector<unsigned int> &vertex_touch_count,
98 const std::vector<std::list<
100 unsigned int>>> &vertex_to_cell,
101 const std::vector<types::global_dof_index>
102 &coarse_cell_to_p4est_tree_permutation,
103 const bool set_vertex_info,
119 if (set_vertex_info ==
true)
120 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
122 connectivity->vertices[3 * v] =
triangulation.get_vertices()[v][0];
123 connectivity->vertices[3 * v + 1] =
125 connectivity->vertices[3 * v + 2] =
137 for (; cell != endc; ++cell)
139 const unsigned int index =
140 coarse_cell_to_p4est_tree_permutation[cell->index()];
144 if (set_vertex_info ==
true)
147 v] = cell->vertex_index(v);
150 v] = cell->vertex_index(v);
156 if (cell->face(f)->at_boundary() == false)
159 coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->
index()];
163 coarse_cell_to_p4est_tree_permutation[cell->index()];
168 if (cell->face(f)->at_boundary() == false)
174 connectivity->tree_to_face
176 cell->neighbor_of_neighbor(f);
195 connectivity->tree_to_face[
index * 6 + f] =
196 cell->neighbor_of_neighbor(f);
198 unsigned int face_idx_list[2] = {
199 f, cell->neighbor_of_neighbor(f)};
201 cell_list[2] = {cell, cell->neighbor(f)};
202 unsigned int smaller_idx = 0;
204 if (f > cell->neighbor_of_neighbor(f))
207 unsigned int larger_idx = (smaller_idx + 1) % 2;
215 unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
217 face_idx_list[smaller_idx],
219 cell_list[smaller_idx]->face_orientation(
220 face_idx_list[smaller_idx]),
221 cell_list[smaller_idx]->face_flip(
222 face_idx_list[smaller_idx]),
223 cell_list[smaller_idx]->face_rotation(
224 face_idx_list[smaller_idx])));
228 for (
unsigned int i = 0;
229 i < GeometryInfo<dim>::vertices_per_face;
233 cell_list[larger_idx]->vertex_index(
235 face_idx_list[larger_idx], i));
244 connectivity->tree_to_face[
index * 6 + f] += 6 * v;
258 connectivity->ctt_offset[0] = 0;
259 std::partial_sum(vertex_touch_count.begin(),
260 vertex_touch_count.end(),
261 &connectivity->ctt_offset[1]);
264 std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
269 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
271 Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
275 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
276 unsigned int>>::const_iterator
p =
277 vertex_to_cell[v].begin();
278 for (
unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++
p)
280 connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
281 coarse_cell_to_p4est_tree_permutation[
p->first->index()];
282 connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
290 template <
int dim,
int spacedim>
296 Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
298 return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
299 (coarse_grid_cell <= parallel_forest->last_local_tree));
303 template <
int dim,
int spacedim>
305 delete_all_children_and_self(
308 if (cell->has_children())
309 for (
unsigned int c = 0; c < cell->n_children(); ++c)
310 delete_all_children_and_self<dim, spacedim>(cell->child(c));
312 cell->set_coarsen_flag();
317 template <
int dim,
int spacedim>
322 if (cell->has_children())
323 for (
unsigned int c = 0; c < cell->n_children(); ++c)
324 delete_all_children_and_self<dim, spacedim>(cell->child(c));
328 template <
int dim,
int spacedim>
330 determine_level_subdomain_id_recursively(
337 const std::vector<std::vector<bool>> &marked_vertices)
346 if (marked_vertices[dealii_cell->level()]
347 [dealii_cell->vertex_index(v)])
366 if (!used && dealii_cell->is_active() &&
367 dealii_cell->is_artificial() ==
false &&
372 if (marked_vertices[dealii_cell->level() + 1]
373 [dealii_cell->vertex_index(v)])
382 if (!used && dealii_cell->is_active() &&
383 dealii_cell->is_artificial() ==
false && dealii_cell->level() > 0)
387 if (marked_vertices[dealii_cell->level() - 1]
388 [dealii_cell->vertex_index(v)])
399 &forest, tree_index, &p4est_cell, my_subdomain);
400 Assert((owner != -2) && (owner != -1),
402 dealii_cell->set_level_subdomain_id(owner);
406 if (dealii_cell->has_children())
410 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
415 P4EST_QUADRANT_INIT(&p4est_child[c]);
418 P8EST_QUADRANT_INIT(&p4est_child[c]);
428 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
431 determine_level_subdomain_id_recursively<dim, spacedim>(
434 dealii_cell->child(c),
444 template <
int dim,
int spacedim>
446 match_tree_recursively(
454 if (sc_array_bsearch(
const_cast<sc_array_t *
>(&tree.quadrants),
460 delete_all_children<dim, spacedim>(dealii_cell);
461 if (dealii_cell->is_active())
462 dealii_cell->set_subdomain_id(my_subdomain);
471 if (dealii_cell->is_active())
472 dealii_cell->set_refine_flag();
477 for (
unsigned int c = 0;
478 c < GeometryInfo<dim>::max_children_per_cell;
483 P4EST_QUADRANT_INIT(&p4est_child[c]);
486 P8EST_QUADRANT_INIT(&p4est_child[c]);
496 for (
unsigned int c = 0;
497 c < GeometryInfo<dim>::max_children_per_cell;
502 &p4est_child[c]) ==
false)
508 delete_all_children<dim, spacedim>(dealii_cell->child(c));
509 dealii_cell->child(c)->recursively_set_subdomain_id(
516 match_tree_recursively<dim, spacedim>(tree,
517 dealii_cell->child(c),
527 template <
int dim,
int spacedim>
530 const ::Triangulation<dim, spacedim> *tria,
531 unsigned int dealii_index,
535 const int l = ghost_quadrant.level;
537 for (
int i = 0; i <
l; ++i)
542 if (cell->is_active())
544 cell->clear_coarsen_flag();
545 cell->set_refine_flag();
552 dealii_index = cell->child_index(child_id);
558 if (cell->has_children())
559 delete_all_children<dim, spacedim>(cell);
562 cell->clear_coarsen_flag();
563 cell->set_subdomain_id(ghost_owner);
566 template <
int dim,
int spacedim>
578 if (cell->refine_flag_set())
580 if (cell->coarsen_flag_set())
587 const std::uint8_t &flag) ->
void {
588 cell->clear_coarsen_flag();
589 cell->clear_refine_flag();
591 cell->set_refine_flag();
593 cell->set_coarsen_flag();
601# ifdef P4EST_SEARCH_LOCAL
603 class PartitionSearch
611 PartitionSearch(
const PartitionSearch<dim> &other) =
delete;
613 PartitionSearch<dim> &
614 operator=(
const PartitionSearch<dim> &other) =
delete;
670 quad_length_on_level);
673 initialize_mapping();
676 map_real_to_unit_cell(
const Point<dim> &p)
const;
679 is_in_this_quadrant(
const Point<dim> &p)
const;
682 std::vector<Point<dim>> cell_vertices;
690 bool are_vertices_initialized;
692 bool is_reference_mapping_initialized;
698 QuadrantData quadrant_data;
705 PartitionSearch<dim>::local_quadrant_fn(
719 PartitionSearch<dim> *this_object =
720 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
724 quad_length_on_level =
726 (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
730 this_object->quadrant_data.set_cell_vertices(forest,
733 quad_length_on_level);
736 this_object->quadrant_data.initialize_mapping();
746 PartitionSearch<dim>::local_point_fn(
759 PartitionSearch<dim> *this_object =
760 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
763 double *this_point_dptr =
static_cast<double *
>(
point);
766 (dim == 2 ?
Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
767 Point<dim>(this_point_dptr[0],
769 this_point_dptr[2]));
772 const bool is_in_this_quadrant =
773 this_object->quadrant_data.is_in_this_quadrant(this_point);
777 if (!is_in_this_quadrant)
786 if (rank_begin < rank_end)
794 this_point_dptr[dim] =
static_cast<double>(rank_begin);
804 PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
807 const Point<dim> p_ref = map_real_to_unit_cell(p);
816 PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
819 Assert(is_reference_mapping_initialized,
821 "Cell vertices and mapping coefficients must be fully "
822 "initialized before transforming a point to the unit cell."));
828 for (
unsigned int alpha = 0;
829 alpha < GeometryInfo<dim>::vertices_per_cell;
835 p_out += (quadrant_mapping_matrix(alpha, 0) +
836 quadrant_mapping_matrix(alpha, 1) *
p(0) +
837 quadrant_mapping_matrix(alpha, 2) *
p(1) +
838 quadrant_mapping_matrix(alpha, 3) *
p(0) *
p(1)) *
844 for (
unsigned int alpha = 0;
845 alpha < GeometryInfo<dim>::vertices_per_cell;
851 p_out += (quadrant_mapping_matrix(alpha, 0) +
852 quadrant_mapping_matrix(alpha, 1) *
p(0) +
853 quadrant_mapping_matrix(alpha, 2) *
p(1) +
854 quadrant_mapping_matrix(alpha, 3) *
p(2) +
855 quadrant_mapping_matrix(alpha, 4) *
p(0) *
p(1) +
856 quadrant_mapping_matrix(alpha, 5) *
p(1) *
p(2) +
857 quadrant_mapping_matrix(alpha, 6) *
p(0) *
p(2) +
858 quadrant_mapping_matrix(alpha, 7) *
p(0) *
p(1) *
p(2)) *
868 PartitionSearch<dim>::QuadrantData::QuadrantData()
870 , quadrant_mapping_matrix(
GeometryInfo<dim>::vertices_per_cell,
872 , are_vertices_initialized(false)
873 , is_reference_mapping_initialized(false)
880 PartitionSearch<dim>::QuadrantData::initialize_mapping()
883 are_vertices_initialized,
885 "Cell vertices must be initialized before the cell mapping can be filled."));
892 for (
unsigned int alpha = 0;
893 alpha < GeometryInfo<dim>::vertices_per_cell;
897 point_matrix(0, alpha) = 1;
898 point_matrix(1, alpha) = cell_vertices[alpha](0);
899 point_matrix(2, alpha) = cell_vertices[alpha](1);
900 point_matrix(3, alpha) =
901 cell_vertices[alpha](0) * cell_vertices[alpha](1);
908 quadrant_mapping_matrix.invert(point_matrix);
912 for (
unsigned int alpha = 0;
913 alpha < GeometryInfo<dim>::vertices_per_cell;
917 point_matrix(0, alpha) = 1;
918 point_matrix(1, alpha) = cell_vertices[alpha](0);
919 point_matrix(2, alpha) = cell_vertices[alpha](1);
920 point_matrix(3, alpha) = cell_vertices[alpha](2);
921 point_matrix(4, alpha) =
922 cell_vertices[alpha](0) * cell_vertices[alpha](1);
923 point_matrix(5, alpha) =
924 cell_vertices[alpha](1) * cell_vertices[alpha](2);
925 point_matrix(6, alpha) =
926 cell_vertices[alpha](0) * cell_vertices[alpha](2);
927 point_matrix(7, alpha) = cell_vertices[alpha](0) *
928 cell_vertices[alpha](1) *
929 cell_vertices[alpha](2);
936 quadrant_mapping_matrix.invert(point_matrix);
939 is_reference_mapping_initialized =
true;
946 PartitionSearch<2>::QuadrantData::set_cell_vertices(
951 quad_length_on_level)
953 constexpr unsigned int dim = 2;
957 double corner_point[dim + 1] = {0};
960 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
962 for (
unsigned int d = 0;
d < dim; ++
d)
964 cell_vertices[vertex_index](
d) = corner_point[d];
974 unsigned int vertex_index = 0;
976 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
979 copy_vertex(vertex_index);
986 forest->connectivity,
988 quadrant->x + quad_length_on_level,
993 copy_vertex(vertex_index);
1000 forest->connectivity,
1003 quadrant->y + quad_length_on_level,
1007 copy_vertex(vertex_index);
1014 forest->connectivity,
1016 quadrant->x + quad_length_on_level,
1017 quadrant->y + quad_length_on_level,
1021 copy_vertex(vertex_index);
1023 are_vertices_initialized =
true;
1030 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1035 quad_length_on_level)
1037 constexpr unsigned int dim = 3;
1039 double corner_point[dim] = {0};
1042 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1044 for (
unsigned int d = 0;
d < dim; ++
d)
1046 cell_vertices[vertex_index](
d) = corner_point[d];
1048 corner_point[
d] = 0;
1056 unsigned int vertex_index = 0;
1058 forest->connectivity,
1066 copy_vertex(vertex_index);
1074 forest->connectivity,
1076 quadrant->x + quad_length_on_level,
1082 copy_vertex(vertex_index);
1089 forest->connectivity,
1092 quadrant->y + quad_length_on_level,
1097 copy_vertex(vertex_index);
1104 forest->connectivity,
1106 quadrant->x + quad_length_on_level,
1107 quadrant->y + quad_length_on_level,
1112 copy_vertex(vertex_index);
1119 forest->connectivity,
1123 quadrant->z + quad_length_on_level,
1127 copy_vertex(vertex_index);
1134 forest->connectivity,
1136 quadrant->x + quad_length_on_level,
1138 quadrant->z + quad_length_on_level,
1142 copy_vertex(vertex_index);
1149 forest->connectivity,
1152 quadrant->y + quad_length_on_level,
1153 quadrant->z + quad_length_on_level,
1157 copy_vertex(vertex_index);
1164 forest->connectivity,
1166 quadrant->x + quad_length_on_level,
1167 quadrant->y + quad_length_on_level,
1168 quadrant->z + quad_length_on_level,
1172 copy_vertex(vertex_index);
1175 are_vertices_initialized =
true;
1185 template <
int dim,
int spacedim>
1186 class RefineAndCoarsenList
1190 const std::vector<types::global_dof_index>
1191 &p4est_tree_to_coarse_cell_permutation,
1219 pointers_are_at_end()
const;
1222 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1223 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1224 const_iterator current_refine_pointer;
1226 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1227 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1228 const_iterator current_coarsen_pointer;
1239 template <
int dim,
int spacedim>
1241 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1243 return ((current_refine_pointer == refine_list.end()) &&
1244 (current_coarsen_pointer == coarsen_list.end()));
1249 template <
int dim,
int spacedim>
1250 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1252 const std::vector<types::global_dof_index>
1253 &p4est_tree_to_coarse_cell_permutation,
1257 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1258 for (
const auto &cell :
triangulation.active_cell_iterators())
1261 if (cell->subdomain_id() != my_subdomain)
1264 if (cell->refine_flag_set())
1266 else if (cell->coarsen_flag_set())
1270 refine_list.reserve(n_refine_flags);
1271 coarsen_list.reserve(n_coarsen_flags);
1283 unsigned int coarse_cell_index =
1284 p4est_tree_to_coarse_cell_permutation[c];
1293 p4est_cell.p.which_tree = c;
1294 build_lists(cell, p4est_cell, my_subdomain);
1302 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1303 Assert(refine_list[i].
p.which_tree >= refine_list[i - 1].p.which_tree,
1305 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1306 Assert(coarsen_list[i].
p.which_tree >= coarsen_list[i - 1].p.which_tree,
1309 current_refine_pointer = refine_list.begin();
1310 current_coarsen_pointer = coarsen_list.begin();
1315 template <
int dim,
int spacedim>
1317 RefineAndCoarsenList<dim, spacedim>::build_lists(
1322 if (cell->is_active())
1324 if (cell->subdomain_id() == my_subdomain)
1326 if (cell->refine_flag_set())
1327 refine_list.push_back(p4est_cell);
1328 else if (cell->coarsen_flag_set())
1329 coarsen_list.push_back(p4est_cell);
1336 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1341 P4EST_QUADRANT_INIT(&p4est_child[c]);
1344 P8EST_QUADRANT_INIT(&p4est_child[c]);
1351 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1354 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1355 build_lists(cell->child(c), p4est_child[c], my_subdomain);
1361 template <
int dim,
int spacedim>
1363 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1368 RefineAndCoarsenList<dim, spacedim> *this_object =
1369 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1370 forest->user_pointer);
1374 if (this_object->current_refine_pointer == this_object->refine_list.end())
1377 Assert(coarse_cell_index <=
1378 this_object->current_refine_pointer->p.which_tree,
1383 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1387 Assert(coarse_cell_index <=
1388 this_object->current_refine_pointer->p.which_tree,
1394 quadrant, &*this_object->current_refine_pointer) <= 0,
1399 quadrant, &*this_object->current_refine_pointer))
1401 ++this_object->current_refine_pointer;
1411 template <
int dim,
int spacedim>
1413 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1418 RefineAndCoarsenList<dim, spacedim> *this_object =
1419 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1420 forest->user_pointer);
1424 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1427 Assert(coarse_cell_index <=
1428 this_object->current_coarsen_pointer->p.which_tree,
1433 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1437 Assert(coarse_cell_index <=
1438 this_object->current_coarsen_pointer->p.which_tree,
1444 children[0], &*this_object->current_coarsen_pointer) <= 0,
1450 children[0], &*this_object->current_coarsen_pointer))
1453 ++this_object->current_coarsen_pointer;
1457 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1461 children[c], &*this_object->current_coarsen_pointer),
1463 ++this_object->current_coarsen_pointer;
1481 template <
int dim,
int spacedim>
1482 class PartitionWeights
1490 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1505 std::vector<unsigned int> cell_weights_list;
1506 std::vector<unsigned int>::const_iterator current_pointer;
1510 template <
int dim,
int spacedim>
1511 PartitionWeights<dim, spacedim>::PartitionWeights(
1512 const std::vector<unsigned int> &cell_weights)
1513 : cell_weights_list(cell_weights)
1517 current_pointer = cell_weights_list.begin();
1521 template <
int dim,
int spacedim>
1523 PartitionWeights<dim, spacedim>::cell_weight(
1532 PartitionWeights<dim, spacedim> *this_object =
1533 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1535 Assert(this_object->current_pointer >=
1536 this_object->cell_weights_list.begin(),
1538 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1544 const unsigned int weight = *this_object->current_pointer;
1545 ++this_object->current_pointer;
1547 Assert(weight <
static_cast<unsigned int>(std::numeric_limits<int>::max()),
1548 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1549 "weights for cells. The weight provided here exceeds "
1550 "the maximum value represented as a 'signed int'."));
1551 return static_cast<int>(weight);
1554 template <
int dim,
int spacedim>
1555 using cell_relation_t =
typename std::pair<
1556 typename ::Triangulation<dim, spacedim>::cell_iterator,
1568 template <
int dim,
int spacedim>
1570 add_single_cell_relation(
1571 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1572 const typename ::internal::p4est::types<dim>::tree &tree,
1573 const unsigned int idx,
1577 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1583 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1597 template <
int dim,
int spacedim>
1599 update_cell_relations_recursively(
1600 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1601 const typename ::internal::p4est::types<dim>::tree &tree,
1603 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1606 const int idx = sc_array_bsearch(
1607 const_cast<sc_array_t *
>(&tree.quadrants),
1612 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1614 &p4est_cell) ==
false))
1619 const bool p4est_has_children = (idx == -1);
1620 if (p4est_has_children && dealii_cell->has_children())
1623 typename ::internal::p4est::types<dim>::quadrant
1626 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1631 P4EST_QUADRANT_INIT(&p4est_child[c]);
1634 P8EST_QUADRANT_INIT(&p4est_child[c]);
1641 &p4est_cell, p4est_child);
1643 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1646 update_cell_relations_recursively<dim, spacedim>(
1647 cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1650 else if (!p4est_has_children && !dealii_cell->has_children())
1654 add_single_cell_relation<dim, spacedim>(
1657 else if (p4est_has_children)
1665 typename ::internal::p4est::types<dim>::quadrant
1667 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1672 P4EST_QUADRANT_INIT(&p4est_child[c]);
1675 P8EST_QUADRANT_INIT(&p4est_child[c]);
1682 &p4est_cell, p4est_child);
1690 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1693 child_idx = sc_array_bsearch(
1694 const_cast<sc_array_t *
>(&tree.quadrants),
1701 add_single_cell_relation<dim, spacedim>(
1702 cell_rel, tree, child_idx, dealii_cell, cell_status);
1710 add_single_cell_relation<dim, spacedim>(
1724 namespace distributed
1727 template <
int dim,
int spacedim>
1739 (settings & construct_multigrid_hierarchy) ?
1743 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1746 , settings(settings)
1747 , triangulation_has_content(false)
1748 , connectivity(
nullptr)
1749 , parallel_forest(
nullptr)
1751 parallel_ghost =
nullptr;
1756 template <
int dim,
int spacedim>
1777 template <
int dim,
int spacedim>
1790 const typename ::Triangulation<dim, spacedim>::DistortedCellList
1799 this->all_reference_cells_are_hyper_cube(),
1801 "The class parallel::distributed::Triangulation only supports meshes "
1802 "consisting only of hypercube-like cells."));
1807 triangulation_has_content =
true;
1809 setup_coarse_cell_to_p4est_tree_permutation();
1811 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1815 copy_local_forest_to_triangulation();
1824 this->update_periodic_face_map();
1825 this->update_number_cache();
1830 template <
int dim,
int spacedim>
1836 (void)construction_data;
1843 template <
int dim,
int spacedim>
1847 triangulation_has_content =
false;
1849 if (parallel_ghost !=
nullptr)
1853 parallel_ghost =
nullptr;
1856 if (parallel_forest !=
nullptr)
1859 parallel_forest =
nullptr;
1862 if (connectivity !=
nullptr)
1866 connectivity =
nullptr;
1869 coarse_cell_to_p4est_tree_permutation.resize(0);
1870 p4est_tree_to_coarse_cell_permutation.resize(0);
1874 this->update_number_cache();
1879 template <
int dim,
int spacedim>
1890 template <
int dim,
int spacedim>
1901 template <
int dim,
int spacedim>
1907 *previous_global_first_quadrant)
1909 Assert(this->data_serializer.sizes_fixed_cumulative.size() > 0,
1913 this->data_serializer.dest_data_fixed.resize(
1914 parallel_forest->local_num_quadrants *
1915 this->data_serializer.sizes_fixed_cumulative.back());
1918 typename ::internal::p4est::types<dim>::transfer_context
1922 parallel_forest->global_first_quadrant,
1923 previous_global_first_quadrant,
1924 parallel_forest->mpicomm,
1926 this->data_serializer.dest_data_fixed.data(),
1927 this->data_serializer.src_data_fixed.data(),
1928 this->data_serializer.sizes_fixed_cumulative.back());
1930 if (this->data_serializer.variable_size_data_stored)
1933 this->data_serializer.dest_sizes_variable.resize(
1934 parallel_forest->local_num_quadrants);
1939 parallel_forest->global_first_quadrant,
1940 previous_global_first_quadrant,
1941 parallel_forest->mpicomm,
1943 this->data_serializer.dest_sizes_variable.data(),
1944 this->data_serializer.src_sizes_variable.data(),
1945 sizeof(
unsigned int));
1951 this->data_serializer.src_data_fixed.clear();
1952 this->data_serializer.src_data_fixed.shrink_to_fit();
1954 if (this->data_serializer.variable_size_data_stored)
1957 this->data_serializer.dest_data_variable.resize(
1958 std::accumulate(this->data_serializer.dest_sizes_variable.begin(),
1959 this->data_serializer.dest_sizes_variable.end(),
1960 std::vector<int>::size_type(0)));
1962# if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1969 if (this->data_serializer.src_sizes_variable.empty())
1970 this->data_serializer.src_sizes_variable.resize(1);
1971 if (this->data_serializer.dest_sizes_variable.empty())
1972 this->data_serializer.dest_sizes_variable.resize(1);
1977 parallel_forest->global_first_quadrant,
1978 previous_global_first_quadrant,
1979 parallel_forest->mpicomm,
1981 this->data_serializer.dest_data_variable.data(),
1982 this->data_serializer.dest_sizes_variable.data(),
1983 this->data_serializer.src_data_variable.data(),
1984 this->data_serializer.src_sizes_variable.data());
1987 this->data_serializer.src_sizes_variable.clear();
1988 this->data_serializer.src_sizes_variable.shrink_to_fit();
1989 this->data_serializer.src_data_variable.clear();
1990 this->data_serializer.src_data_variable.shrink_to_fit();
1996 template <
int dim,
int spacedim>
1999 spacedim>::setup_coarse_cell_to_p4est_tree_permutation()
2004 coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
2006 cell_connectivity, coarse_cell_to_p4est_tree_permutation);
2008 p4est_tree_to_coarse_cell_permutation =
2014 template <
int dim,
int spacedim>
2017 const
std::
string &file_basename)
const
2019 Assert(parallel_forest !=
nullptr,
2020 ExcMessage(
"Can't produce output when no forest is created yet."));
2024 "To use this function the triangulation's flag "
2025 "Settings::communicate_vertices_to_p4est must be set."));
2028 parallel_forest,
nullptr, file_basename.c_str());
2033 template <
int dim,
int spacedim>
2038 this->cell_attached_data.n_attached_deserialize == 0,
2040 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2041 Assert(this->n_cells() > 0,
2042 ExcMessage(
"Can not save() an empty Triangulation."));
2048 this->signals.pre_distributed_save();
2050 if (this->my_subdomain == 0)
2052 std::string fname = std::string(filename) +
".info";
2053 std::ofstream f(fname);
2054 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2058 << this->cell_attached_data.pack_callbacks_fixed.size() <<
" "
2059 << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
2060 << this->n_cells(0) << std::endl;
2064 for (
const auto &cell_rel : this->local_cell_relations)
2067 Assert((cell_rel.second ==
2073 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2074 parallel_forest->global_num_quadrants,
2082 this->signals.post_distributed_save();
2087 template <
int dim,
int spacedim>
2092 this->n_cells() > 0,
2094 "load() only works if the Triangulation already contains a coarse mesh!"));
2096 this->n_levels() == 1,
2098 "Triangulation may only contain coarse cells when calling load()."));
2104 this->signals.pre_distributed_load();
2106 if (parallel_ghost !=
nullptr)
2110 parallel_ghost =
nullptr;
2113 parallel_forest =
nullptr;
2116 connectivity =
nullptr;
2118 unsigned int version, numcpus, attached_count_fixed,
2119 attached_count_variable, n_coarse_cells;
2121 std::string fname = std::string(filename) +
".info";
2122 std::ifstream f(fname);
2124 std::string firstline;
2125 getline(f, firstline);
2126 f >> version >> numcpus >> attached_count_fixed >>
2127 attached_count_variable >> n_coarse_cells;
2131 ExcMessage(
"Incompatible version found in .info file."));
2132 Assert(this->n_cells(0) == n_coarse_cells,
2133 ExcMessage(
"Number of coarse cells differ!"));
2137 this->cell_attached_data.n_attached_data_sets = 0;
2138 this->cell_attached_data.n_attached_deserialize =
2139 attached_count_fixed + attached_count_variable;
2143 this->mpi_communicator,
2161 copy_local_forest_to_triangulation();
2171 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2172 parallel_forest->global_num_quadrants,
2173 parallel_forest->local_num_quadrants,
2175 attached_count_fixed,
2176 attached_count_variable);
2179 this->signals.post_distributed_load();
2181 this->update_periodic_face_map();
2182 this->update_number_cache();
2187 template <
int dim,
int spacedim>
2190 const
bool autopartition)
2192 (void)autopartition;
2198 template <
int dim,
int spacedim>
2201 const typename ::
internal::p4est::
types<dim>::forest *forest)
2203 Assert(this->n_cells() > 0,
2205 "load() only works if the Triangulation already contains "
2207 Assert(this->n_cells() == forest->trees->elem_count,
2209 "Coarse mesh of the Triangulation does not match the one "
2210 "of the provided forest!"));
2213 if (parallel_ghost !=
nullptr)
2217 parallel_ghost =
nullptr;
2220 parallel_forest =
nullptr;
2226 typename ::internal::p4est::types<dim>::forest *temp =
2227 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2231 parallel_forest->connectivity = connectivity;
2232 parallel_forest->user_pointer =
this;
2236 copy_local_forest_to_triangulation();
2245 this->update_periodic_face_map();
2246 this->update_number_cache();
2251 template <
int dim,
int spacedim>
2255 Assert(parallel_forest !=
nullptr,
2257 "Can't produce a check sum when no forest is created yet."));
2258 return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2263 template <
int dim,
int spacedim>
2265 const typename ::internal::p4est::types<dim>::forest
2268 Assert(parallel_forest !=
nullptr,
2269 ExcMessage(
"The forest has not been allocated yet."));
2270 return parallel_forest;
2275 template <
int dim,
int spacedim>
2277 typename ::internal::p4est::types<dim>::tree
2279 const int dealii_coarse_cell_index)
const
2281 const unsigned int tree_index =
2282 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2283 typename ::internal::p4est::types<dim>::tree *tree =
2284 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2285 sc_array_index(parallel_forest->trees, tree_index));
2299 std::integral_constant<int, 2>)
2301 const unsigned int dim = 2, spacedim = 2;
2309 std::vector<unsigned int> vertex_touch_count;
2311 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2314 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2315 const ::internal::p4est::types<2>::locidx num_vtt =
2316 std::accumulate(vertex_touch_count.begin(),
2317 vertex_touch_count.end(),
2323 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2326 (set_vertex_info ==
true ? this->n_vertices() : 0),
2331 set_vertex_and_cell_info(*
this,
2334 coarse_cell_to_p4est_tree_permutation,
2338 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2343 this->mpi_communicator,
2360 std::integral_constant<int, 2>)
2362 const unsigned int dim = 2, spacedim = 3;
2370 std::vector<unsigned int> vertex_touch_count;
2372 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2375 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2376 const ::internal::p4est::types<2>::locidx num_vtt =
2377 std::accumulate(vertex_touch_count.begin(),
2378 vertex_touch_count.end(),
2384 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2387 (set_vertex_info ==
true ? this->n_vertices() : 0),
2392 set_vertex_and_cell_info(*
this,
2395 coarse_cell_to_p4est_tree_permutation,
2399 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2404 this->mpi_communicator,
2419 std::integral_constant<int, 3>)
2421 const int dim = 3, spacedim = 3;
2429 std::vector<unsigned int> vertex_touch_count;
2430 std::vector<std::list<
2431 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2433 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2434 const ::internal::p4est::types<2>::locidx num_vtt =
2435 std::accumulate(vertex_touch_count.begin(),
2436 vertex_touch_count.end(),
2439 std::vector<unsigned int> edge_touch_count;
2440 std::vector<std::list<
2441 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2443 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2444 const ::internal::p4est::types<2>::locidx num_ett =
2445 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2448 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2451 (set_vertex_info ==
true ? this->n_vertices() : 0),
2453 this->n_active_lines(),
2458 set_vertex_and_cell_info(*
this,
2461 coarse_cell_to_p4est_tree_permutation,
2490 const unsigned int deal_to_p4est_line_index[12] = {
2491 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2494 this->begin_active();
2495 cell != this->
end();
2498 const unsigned int index =
2499 coarse_cell_to_p4est_tree_permutation[cell->index()];
2500 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2502 deal_to_p4est_line_index[e]] =
2503 cell->line(e)->index();
2508 connectivity->ett_offset[0] = 0;
2509 std::partial_sum(edge_touch_count.begin(),
2510 edge_touch_count.end(),
2511 &connectivity->ett_offset[1]);
2513 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2516 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2518 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2522 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2523 unsigned int>>::const_iterator
p =
2524 edge_to_cell[v].begin();
2525 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++
p)
2527 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2528 coarse_cell_to_p4est_tree_permutation[
p->first->index()];
2529 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2530 deal_to_p4est_line_index[
p->second];
2534 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2539 this->mpi_communicator,
2556 template <
int dim,
int spacedim>
2558 enforce_mesh_balance_over_periodic_boundaries(
2564 std::vector<bool> flags_before[2];
2568 std::vector<unsigned int> topological_vertex_numbering(
2570 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2571 topological_vertex_numbering[i] = i;
2589 using cell_iterator =
2591 typename std::map<std::pair<cell_iterator, unsigned int>,
2592 std::pair<std::pair<cell_iterator, unsigned int>,
2593 unsigned char>>::const_iterator
it;
2598 const cell_iterator &cell_1 =
it->first.first;
2599 const unsigned int face_no_1 =
it->first.second;
2600 const cell_iterator &cell_2 =
it->second.first.first;
2601 const unsigned int face_no_2 =
it->second.first.second;
2602 const unsigned char combined_orientation =
it->second.second;
2603 const auto [orientation, rotation, flip] =
2606 if (cell_1->level() == cell_2->level())
2608 for (
unsigned int v = 0;
2614 const unsigned int vface0 =
2616 v, orientation, flip, rotation);
2617 const unsigned int vi0 =
2618 topological_vertex_numbering[cell_1->face(face_no_1)
2619 ->vertex_index(vface0)];
2620 const unsigned int vi1 =
2621 topological_vertex_numbering[cell_2->face(face_no_2)
2623 const unsigned int min_index =
std::min(vi0, vi1);
2624 topological_vertex_numbering[cell_1->face(face_no_1)
2625 ->vertex_index(vface0)] =
2626 topological_vertex_numbering[cell_2->face(face_no_2)
2627 ->vertex_index(v)] =
2634 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2636 const unsigned int j = topological_vertex_numbering[i];
2644 bool continue_iterating =
true;
2646 while (continue_iterating)
2650 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2654 for (; cell != endc; ++cell)
2656 if (cell->refine_flag_set())
2657 for (
const unsigned int vertex :
2659 vertex_level[topological_vertex_numbering
2660 [cell->vertex_index(vertex)]] =
2661 std::
max(vertex_level[topological_vertex_numbering
2662 [cell->vertex_index(vertex)]],
2664 else if (!cell->coarsen_flag_set())
2665 for (
const unsigned int vertex :
2667 vertex_level[topological_vertex_numbering
2668 [cell->vertex_index(vertex)]] =
2669 std::
max(vertex_level[topological_vertex_numbering
2670 [cell->vertex_index(vertex)]],
2681 for (
const unsigned int vertex :
2683 vertex_level[topological_vertex_numbering
2684 [cell->vertex_index(vertex)]] =
2685 std::
max(vertex_level[topological_vertex_numbering
2686 [cell->vertex_index(vertex)]],
2691 continue_iterating =
false;
2703 if (cell->refine_flag_set() ==
false)
2705 for (
const unsigned int vertex :
2707 if (vertex_level[topological_vertex_numbering
2708 [cell->vertex_index(vertex)]] >=
2712 cell->clear_coarsen_flag();
2717 if (vertex_level[topological_vertex_numbering
2718 [cell->vertex_index(vertex)]] >
2721 cell->set_refine_flag();
2722 continue_iterating =
true;
2724 for (
const unsigned int v :
2726 vertex_level[topological_vertex_numbering
2727 [cell->vertex_index(v)]] =
2729 vertex_level[topological_vertex_numbering
2730 [cell->vertex_index(v)]],
2741 for (
const auto &cell :
tria.cell_iterators())
2744 if (cell->is_active())
2747 const unsigned int n_children = cell->n_children();
2748 unsigned int flagged_children = 0;
2749 for (
unsigned int child = 0; child < n_children; ++child)
2750 if (cell->child(child)->is_active() &&
2751 cell->child(child)->coarsen_flag_set())
2756 if (flagged_children < n_children)
2757 for (
unsigned int child = 0; child < n_children; ++child)
2758 if (cell->child(child)->is_active())
2759 cell->child(child)->clear_coarsen_flag();
2762 std::vector<bool> flags_after[2];
2765 return ((flags_before[0] != flags_after[0]) ||
2766 (flags_before[1] != flags_after[1]));
2772 template <
int dim,
int spacedim>
2779 exchange_refinement_flags(*
this);
2790 template <
int dim,
int spacedim>
2811 if (settings & construct_multigrid_hierarchy)
2814 spacedim>::limit_level_difference_at_vertices;
2818 bool mesh_changed =
false;
2826 if (settings & mesh_reconstruction_after_repartitioning)
2827 while (this->n_levels() > 1)
2834 for (
const auto &cell :
2835 this->active_cell_iterators_on_level(this->n_levels() - 1))
2837 cell->set_coarsen_flag();
2855 if (parallel_ghost !=
nullptr)
2859 parallel_ghost =
nullptr;
2863 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2864 P4EST_CONNECT_CORNER) :
2865 typename ::internal::p4est::types<dim>::balance_type(
2866 P8EST_CONNECT_CORNER)));
2873 for (
const auto &cell : this->cell_iterators_on_level(0))
2878 for (
const auto &cell : this->cell_iterators_on_level(0))
2883 if (tree_exists_locally<dim, spacedim>(
2885 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2888 delete_all_children<dim, spacedim>(cell);
2889 if (cell->is_active())
2898 typename ::internal::p4est::types<dim>::quadrant
2900 typename ::internal::p4est::types<dim>::tree *tree =
2901 init_tree(cell->index());
2903 ::internal::p4est::init_coarse_quadrant<dim>(
2906 match_tree_recursively<dim, spacedim>(*tree,
2910 this->my_subdomain);
2917 typename ::internal::p4est::types<dim>::quadrant *quadr;
2919 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2921 for (
unsigned int g_idx = 0;
2922 g_idx < parallel_ghost->ghosts.elem_count;
2925 while (g_idx >=
static_cast<unsigned int>(
2926 parallel_ghost->proc_offsets[ghost_owner + 1]))
2928 while (g_idx >=
static_cast<unsigned int>(
2929 parallel_ghost->tree_offsets[ghost_tree + 1]))
2932 quadr =
static_cast<
2933 typename ::internal::p4est::types<dim>::quadrant *
>(
2934 sc_array_index(¶llel_ghost->ghosts, g_idx));
2936 unsigned int coarse_cell_index =
2937 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2939 match_quadrant<dim, spacedim>(
this,
2954 bool mesh_changed =
true;
2955 unsigned int loop_counter = 0;
2962 this->update_periodic_face_map();
2965 enforce_mesh_balance_over_periodic_boundaries(*
this);
2977 "parallel::distributed::Triangulation::copy_local_forest_to_triangulation() "
2978 "for periodic boundaries detected. Aborting."));
2980 while (mesh_changed);
2985 std::any_of(this->begin_active(),
2988 return cell.refine_flag_set() ||
2989 cell.coarsen_flag_set();
3007 while (mesh_changed);
3011 unsigned int num_ghosts = 0;
3013 for (
const auto &cell : this->active_cell_iterators())
3015 if (cell->subdomain_id() != this->my_subdomain &&
3020 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
3031 if (settings & construct_multigrid_hierarchy)
3037 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3040 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3042 if ((cell->is_active() &&
3043 cell->subdomain_id() ==
3044 this->locally_owned_subdomain()) ||
3045 (cell->has_children() &&
3046 cell->child(0)->level_subdomain_id() ==
3047 this->locally_owned_subdomain()))
3048 cell->set_level_subdomain_id(
3049 this->locally_owned_subdomain());
3053 cell->set_level_subdomain_id(
3061 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3062 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3063 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3065 for (
const auto &cell : this->cell_iterators_on_level(0))
3067 typename ::internal::p4est::types<dim>::quadrant
3069 const unsigned int tree_index =
3070 coarse_cell_to_p4est_tree_permutation[cell->index()];
3071 typename ::internal::p4est::types<dim>::tree *tree =
3072 init_tree(cell->index());
3074 ::internal::p4est::init_coarse_quadrant<dim>(
3077 determine_level_subdomain_id_recursively<dim, spacedim>(
3088 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3091 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3093 if (cell->has_children())
3094 for (
unsigned int c = 0;
3095 c < GeometryInfo<dim>::max_children_per_cell;
3098 if (cell->child(c)->level_subdomain_id() ==
3099 this->locally_owned_subdomain())
3104 cell->child(0)->level_subdomain_id();
3108 cell->set_level_subdomain_id(mark);
3124 const unsigned int total_local_cells = this->n_active_cells();
3125 (void)total_local_cells;
3129 Assert(
static_cast<unsigned int>(
3130 parallel_forest->local_num_quadrants) == total_local_cells,
3135 Assert(
static_cast<unsigned int>(
3136 parallel_forest->local_num_quadrants) <= total_local_cells,
3142 unsigned int n_owned = 0;
3143 for (
const auto &cell : this->active_cell_iterators())
3145 if (cell->subdomain_id() == this->my_subdomain)
3149 Assert(
static_cast<unsigned int>(
3150 parallel_forest->local_num_quadrants) == n_owned,
3155 this->smooth_grid = save_smooth;
3161 update_cell_relations();
3166 template <
int dim,
int spacedim>
3172 std::vector<Point<dim>> point{p};
3173 std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
3180 template <
int dim,
int spacedim>
3183 find_point_owner_rank(const
std::vector<
Point<dim>> &points)
3185# ifndef P4EST_SEARCH_LOCAL
3190 "This function is only available if p4est is version 2.2 and higher."));
3192 return std::vector<unsigned int>(1,
3196 AssertThrow(this->are_vertices_communicated_to_p4est(),
3198 "Vertices need to be communicated to p4est to use this "
3199 "function. This must explicitly be turned on in the "
3200 "settings of the triangulation's constructor."));
3203 for (
const auto &manifold_id : this->get_manifold_ids())
3208 "This function can only be used if the triangulation "
3209 "has no other manifold than a Cartesian (flat) manifold attached."));
3213 PartitionSearch<dim> partition_search;
3219 parallel_forest->user_pointer = &partition_search;
3225 sc_array_t *point_sc_array;
3229 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3232 for (
size_t i = 0; i < points.size(); ++i)
3237 double *this_sc_point =
3238 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3240 for (
unsigned int d = 0; d < dim; ++d)
3242 this_sc_point[d] = p(d);
3244 this_sc_point[dim] = -1.0;
3250 static_cast<int>(
false),
3251 &PartitionSearch<dim>::local_quadrant_fn,
3252 &PartitionSearch<dim>::local_point_fn,
3256 std::vector<types::subdomain_id> owner_rank(
3260 for (
size_t i = 0; i < points.size(); ++i)
3263 double *this_sc_point =
3264 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3265 Assert(this_sc_point[dim] >= 0. || this_sc_point[dim] == -1.,
3267 if (this_sc_point[dim] < 0.)
3275 parallel_forest->user_pointer =
this;
3278 sc_array_destroy_null(&point_sc_array);
3286 template <
int dim,
int spacedim>
3292 for (
const auto &cell : this->active_cell_iterators())
3293 if (cell->is_locally_owned() && cell->refine_flag_set())
3294 Assert(cell->refine_flag_set() ==
3297 "This class does not support anisotropic refinement"));
3302 if (this->n_levels() ==
3306 cell = this->begin_active(
3314 !(cell->refine_flag_set()),
3316 "Fatal Error: maximum refinement level of p4est reached."));
3320 this->prepare_coarsening_and_refinement();
3323 this->signals.pre_distributed_refinement();
3328 for (
const auto &cell : this->active_cell_iterators())
3329 if (cell->is_ghost() || cell->is_artificial())
3331 cell->clear_refine_flag();
3332 cell->clear_coarsen_flag();
3338 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3339 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3346 parallel_forest->user_pointer = &refine_and_coarsen_list;
3348 if (parallel_ghost !=
nullptr)
3352 parallel_ghost =
nullptr;
3357 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3362 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3369 parallel_forest->user_pointer =
this;
3375 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3376 P4EST_CONNECT_FULL) :
3377 typename ::internal::p4est::types<dim>::balance_type(
3378 P8EST_CONNECT_FULL)),
3383 update_cell_relations();
3387 this->signals.post_p4est_refinement();
3391 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3392 previous_global_first_quadrant;
3394 if (this->cell_attached_data.n_attached_data_sets > 0)
3396 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3397 std::memcpy(previous_global_first_quadrant.data(),
3398 parallel_forest->global_first_quadrant,
3400 typename ::internal::p4est::types<dim>::gloidx) *
3401 (parallel_forest->mpisize + 1));
3404 if (!(settings & no_automatic_repartitioning))
3408 if (this->signals.weight.empty())
3416 const std::vector<unsigned int> cell_weights = get_cell_weights();
3422 this->mpi_communicator) > 0,
3424 "The global sum of weights over all active cells "
3425 "is zero. Please verify how you generate weights."));
3427 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3432 parallel_forest->user_pointer = &partition_weights;
3438 &PartitionWeights<dim, spacedim>::cell_weight);
3442 parallel_forest, 0,
nullptr,
nullptr);
3444 parallel_forest->user_pointer =
this;
3449 if (this->cell_attached_data.n_attached_data_sets > 0)
3451 this->data_serializer.pack_data(
3452 this->local_cell_relations,
3453 this->cell_attached_data.pack_callbacks_fixed,
3454 this->cell_attached_data.pack_callbacks_variable,
3455 this->get_communicator());
3461 for (
const auto &cell : this->active_cell_iterators())
3463 cell->clear_refine_flag();
3464 cell->clear_coarsen_flag();
3469 copy_local_forest_to_triangulation();
3479 if (this->cell_attached_data.n_attached_data_sets > 0)
3481 this->execute_transfer(parallel_forest,
3482 previous_global_first_quadrant.data());
3485 this->data_serializer.unpack_cell_status(this->local_cell_relations);
3504 if (settings & construct_multigrid_hierarchy)
3506 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3508 std::vector<bool> active_verts =
3509 this->mark_locally_active_vertices_on_level(lvl);
3511 const unsigned int maybe_coarser_lvl =
3512 (lvl > 0) ? (lvl - 1) : lvl;
3514 cell = this->
begin(maybe_coarser_lvl),
3515 endc = this->
end(lvl);
3516 for (; cell != endc; ++cell)
3517 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3519 const bool is_level_artificial =
3520 (cell->level_subdomain_id() ==
3522 bool need_to_know =
false;
3523 for (
const unsigned int vertex :
3525 if (active_verts[cell->vertex_index(vertex)])
3527 need_to_know =
true;
3532 !need_to_know || !is_level_artificial,
3534 "Internal error: the owner of cell" +
3535 cell->id().to_string() +
3536 " is unknown even though it is needed for geometric multigrid."));
3542 this->update_periodic_face_map();
3543 this->update_number_cache();
3546 this->signals.post_distributed_refinement();
3551 template <
int dim,
int spacedim>
3556 for (
const auto &cell : this->active_cell_iterators())
3557 if (cell->is_locally_owned())
3559 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3561 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3565 this->signals.pre_distributed_repartition();
3569 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3570 previous_global_first_quadrant;
3572 if (this->cell_attached_data.n_attached_data_sets > 0)
3574 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3575 std::memcpy(previous_global_first_quadrant.data(),
3576 parallel_forest->global_first_quadrant,
3578 typename ::internal::p4est::types<dim>::gloidx) *
3579 (parallel_forest->mpisize + 1));
3582 if (this->signals.weight.empty())
3594 const std::vector<unsigned int> cell_weights = get_cell_weights();
3600 this->mpi_communicator) > 0,
3602 "The global sum of weights over all active cells "
3603 "is zero. Please verify how you generate weights."));
3605 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3610 parallel_forest->user_pointer = &partition_weights;
3616 &PartitionWeights<dim, spacedim>::cell_weight);
3619 parallel_forest->user_pointer =
this;
3623 if (this->cell_attached_data.n_attached_data_sets > 0)
3625 this->data_serializer.pack_data(
3626 this->local_cell_relations,
3627 this->cell_attached_data.pack_callbacks_fixed,
3628 this->cell_attached_data.pack_callbacks_variable,
3629 this->get_communicator());
3634 copy_local_forest_to_triangulation();
3644 if (this->cell_attached_data.n_attached_data_sets > 0)
3646 this->execute_transfer(parallel_forest,
3647 previous_global_first_quadrant.data());
3650 this->update_periodic_face_map();
3653 this->update_number_cache();
3656 this->signals.post_distributed_repartition();
3661 template <
int dim,
int spacedim>
3663 const std::vector<types::global_dof_index>
3667 return p4est_tree_to_coarse_cell_permutation;
3672 template <
int dim,
int spacedim>
3674 const std::vector<types::global_dof_index>
3678 return coarse_cell_to_p4est_tree_permutation;
3683 template <
int dim,
int spacedim>
3686 mark_locally_active_vertices_on_level(const
int level)
const
3690 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3691 for (
const auto &cell : this->cell_iterators_on_level(
level))
3692 if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3694 marked_vertices[cell->vertex_index(v)] =
true;
3711 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3712 for (
const auto &it : this->get_periodic_face_map())
3715 const unsigned int face_no_1 = it.first.second;
3717 const unsigned int face_no_2 = it.second.first.second;
3718 const unsigned char combined_orientation = it.second.second;
3719 const auto [orientation, rotation, flip] =
3722 if (cell_1->level() ==
level && cell_2->level() ==
level)
3724 for (
unsigned int v = 0;
3730 const unsigned int vface0 =
3732 v, orientation, flip, rotation);
3733 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3735 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3737 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3739 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3745 return marked_vertices;
3750 template <
int dim,
int spacedim>
3753 coarse_cell_id_to_coarse_cell_index(
3754 const
types::coarse_cell_id coarse_cell_id)
const
3756 return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3761 template <
int dim,
int spacedim>
3765 const unsigned int coarse_cell_index)
const
3767 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3772 template <
int dim,
int spacedim>
3776 &periodicity_vector)
3778 Assert(triangulation_has_content ==
true,
3780 Assert(this->n_levels() == 1,
3781 ExcMessage(
"The triangulation is refined!"));
3788 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
3790 for (
const auto &face_pair : periodicity_vector)
3794 const unsigned int face_left = face_pair.face_idx[0];
3795 const unsigned int face_right = face_pair.face_idx[1];
3798 const unsigned int tree_left =
3799 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3800 const unsigned int tree_right =
3801 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3809 unsigned int p4est_orientation = 0;
3814 face_pair.orientation ==
3821 const unsigned int face_idx_list[] = {face_left, face_right};
3822 const cell_iterator cell_list[] = {first_cell, second_cell};
3823 unsigned int lower_idx, higher_idx;
3824 unsigned char orientation;
3825 if (face_left <= face_right)
3830 face_reference_cell.get_inverse_combined_orientation(
3831 face_pair.orientation);
3837 orientation = face_pair.orientation;
3842 unsigned int first_p4est_idx_on_cell =
3843 p8est_face_corners[face_idx_list[lower_idx]][0];
3844 unsigned int first_dealii_idx_on_face =
3846 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3849 const unsigned int first_dealii_idx_on_cell =
3851 face_idx_list[lower_idx],
3853 cell_list[lower_idx]->face_orientation(
3854 face_idx_list[lower_idx]),
3855 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3856 cell_list[lower_idx]->face_rotation(
3857 face_idx_list[lower_idx]));
3858 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3860 first_dealii_idx_on_face = i;
3868 const unsigned int second_dealii_idx_on_face =
3869 reference_cell.standard_to_real_face_vertex(
3870 first_dealii_idx_on_face,
3871 face_idx_list[lower_idx],
3873 const unsigned int second_dealii_idx_on_cell =
3874 reference_cell.face_to_cell_vertices(
3875 face_idx_list[higher_idx],
3876 second_dealii_idx_on_face,
3877 cell_list[higher_idx]->combined_face_orientation(
3878 face_idx_list[higher_idx]));
3880 const unsigned int second_p4est_idx_on_face =
3881 p8est_corner_face_corners[second_dealii_idx_on_cell]
3882 [face_idx_list[higher_idx]];
3883 p4est_orientation = second_p4est_idx_on_face;
3903 this->mpi_communicator,
3914 copy_local_forest_to_triangulation();
3925 this->update_number_cache();
3930 template <
int dim,
int spacedim>
3941 this->cell_attached_data.n_attached_data_sets) +
3948 coarse_cell_to_p4est_tree_permutation) +
3950 p4est_tree_to_coarse_cell_permutation) +
3951 memory_consumption_p4est();
3958 template <
int dim,
int spacedim>
3962 return ::internal::p4est::functions<dim>::forest_memory_used(
3970 template <
int dim,
int spacedim>
3977 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
3979 (other_tria.n_global_levels() == 1),
3990 const typename ::Triangulation<dim, spacedim>::DistortedCellList
3998 if (const ::parallel::distributed::Triangulation<dim, spacedim>
3999 *other_distributed =
4000 dynamic_cast<const ::parallel::distributed::
4001 Triangulation<dim, spacedim> *
>(&other_tria))
4004 settings = other_distributed->settings;
4005 triangulation_has_content =
4006 other_distributed->triangulation_has_content;
4007 coarse_cell_to_p4est_tree_permutation =
4008 other_distributed->coarse_cell_to_p4est_tree_permutation;
4009 p4est_tree_to_coarse_cell_permutation =
4010 other_distributed->p4est_tree_to_coarse_cell_permutation;
4013 typename ::internal::p4est::types<dim>::connectivity
4014 *temp_connectivity =
const_cast<
4015 typename ::internal::p4est::types<dim>::connectivity *
>(
4016 other_distributed->connectivity);
4018 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
4021 typename ::internal::p4est::types<dim>::forest *temp_forest =
4022 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
4023 other_distributed->parallel_forest);
4027 parallel_forest->connectivity = connectivity;
4028 parallel_forest->user_pointer =
this;
4032 triangulation_has_content =
true;
4033 setup_coarse_cell_to_p4est_tree_permutation();
4034 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
4039 copy_local_forest_to_triangulation();
4048 this->update_periodic_face_map();
4049 this->update_number_cache();
4054 template <
int dim,
int spacedim>
4059 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
4060 this->local_cell_relations.shrink_to_fit();
4063 for (
const auto &cell : this->cell_iterators_on_level(0))
4066 if (tree_exists_locally<dim, spacedim>(
4068 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4072 typename ::internal::p4est::types<dim>::quadrant
4074 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4077 typename ::internal::p4est::types<dim>::tree *tree =
4078 init_tree(cell->index());
4080 update_cell_relations_recursively<dim, spacedim>(
4081 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4087 template <
int dim,
int spacedim>
4094 Assert(this->local_cell_relations.size() ==
4095 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4100 std::vector<unsigned int> weights;
4101 weights.reserve(this->local_cell_relations.size());
4109 for (
const auto &cell_rel : this->local_cell_relations)
4111 const auto &cell_it = cell_rel.first;
4112 const auto &cell_status = cell_rel.second;
4114 weights.push_back(this->signals.weight(cell_it, cell_status));
4122 template <
int spacedim>
4138 template <
int spacedim>
4147 template <
int spacedim>
4149 const std::vector<types::global_dof_index>
4153 static std::vector<types::global_dof_index> a;
4159 template <
int spacedim>
4161 std::map<
unsigned int,
4164 const unsigned int )
const
4168 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4173 template <
int spacedim>
4176 mark_locally_active_vertices_on_level(const
unsigned int)
const
4179 return std::vector<bool>();
4184 template <
int spacedim>
4187 coarse_cell_id_to_coarse_cell_index(const
types::coarse_cell_id)
const
4195 template <
int spacedim>
4199 const unsigned int)
const
4207 template <
int spacedim>
4216 template <
int spacedim>
4225 template <
int spacedim>
4234 template <
int spacedim>
4244 template <
int spacedim>
4254 template <
int spacedim>
4271 namespace distributed
4273 template <
int dim,
int spacedim>
4281#ifdef DEAL_II_WITH_P4EST
4291 const auto &cell = pair.first;
4292 const auto &status = pair.second;
4298 cell->clear_refine_flag();
4299 cell->clear_coarsen_flag();
4304 cell->clear_coarsen_flag();
4305 cell->set_refine_flag();
4310 for (
const auto &child : cell->child_iterators())
4312 child->clear_refine_flag();
4313 child->set_coarsen_flag();
4332 template <
int dim,
int spacedim>
4335#ifdef DEAL_II_WITH_P4EST
4336 if (distributed_tria)
4339 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4340 distributed_tria->load_refine_flags(saved_refine_flags);
4344 (void)distributed_tria;
@ children_will_be_coarsened
static constexpr unsigned char default_combined_face_orientation()
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, unsigned char > > & get_periodic_face_map() const
active_cell_iterator last_active() const
cell_iterator end() const
void save_refine_flags(std::ostream &out) const
unsigned int n_vertices() const
void save_coarsen_flags(std::ostream &out) const
active_cell_iterator begin_active(const unsigned int level=0) const
virtual void clear() override
virtual std::size_t memory_consumption() const override
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
TemporarilyMatchRefineFlags(::Triangulation< dim, spacedim > &tria)
~TemporarilyMatchRefineFlags()
const SmartPointer< ::parallel::distributed::Triangulation< dim, spacedim > > distributed_tria
std::vector< bool > saved_refine_flags
std::vector< bool > saved_coarsen_flags
virtual void execute_coarsening_and_refinement() override
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator > > &) override
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const override
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
typename::internal::p4est::types< dim >::tree * init_tree(const int dealii_coarse_cell_index) const
types::subdomain_id find_point_owner_rank(const Point< dim > &p)
virtual bool prepare_coarsening_and_refinement() override
const ::internal::p4est::types< dim >::forest * get_p4est() const
const std::vector< types::global_dof_index > & get_coarse_cell_to_p4est_tree_permutation() const
void copy_new_triangulation_to_p4est(std::integral_constant< int, 2 >)
virtual void clear() override
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata) 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()
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)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
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)
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
const Iterator const std_cxx20::type_identity_t< Iterator > & end
::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)