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,
45 namespace TriangulationImplementation {
…}
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 num_vtt = std::accumulate(vertex_touch_count.begin(),
318 vertex_touch_count.end(),
323 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
325 Assert(vertex_to_cell[v].
size() == vertex_touch_count[v],
329 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
330 unsigned int>>::const_iterator p =
331 vertex_to_cell[v].begin();
332 for (
unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
334 connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
335 coarse_cell_to_p4est_tree_permutation[p->first->index()];
336 connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
344 template <
int dim,
int spacedim>
350 Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
352 return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
353 (coarse_grid_cell <= parallel_forest->last_local_tree));
357 template <
int dim,
int spacedim>
359 delete_all_children_and_self(
362 if (cell->has_children())
363 for (
unsigned int c = 0; c < cell->n_children(); ++c)
364 delete_all_children_and_self<dim, spacedim>(cell->child(c));
366 cell->set_coarsen_flag();
371 template <
int dim,
int spacedim>
376 if (cell->has_children())
377 for (
unsigned int c = 0; c < cell->n_children(); ++c)
378 delete_all_children_and_self<dim, spacedim>(cell->child(c));
382 template <
int dim,
int spacedim>
384 determine_level_subdomain_id_recursively(
391 const std::vector<std::vector<bool>> &marked_vertices)
400 if (marked_vertices[dealii_cell->level()]
401 [dealii_cell->vertex_index(v)])
420 if (!used && dealii_cell->is_active() &&
421 dealii_cell->is_artificial() ==
false &&
422 dealii_cell->level() + 1 <
static_cast<int>(marked_vertices.size()))
426 if (marked_vertices[dealii_cell->level() + 1]
427 [dealii_cell->vertex_index(v)])
436 if (!used && dealii_cell->is_active() &&
437 dealii_cell->is_artificial() ==
false && dealii_cell->level() > 0)
441 if (marked_vertices[dealii_cell->level() - 1]
442 [dealii_cell->vertex_index(v)])
453 &forest, tree_index, &p4est_cell, my_subdomain);
454 Assert((owner != -2) && (owner != -1),
456 dealii_cell->set_level_subdomain_id(owner);
460 if (dealii_cell->has_children())
464 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
469 P4EST_QUADRANT_INIT(&p4est_child[c]);
472 P8EST_QUADRANT_INIT(&p4est_child[c]);
482 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
485 determine_level_subdomain_id_recursively<dim, spacedim>(
488 dealii_cell->child(c),
498 template <
int dim,
int spacedim>
500 match_tree_recursively(
508 if (sc_array_bsearch(
const_cast<sc_array_t *
>(&tree.quadrants),
514 delete_all_children<dim, spacedim>(dealii_cell);
515 if (dealii_cell->is_active())
516 dealii_cell->set_subdomain_id(my_subdomain);
525 if (dealii_cell->is_active())
526 dealii_cell->set_refine_flag();
531 for (
unsigned int c = 0;
532 c < GeometryInfo<dim>::max_children_per_cell;
537 P4EST_QUADRANT_INIT(&p4est_child[c]);
540 P8EST_QUADRANT_INIT(&p4est_child[c]);
550 for (
unsigned int c = 0;
551 c < GeometryInfo<dim>::max_children_per_cell;
556 &p4est_child[c]) ==
false)
562 delete_all_children<dim, spacedim>(dealii_cell->child(c));
563 dealii_cell->child(c)->recursively_set_subdomain_id(
570 match_tree_recursively<dim, spacedim>(tree,
571 dealii_cell->child(c),
581 template <
int dim,
int spacedim>
584 const ::Triangulation<dim, spacedim> *tria,
585 unsigned int dealii_index,
589 const int l = ghost_quadrant.level;
591 for (
int i = 0; i <
l; ++i)
596 if (cell->is_active())
598 cell->clear_coarsen_flag();
599 cell->set_refine_flag();
606 dealii_index = cell->child_index(child_id);
612 if (cell->has_children())
613 delete_all_children<dim, spacedim>(cell);
616 cell->clear_coarsen_flag();
617 cell->set_subdomain_id(ghost_owner);
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(
737 PartitionSearch<dim> *this_object =
738 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
742 quad_length_on_level =
744 (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
748 this_object->quadrant_data.set_cell_vertices(forest,
751 quad_length_on_level);
754 this_object->quadrant_data.initialize_mapping();
764 PartitionSearch<dim>::local_point_fn(
777 PartitionSearch<dim> *this_object =
778 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
781 double *this_point_dptr =
static_cast<double *
>(
point);
784 (dim == 2 ?
Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
785 Point<dim>(this_point_dptr[0],
787 this_point_dptr[2]));
790 const bool is_in_this_quadrant =
791 this_object->quadrant_data.is_in_this_quadrant(this_point);
795 if (!is_in_this_quadrant)
804 if (rank_begin < rank_end)
812 this_point_dptr[dim] =
static_cast<double>(rank_begin);
822 PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
825 const Point<dim> p_ref = map_real_to_unit_cell(p);
834 PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
837 Assert(is_reference_mapping_initialized,
839 "Cell vertices and mapping coefficients must be fully "
840 "initialized before transforming a point to the unit cell."));
846 for (
unsigned int alpha = 0;
847 alpha < GeometryInfo<dim>::vertices_per_cell;
853 p_out += (quadrant_mapping_matrix(alpha, 0) +
854 quadrant_mapping_matrix(alpha, 1) * p(0) +
855 quadrant_mapping_matrix(alpha, 2) * p(1) +
856 quadrant_mapping_matrix(alpha, 3) * p(0) * p(1)) *
862 for (
unsigned int alpha = 0;
863 alpha < GeometryInfo<dim>::vertices_per_cell;
869 p_out += (quadrant_mapping_matrix(alpha, 0) +
870 quadrant_mapping_matrix(alpha, 1) * p(0) +
871 quadrant_mapping_matrix(alpha, 2) * p(1) +
872 quadrant_mapping_matrix(alpha, 3) * p(2) +
873 quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
874 quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
875 quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
876 quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
886 PartitionSearch<dim>::QuadrantData::QuadrantData()
888 , quadrant_mapping_matrix(
GeometryInfo<dim>::vertices_per_cell,
890 , are_vertices_initialized(false)
891 , is_reference_mapping_initialized(false)
898 PartitionSearch<dim>::QuadrantData::initialize_mapping()
901 are_vertices_initialized,
903 "Cell vertices must be initialized before the cell mapping can be filled."));
910 for (
unsigned int alpha = 0;
911 alpha < GeometryInfo<dim>::vertices_per_cell;
915 point_matrix(0, alpha) = 1;
916 point_matrix(1, alpha) = cell_vertices[alpha](0);
917 point_matrix(2, alpha) = cell_vertices[alpha](1);
918 point_matrix(3, alpha) =
919 cell_vertices[alpha](0) * cell_vertices[alpha](1);
926 quadrant_mapping_matrix.invert(point_matrix);
930 for (
unsigned int alpha = 0;
931 alpha < GeometryInfo<dim>::vertices_per_cell;
935 point_matrix(0, alpha) = 1;
936 point_matrix(1, alpha) = cell_vertices[alpha](0);
937 point_matrix(2, alpha) = cell_vertices[alpha](1);
938 point_matrix(3, alpha) = cell_vertices[alpha](2);
939 point_matrix(4, alpha) =
940 cell_vertices[alpha](0) * cell_vertices[alpha](1);
941 point_matrix(5, alpha) =
942 cell_vertices[alpha](1) * cell_vertices[alpha](2);
943 point_matrix(6, alpha) =
944 cell_vertices[alpha](0) * cell_vertices[alpha](2);
945 point_matrix(7, alpha) = cell_vertices[alpha](0) *
946 cell_vertices[alpha](1) *
947 cell_vertices[alpha](2);
954 quadrant_mapping_matrix.invert(point_matrix);
957 is_reference_mapping_initialized =
true;
964 PartitionSearch<2>::QuadrantData::set_cell_vertices(
969 quad_length_on_level)
971 constexpr unsigned int dim = 2;
975 double corner_point[dim + 1] = {0};
978 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
980 for (
unsigned int d = 0;
d < dim; ++
d)
982 cell_vertices[vertex_index](
d) = corner_point[d];
992 unsigned int vertex_index = 0;
994 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
997 copy_vertex(vertex_index);
1004 forest->connectivity,
1006 quadrant->x + quad_length_on_level,
1011 copy_vertex(vertex_index);
1018 forest->connectivity,
1021 quadrant->y + quad_length_on_level,
1025 copy_vertex(vertex_index);
1032 forest->connectivity,
1034 quadrant->x + quad_length_on_level,
1035 quadrant->y + quad_length_on_level,
1039 copy_vertex(vertex_index);
1041 are_vertices_initialized =
true;
1048 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1053 quad_length_on_level)
1055 constexpr unsigned int dim = 3;
1057 double corner_point[dim] = {0};
1060 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1062 for (
unsigned int d = 0;
d < dim; ++
d)
1064 cell_vertices[vertex_index](
d) = corner_point[d];
1066 corner_point[
d] = 0;
1074 unsigned int vertex_index = 0;
1076 forest->connectivity,
1084 copy_vertex(vertex_index);
1092 forest->connectivity,
1094 quadrant->x + quad_length_on_level,
1100 copy_vertex(vertex_index);
1107 forest->connectivity,
1110 quadrant->y + quad_length_on_level,
1115 copy_vertex(vertex_index);
1122 forest->connectivity,
1124 quadrant->x + quad_length_on_level,
1125 quadrant->y + quad_length_on_level,
1130 copy_vertex(vertex_index);
1137 forest->connectivity,
1141 quadrant->z + quad_length_on_level,
1145 copy_vertex(vertex_index);
1152 forest->connectivity,
1154 quadrant->x + quad_length_on_level,
1156 quadrant->z + quad_length_on_level,
1160 copy_vertex(vertex_index);
1167 forest->connectivity,
1170 quadrant->y + quad_length_on_level,
1171 quadrant->z + quad_length_on_level,
1175 copy_vertex(vertex_index);
1182 forest->connectivity,
1184 quadrant->x + quad_length_on_level,
1185 quadrant->y + quad_length_on_level,
1186 quadrant->z + quad_length_on_level,
1190 copy_vertex(vertex_index);
1193 are_vertices_initialized =
true;
1203 template <
int dim,
int spacedim>
1204 class RefineAndCoarsenList
1208 const std::vector<types::global_dof_index>
1209 &p4est_tree_to_coarse_cell_permutation,
1237 pointers_are_at_end()
const;
1240 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1241 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1242 const_iterator current_refine_pointer;
1244 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1245 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1246 const_iterator current_coarsen_pointer;
1257 template <
int dim,
int spacedim>
1259 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1261 return ((current_refine_pointer == refine_list.end()) &&
1262 (current_coarsen_pointer == coarsen_list.end()));
1267 template <
int dim,
int spacedim>
1268 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1270 const std::vector<types::global_dof_index>
1271 &p4est_tree_to_coarse_cell_permutation,
1275 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1276 for (
const auto &cell :
triangulation.active_cell_iterators())
1279 if (cell->subdomain_id() != my_subdomain)
1282 if (cell->refine_flag_set())
1284 else if (cell->coarsen_flag_set())
1288 refine_list.reserve(n_refine_flags);
1289 coarsen_list.reserve(n_coarsen_flags);
1301 unsigned int coarse_cell_index =
1302 p4est_tree_to_coarse_cell_permutation[c];
1311 p4est_cell.p.which_tree = c;
1312 build_lists(cell, p4est_cell, my_subdomain);
1320 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1321 Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1323 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1324 Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1327 current_refine_pointer = refine_list.begin();
1328 current_coarsen_pointer = coarsen_list.begin();
1333 template <
int dim,
int spacedim>
1335 RefineAndCoarsenList<dim, spacedim>::build_lists(
1340 if (cell->is_active())
1342 if (cell->subdomain_id() == my_subdomain)
1344 if (cell->refine_flag_set())
1345 refine_list.push_back(p4est_cell);
1346 else if (cell->coarsen_flag_set())
1347 coarsen_list.push_back(p4est_cell);
1354 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1359 P4EST_QUADRANT_INIT(&p4est_child[c]);
1362 P8EST_QUADRANT_INIT(&p4est_child[c]);
1369 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1372 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1373 build_lists(cell->child(c), p4est_child[c], my_subdomain);
1379 template <
int dim,
int spacedim>
1381 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1386 RefineAndCoarsenList<dim, spacedim> *this_object =
1387 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1388 forest->user_pointer);
1392 if (this_object->current_refine_pointer == this_object->refine_list.end())
1395 Assert(coarse_cell_index <=
1396 this_object->current_refine_pointer->p.which_tree,
1401 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1405 Assert(coarse_cell_index <=
1406 this_object->current_refine_pointer->p.which_tree,
1412 quadrant, &*this_object->current_refine_pointer) <= 0,
1417 quadrant, &*this_object->current_refine_pointer))
1419 ++this_object->current_refine_pointer;
1429 template <
int dim,
int spacedim>
1431 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1436 RefineAndCoarsenList<dim, spacedim> *this_object =
1437 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1438 forest->user_pointer);
1442 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1445 Assert(coarse_cell_index <=
1446 this_object->current_coarsen_pointer->p.which_tree,
1451 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1455 Assert(coarse_cell_index <=
1456 this_object->current_coarsen_pointer->p.which_tree,
1462 children[0], &*this_object->current_coarsen_pointer) <= 0,
1468 children[0], &*this_object->current_coarsen_pointer))
1471 ++this_object->current_coarsen_pointer;
1475 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1479 children[c], &*this_object->current_coarsen_pointer),
1481 ++this_object->current_coarsen_pointer;
1499 template <
int dim,
int spacedim>
1500 class PartitionWeights
1508 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1523 std::vector<unsigned int> cell_weights_list;
1524 std::vector<unsigned int>::const_iterator current_pointer;
1528 template <
int dim,
int spacedim>
1529 PartitionWeights<dim, spacedim>::PartitionWeights(
1530 const std::vector<unsigned int> &cell_weights)
1531 : cell_weights_list(cell_weights)
1535 current_pointer = cell_weights_list.begin();
1539 template <
int dim,
int spacedim>
1541 PartitionWeights<dim, spacedim>::cell_weight(
1550 PartitionWeights<dim, spacedim> *this_object =
1551 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1553 Assert(this_object->current_pointer >=
1554 this_object->cell_weights_list.begin(),
1556 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1562 const unsigned int weight = *this_object->current_pointer;
1563 ++this_object->current_pointer;
1565 Assert(weight <
static_cast<unsigned int>(std::numeric_limits<int>::max()),
1566 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1567 "weights for cells. The weight provided here exceeds "
1568 "the maximum value represented as a 'signed int'."));
1569 return static_cast<int>(weight);
1572 template <
int dim,
int spacedim>
1573 using cell_relation_t =
typename std::pair<
1574 typename ::Triangulation<dim, spacedim>::cell_iterator,
1586 template <
int dim,
int spacedim>
1588 add_single_cell_relation(
1589 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1590 const typename ::internal::p4est::types<dim>::tree &tree,
1591 const unsigned int idx,
1595 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1601 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1615 template <
int dim,
int spacedim>
1617 update_cell_relations_recursively(
1618 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1619 const typename ::internal::p4est::types<dim>::tree &tree,
1621 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1624 const int idx = sc_array_bsearch(
1625 const_cast<sc_array_t *
>(&tree.quadrants),
1630 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1632 &p4est_cell) ==
false))
1637 const bool p4est_has_children = (idx == -1);
1638 if (p4est_has_children && dealii_cell->has_children())
1641 typename ::internal::p4est::types<dim>::quadrant
1644 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1649 P4EST_QUADRANT_INIT(&p4est_child[c]);
1652 P8EST_QUADRANT_INIT(&p4est_child[c]);
1659 &p4est_cell, p4est_child);
1661 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1664 update_cell_relations_recursively<dim, spacedim>(
1665 cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1668 else if (!p4est_has_children && !dealii_cell->has_children())
1672 add_single_cell_relation<dim, spacedim>(
1675 else if (p4est_has_children)
1683 typename ::internal::p4est::types<dim>::quadrant
1685 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1690 P4EST_QUADRANT_INIT(&p4est_child[c]);
1693 P8EST_QUADRANT_INIT(&p4est_child[c]);
1700 &p4est_cell, p4est_child);
1708 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1711 child_idx = sc_array_bsearch(
1712 const_cast<sc_array_t *
>(&tree.quadrants),
1719 add_single_cell_relation<dim, spacedim>(
1720 cell_rel, tree, child_idx, dealii_cell, cell_status);
1728 add_single_cell_relation<dim, spacedim>(
1742 namespace distributed
1745 template <
int dim,
int spacedim>
1757 (settings & construct_multigrid_hierarchy) ?
1761 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1764 , settings(settings)
1765 , triangulation_has_content(false)
1766 , connectivity(
nullptr)
1767 , parallel_forest(
nullptr)
1769 parallel_ghost =
nullptr;
1774 template <
int dim,
int spacedim>
1796 template <
int dim,
int spacedim>
1799 const
std::vector<
Point<spacedim>> &vertices,
1806 vertices, cells, subcelldata);
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>
1860 template <
int dim,
int spacedim>
1864 triangulation_has_content =
false;
1866 if (parallel_ghost !=
nullptr)
1870 parallel_ghost =
nullptr;
1873 if (parallel_forest !=
nullptr)
1876 parallel_forest =
nullptr;
1879 if (connectivity !=
nullptr)
1883 connectivity =
nullptr;
1886 coarse_cell_to_p4est_tree_permutation.resize(0);
1887 p4est_tree_to_coarse_cell_permutation.resize(0);
1891 this->update_number_cache();
1896 template <
int dim,
int spacedim>
1907 template <
int dim,
int spacedim>
1918 template <
int dim,
int spacedim>
1924 *previous_global_first_quadrant)
1926 Assert(this->data_serializer.sizes_fixed_cumulative.size() > 0,
1930 this->data_serializer.dest_data_fixed.resize(
1931 parallel_forest->local_num_quadrants *
1932 this->data_serializer.sizes_fixed_cumulative.back());
1935 typename ::internal::p4est::types<dim>::transfer_context
1939 parallel_forest->global_first_quadrant,
1940 previous_global_first_quadrant,
1941 parallel_forest->mpicomm,
1943 this->data_serializer.dest_data_fixed.data(),
1944 this->data_serializer.src_data_fixed.data(),
1945 this->data_serializer.sizes_fixed_cumulative.back());
1947 if (this->data_serializer.variable_size_data_stored)
1950 this->data_serializer.dest_sizes_variable.resize(
1951 parallel_forest->local_num_quadrants);
1956 parallel_forest->global_first_quadrant,
1957 previous_global_first_quadrant,
1958 parallel_forest->mpicomm,
1960 this->data_serializer.dest_sizes_variable.data(),
1961 this->data_serializer.src_sizes_variable.data(),
1962 sizeof(
unsigned int));
1968 this->data_serializer.src_data_fixed.clear();
1969 this->data_serializer.src_data_fixed.shrink_to_fit();
1971 if (this->data_serializer.variable_size_data_stored)
1974 this->data_serializer.dest_data_variable.resize(
1975 std::accumulate(this->data_serializer.dest_sizes_variable.begin(),
1976 this->data_serializer.dest_sizes_variable.end(),
1977 std::vector<int>::size_type(0)));
1981 parallel_forest->global_first_quadrant,
1982 previous_global_first_quadrant,
1983 parallel_forest->mpicomm,
1985 this->data_serializer.dest_data_variable.data(),
1986 this->data_serializer.dest_sizes_variable.data(),
1987 this->data_serializer.src_data_variable.data(),
1988 this->data_serializer.src_sizes_variable.data());
1991 this->data_serializer.src_sizes_variable.clear();
1992 this->data_serializer.src_sizes_variable.shrink_to_fit();
1993 this->data_serializer.src_data_variable.clear();
1994 this->data_serializer.src_data_variable.shrink_to_fit();
2000 template <
int dim,
int spacedim>
2003 spacedim>::setup_coarse_cell_to_p4est_tree_permutation()
2008 coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
2010 cell_connectivity, coarse_cell_to_p4est_tree_permutation);
2012 p4est_tree_to_coarse_cell_permutation =
2003 spacedim>::setup_coarse_cell_to_p4est_tree_permutation() {
…}
2018 template <
int dim,
int spacedim>
2021 const
std::
string &file_basename)
const
2023 Assert(parallel_forest !=
nullptr,
2024 ExcMessage(
"Can't produce output when no forest is created yet."));
2028 "To use this function the triangulation's flag "
2029 "Settings::communicate_vertices_to_p4est must be set."));
2032 parallel_forest,
nullptr, file_basename.c_str());
2037 template <
int dim,
int spacedim>
2040 const
std::
string &file_basename)
const
2043 this->cell_attached_data.n_attached_deserialize == 0,
2045 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2046 Assert(this->n_cells() > 0,
2047 ExcMessage(
"Can not save() an empty Triangulation."));
2053 this->signals.pre_distributed_save();
2055 if (this->my_subdomain == 0)
2057 std::string fname = file_basename +
".info";
2058 std::ofstream f(fname);
2059 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2063 << this->cell_attached_data.pack_callbacks_fixed.size() <<
" "
2064 << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
2065 << this->n_cells(0) << std::endl;
2069 for ([[maybe_unused]]
const auto &cell_rel : this->local_cell_relations)
2071 Assert((cell_rel.second ==
2077 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2078 parallel_forest->global_num_quadrants,
2086 this->signals.post_distributed_save();
2091 template <
int dim,
int spacedim>
2096 this->n_cells() > 0,
2098 "load() only works if the Triangulation already contains a coarse mesh!"));
2100 this->n_levels() == 1,
2102 "Triangulation may only contain coarse cells when calling load()."));
2108 this->signals.pre_distributed_load();
2110 if (parallel_ghost !=
nullptr)
2114 parallel_ghost =
nullptr;
2117 parallel_forest =
nullptr;
2120 connectivity =
nullptr;
2122 unsigned int version, numcpus, attached_count_fixed,
2123 attached_count_variable, n_coarse_cells;
2125 std::string fname = std::string(file_basename) +
".info";
2126 std::ifstream f(fname);
2128 std::string firstline;
2129 getline(f, firstline);
2130 f >> version >> numcpus >> attached_count_fixed >>
2131 attached_count_variable >> n_coarse_cells;
2135 ExcMessage(
"Incompatible version found in .info file."));
2136 Assert(this->n_cells(0) == n_coarse_cells,
2137 ExcMessage(
"Number of coarse cells differ!"));
2141 this->cell_attached_data.n_attached_data_sets = 0;
2142 this->cell_attached_data.n_attached_deserialize =
2143 attached_count_fixed + attached_count_variable;
2146 file_basename.c_str(),
2147 this->mpi_communicator,
2165 copy_local_forest_to_triangulation();
2175 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2176 parallel_forest->global_num_quadrants,
2177 parallel_forest->local_num_quadrants,
2179 attached_count_fixed,
2180 attached_count_variable);
2183 this->signals.post_distributed_load();
2185 this->update_periodic_face_map();
2186 this->update_number_cache();
2191 template <
int dim,
int spacedim>
2194 const typename ::
internal::p4est::
types<dim>::forest *forest)
2196 Assert(this->n_cells() > 0,
2198 "load() only works if the Triangulation already contains "
2200 Assert(this->n_cells() == forest->trees->elem_count,
2202 "Coarse mesh of the Triangulation does not match the one "
2203 "of the provided forest!"));
2206 if (parallel_ghost !=
nullptr)
2210 parallel_ghost =
nullptr;
2213 parallel_forest =
nullptr;
2219 typename ::internal::p4est::types<dim>::forest *temp =
2220 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2224 parallel_forest->connectivity = connectivity;
2225 parallel_forest->user_pointer =
this;
2229 copy_local_forest_to_triangulation();
2238 this->update_periodic_face_map();
2239 this->update_number_cache();
2244 template <
int dim,
int spacedim>
2248 Assert(parallel_forest !=
nullptr,
2250 "Can't produce a check sum when no forest is created yet."));
2255# if !DEAL_II_P4EST_VERSION_GTE(2, 8, 6, 0)
2271 template <
int dim,
int spacedim>
2273 const typename ::internal::p4est::types<dim>::forest
2276 Assert(parallel_forest !=
nullptr,
2277 ExcMessage(
"The forest has not been allocated yet."));
2278 return parallel_forest;
2283 template <
int dim,
int spacedim>
2285 typename ::internal::p4est::types<dim>::tree
2287 const int dealii_coarse_cell_index)
const
2289 const unsigned int tree_index =
2290 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2291 typename ::internal::p4est::types<dim>::tree *tree =
2292 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2293 sc_array_index(parallel_forest->trees, tree_index));
2307 std::integral_constant<int, 2>)
2309 const unsigned int dim = 2, spacedim = 2;
2317 std::vector<unsigned int> vertex_touch_count;
2319 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2322 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2323 const ::internal::p4est::types<2>::locidx num_vtt =
2324 std::accumulate(vertex_touch_count.begin(),
2325 vertex_touch_count.end(),
2331 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2334 (set_vertex_info ==
true ? this->n_vertices() : 0),
2339 set_vertex_and_cell_info(*
this,
2342 coarse_cell_to_p4est_tree_permutation,
2346 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2351 this->mpi_communicator,
2368 std::integral_constant<int, 2>)
2370 const unsigned int dim = 2, spacedim = 3;
2378 std::vector<unsigned int> vertex_touch_count;
2380 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2383 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2384 const ::internal::p4est::types<2>::locidx num_vtt =
2385 std::accumulate(vertex_touch_count.begin(),
2386 vertex_touch_count.end(),
2392 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2395 (set_vertex_info ==
true ? this->n_vertices() : 0),
2400 set_vertex_and_cell_info(*
this,
2403 coarse_cell_to_p4est_tree_permutation,
2407 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2412 this->mpi_communicator,
2427 std::integral_constant<int, 3>)
2429 const int dim = 3, spacedim = 3;
2437 std::vector<unsigned int> vertex_touch_count;
2438 std::vector<std::list<
2439 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2441 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2442 const ::internal::p4est::types<2>::locidx num_vtt =
2443 std::accumulate(vertex_touch_count.begin(),
2444 vertex_touch_count.end(),
2447 std::vector<unsigned int> edge_touch_count;
2448 std::vector<std::list<
2449 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2451 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2452 const ::internal::p4est::types<2>::locidx num_ett =
2453 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2456 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2459 (set_vertex_info ==
true ? this->n_vertices() : 0),
2461 this->n_active_lines(),
2466 set_vertex_and_cell_info(*
this,
2469 coarse_cell_to_p4est_tree_permutation,
2498 const unsigned int deal_to_p4est_line_index[12] = {
2499 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2502 this->begin_active();
2503 cell != this->
end();
2506 const unsigned int index =
2507 coarse_cell_to_p4est_tree_permutation[cell->index()];
2508 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2510 deal_to_p4est_line_index[e]] =
2511 cell->line(e)->index();
2516 connectivity->ett_offset[0] = 0;
2517 std::partial_sum(edge_touch_count.begin(),
2518 edge_touch_count.end(),
2519 &connectivity->ett_offset[1]);
2521 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2524 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2526 Assert(edge_to_cell[v].
size() == edge_touch_count[v],
2530 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2531 unsigned int>>::const_iterator p =
2532 edge_to_cell[v].begin();
2533 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2535 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2536 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2537 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2538 deal_to_p4est_line_index[p->second];
2542 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2547 this->mpi_communicator,
2564 template <
int dim,
int spacedim>
2566 enforce_mesh_balance_over_periodic_boundaries(
2572 std::vector<bool> flags_before[2];
2576 std::vector<unsigned int> topological_vertex_numbering(
2578 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2579 topological_vertex_numbering[i] = i;
2597 using cell_iterator =
2599 for (
const auto &it : tria.get_periodic_face_map())
2601 const cell_iterator &cell_1 = it.first.first;
2602 const unsigned int face_no_1 = it.first.second;
2603 const cell_iterator &cell_2 = it.second.first.first;
2604 const unsigned int face_no_2 = it.second.first.second;
2605 const auto combined_orientation = it.second.second;
2607 if (cell_1->level() == cell_2->level())
2609 for (
const unsigned int v :
2614 const unsigned int vface1 =
2615 cell_1->reference_cell().standard_to_real_face_vertex(
2616 v, face_no_1, combined_orientation);
2617 const unsigned int vi1 =
2618 topological_vertex_numbering[cell_1->face(face_no_1)
2619 ->vertex_index(vface1)];
2620 const unsigned int vi2 =
2621 topological_vertex_numbering[cell_2->face(face_no_2)
2623 const unsigned int min_index =
std::min(vi1, vi2);
2624 topological_vertex_numbering[cell_1->face(face_no_1)
2625 ->vertex_index(vface1)] =
2626 topological_vertex_numbering[cell_2->face(face_no_2)
2627 ->vertex_index(v)] =
2636 for (
unsigned int i = 0; i < topological_vertex_numbering.size();
2639 const unsigned int j = topological_vertex_numbering[i];
2640 Assert(j == i || topological_vertex_numbering[j] == j,
2642 "Got inconclusive constraints with chain: " +
2643 std::to_string(i) +
" vs " + std::to_string(j) +
2644 " which should be equal to " +
2645 std::to_string(topological_vertex_numbering[j])));
2652 bool continue_iterating =
true;
2653 std::vector<int> vertex_level(tria.
n_vertices());
2654 while (continue_iterating)
2658 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2662 for (; cell != endc; ++cell)
2664 if (cell->refine_flag_set())
2665 for (
const unsigned int vertex :
2667 vertex_level[topological_vertex_numbering
2668 [cell->vertex_index(vertex)]] =
2669 std::
max(vertex_level[topological_vertex_numbering
2670 [cell->vertex_index(vertex)]],
2672 else if (!cell->coarsen_flag_set())
2673 for (
const unsigned int vertex :
2675 vertex_level[topological_vertex_numbering
2676 [cell->vertex_index(vertex)]] =
2677 std::
max(vertex_level[topological_vertex_numbering
2678 [cell->vertex_index(vertex)]],
2689 for (
const unsigned int vertex :
2691 vertex_level[topological_vertex_numbering
2692 [cell->vertex_index(vertex)]] =
2693 std::
max(vertex_level[topological_vertex_numbering
2694 [cell->vertex_index(vertex)]],
2699 continue_iterating =
false;
2710 for (cell = tria.
last_active(); cell != endc; --cell)
2711 if (cell->refine_flag_set() ==
false)
2713 for (
const unsigned int vertex :
2715 if (vertex_level[topological_vertex_numbering
2716 [cell->vertex_index(vertex)]] >=
2720 cell->clear_coarsen_flag();
2725 if (vertex_level[topological_vertex_numbering
2726 [cell->vertex_index(vertex)]] >
2729 cell->set_refine_flag();
2730 continue_iterating =
true;
2732 for (
const unsigned int v :
2734 vertex_level[topological_vertex_numbering
2735 [cell->vertex_index(v)]] =
2737 vertex_level[topological_vertex_numbering
2738 [cell->vertex_index(v)]],
2749 for (
const auto &cell : tria.cell_iterators())
2752 if (cell->is_active())
2755 const unsigned int n_children = cell->n_children();
2756 unsigned int flagged_children = 0;
2757 for (
unsigned int child = 0; child < n_children; ++child)
2758 if (cell->child(child)->is_active() &&
2759 cell->child(child)->coarsen_flag_set())
2764 if (flagged_children < n_children)
2765 for (
unsigned int child = 0; child < n_children; ++child)
2766 if (cell->child(child)->is_active())
2767 cell->child(child)->clear_coarsen_flag();
2770 std::vector<bool> flags_after[2];
2773 return ((flags_before[0] != flags_after[0]) ||
2774 (flags_before[1] != flags_after[1]));
2780 template <
int dim,
int spacedim>
2799 template <
int dim,
int spacedim>
2820 if (settings & construct_multigrid_hierarchy)
2823 spacedim>::limit_level_difference_at_vertices;
2827 bool mesh_changed =
false;
2835 if (settings & mesh_reconstruction_after_repartitioning)
2836 while (this->n_levels() > 1)
2843 for (
const auto &cell :
2844 this->active_cell_iterators_on_level(this->n_levels() - 1))
2846 cell->set_coarsen_flag();
2864 if (parallel_ghost !=
nullptr)
2868 parallel_ghost =
nullptr;
2872 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2873 P4EST_CONNECT_CORNER) :
2874 typename ::internal::p4est::types<dim>::balance_type(
2875 P8EST_CONNECT_CORNER)));
2882 for (
const auto &cell : this->cell_iterators_on_level(0))
2887 for (
const auto &cell : this->cell_iterators_on_level(0))
2892 if (tree_exists_locally<dim, spacedim>(
2894 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2897 delete_all_children<dim, spacedim>(cell);
2898 if (cell->is_active())
2907 typename ::internal::p4est::types<dim>::quadrant
2909 typename ::internal::p4est::types<dim>::tree *tree =
2910 init_tree(cell->index());
2912 ::internal::p4est::init_coarse_quadrant<dim>(
2915 match_tree_recursively<dim, spacedim>(*tree,
2919 this->my_subdomain);
2926 typename ::internal::p4est::types<dim>::quadrant *quadr;
2928 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2930 for (
unsigned int g_idx = 0;
2931 g_idx < parallel_ghost->ghosts.elem_count;
2934 while (g_idx >=
static_cast<unsigned int>(
2935 parallel_ghost->proc_offsets[ghost_owner + 1]))
2937 while (g_idx >=
static_cast<unsigned int>(
2938 parallel_ghost->tree_offsets[ghost_tree + 1]))
2941 quadr =
static_cast<
2942 typename ::internal::p4est::types<dim>::quadrant *
>(
2943 sc_array_index(¶llel_ghost->ghosts, g_idx));
2945 unsigned int coarse_cell_index =
2946 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2948 match_quadrant<dim, spacedim>(
this,
2963 bool mesh_changed =
true;
2964 unsigned int loop_counter = 0;
2971 this->update_periodic_face_map();
2974 enforce_mesh_balance_over_periodic_boundaries(*
this);
2986 "parallel::distributed::Triangulation::copy_local_forest_to_triangulation() "
2987 "for periodic boundaries detected. Aborting."));
2989 while (mesh_changed);
2994 std::any_of(this->begin_active(),
2997 return cell.refine_flag_set() ||
2998 cell.coarsen_flag_set();
3016 while (mesh_changed);
3021 unsigned int num_ghosts = 0;
3023 for (
const auto &cell : this->active_cell_iterators())
3025 if (cell->subdomain_id() != this->my_subdomain &&
3030 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
3041 if (settings & construct_multigrid_hierarchy)
3047 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3050 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3052 if ((cell->is_active() &&
3053 cell->subdomain_id() ==
3054 this->locally_owned_subdomain()) ||
3055 (cell->has_children() &&
3056 cell->child(0)->level_subdomain_id() ==
3057 this->locally_owned_subdomain()))
3058 cell->set_level_subdomain_id(
3059 this->locally_owned_subdomain());
3063 cell->set_level_subdomain_id(
3071 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3072 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3073 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3075 for (
const auto &cell : this->cell_iterators_on_level(0))
3077 typename ::internal::p4est::types<dim>::quadrant
3079 const unsigned int tree_index =
3080 coarse_cell_to_p4est_tree_permutation[cell->index()];
3081 typename ::internal::p4est::types<dim>::tree *tree =
3082 init_tree(cell->index());
3084 ::internal::p4est::init_coarse_quadrant<dim>(
3087 determine_level_subdomain_id_recursively<dim, spacedim>(
3098 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3101 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3103 if (cell->has_children())
3104 for (
unsigned int c = 0;
3105 c < GeometryInfo<dim>::max_children_per_cell;
3108 if (cell->child(c)->level_subdomain_id() ==
3109 this->locally_owned_subdomain())
3114 cell->child(0)->level_subdomain_id();
3118 cell->set_level_subdomain_id(mark);
3137 const unsigned int total_local_cells = this->n_active_cells();
3142 Assert(
static_cast<unsigned int>(
3143 parallel_forest->local_num_quadrants) ==
3149 Assert(
static_cast<unsigned int>(
3150 parallel_forest->local_num_quadrants) <=
3156 unsigned int n_owned = 0;
3157 for (
const auto &cell : this->active_cell_iterators())
3159 if (cell->subdomain_id() == this->my_subdomain)
3163 Assert(
static_cast<unsigned int>(
3164 parallel_forest->local_num_quadrants) == n_owned,
3169 this->smooth_grid = save_smooth;
3175 update_cell_relations();
3180 template <
int dim,
int spacedim>
3186 std::vector<Point<dim>> point{p};
3187 std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
3194 template <
int dim,
int spacedim>
3197 find_point_owner_rank(const
std::vector<
Point<dim>> &points)
3200 AssertThrow(this->are_vertices_communicated_to_p4est(),
3202 "Vertices need to be communicated to p4est to use this "
3203 "function. This must explicitly be turned on in the "
3204 "settings of the triangulation's constructor."));
3207 for (
const auto &manifold_id : this->get_manifold_ids())
3212 "This function can only be used if the triangulation "
3213 "has no other manifold than a Cartesian (flat) manifold attached."));
3217 PartitionSearch<dim> partition_search;
3223 parallel_forest->user_pointer = &partition_search;
3229 sc_array_t *point_sc_array;
3233 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3236 for (
size_t i = 0; i < points.size(); ++i)
3241 double *this_sc_point =
3242 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3244 for (
unsigned int d = 0; d < dim; ++d)
3246 this_sc_point[d] = p(d);
3248 this_sc_point[dim] = -1.0;
3254 static_cast<int>(
false),
3255 &PartitionSearch<dim>::local_quadrant_fn,
3256 &PartitionSearch<dim>::local_point_fn,
3260 std::vector<types::subdomain_id> owner_rank(
3264 for (
size_t i = 0; i < points.size(); ++i)
3267 double *this_sc_point =
3268 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3269 Assert(this_sc_point[dim] >= 0. || this_sc_point[dim] == -1.,
3271 if (this_sc_point[dim] < 0.)
3279 parallel_forest->user_pointer =
this;
3282 sc_array_destroy_null(&point_sc_array);
3289 template <
int dim,
int spacedim>
3296 for (
const auto &cell : this->active_cell_iterators())
3297 if (cell->is_locally_owned() && cell->refine_flag_set())
3298 Assert(cell->refine_flag_set() ==
3301 "This class does not support anisotropic refinement"));
3306 if (this->n_levels() ==
3310 cell = this->begin_active(
3318 !(cell->refine_flag_set()),
3320 "Fatal Error: maximum refinement level of p4est reached."));
3324 this->prepare_coarsening_and_refinement();
3327 this->signals.pre_distributed_refinement();
3332 for (
const auto &cell : this->active_cell_iterators())
3333 if (cell->is_ghost() || cell->is_artificial())
3335 cell->clear_refine_flag();
3336 cell->clear_coarsen_flag();
3342 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3343 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3350 parallel_forest->user_pointer = &refine_and_coarsen_list;
3352 if (parallel_ghost !=
nullptr)
3356 parallel_ghost =
nullptr;
3361 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3366 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3373 parallel_forest->user_pointer =
this;
3379 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3380 P4EST_CONNECT_FULL) :
3381 typename ::internal::p4est::types<dim>::balance_type(
3382 P8EST_CONNECT_FULL)),
3387 update_cell_relations();
3391 this->signals.post_p4est_refinement();
3395 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3396 previous_global_first_quadrant;
3398 if (this->cell_attached_data.n_attached_data_sets > 0)
3400 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3401 std::memcpy(previous_global_first_quadrant.data(),
3402 parallel_forest->global_first_quadrant,
3404 typename ::internal::p4est::types<dim>::gloidx) *
3405 (parallel_forest->mpisize + 1));
3408 if (!(settings & no_automatic_repartitioning))
3412 if (this->signals.weight.empty())
3420 const std::vector<unsigned int> cell_weights = get_cell_weights();
3426 this->mpi_communicator) > 0,
3428 "The global sum of weights over all active cells "
3429 "is zero. Please verify how you generate weights."));
3431 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3436 parallel_forest->user_pointer = &partition_weights;
3442 &PartitionWeights<dim, spacedim>::cell_weight);
3446 parallel_forest, 0,
nullptr,
nullptr);
3448 parallel_forest->user_pointer =
this;
3453 if (this->cell_attached_data.n_attached_data_sets > 0)
3455 this->data_serializer.pack_data(
3456 this->local_cell_relations,
3457 this->cell_attached_data.pack_callbacks_fixed,
3458 this->cell_attached_data.pack_callbacks_variable,
3459 this->get_mpi_communicator());
3465 for (
const auto &cell : this->active_cell_iterators())
3467 cell->clear_refine_flag();
3468 cell->clear_coarsen_flag();
3473 copy_local_forest_to_triangulation();
3483 if (this->cell_attached_data.n_attached_data_sets > 0)
3485 this->execute_transfer(parallel_forest,
3486 previous_global_first_quadrant.data());
3489 this->data_serializer.unpack_cell_status(this->local_cell_relations);
3511 if (settings & construct_multigrid_hierarchy)
3513 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3515 std::vector<bool> active_verts =
3516 this->mark_locally_active_vertices_on_level(lvl);
3518 const unsigned int maybe_coarser_lvl =
3519 (lvl > 0) ? (lvl - 1) : lvl;
3521 cell = this->begin(maybe_coarser_lvl),
3522 endc = this->end(lvl);
3523 for (; cell != endc; ++cell)
3524 if (cell->level() ==
static_cast<int>(lvl) ||
3527 const bool is_level_artificial =
3528 (cell->level_subdomain_id() ==
3530 bool need_to_know =
false;
3531 for (
const unsigned int vertex :
3533 if (active_verts[cell->vertex_index(vertex)])
3535 need_to_know =
true;
3540 !need_to_know || !is_level_artificial,
3542 "Internal error: the owner of cell" +
3543 cell->id().to_string() +
3544 " is unknown even though it is needed for geometric multigrid."));
3550 this->update_periodic_face_map();
3551 this->update_number_cache();
3554 this->signals.post_distributed_refinement();
3559 template <
int dim,
int spacedim>
3565 for (
const auto &cell : this->active_cell_iterators())
3566 if (cell->is_locally_owned())
3568 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3570 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3574 this->signals.pre_distributed_repartition();
3578 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3579 previous_global_first_quadrant;
3581 if (this->cell_attached_data.n_attached_data_sets > 0)
3583 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3584 std::memcpy(previous_global_first_quadrant.data(),
3585 parallel_forest->global_first_quadrant,
3587 typename ::internal::p4est::types<dim>::gloidx) *
3588 (parallel_forest->mpisize + 1));
3591 if (this->signals.weight.empty())
3603 const std::vector<unsigned int> cell_weights = get_cell_weights();
3609 this->mpi_communicator) > 0,
3611 "The global sum of weights over all active cells "
3612 "is zero. Please verify how you generate weights."));
3614 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3619 parallel_forest->user_pointer = &partition_weights;
3625 &PartitionWeights<dim, spacedim>::cell_weight);
3628 parallel_forest->user_pointer =
this;
3632 if (this->cell_attached_data.n_attached_data_sets > 0)
3634 this->data_serializer.pack_data(
3635 this->local_cell_relations,
3636 this->cell_attached_data.pack_callbacks_fixed,
3637 this->cell_attached_data.pack_callbacks_variable,
3638 this->get_mpi_communicator());
3643 copy_local_forest_to_triangulation();
3653 if (this->cell_attached_data.n_attached_data_sets > 0)
3655 this->execute_transfer(parallel_forest,
3656 previous_global_first_quadrant.data());
3659 this->update_periodic_face_map();
3662 this->update_number_cache();
3665 this->signals.post_distributed_repartition();
3670 template <
int dim,
int spacedim>
3672 const std::vector<types::global_dof_index>
3676 return p4est_tree_to_coarse_cell_permutation;
3681 template <
int dim,
int spacedim>
3683 const std::vector<types::global_dof_index>
3687 return coarse_cell_to_p4est_tree_permutation;
3692 template <
int dim,
int spacedim>
3695 mark_locally_active_vertices_on_level(const
int level)
const
3699 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3700 for (
const auto &cell : this->cell_iterators_on_level(
level))
3701 if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3703 marked_vertices[cell->vertex_index(v)] =
true;
3720 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3721 for (
const auto &it : this->get_periodic_face_map())
3724 const unsigned int face_no_1 = it.first.second;
3726 const unsigned int face_no_2 = it.second.first.second;
3727 const auto combined_orientation = it.second.second;
3728 const auto [orientation, rotation, flip] =
3731 if (cell_1->level() ==
level && cell_2->level() ==
level)
3733 for (
unsigned int v = 0;
3739 const unsigned int vface0 =
3741 v, orientation, flip, rotation);
3742 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3744 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3746 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3748 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3754 return marked_vertices;
3759 template <
int dim,
int spacedim>
3762 coarse_cell_id_to_coarse_cell_index(
3763 const
types::coarse_cell_id coarse_cell_id)
const
3765 return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3770 template <
int dim,
int spacedim>
3774 const unsigned int coarse_cell_index)
const
3776 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3781 template <
int dim,
int spacedim>
3785 &periodicity_vector)
3787 Assert(triangulation_has_content ==
true,
3789 Assert(this->n_levels() == 1,
3790 ExcMessage(
"The triangulation is refined!"));
3797 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
3799 for (
const auto &face_pair : periodicity_vector)
3803 const unsigned int face_left = face_pair.face_idx[0];
3804 const unsigned int face_right = face_pair.face_idx[1];
3807 const unsigned int tree_left =
3808 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3809 const unsigned int tree_right =
3810 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3818 unsigned int p4est_orientation = 0;
3822 p4est_orientation = face_pair.orientation ==
3829 const unsigned int face_idx_list[] = {face_left, face_right};
3830 const cell_iterator cell_list[] = {first_cell, second_cell};
3831 unsigned int lower_idx, higher_idx;
3833 if (face_left <= face_right)
3838 face_reference_cell.get_inverse_combined_orientation(
3839 face_pair.orientation);
3845 orientation = face_pair.orientation;
3850 unsigned int first_p4est_idx_on_cell =
3851 p8est_face_corners[face_idx_list[lower_idx]][0];
3852 unsigned int first_dealii_idx_on_face =
3854 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3857 const unsigned int first_dealii_idx_on_cell =
3859 face_idx_list[lower_idx],
3861 cell_list[lower_idx]->face_orientation(
3862 face_idx_list[lower_idx]),
3863 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3864 cell_list[lower_idx]->face_rotation(
3865 face_idx_list[lower_idx]));
3866 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3868 first_dealii_idx_on_face = i;
3876 const unsigned int second_dealii_idx_on_face =
3877 reference_cell.standard_to_real_face_vertex(
3878 first_dealii_idx_on_face,
3879 face_idx_list[lower_idx],
3881 const unsigned int second_dealii_idx_on_cell =
3882 reference_cell.face_to_cell_vertices(
3883 face_idx_list[higher_idx],
3884 second_dealii_idx_on_face,
3885 cell_list[higher_idx]->combined_face_orientation(
3886 face_idx_list[higher_idx]));
3888 const unsigned int second_p4est_idx_on_face =
3889 p8est_corner_face_corners[second_dealii_idx_on_cell]
3890 [face_idx_list[higher_idx]];
3891 p4est_orientation = second_p4est_idx_on_face;
3911 this->mpi_communicator,
3922 copy_local_forest_to_triangulation();
3933 this->update_number_cache();
3938 template <
int dim,
int spacedim>
3949 this->cell_attached_data.n_attached_data_sets) +
3956 coarse_cell_to_p4est_tree_permutation) +
3958 p4est_tree_to_coarse_cell_permutation) +
3959 memory_consumption_p4est();
3966 template <
int dim,
int spacedim>
3970 return ::internal::p4est::functions<dim>::forest_memory_used(
3978 template <
int dim,
int spacedim>
3985 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
3987 (other_tria.n_global_levels() == 1),
3998 const typename ::Triangulation<dim, spacedim>::DistortedCellList
4006 if (const ::parallel::distributed::Triangulation<dim, spacedim>
4007 *other_distributed =
4008 dynamic_cast<const ::parallel::distributed::
4009 Triangulation<dim, spacedim> *
>(&other_tria))
4012 settings = other_distributed->settings;
4013 triangulation_has_content =
4014 other_distributed->triangulation_has_content;
4015 coarse_cell_to_p4est_tree_permutation =
4016 other_distributed->coarse_cell_to_p4est_tree_permutation;
4017 p4est_tree_to_coarse_cell_permutation =
4018 other_distributed->p4est_tree_to_coarse_cell_permutation;
4021 typename ::internal::p4est::types<dim>::connectivity
4022 *temp_connectivity =
const_cast<
4023 typename ::internal::p4est::types<dim>::connectivity *
>(
4024 other_distributed->connectivity);
4026 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
4029 typename ::internal::p4est::types<dim>::forest *temp_forest =
4030 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
4031 other_distributed->parallel_forest);
4035 parallel_forest->connectivity = connectivity;
4036 parallel_forest->user_pointer =
this;
4040 triangulation_has_content =
true;
4041 setup_coarse_cell_to_p4est_tree_permutation();
4042 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
4047 copy_local_forest_to_triangulation();
4056 this->update_periodic_face_map();
4057 this->update_number_cache();
4062 template <
int dim,
int spacedim>
4067 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
4068 this->local_cell_relations.shrink_to_fit();
4071 for (
const auto &cell : this->cell_iterators_on_level(0))
4074 if (tree_exists_locally<dim, spacedim>(
4076 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4080 typename ::internal::p4est::types<dim>::quadrant
4082 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4085 typename ::internal::p4est::types<dim>::tree *tree =
4086 init_tree(cell->index());
4088 update_cell_relations_recursively<dim, spacedim>(
4089 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4095 template <
int dim,
int spacedim>
4102 Assert(this->local_cell_relations.size() ==
4103 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4108 std::vector<unsigned int> weights;
4109 weights.reserve(this->local_cell_relations.size());
4117 for (
const auto &cell_rel : this->local_cell_relations)
4119 const auto &cell_it = cell_rel.first;
4120 const auto &cell_status = cell_rel.second;
4122 weights.push_back(this->signals.weight(cell_it, cell_status));
4130 template <
int spacedim>
4146 template <
int spacedim>
4155 template <
int spacedim>
4157 const std::vector<types::global_dof_index>
4161 static std::vector<types::global_dof_index> a;
4167 template <
int spacedim>
4169 std::map<
unsigned int,
4172 const unsigned int )
const
4176 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4181 template <
int spacedim>
4184 mark_locally_active_vertices_on_level(const
unsigned int)
const
4187 return std::vector<bool>();
4192 template <
int spacedim>
4195 coarse_cell_id_to_coarse_cell_index(const
types::coarse_cell_id)
const
4203 template <
int spacedim>
4207 const unsigned int)
const
4215 template <
int spacedim>
4224 template <
int spacedim>
4233 template <
int spacedim>
4243 template <
int spacedim>
4253 template <
int spacedim>
4270 namespace distributed
4272 template <
int dim,
int spacedim>
4280#ifdef DEAL_II_WITH_P4EST
4290 const auto &cell = pair.first;
4291 const auto &status = pair.second;
4297 cell->clear_refine_flag();
4298 cell->clear_coarsen_flag();
4303 cell->clear_coarsen_flag();
4304 cell->set_refine_flag();
4309 for (
const auto &child : cell->child_iterators())
4311 child->clear_refine_flag();
4312 child->set_coarsen_flag();
4331 template <
int dim,
int spacedim>
4334#ifdef DEAL_II_WITH_P4EST
4335 if (distributed_tria)
4338 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4339 distributed_tria->load_refine_flags(saved_refine_flags);
4343 (void)distributed_tria;
4352#include "distributed/tria.inst"
@ children_will_be_coarsened
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator > > &)
active_cell_iterator last_active() const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
cell_iterator end() const
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const
virtual void execute_coarsening_and_refinement()
virtual bool prepare_coarsening_and_refinement()
void save_refine_flags(std::ostream &out) const
unsigned int n_vertices() const
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, types::geometric_orientation > > & get_periodic_face_map() const
void save_coarsen_flags(std::ostream &out) const
active_cell_iterator begin_active(const unsigned int level=0) const
virtual std::size_t memory_consumption() const override
virtual void clear() override
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
~TemporarilyMatchRefineFlags()
const ObserverPointer< ::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
constexpr bool running_in_debug_mode()
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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
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)
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)
T broadcast(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
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 types::geometric_orientation combined_face_orientation)
constexpr unsigned int invalid_unsigned_int
constexpr types::manifold_id flat_manifold_id
constexpr types::subdomain_id artificial_subdomain_id
constexpr types::subdomain_id invalid_subdomain_id
constexpr types::geometric_orientation default_geometric_orientation
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned char geometric_orientation
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)