40 #ifdef DEAL_II_WITH_P4EST
44 template <
int dim,
int spacedim>
46 get_vertex_to_cell_mappings(
48 std::vector<unsigned int> &vertex_touch_count,
49 std::vector<std::list<
51 unsigned int>>> &vertex_to_cell)
56 for (
const auto &cell :
triangulation.active_cell_iterators())
59 ++vertex_touch_count[cell->vertex_index(v)];
60 vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
66 template <
int dim,
int spacedim>
68 get_edge_to_cell_mappings(
70 std::vector<unsigned int> &edge_touch_count,
71 std::vector<std::list<
73 unsigned int>>> &edge_to_cell)
80 for (
const auto &cell :
triangulation.active_cell_iterators())
81 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
83 ++edge_touch_count[cell->line(
l)->index()];
84 edge_to_cell[cell->line(
l)->index()].emplace_back(cell,
l);
94 template <
int dim,
int spacedim>
96 set_vertex_and_cell_info(
98 const std::vector<unsigned int> &vertex_touch_count,
99 const std::vector<std::list<
101 unsigned int>>> &vertex_to_cell,
102 const std::vector<types::global_dof_index>
103 &coarse_cell_to_p4est_tree_permutation,
104 const bool set_vertex_info,
120 if (set_vertex_info ==
true)
121 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
123 connectivity->vertices[3 * v] =
triangulation.get_vertices()[v][0];
124 connectivity->vertices[3 * v + 1] =
126 connectivity->vertices[3 * v + 2] =
138 for (; cell != endc; ++cell)
140 const unsigned int index =
141 coarse_cell_to_p4est_tree_permutation[cell->index()];
145 if (set_vertex_info ==
true)
148 v] = cell->vertex_index(v);
151 v] = cell->vertex_index(v);
157 if (cell->face(f)->at_boundary() ==
false)
160 coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->index()];
164 coarse_cell_to_p4est_tree_permutation[cell->index()];
169 if (cell->face(f)->at_boundary() ==
false)
175 connectivity->tree_to_face
177 cell->neighbor_of_neighbor(f);
196 connectivity->tree_to_face[
index * 6 + f] =
197 cell->neighbor_of_neighbor(f);
199 unsigned int face_idx_list[2] = {
200 f, cell->neighbor_of_neighbor(f)};
202 cell_list[2] = {cell, cell->neighbor(f)};
203 unsigned int smaller_idx = 0;
205 if (f > cell->neighbor_of_neighbor(f))
208 unsigned int larger_idx = (smaller_idx + 1) % 2;
216 unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
218 face_idx_list[smaller_idx],
220 cell_list[smaller_idx]->face_orientation(
221 face_idx_list[smaller_idx]),
222 cell_list[smaller_idx]->face_flip(
223 face_idx_list[smaller_idx]),
224 cell_list[smaller_idx]->face_rotation(
225 face_idx_list[smaller_idx])));
229 for (
unsigned int i = 0;
230 i < GeometryInfo<dim>::vertices_per_face;
234 cell_list[larger_idx]->vertex_index(
236 face_idx_list[larger_idx], i));
245 connectivity->tree_to_face[
index * 6 + f] += 6 * v;
259 connectivity->ctt_offset[0] = 0;
260 std::partial_sum(vertex_touch_count.begin(),
261 vertex_touch_count.end(),
262 &connectivity->ctt_offset[1]);
265 std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
270 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
272 Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
276 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
277 unsigned int>>::const_iterator p =
278 vertex_to_cell[v].begin();
279 for (
unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
281 connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
282 coarse_cell_to_p4est_tree_permutation[p->first->index()];
283 connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
291 template <
int dim,
int spacedim>
297 Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
299 return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
300 (coarse_grid_cell <= parallel_forest->last_local_tree));
304 template <
int dim,
int spacedim>
306 delete_all_children_and_self(
309 if (cell->has_children())
310 for (
unsigned int c = 0; c < cell->n_children(); ++c)
311 delete_all_children_and_self<dim, spacedim>(cell->child(c));
313 cell->set_coarsen_flag();
318 template <
int dim,
int spacedim>
323 if (cell->has_children())
324 for (
unsigned int c = 0; c < cell->n_children(); ++c)
325 delete_all_children_and_self<dim, spacedim>(cell->child(c));
329 template <
int dim,
int spacedim>
331 determine_level_subdomain_id_recursively(
338 const std::vector<std::vector<bool>> &marked_vertices)
347 if (marked_vertices[dealii_cell->level()]
348 [dealii_cell->vertex_index(v)])
367 if (!used && dealii_cell->is_active() &&
368 dealii_cell->is_artificial() ==
false &&
369 dealii_cell->level() + 1 <
static_cast<int>(marked_vertices.size()))
373 if (marked_vertices[dealii_cell->level() + 1]
374 [dealii_cell->vertex_index(v)])
383 if (!used && dealii_cell->is_active() &&
384 dealii_cell->is_artificial() ==
false && dealii_cell->level() > 0)
388 if (marked_vertices[dealii_cell->level() - 1]
389 [dealii_cell->vertex_index(v)])
400 &forest, tree_index, &p4est_cell, my_subdomain);
401 Assert((owner != -2) && (owner != -1),
403 dealii_cell->set_level_subdomain_id(owner);
407 if (dealii_cell->has_children())
411 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
416 P4EST_QUADRANT_INIT(&p4est_child[c]);
419 P8EST_QUADRANT_INIT(&p4est_child[c]);
429 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
432 determine_level_subdomain_id_recursively<dim, spacedim>(
435 dealii_cell->child(c),
445 template <
int dim,
int spacedim>
447 match_tree_recursively(
455 if (sc_array_bsearch(
const_cast<sc_array_t *
>(&tree.quadrants),
461 delete_all_children<dim, spacedim>(dealii_cell);
462 if (dealii_cell->is_active())
463 dealii_cell->set_subdomain_id(my_subdomain);
472 if (dealii_cell->is_active())
473 dealii_cell->set_refine_flag();
478 for (
unsigned int c = 0;
479 c < GeometryInfo<dim>::max_children_per_cell;
484 P4EST_QUADRANT_INIT(&p4est_child[c]);
487 P8EST_QUADRANT_INIT(&p4est_child[c]);
497 for (
unsigned int c = 0;
498 c < GeometryInfo<dim>::max_children_per_cell;
503 &p4est_child[c]) ==
false)
509 delete_all_children<dim, spacedim>(dealii_cell->child(c));
510 dealii_cell->child(c)->recursively_set_subdomain_id(
517 match_tree_recursively<dim, spacedim>(tree,
518 dealii_cell->child(c),
528 template <
int dim,
int spacedim>
531 const ::Triangulation<dim, spacedim> *
tria,
532 unsigned int dealii_index,
536 const int l = ghost_quadrant.level;
538 for (
int i = 0; i <
l; ++i)
543 if (cell->is_active())
545 cell->clear_coarsen_flag();
546 cell->set_refine_flag();
553 dealii_index = cell->child_index(child_id);
559 if (cell->has_children())
560 delete_all_children<dim, spacedim>(cell);
563 cell->clear_coarsen_flag();
564 cell->set_subdomain_id(ghost_owner);
569 # ifdef P4EST_SEARCH_LOCAL
571 class PartitionSearch
579 PartitionSearch(
const PartitionSearch<dim> &other) =
delete;
581 PartitionSearch<dim> &
582 operator=(
const PartitionSearch<dim> &other) =
delete;
638 quad_length_on_level);
641 initialize_mapping();
644 map_real_to_unit_cell(
const Point<dim> &p)
const;
647 is_in_this_quadrant(
const Point<dim> &p)
const;
650 std::vector<Point<dim>> cell_vertices;
658 bool are_vertices_initialized;
660 bool is_reference_mapping_initialized;
666 QuadrantData quadrant_data;
673 PartitionSearch<dim>::local_quadrant_fn(
687 PartitionSearch<dim> *this_object =
688 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
692 quad_length_on_level =
694 (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
698 this_object->quadrant_data.set_cell_vertices(forest,
701 quad_length_on_level);
704 this_object->quadrant_data.initialize_mapping();
714 PartitionSearch<dim>::local_point_fn(
727 PartitionSearch<dim> *this_object =
728 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
731 double *this_point_dptr =
static_cast<double *
>(
point);
734 (dim == 2 ?
Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
735 Point<dim>(this_point_dptr[0],
737 this_point_dptr[2]));
740 const bool is_in_this_quadrant =
741 this_object->quadrant_data.is_in_this_quadrant(this_point);
745 if (!is_in_this_quadrant)
754 if (rank_begin < rank_end)
762 this_point_dptr[dim] =
static_cast<double>(rank_begin);
772 PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
775 const Point<dim> p_ref = map_real_to_unit_cell(p);
784 PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
787 Assert(is_reference_mapping_initialized,
789 "Cell vertices and mapping coefficients must be fully "
790 "initialized before transforming a point to the unit cell."));
796 for (
unsigned int alpha = 0;
797 alpha < GeometryInfo<dim>::vertices_per_cell;
803 p_out += (quadrant_mapping_matrix(alpha, 0) +
804 quadrant_mapping_matrix(alpha, 1) * p(0) +
805 quadrant_mapping_matrix(alpha, 2) * p(1) +
806 quadrant_mapping_matrix(alpha, 3) * p(0) * p(1)) *
812 for (
unsigned int alpha = 0;
813 alpha < GeometryInfo<dim>::vertices_per_cell;
819 p_out += (quadrant_mapping_matrix(alpha, 0) +
820 quadrant_mapping_matrix(alpha, 1) * p(0) +
821 quadrant_mapping_matrix(alpha, 2) * p(1) +
822 quadrant_mapping_matrix(alpha, 3) * p(2) +
823 quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
824 quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
825 quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
826 quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
836 PartitionSearch<dim>::QuadrantData::QuadrantData()
838 , quadrant_mapping_matrix(
GeometryInfo<dim>::vertices_per_cell,
840 , are_vertices_initialized(false)
841 , is_reference_mapping_initialized(false)
848 PartitionSearch<dim>::QuadrantData::initialize_mapping()
851 are_vertices_initialized,
853 "Cell vertices must be initialized before the cell mapping can be filled."));
860 for (
unsigned int alpha = 0;
861 alpha < GeometryInfo<dim>::vertices_per_cell;
865 point_matrix(0, alpha) = 1;
866 point_matrix(1, alpha) = cell_vertices[alpha](0);
867 point_matrix(2, alpha) = cell_vertices[alpha](1);
868 point_matrix(3, alpha) =
869 cell_vertices[alpha](0) * cell_vertices[alpha](1);
876 quadrant_mapping_matrix.invert(point_matrix);
880 for (
unsigned int alpha = 0;
881 alpha < GeometryInfo<dim>::vertices_per_cell;
885 point_matrix(0, alpha) = 1;
886 point_matrix(1, alpha) = cell_vertices[alpha](0);
887 point_matrix(2, alpha) = cell_vertices[alpha](1);
888 point_matrix(3, alpha) = cell_vertices[alpha](2);
889 point_matrix(4, alpha) =
890 cell_vertices[alpha](0) * cell_vertices[alpha](1);
891 point_matrix(5, alpha) =
892 cell_vertices[alpha](1) * cell_vertices[alpha](2);
893 point_matrix(6, alpha) =
894 cell_vertices[alpha](0) * cell_vertices[alpha](2);
895 point_matrix(7, alpha) = cell_vertices[alpha](0) *
896 cell_vertices[alpha](1) *
897 cell_vertices[alpha](2);
904 quadrant_mapping_matrix.invert(point_matrix);
907 is_reference_mapping_initialized =
true;
914 PartitionSearch<2>::QuadrantData::set_cell_vertices(
919 quad_length_on_level)
921 constexpr
unsigned int dim = 2;
925 double corner_point[dim + 1] = {0};
928 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
930 for (
unsigned int d = 0;
d < dim; ++
d)
932 cell_vertices[vertex_index](
d) = corner_point[
d];
942 unsigned int vertex_index = 0;
944 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
947 copy_vertex(vertex_index);
954 forest->connectivity,
956 quadrant->x + quad_length_on_level,
961 copy_vertex(vertex_index);
968 forest->connectivity,
971 quadrant->y + quad_length_on_level,
975 copy_vertex(vertex_index);
982 forest->connectivity,
984 quadrant->x + quad_length_on_level,
985 quadrant->y + quad_length_on_level,
989 copy_vertex(vertex_index);
991 are_vertices_initialized =
true;
998 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1003 quad_length_on_level)
1005 constexpr
unsigned int dim = 3;
1007 double corner_point[dim] = {0};
1010 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1012 for (
unsigned int d = 0;
d < dim; ++
d)
1014 cell_vertices[vertex_index](
d) = corner_point[
d];
1016 corner_point[
d] = 0;
1024 unsigned int vertex_index = 0;
1026 forest->connectivity,
1034 copy_vertex(vertex_index);
1042 forest->connectivity,
1044 quadrant->x + quad_length_on_level,
1050 copy_vertex(vertex_index);
1057 forest->connectivity,
1060 quadrant->y + quad_length_on_level,
1065 copy_vertex(vertex_index);
1072 forest->connectivity,
1074 quadrant->x + quad_length_on_level,
1075 quadrant->y + quad_length_on_level,
1080 copy_vertex(vertex_index);
1087 forest->connectivity,
1091 quadrant->z + quad_length_on_level,
1095 copy_vertex(vertex_index);
1102 forest->connectivity,
1104 quadrant->x + quad_length_on_level,
1106 quadrant->z + quad_length_on_level,
1110 copy_vertex(vertex_index);
1117 forest->connectivity,
1120 quadrant->y + quad_length_on_level,
1121 quadrant->z + quad_length_on_level,
1125 copy_vertex(vertex_index);
1132 forest->connectivity,
1134 quadrant->x + quad_length_on_level,
1135 quadrant->y + quad_length_on_level,
1136 quadrant->z + quad_length_on_level,
1140 copy_vertex(vertex_index);
1143 are_vertices_initialized =
true;
1153 template <
int dim,
int spacedim>
1154 class RefineAndCoarsenList
1158 const std::vector<types::global_dof_index>
1159 &p4est_tree_to_coarse_cell_permutation,
1187 pointers_are_at_end()
const;
1190 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1191 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1192 const_iterator current_refine_pointer;
1194 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1195 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1196 const_iterator current_coarsen_pointer;
1207 template <
int dim,
int spacedim>
1209 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1211 return ((current_refine_pointer == refine_list.end()) &&
1212 (current_coarsen_pointer == coarsen_list.end()));
1217 template <
int dim,
int spacedim>
1218 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1220 const std::vector<types::global_dof_index>
1221 &p4est_tree_to_coarse_cell_permutation,
1225 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1226 for (
const auto &cell :
triangulation.active_cell_iterators())
1229 if (cell->subdomain_id() != my_subdomain)
1232 if (cell->refine_flag_set())
1234 else if (cell->coarsen_flag_set())
1238 refine_list.reserve(n_refine_flags);
1239 coarsen_list.reserve(n_coarsen_flags);
1251 unsigned int coarse_cell_index =
1252 p4est_tree_to_coarse_cell_permutation[c];
1261 p4est_cell.p.which_tree = c;
1262 build_lists(cell, p4est_cell, my_subdomain);
1270 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1271 Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1273 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1274 Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1277 current_refine_pointer = refine_list.begin();
1278 current_coarsen_pointer = coarsen_list.begin();
1283 template <
int dim,
int spacedim>
1285 RefineAndCoarsenList<dim, spacedim>::build_lists(
1290 if (cell->is_active())
1292 if (cell->subdomain_id() == my_subdomain)
1294 if (cell->refine_flag_set())
1295 refine_list.push_back(p4est_cell);
1296 else if (cell->coarsen_flag_set())
1297 coarsen_list.push_back(p4est_cell);
1304 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1309 P4EST_QUADRANT_INIT(&p4est_child[c]);
1312 P8EST_QUADRANT_INIT(&p4est_child[c]);
1319 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1322 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1323 build_lists(cell->child(c), p4est_child[c], my_subdomain);
1329 template <
int dim,
int spacedim>
1331 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1336 RefineAndCoarsenList<dim, spacedim> *this_object =
1337 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1338 forest->user_pointer);
1342 if (this_object->current_refine_pointer == this_object->refine_list.end())
1345 Assert(coarse_cell_index <=
1346 this_object->current_refine_pointer->p.which_tree,
1351 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1355 Assert(coarse_cell_index <=
1356 this_object->current_refine_pointer->p.which_tree,
1362 quadrant, &*this_object->current_refine_pointer) <= 0,
1367 quadrant, &*this_object->current_refine_pointer))
1369 ++this_object->current_refine_pointer;
1379 template <
int dim,
int spacedim>
1381 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1386 RefineAndCoarsenList<dim, spacedim> *this_object =
1387 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1388 forest->user_pointer);
1392 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1395 Assert(coarse_cell_index <=
1396 this_object->current_coarsen_pointer->p.which_tree,
1401 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1405 Assert(coarse_cell_index <=
1406 this_object->current_coarsen_pointer->p.which_tree,
1412 children[0], &*this_object->current_coarsen_pointer) <= 0,
1418 children[0], &*this_object->current_coarsen_pointer))
1421 ++this_object->current_coarsen_pointer;
1425 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1429 children[c], &*this_object->current_coarsen_pointer),
1431 ++this_object->current_coarsen_pointer;
1449 template <
int dim,
int spacedim>
1450 class PartitionWeights
1458 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1473 std::vector<unsigned int> cell_weights_list;
1474 std::vector<unsigned int>::const_iterator current_pointer;
1478 template <
int dim,
int spacedim>
1479 PartitionWeights<dim, spacedim>::PartitionWeights(
1480 const std::vector<unsigned int> &cell_weights)
1481 : cell_weights_list(cell_weights)
1485 current_pointer = cell_weights_list.begin();
1489 template <
int dim,
int spacedim>
1491 PartitionWeights<dim, spacedim>::cell_weight(
1500 PartitionWeights<dim, spacedim> *this_object =
1501 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1503 Assert(this_object->current_pointer >=
1504 this_object->cell_weights_list.begin(),
1506 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1512 const unsigned int weight = *this_object->current_pointer;
1513 ++this_object->current_pointer;
1516 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1517 "weights for cells. The weight provided here exceeds "
1518 "the maximum value represented as a 'signed int'."));
1519 return static_cast<int>(weight);
1522 template <
int dim,
int spacedim>
1523 using cell_relation_t =
typename std::pair<
1524 typename ::Triangulation<dim, spacedim>::cell_iterator,
1536 template <
int dim,
int spacedim>
1538 add_single_cell_relation(
1539 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1540 const typename ::internal::p4est::types<dim>::tree &tree,
1541 const unsigned int idx,
1545 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1551 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1565 template <
int dim,
int spacedim>
1567 update_cell_relations_recursively(
1568 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1569 const typename ::internal::p4est::types<dim>::tree &tree,
1571 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1574 const int idx = sc_array_bsearch(
1575 const_cast<sc_array_t *
>(&tree.quadrants),
1580 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1582 &p4est_cell) ==
false))
1587 const bool p4est_has_children = (idx == -1);
1588 if (p4est_has_children && dealii_cell->has_children())
1591 typename ::internal::p4est::types<dim>::quadrant
1594 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1599 P4EST_QUADRANT_INIT(&p4est_child[c]);
1602 P8EST_QUADRANT_INIT(&p4est_child[c]);
1609 &p4est_cell, p4est_child);
1611 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1614 update_cell_relations_recursively<dim, spacedim>(
1615 cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1618 else if (!p4est_has_children && !dealii_cell->has_children())
1622 add_single_cell_relation<dim, spacedim>(
1625 else if (p4est_has_children)
1633 typename ::internal::p4est::types<dim>::quadrant
1635 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1640 P4EST_QUADRANT_INIT(&p4est_child[c]);
1643 P8EST_QUADRANT_INIT(&p4est_child[c]);
1650 &p4est_cell, p4est_child);
1658 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1661 child_idx = sc_array_bsearch(
1662 const_cast<sc_array_t *
>(&tree.quadrants),
1669 add_single_cell_relation<dim, spacedim>(
1670 cell_rel, tree, child_idx, dealii_cell, cell_status);
1678 add_single_cell_relation<dim, spacedim>(
1692 namespace distributed
1695 template <
int dim,
int spacedim>
1698 const MPI_Comm mpi_communicator,
1711 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1715 , triangulation_has_content(false)
1716 , connectivity(
nullptr)
1717 , parallel_forest(
nullptr)
1719 parallel_ghost =
nullptr;
1724 template <
int dim,
int spacedim>
1745 template <
int dim,
int spacedim>
1749 const std::vector<
CellData<dim>> &cells,
1758 const typename ::Triangulation<dim, spacedim>::DistortedCellList
1767 this->all_reference_cells_are_hyper_cube(),
1769 "The class parallel::distributed::Triangulation only supports meshes "
1770 "consisting only of hypercube-like cells."));
1775 triangulation_has_content =
true;
1777 setup_coarse_cell_to_p4est_tree_permutation();
1779 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1783 copy_local_forest_to_triangulation();
1792 this->update_periodic_face_map();
1793 this->update_number_cache();
1798 template <
int dim,
int spacedim>
1804 (void)construction_data;
1811 template <
int dim,
int spacedim>
1815 triangulation_has_content =
false;
1817 if (parallel_ghost !=
nullptr)
1821 parallel_ghost =
nullptr;
1824 if (parallel_forest !=
nullptr)
1827 parallel_forest =
nullptr;
1830 if (connectivity !=
nullptr)
1834 connectivity =
nullptr;
1837 coarse_cell_to_p4est_tree_permutation.resize(0);
1838 p4est_tree_to_coarse_cell_permutation.resize(0);
1842 this->update_number_cache();
1847 template <
int dim,
int spacedim>
1858 template <
int dim,
int spacedim>
1869 template <
int dim,
int spacedim>
1875 *previous_global_first_quadrant)
1877 Assert(this->data_serializer.sizes_fixed_cumulative.size() > 0,
1881 this->data_serializer.dest_data_fixed.resize(
1882 parallel_forest->local_num_quadrants *
1883 this->data_serializer.sizes_fixed_cumulative.back());
1886 typename ::internal::p4est::types<dim>::transfer_context
1890 parallel_forest->global_first_quadrant,
1891 previous_global_first_quadrant,
1892 parallel_forest->mpicomm,
1894 this->data_serializer.dest_data_fixed.data(),
1895 this->data_serializer.src_data_fixed.data(),
1896 this->data_serializer.sizes_fixed_cumulative.back());
1898 if (this->data_serializer.variable_size_data_stored)
1901 this->data_serializer.dest_sizes_variable.resize(
1902 parallel_forest->local_num_quadrants);
1907 parallel_forest->global_first_quadrant,
1908 previous_global_first_quadrant,
1909 parallel_forest->mpicomm,
1911 this->data_serializer.dest_sizes_variable.data(),
1912 this->data_serializer.src_sizes_variable.data(),
1913 sizeof(
unsigned int));
1919 this->data_serializer.src_data_fixed.clear();
1920 this->data_serializer.src_data_fixed.shrink_to_fit();
1922 if (this->data_serializer.variable_size_data_stored)
1925 this->data_serializer.dest_data_variable.resize(
1926 std::accumulate(this->data_serializer.dest_sizes_variable.begin(),
1927 this->data_serializer.dest_sizes_variable.end(),
1930 # if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1937 if (this->data_serializer.src_sizes_variable.empty())
1938 this->data_serializer.src_sizes_variable.resize(1);
1939 if (this->data_serializer.dest_sizes_variable.empty())
1940 this->data_serializer.dest_sizes_variable.resize(1);
1945 parallel_forest->global_first_quadrant,
1946 previous_global_first_quadrant,
1947 parallel_forest->mpicomm,
1949 this->data_serializer.dest_data_variable.data(),
1950 this->data_serializer.dest_sizes_variable.data(),
1951 this->data_serializer.src_data_variable.data(),
1952 this->data_serializer.src_sizes_variable.data());
1955 this->data_serializer.src_sizes_variable.clear();
1956 this->data_serializer.src_sizes_variable.shrink_to_fit();
1957 this->data_serializer.src_data_variable.clear();
1958 this->data_serializer.src_data_variable.shrink_to_fit();
1964 template <
int dim,
int spacedim>
1967 spacedim>::setup_coarse_cell_to_p4est_tree_permutation()
1972 coarse_cell_to_p4est_tree_permutation.resize(this->
n_cells(0));
1974 cell_connectivity, coarse_cell_to_p4est_tree_permutation);
1976 p4est_tree_to_coarse_cell_permutation =
1982 template <
int dim,
int spacedim>
1985 const std::
string &file_basename)
const
1987 Assert(parallel_forest !=
nullptr,
1988 ExcMessage(
"Can't produce output when no forest is created yet."));
1992 "To use this function the triangulation's flag "
1993 "Settings::communicate_vertices_to_p4est must be set."));
1996 parallel_forest,
nullptr, file_basename.c_str());
2001 template <
int dim,
int spacedim>
2006 this->cell_attached_data.n_attached_deserialize == 0,
2008 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2010 ExcMessage(
"Can not save() an empty Triangulation."));
2016 this->signals.pre_distributed_save();
2018 if (this->my_subdomain == 0)
2020 std::string fname = std::string(filename) +
".info";
2021 std::ofstream f(fname);
2022 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2026 << this->cell_attached_data.pack_callbacks_fixed.size() <<
" "
2027 << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
2028 << this->
n_cells(0) << std::endl;
2032 for (
const auto &cell_rel : this->local_cell_relations)
2035 Assert((cell_rel.second ==
2041 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2042 parallel_forest->global_num_quadrants,
2050 this->signals.post_distributed_save();
2055 template <
int dim,
int spacedim>
2062 "load() only works if the Triangulation already contains a coarse mesh!"));
2064 this->n_levels() == 1,
2066 "Triangulation may only contain coarse cells when calling load()."));
2072 this->signals.pre_distributed_load();
2074 if (parallel_ghost !=
nullptr)
2078 parallel_ghost =
nullptr;
2081 parallel_forest =
nullptr;
2084 connectivity =
nullptr;
2086 unsigned int version, numcpus, attached_count_fixed,
2087 attached_count_variable, n_coarse_cells;
2089 std::string fname = std::string(filename) +
".info";
2090 std::ifstream f(fname);
2092 std::string firstline;
2093 getline(f, firstline);
2094 f >> version >> numcpus >> attached_count_fixed >>
2095 attached_count_variable >> n_coarse_cells;
2099 ExcMessage(
"Incompatible version found in .info file."));
2101 ExcMessage(
"Number of coarse cells differ!"));
2105 this->cell_attached_data.n_attached_data_sets = 0;
2106 this->cell_attached_data.n_attached_deserialize =
2107 attached_count_fixed + attached_count_variable;
2111 this->mpi_communicator,
2129 copy_local_forest_to_triangulation();
2139 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2140 parallel_forest->global_num_quadrants,
2141 parallel_forest->local_num_quadrants,
2143 attached_count_fixed,
2144 attached_count_variable);
2147 this->signals.post_distributed_load();
2149 this->update_periodic_face_map();
2150 this->update_number_cache();
2155 template <
int dim,
int spacedim>
2158 const
bool autopartition)
2160 (void)autopartition;
2166 template <
int dim,
int spacedim>
2169 const typename ::
internal::p4est::
types<dim>::forest *forest)
2173 "load() only works if the Triangulation already contains "
2177 "Coarse mesh of the Triangulation does not match the one "
2178 "of the provided forest!"));
2181 if (parallel_ghost !=
nullptr)
2185 parallel_ghost =
nullptr;
2188 parallel_forest =
nullptr;
2194 typename ::internal::p4est::types<dim>::forest *temp =
2195 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2199 parallel_forest->connectivity = connectivity;
2200 parallel_forest->user_pointer =
this;
2204 copy_local_forest_to_triangulation();
2213 this->update_periodic_face_map();
2214 this->update_number_cache();
2219 template <
int dim,
int spacedim>
2223 Assert(parallel_forest !=
nullptr,
2225 "Can't produce a check sum when no forest is created yet."));
2226 return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2231 template <
int dim,
int spacedim>
2233 const typename ::internal::p4est::types<dim>::forest
2236 Assert(parallel_forest !=
nullptr,
2237 ExcMessage(
"The forest has not been allocated yet."));
2238 return parallel_forest;
2243 template <
int dim,
int spacedim>
2245 typename ::internal::p4est::types<dim>::tree
2247 const int dealii_coarse_cell_index)
const
2249 const unsigned int tree_index =
2250 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2251 typename ::internal::p4est::types<dim>::tree *tree =
2252 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2253 sc_array_index(parallel_forest->trees, tree_index));
2267 std::integral_constant<int, 2>)
2269 const unsigned int dim = 2, spacedim = 2;
2277 std::vector<unsigned int> vertex_touch_count;
2279 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2282 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2283 const ::internal::p4est::types<2>::locidx num_vtt =
2284 std::accumulate(vertex_touch_count.begin(),
2285 vertex_touch_count.end(),
2291 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2294 (set_vertex_info ==
true ? this->n_vertices() : 0),
2299 set_vertex_and_cell_info(*
this,
2302 coarse_cell_to_p4est_tree_permutation,
2306 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2311 this->mpi_communicator,
2328 std::integral_constant<int, 2>)
2330 const unsigned int dim = 2, spacedim = 3;
2338 std::vector<unsigned int> vertex_touch_count;
2340 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2343 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2344 const ::internal::p4est::types<2>::locidx num_vtt =
2345 std::accumulate(vertex_touch_count.begin(),
2346 vertex_touch_count.end(),
2352 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2355 (set_vertex_info ==
true ? this->n_vertices() : 0),
2360 set_vertex_and_cell_info(*
this,
2363 coarse_cell_to_p4est_tree_permutation,
2367 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2372 this->mpi_communicator,
2387 std::integral_constant<int, 3>)
2389 const int dim = 3, spacedim = 3;
2397 std::vector<unsigned int> vertex_touch_count;
2398 std::vector<std::list<
2399 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2401 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2402 const ::internal::p4est::types<2>::locidx num_vtt =
2403 std::accumulate(vertex_touch_count.begin(),
2404 vertex_touch_count.end(),
2407 std::vector<unsigned int> edge_touch_count;
2408 std::vector<std::list<
2409 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2411 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2412 const ::internal::p4est::types<2>::locidx num_ett =
2413 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2416 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2419 (set_vertex_info ==
true ? this->n_vertices() : 0),
2421 this->n_active_lines(),
2426 set_vertex_and_cell_info(*
this,
2429 coarse_cell_to_p4est_tree_permutation,
2458 const unsigned int deal_to_p4est_line_index[12] = {
2459 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2462 this->begin_active();
2463 cell != this->
end();
2466 const unsigned int index =
2467 coarse_cell_to_p4est_tree_permutation[cell->index()];
2468 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2470 deal_to_p4est_line_index[
e]] =
2471 cell->line(
e)->index();
2476 connectivity->ett_offset[0] = 0;
2477 std::partial_sum(edge_touch_count.begin(),
2478 edge_touch_count.end(),
2479 &connectivity->ett_offset[1]);
2481 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2484 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2486 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2490 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2491 unsigned int>>::const_iterator p =
2492 edge_to_cell[v].begin();
2493 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2495 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2496 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2497 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2498 deal_to_p4est_line_index[p->second];
2502 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2507 this->mpi_communicator,
2524 template <
int dim,
int spacedim>
2526 enforce_mesh_balance_over_periodic_boundaries(
2532 std::vector<bool> flags_before[2];
2536 std::vector<unsigned int> topological_vertex_numbering(
2538 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2539 topological_vertex_numbering[i] = i;
2557 using cell_iterator =
2559 typename std::map<std::pair<cell_iterator, unsigned int>,
2560 std::pair<std::pair<cell_iterator, unsigned int>,
2561 std::bitset<3>>>::const_iterator it;
2566 const cell_iterator &cell_1 = it->first.first;
2567 const unsigned int face_no_1 = it->first.second;
2568 const cell_iterator &cell_2 = it->second.first.first;
2569 const unsigned int face_no_2 = it->second.first.second;
2570 const std::bitset<3> face_orientation = it->second.second;
2572 if (cell_1->level() == cell_2->level())
2574 for (
unsigned int v = 0;
2580 const unsigned int vface0 =
2583 face_orientation[0],
2584 face_orientation[1],
2585 face_orientation[2]);
2586 const unsigned int vi0 =
2587 topological_vertex_numbering[cell_1->face(face_no_1)
2588 ->vertex_index(vface0)];
2589 const unsigned int vi1 =
2590 topological_vertex_numbering[cell_2->face(face_no_2)
2592 const unsigned int min_index =
std::min(vi0, vi1);
2593 topological_vertex_numbering[cell_1->face(face_no_1)
2594 ->vertex_index(vface0)] =
2595 topological_vertex_numbering[cell_2->face(face_no_2)
2596 ->vertex_index(v)] =
2603 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2605 const unsigned int j = topological_vertex_numbering[i];
2613 bool continue_iterating =
true;
2615 while (continue_iterating)
2619 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2623 for (; cell != endc; ++cell)
2625 if (cell->refine_flag_set())
2626 for (
const unsigned int vertex :
2628 vertex_level[topological_vertex_numbering
2629 [cell->vertex_index(vertex)]] =
2630 std::max(vertex_level[topological_vertex_numbering
2631 [cell->vertex_index(vertex)]],
2633 else if (!cell->coarsen_flag_set())
2634 for (
const unsigned int vertex :
2636 vertex_level[topological_vertex_numbering
2637 [cell->vertex_index(vertex)]] =
2638 std::max(vertex_level[topological_vertex_numbering
2639 [cell->vertex_index(vertex)]],
2650 for (
const unsigned int vertex :
2652 vertex_level[topological_vertex_numbering
2653 [cell->vertex_index(vertex)]] =
2654 std::max(vertex_level[topological_vertex_numbering
2655 [cell->vertex_index(vertex)]],
2660 continue_iterating =
false;
2672 if (cell->refine_flag_set() ==
false)
2674 for (
const unsigned int vertex :
2676 if (vertex_level[topological_vertex_numbering
2677 [cell->vertex_index(vertex)]] >=
2681 cell->clear_coarsen_flag();
2686 if (vertex_level[topological_vertex_numbering
2687 [cell->vertex_index(vertex)]] >
2690 cell->set_refine_flag();
2691 continue_iterating =
true;
2693 for (
const unsigned int v :
2695 vertex_level[topological_vertex_numbering
2696 [cell->vertex_index(v)]] =
2698 vertex_level[topological_vertex_numbering
2699 [cell->vertex_index(v)]],
2713 if (cell->is_active())
2716 const unsigned int n_children = cell->n_children();
2717 unsigned int flagged_children = 0;
2718 for (
unsigned int child = 0; child < n_children; ++child)
2719 if (cell->child(child)->is_active() &&
2720 cell->child(child)->coarsen_flag_set())
2725 if (flagged_children < n_children)
2726 for (
unsigned int child = 0; child < n_children; ++child)
2727 if (cell->child(child)->is_active())
2728 cell->child(child)->clear_coarsen_flag();
2731 std::vector<bool> flags_after[2];
2734 return ((flags_before[0] != flags_after[0]) ||
2735 (flags_before[1] != flags_after[1]));
2741 template <
int dim,
int spacedim>
2745 bool mesh_changed =
false;
2746 unsigned int loop_counter = 0;
2747 unsigned int n_changes = 0;
2752 this->update_periodic_face_map();
2754 mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*
this);
2755 n_changes += mesh_changed;
2766 "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2767 "for periodic boundaries detected. Aborting."));
2769 while (mesh_changed);
2772 return n_changes > 0;
2777 template <
int dim,
int spacedim>
2801 spacedim>::limit_level_difference_at_vertices;
2805 bool mesh_changed =
false;
2813 if (
settings & mesh_reconstruction_after_repartitioning)
2814 while (this->n_levels() > 1)
2821 for (
const auto &cell :
2822 this->active_cell_iterators_on_level(this->n_levels() - 1))
2824 cell->set_coarsen_flag();
2842 if (parallel_ghost !=
nullptr)
2846 parallel_ghost =
nullptr;
2850 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2851 P4EST_CONNECT_CORNER) :
2852 typename ::internal::p4est::types<dim>::balance_type(
2853 P8EST_CONNECT_CORNER)));
2860 for (
const auto &cell : this->cell_iterators_on_level(0))
2865 for (
const auto &cell : this->cell_iterators_on_level(0))
2870 if (tree_exists_locally<dim, spacedim>(
2872 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2875 delete_all_children<dim, spacedim>(cell);
2876 if (cell->is_active())
2885 typename ::internal::p4est::types<dim>::quadrant
2887 typename ::internal::p4est::types<dim>::tree *tree =
2888 init_tree(cell->index());
2890 ::internal::p4est::init_coarse_quadrant<dim>(
2893 match_tree_recursively<dim, spacedim>(*tree,
2897 this->my_subdomain);
2904 typename ::internal::p4est::types<dim>::quadrant *quadr;
2906 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2908 for (
unsigned int g_idx = 0;
2909 g_idx < parallel_ghost->ghosts.elem_count;
2912 while (g_idx >=
static_cast<unsigned int>(
2913 parallel_ghost->proc_offsets[ghost_owner + 1]))
2915 while (g_idx >=
static_cast<unsigned int>(
2916 parallel_ghost->tree_offsets[ghost_tree + 1]))
2919 quadr =
static_cast<
2920 typename ::internal::p4est::types<dim>::quadrant *
>(
2921 sc_array_index(¶llel_ghost->ghosts, g_idx));
2923 unsigned int coarse_cell_index =
2924 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2926 match_quadrant<dim, spacedim>(
this,
2933 this->prepare_coarsening_and_refinement();
2937 std::any_of(this->begin_active(),
2940 return cell.refine_flag_set() ||
2941 cell.coarsen_flag_set();
2959 while (mesh_changed);
2963 unsigned int num_ghosts = 0;
2965 for (
const auto &cell : this->active_cell_iterators())
2967 if (cell->subdomain_id() != this->my_subdomain &&
2972 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2989 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
2992 for (
const auto &cell : this->cell_iterators_on_level(lvl))
2994 if ((cell->is_active() &&
2995 cell->subdomain_id() ==
2997 (cell->has_children() &&
2998 cell->child(0)->level_subdomain_id() ==
3000 cell->set_level_subdomain_id(
3001 this->locally_owned_subdomain());
3005 cell->set_level_subdomain_id(
3013 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3014 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3015 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3017 for (
const auto &cell : this->cell_iterators_on_level(0))
3019 typename ::internal::p4est::types<dim>::quadrant
3021 const unsigned int tree_index =
3022 coarse_cell_to_p4est_tree_permutation[cell->index()];
3023 typename ::internal::p4est::types<dim>::tree *tree =
3024 init_tree(cell->index());
3026 ::internal::p4est::init_coarse_quadrant<dim>(
3029 determine_level_subdomain_id_recursively<dim, spacedim>(
3040 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3043 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3045 if (cell->has_children())
3046 for (
unsigned int c = 0;
3047 c < GeometryInfo<dim>::max_children_per_cell;
3050 if (cell->child(c)->level_subdomain_id() ==
3056 cell->child(0)->level_subdomain_id();
3060 cell->set_level_subdomain_id(mark);
3077 (void)total_local_cells;
3081 Assert(
static_cast<unsigned int>(
3082 parallel_forest->local_num_quadrants) == total_local_cells,
3087 Assert(
static_cast<unsigned int>(
3088 parallel_forest->local_num_quadrants) <= total_local_cells,
3094 unsigned int n_owned = 0;
3095 for (
const auto &cell : this->active_cell_iterators())
3097 if (cell->subdomain_id() == this->my_subdomain)
3101 Assert(
static_cast<unsigned int>(
3102 parallel_forest->local_num_quadrants) == n_owned,
3107 this->smooth_grid = save_smooth;
3113 update_cell_relations();
3118 template <
int dim,
int spacedim>
3124 std::vector<Point<dim>>
point{p};
3125 std::vector<types::subdomain_id> owner = find_point_owner_rank(
point);
3132 template <
int dim,
int spacedim>
3135 find_point_owner_rank(const std::vector<
Point<dim>> &points)
3137 # ifndef P4EST_SEARCH_LOCAL
3142 "This function is only available if p4est is version 2.2 and higher."));
3144 return std::vector<unsigned int>(1,
3148 AssertThrow(this->are_vertices_communicated_to_p4est(),
3150 "Vertices need to be communicated to p4est to use this "
3151 "function. This must explicitly be turned on in the "
3152 "settings of the triangulation's constructor."));
3155 for (
const auto &
manifold_id : this->get_manifold_ids())
3160 "This function can only be used if the triangulation "
3161 "has no other manifold than a Cartesian (flat) manifold attached."));
3165 PartitionSearch<dim> partition_search;
3171 parallel_forest->user_pointer = &partition_search;
3177 sc_array_t *point_sc_array;
3181 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3184 for (
size_t i = 0; i < points.size(); ++i)
3189 double *this_sc_point =
3190 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3192 for (
unsigned int d = 0;
d < dim; ++
d)
3194 this_sc_point[
d] = p(
d);
3196 this_sc_point[dim] = -1.0;
3202 static_cast<int>(
false),
3203 &PartitionSearch<dim>::local_quadrant_fn,
3204 &PartitionSearch<dim>::local_point_fn,
3208 std::vector<types::subdomain_id> owner_rank(
3212 for (
size_t i = 0; i < points.size(); ++i)
3215 double *this_sc_point =
3216 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3221 parallel_forest->user_pointer =
this;
3224 sc_array_destroy_null(&point_sc_array);
3232 template <
int dim,
int spacedim>
3238 for (
const auto &cell : this->active_cell_iterators())
3239 if (cell->is_locally_owned() && cell->refine_flag_set())
3240 Assert(cell->refine_flag_set() ==
3243 "This class does not support anisotropic refinement"));
3248 if (this->n_levels() ==
3252 cell = this->begin_active(
3260 !(cell->refine_flag_set()),
3262 "Fatal Error: maximum refinement level of p4est reached."));
3266 this->prepare_coarsening_and_refinement();
3269 this->signals.pre_distributed_refinement();
3274 for (
const auto &cell : this->active_cell_iterators())
3275 if (cell->is_ghost() || cell->is_artificial())
3277 cell->clear_refine_flag();
3278 cell->clear_coarsen_flag();
3284 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3285 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3292 parallel_forest->user_pointer = &refine_and_coarsen_list;
3294 if (parallel_ghost !=
nullptr)
3298 parallel_ghost =
nullptr;
3303 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3308 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3315 parallel_forest->user_pointer =
this;
3321 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3322 P4EST_CONNECT_FULL) :
3323 typename ::internal::p4est::types<dim>::balance_type(
3324 P8EST_CONNECT_FULL)),
3329 update_cell_relations();
3333 this->signals.post_p4est_refinement();
3337 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3338 previous_global_first_quadrant;
3340 if (this->cell_attached_data.n_attached_data_sets > 0)
3342 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3343 std::memcpy(previous_global_first_quadrant.data(),
3344 parallel_forest->global_first_quadrant,
3346 typename ::internal::p4est::types<dim>::gloidx) *
3347 (parallel_forest->mpisize + 1));
3350 if (!(
settings & no_automatic_repartitioning))
3354 if (this->signals.weight.empty())
3362 const std::vector<unsigned int> cell_weights = get_cell_weights();
3368 this->mpi_communicator) > 0,
3370 "The global sum of weights over all active cells "
3371 "is zero. Please verify how you generate weights."));
3373 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3378 parallel_forest->user_pointer = &partition_weights;
3384 &PartitionWeights<dim, spacedim>::cell_weight);
3388 parallel_forest, 0,
nullptr,
nullptr);
3390 parallel_forest->user_pointer =
this;
3395 if (this->cell_attached_data.n_attached_data_sets > 0)
3397 this->data_serializer.pack_data(
3398 this->local_cell_relations,
3399 this->cell_attached_data.pack_callbacks_fixed,
3400 this->cell_attached_data.pack_callbacks_variable,
3401 this->get_communicator());
3407 for (
const auto &cell : this->active_cell_iterators())
3409 cell->clear_refine_flag();
3410 cell->clear_coarsen_flag();
3415 copy_local_forest_to_triangulation();
3425 if (this->cell_attached_data.n_attached_data_sets > 0)
3427 this->execute_transfer(parallel_forest,
3428 previous_global_first_quadrant.data());
3431 this->data_serializer.unpack_cell_status(this->local_cell_relations);
3452 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3454 std::vector<bool> active_verts =
3455 this->mark_locally_active_vertices_on_level(lvl);
3457 const unsigned int maybe_coarser_lvl =
3458 (lvl > 0) ? (lvl - 1) : lvl;
3460 cell = this->
begin(maybe_coarser_lvl),
3461 endc = this->
end(lvl);
3462 for (; cell != endc; ++cell)
3463 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3465 const bool is_level_artificial =
3466 (cell->level_subdomain_id() ==
3468 bool need_to_know =
false;
3469 for (
const unsigned int vertex :
3471 if (active_verts[cell->vertex_index(vertex)])
3473 need_to_know =
true;
3478 !need_to_know || !is_level_artificial,
3480 "Internal error: the owner of cell" +
3481 cell->id().to_string() +
3482 " is unknown even though it is needed for geometric multigrid."));
3488 this->update_periodic_face_map();
3489 this->update_number_cache();
3492 this->signals.post_distributed_refinement();
3497 template <
int dim,
int spacedim>
3502 for (
const auto &cell : this->active_cell_iterators())
3503 if (cell->is_locally_owned())
3505 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3507 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3511 this->signals.pre_distributed_repartition();
3515 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3516 previous_global_first_quadrant;
3518 if (this->cell_attached_data.n_attached_data_sets > 0)
3520 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3521 std::memcpy(previous_global_first_quadrant.data(),
3522 parallel_forest->global_first_quadrant,
3524 typename ::internal::p4est::types<dim>::gloidx) *
3525 (parallel_forest->mpisize + 1));
3528 if (this->signals.weight.empty())
3540 const std::vector<unsigned int> cell_weights = get_cell_weights();
3546 this->mpi_communicator) > 0,
3548 "The global sum of weights over all active cells "
3549 "is zero. Please verify how you generate weights."));
3551 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3556 parallel_forest->user_pointer = &partition_weights;
3562 &PartitionWeights<dim, spacedim>::cell_weight);
3565 parallel_forest->user_pointer =
this;
3569 if (this->cell_attached_data.n_attached_data_sets > 0)
3571 this->data_serializer.pack_data(
3572 this->local_cell_relations,
3573 this->cell_attached_data.pack_callbacks_fixed,
3574 this->cell_attached_data.pack_callbacks_variable,
3575 this->get_communicator());
3580 copy_local_forest_to_triangulation();
3590 if (this->cell_attached_data.n_attached_data_sets > 0)
3592 this->execute_transfer(parallel_forest,
3593 previous_global_first_quadrant.data());
3596 this->update_periodic_face_map();
3599 this->update_number_cache();
3602 this->signals.post_distributed_repartition();
3607 template <
int dim,
int spacedim>
3609 const std::vector<types::global_dof_index>
3613 return p4est_tree_to_coarse_cell_permutation;
3618 template <
int dim,
int spacedim>
3620 const std::vector<types::global_dof_index>
3624 return coarse_cell_to_p4est_tree_permutation;
3629 template <
int dim,
int spacedim>
3632 mark_locally_active_vertices_on_level(const
int level)
const
3636 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3637 for (
const auto &cell : this->cell_iterators_on_level(
level))
3640 marked_vertices[cell->vertex_index(v)] =
true;
3657 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3658 for (
const auto &it : this->get_periodic_face_map())
3661 const unsigned int face_no_1 = it.first.second;
3663 const unsigned int face_no_2 = it.second.first.second;
3664 const std::bitset<3> &face_orientation = it.second.second;
3666 if (cell_1->level() ==
level && cell_2->level() ==
level)
3668 for (
unsigned int v = 0;
3674 const unsigned int vface0 =
3677 face_orientation[0],
3678 face_orientation[1],
3679 face_orientation[2]);
3680 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3682 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3684 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3686 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3692 return marked_vertices;
3697 template <
int dim,
int spacedim>
3700 coarse_cell_id_to_coarse_cell_index(
3708 template <
int dim,
int spacedim>
3712 const unsigned int coarse_cell_index)
const
3714 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3719 template <
int dim,
int spacedim>
3723 &periodicity_vector)
3725 Assert(triangulation_has_content ==
true,
3727 Assert(this->n_levels() == 1,
3728 ExcMessage(
"The triangulation is refined!"));
3735 for (
const auto &face_pair : periodicity_vector)
3739 const unsigned int face_left = face_pair.face_idx[0];
3740 const unsigned int face_right = face_pair.face_idx[1];
3743 const unsigned int tree_left =
3744 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3745 const unsigned int tree_right =
3746 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3755 unsigned int p4est_orientation = 0;
3757 p4est_orientation = face_pair.orientation[1];
3760 const unsigned int face_idx_list[] = {face_left, face_right};
3761 const cell_iterator cell_list[] = {first_cell, second_cell};
3762 unsigned int lower_idx, higher_idx;
3763 if (face_left <= face_right)
3776 unsigned int first_p4est_idx_on_cell =
3777 p8est_face_corners[face_idx_list[lower_idx]][0];
3778 unsigned int first_dealii_idx_on_face =
3780 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3783 const unsigned int first_dealii_idx_on_cell =
3785 face_idx_list[lower_idx],
3787 cell_list[lower_idx]->face_orientation(
3788 face_idx_list[lower_idx]),
3789 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3790 cell_list[lower_idx]->face_rotation(
3791 face_idx_list[lower_idx]));
3792 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3794 first_dealii_idx_on_face = i;
3801 constexpr
unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
3809 constexpr
unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
3817 const unsigned int second_dealii_idx_on_face =
3818 lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
3819 [first_dealii_idx_on_face] :
3820 right_to_left[face_pair.orientation.to_ulong()]
3821 [first_dealii_idx_on_face];
3822 const unsigned int second_dealii_idx_on_cell =
3824 face_idx_list[higher_idx],
3825 second_dealii_idx_on_face,
3826 cell_list[higher_idx]->face_orientation(
3827 face_idx_list[higher_idx]),
3828 cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
3829 cell_list[higher_idx]->face_rotation(
3830 face_idx_list[higher_idx]));
3832 const unsigned int second_p4est_idx_on_face =
3833 p8est_corner_face_corners[second_dealii_idx_on_cell]
3834 [face_idx_list[higher_idx]];
3835 p4est_orientation = second_p4est_idx_on_face;
3855 this->mpi_communicator,
3866 copy_local_forest_to_triangulation();
3877 this->update_number_cache();
3882 template <
int dim,
int spacedim>
3893 this->cell_attached_data.n_attached_data_sets) +
3900 coarse_cell_to_p4est_tree_permutation) +
3902 p4est_tree_to_coarse_cell_permutation) +
3903 memory_consumption_p4est();
3910 template <
int dim,
int spacedim>
3914 return ::internal::p4est::functions<dim>::forest_memory_used(
3922 template <
int dim,
int spacedim>
3929 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
3931 (other_tria.n_global_levels() == 1),
3942 const typename ::Triangulation<dim, spacedim>::DistortedCellList
3950 if (const ::parallel::distributed::Triangulation<dim, spacedim>
3951 *other_distributed =
3952 dynamic_cast<const ::parallel::distributed::
3953 Triangulation<dim, spacedim> *
>(&other_tria))
3956 settings = other_distributed->settings;
3957 triangulation_has_content =
3958 other_distributed->triangulation_has_content;
3959 coarse_cell_to_p4est_tree_permutation =
3960 other_distributed->coarse_cell_to_p4est_tree_permutation;
3961 p4est_tree_to_coarse_cell_permutation =
3962 other_distributed->p4est_tree_to_coarse_cell_permutation;
3965 typename ::internal::p4est::types<dim>::connectivity
3966 *temp_connectivity =
const_cast<
3967 typename ::internal::p4est::types<dim>::connectivity *
>(
3968 other_distributed->connectivity);
3970 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
3973 typename ::internal::p4est::types<dim>::forest *temp_forest =
3974 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
3975 other_distributed->parallel_forest);
3979 parallel_forest->connectivity = connectivity;
3980 parallel_forest->user_pointer =
this;
3984 triangulation_has_content =
true;
3985 setup_coarse_cell_to_p4est_tree_permutation();
3986 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
3991 copy_local_forest_to_triangulation();
4000 this->update_periodic_face_map();
4001 this->update_number_cache();
4006 template <
int dim,
int spacedim>
4011 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
4012 this->local_cell_relations.shrink_to_fit();
4015 for (
const auto &cell : this->cell_iterators_on_level(0))
4018 if (tree_exists_locally<dim, spacedim>(
4020 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4024 typename ::internal::p4est::types<dim>::quadrant
4026 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4029 typename ::internal::p4est::types<dim>::tree *tree =
4030 init_tree(cell->index());
4032 update_cell_relations_recursively<dim, spacedim>(
4033 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4039 template <
int dim,
int spacedim>
4046 Assert(this->local_cell_relations.size() ==
4047 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4052 std::vector<unsigned int> weights;
4053 weights.reserve(this->local_cell_relations.size());
4061 for (
const auto &cell_rel : this->local_cell_relations)
4063 const auto &cell_it = cell_rel.first;
4064 const auto &cell_status = cell_rel.second;
4066 weights.push_back(this->signals.weight(cell_it, cell_status));
4074 template <
int spacedim>
4077 const MPI_Comm mpi_communicator,
4090 template <
int spacedim>
4099 template <
int spacedim>
4101 const std::vector<types::global_dof_index>
4105 static std::vector<types::global_dof_index> a;
4111 template <
int spacedim>
4113 std::map<
unsigned int,
4116 const unsigned int )
const
4120 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4125 template <
int spacedim>
4128 mark_locally_active_vertices_on_level(const
unsigned int)
const
4131 return std::vector<bool>();
4136 template <
int spacedim>
4147 template <
int spacedim>
4151 const unsigned int)
const
4159 template <
int spacedim>
4168 template <
int spacedim>
4177 template <
int spacedim>
4186 template <
int spacedim>
4196 template <
int spacedim>
4206 template <
int spacedim>
4223 namespace distributed
4225 template <
int dim,
int spacedim>
4233 #ifdef DEAL_II_WITH_P4EST
4243 const auto &cell = pair.first;
4244 const auto &status = pair.second;
4250 cell->clear_refine_flag();
4251 cell->clear_coarsen_flag();
4256 cell->clear_coarsen_flag();
4257 cell->set_refine_flag();
4262 for (
const auto &child : cell->child_iterators())
4264 child->clear_refine_flag();
4265 child->set_coarsen_flag();
4284 template <
int dim,
int spacedim>
4287 #ifdef DEAL_II_WITH_P4EST
4288 if (distributed_tria)
4291 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4292 distributed_tria->load_refine_flags(saved_refine_flags);
4296 (void)distributed_tria;
4305 #include "tria.inst"
@ children_will_be_coarsened
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
active_cell_iterator last_active() const
virtual types::subdomain_id locally_owned_subdomain() const
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator >> &)
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
virtual bool prepare_coarsening_and_refinement()
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > & get_periodic_face_map() const
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 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)
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
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcIO()
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
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
types::global_dof_index size_type
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 > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
@ construct_multigrid_hierarchy
VectorType::value_type * begin(VectorType &V)
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)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 3 > &c)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 3 > &c)
bool tree_exists_locally(const typename types< dim >::forest *parallel_forest, const typename types< dim >::topidx coarse_grid_cell)
const types::subdomain_id artificial_subdomain_id
const types::subdomain_id invalid_subdomain_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
global_cell_index coarse_cell_id
unsigned int subdomain_id
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 Point< dim > unit_cell_vertex(const unsigned int vertex)
static bool is_inside_unit_cell(const Point< dim > &p)
const TriangulationDescription::Settings settings
const ::Triangulation< dim, spacedim > & tria