42 #ifdef DEAL_II_WITH_P4EST
46 template <
int dim,
int spacedim>
48 get_vertex_to_cell_mappings(
50 std::vector<unsigned int> & vertex_touch_count,
51 std::vector<std::list<
53 unsigned int>>> & vertex_to_cell)
58 for (
const auto &cell :
triangulation.active_cell_iterators())
61 ++vertex_touch_count[cell->vertex_index(v)];
62 vertex_to_cell[cell->vertex_index(v)].emplace_back(cell, v);
68 template <
int dim,
int spacedim>
70 get_edge_to_cell_mappings(
72 std::vector<unsigned int> & edge_touch_count,
73 std::vector<std::list<
75 unsigned int>>> & edge_to_cell)
82 for (
const auto &cell :
triangulation.active_cell_iterators())
83 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
85 ++edge_touch_count[cell->line(
l)->index()];
86 edge_to_cell[cell->line(
l)->index()].emplace_back(cell,
l);
96 template <
int dim,
int spacedim>
98 set_vertex_and_cell_info(
100 const std::vector<unsigned int> & vertex_touch_count,
101 const std::vector<std::list<
103 unsigned int>>> & vertex_to_cell,
104 const std::vector<types::global_dof_index>
105 & coarse_cell_to_p4est_tree_permutation,
106 const bool set_vertex_info,
122 if (set_vertex_info ==
true)
123 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
125 connectivity->vertices[3 * v] =
triangulation.get_vertices()[v][0];
126 connectivity->vertices[3 * v + 1] =
128 connectivity->vertices[3 * v + 2] =
140 for (; cell != endc; ++cell)
142 const unsigned int index =
143 coarse_cell_to_p4est_tree_permutation[cell->index()];
147 if (set_vertex_info ==
true)
150 v] = cell->vertex_index(v);
153 v] = cell->vertex_index(v);
159 if (cell->face(f)->at_boundary() ==
false)
162 coarse_cell_to_p4est_tree_permutation[cell->neighbor(f)->index()];
166 coarse_cell_to_p4est_tree_permutation[cell->index()];
171 if (cell->face(f)->at_boundary() ==
false)
177 connectivity->tree_to_face
179 cell->neighbor_of_neighbor(f);
198 connectivity->tree_to_face[
index * 6 + f] =
199 cell->neighbor_of_neighbor(f);
201 unsigned int face_idx_list[2] = {
202 f, cell->neighbor_of_neighbor(f)};
204 cell_list[2] = {cell, cell->neighbor(f)};
205 unsigned int smaller_idx = 0;
207 if (f > cell->neighbor_of_neighbor(f))
210 unsigned int larger_idx = (smaller_idx + 1) % 2;
218 unsigned int g_idx = cell_list[smaller_idx]->vertex_index(
220 face_idx_list[smaller_idx],
222 cell_list[smaller_idx]->face_orientation(
223 face_idx_list[smaller_idx]),
224 cell_list[smaller_idx]->face_flip(
225 face_idx_list[smaller_idx]),
226 cell_list[smaller_idx]->face_rotation(
227 face_idx_list[smaller_idx])));
231 for (
unsigned int i = 0;
232 i < GeometryInfo<dim>::vertices_per_face;
236 cell_list[larger_idx]->vertex_index(
238 face_idx_list[larger_idx], i));
247 connectivity->tree_to_face[
index * 6 + f] += 6 * v;
261 connectivity->ctt_offset[0] = 0;
262 std::partial_sum(vertex_touch_count.begin(),
263 vertex_touch_count.end(),
264 &connectivity->ctt_offset[1]);
267 std::accumulate(vertex_touch_count.begin(), vertex_touch_count.end(), 0u);
272 for (
unsigned int v = 0; v <
triangulation.n_vertices(); ++v)
274 Assert(vertex_to_cell[v].size() == vertex_touch_count[v],
278 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
279 unsigned int>>::const_iterator p =
280 vertex_to_cell[v].begin();
281 for (
unsigned int c = 0; c < vertex_touch_count[v]; ++c, ++p)
283 connectivity->corner_to_tree[connectivity->ctt_offset[v] + c] =
284 coarse_cell_to_p4est_tree_permutation[p->first->index()];
285 connectivity->corner_to_corner[connectivity->ctt_offset[v] + c] =
293 template <
int dim,
int spacedim>
299 Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees,
301 return ((coarse_grid_cell >= parallel_forest->first_local_tree) &&
302 (coarse_grid_cell <= parallel_forest->last_local_tree));
306 template <
int dim,
int spacedim>
308 delete_all_children_and_self(
311 if (cell->has_children())
312 for (
unsigned int c = 0; c < cell->n_children(); ++c)
313 delete_all_children_and_self<dim, spacedim>(cell->child(c));
315 cell->set_coarsen_flag();
320 template <
int dim,
int spacedim>
325 if (cell->has_children())
326 for (
unsigned int c = 0; c < cell->n_children(); ++c)
327 delete_all_children_and_self<dim, spacedim>(cell->child(c));
331 template <
int dim,
int spacedim>
333 determine_level_subdomain_id_recursively(
340 const std::vector<std::vector<bool>> & marked_vertices)
349 if (marked_vertices[dealii_cell->level()]
350 [dealii_cell->vertex_index(v)])
369 if (!used && dealii_cell->is_active() &&
370 dealii_cell->is_artificial() ==
false &&
371 dealii_cell->level() + 1 <
static_cast<int>(marked_vertices.size()))
375 if (marked_vertices[dealii_cell->level() + 1]
376 [dealii_cell->vertex_index(v)])
385 if (!used && dealii_cell->is_active() &&
386 dealii_cell->is_artificial() ==
false && dealii_cell->level() > 0)
390 if (marked_vertices[dealii_cell->level() - 1]
391 [dealii_cell->vertex_index(v)])
402 &forest, tree_index, &p4est_cell, my_subdomain);
403 Assert((owner != -2) && (owner != -1),
405 dealii_cell->set_level_subdomain_id(owner);
409 if (dealii_cell->has_children())
413 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
418 P4EST_QUADRANT_INIT(&p4est_child[c]);
421 P8EST_QUADRANT_INIT(&p4est_child[c]);
431 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
434 determine_level_subdomain_id_recursively<dim, spacedim>(
437 dealii_cell->child(c),
447 template <
int dim,
int spacedim>
449 match_tree_recursively(
457 if (sc_array_bsearch(
const_cast<sc_array_t *
>(&tree.quadrants),
463 delete_all_children<dim, spacedim>(dealii_cell);
464 if (dealii_cell->is_active())
465 dealii_cell->set_subdomain_id(my_subdomain);
474 if (dealii_cell->is_active())
475 dealii_cell->set_refine_flag();
480 for (
unsigned int c = 0;
481 c < GeometryInfo<dim>::max_children_per_cell;
486 P4EST_QUADRANT_INIT(&p4est_child[c]);
489 P8EST_QUADRANT_INIT(&p4est_child[c]);
499 for (
unsigned int c = 0;
500 c < GeometryInfo<dim>::max_children_per_cell;
505 &p4est_child[c]) ==
false)
511 delete_all_children<dim, spacedim>(dealii_cell->child(c));
512 dealii_cell->child(c)->recursively_set_subdomain_id(
519 match_tree_recursively<dim, spacedim>(tree,
520 dealii_cell->child(c),
530 template <
int dim,
int spacedim>
533 const ::Triangulation<dim, spacedim> *
tria,
534 unsigned int dealii_index,
538 const int l = ghost_quadrant.level;
540 for (
int i = 0; i <
l; ++i)
545 if (cell->is_active())
547 cell->clear_coarsen_flag();
548 cell->set_refine_flag();
555 dealii_index = cell->child_index(child_id);
561 if (cell->has_children())
562 delete_all_children<dim, spacedim>(cell);
565 cell->clear_coarsen_flag();
566 cell->set_subdomain_id(ghost_owner);
571 # ifdef P4EST_SEARCH_LOCAL
573 class PartitionSearch
581 PartitionSearch(
const PartitionSearch<dim> &other) =
delete;
583 PartitionSearch<dim> &
584 operator=(
const PartitionSearch<dim> &other) =
delete;
640 quad_length_on_level);
643 initialize_mapping();
646 map_real_to_unit_cell(
const Point<dim> &p)
const;
649 is_in_this_quadrant(
const Point<dim> &p)
const;
652 std::vector<Point<dim>> cell_vertices;
660 bool are_vertices_initialized;
662 bool is_reference_mapping_initialized;
668 QuadrantData quadrant_data;
675 PartitionSearch<dim>::local_quadrant_fn(
689 PartitionSearch<dim> *this_object =
690 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
694 quad_length_on_level =
696 (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
700 this_object->quadrant_data.set_cell_vertices(forest,
703 quad_length_on_level);
706 this_object->quadrant_data.initialize_mapping();
716 PartitionSearch<dim>::local_point_fn(
729 PartitionSearch<dim> *this_object =
730 reinterpret_cast<PartitionSearch<dim> *
>(forest->user_pointer);
733 double *this_point_dptr =
static_cast<double *
>(
point);
736 (dim == 2 ?
Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
737 Point<dim>(this_point_dptr[0],
739 this_point_dptr[2]));
742 const bool is_in_this_quadrant =
743 this_object->quadrant_data.is_in_this_quadrant(this_point);
747 if (!is_in_this_quadrant)
756 if (rank_begin < rank_end)
764 this_point_dptr[dim] =
static_cast<double>(rank_begin);
774 PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
777 const Point<dim> p_ref = map_real_to_unit_cell(p);
786 PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
789 Assert(is_reference_mapping_initialized,
791 "Cell vertices and mapping coefficients must be fully "
792 "initialized before transforming a point to the unit cell."));
798 for (
unsigned int alpha = 0;
799 alpha < GeometryInfo<dim>::vertices_per_cell;
805 p_out += (quadrant_mapping_matrix(alpha, 0) +
806 quadrant_mapping_matrix(alpha, 1) * p(0) +
807 quadrant_mapping_matrix(alpha, 2) * p(1) +
808 quadrant_mapping_matrix(alpha, 3) * p(0) * p(1)) *
814 for (
unsigned int alpha = 0;
815 alpha < GeometryInfo<dim>::vertices_per_cell;
821 p_out += (quadrant_mapping_matrix(alpha, 0) +
822 quadrant_mapping_matrix(alpha, 1) * p(0) +
823 quadrant_mapping_matrix(alpha, 2) * p(1) +
824 quadrant_mapping_matrix(alpha, 3) * p(2) +
825 quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
826 quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
827 quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
828 quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
838 PartitionSearch<dim>::QuadrantData::QuadrantData()
840 , quadrant_mapping_matrix(
GeometryInfo<dim>::vertices_per_cell,
842 , are_vertices_initialized(false)
843 , is_reference_mapping_initialized(false)
850 PartitionSearch<dim>::QuadrantData::initialize_mapping()
853 are_vertices_initialized,
855 "Cell vertices must be initialized before the cell mapping can be filled."));
862 for (
unsigned int alpha = 0;
863 alpha < GeometryInfo<dim>::vertices_per_cell;
867 point_matrix(0, alpha) = 1;
868 point_matrix(1, alpha) = cell_vertices[alpha](0);
869 point_matrix(2, alpha) = cell_vertices[alpha](1);
870 point_matrix(3, alpha) =
871 cell_vertices[alpha](0) * cell_vertices[alpha](1);
878 quadrant_mapping_matrix.invert(point_matrix);
882 for (
unsigned int alpha = 0;
883 alpha < GeometryInfo<dim>::vertices_per_cell;
887 point_matrix(0, alpha) = 1;
888 point_matrix(1, alpha) = cell_vertices[alpha](0);
889 point_matrix(2, alpha) = cell_vertices[alpha](1);
890 point_matrix(3, alpha) = cell_vertices[alpha](2);
891 point_matrix(4, alpha) =
892 cell_vertices[alpha](0) * cell_vertices[alpha](1);
893 point_matrix(5, alpha) =
894 cell_vertices[alpha](1) * cell_vertices[alpha](2);
895 point_matrix(6, alpha) =
896 cell_vertices[alpha](0) * cell_vertices[alpha](2);
897 point_matrix(7, alpha) = cell_vertices[alpha](0) *
898 cell_vertices[alpha](1) *
899 cell_vertices[alpha](2);
906 quadrant_mapping_matrix.invert(point_matrix);
909 is_reference_mapping_initialized =
true;
916 PartitionSearch<2>::QuadrantData::set_cell_vertices(
921 quad_length_on_level)
923 constexpr
unsigned int dim = 2;
927 double corner_point[dim + 1] = {0};
930 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
932 for (
unsigned int d = 0;
d < dim; ++
d)
934 cell_vertices[vertex_index](
d) = corner_point[
d];
944 unsigned int vertex_index = 0;
946 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
949 copy_vertex(vertex_index);
956 forest->connectivity,
958 quadrant->x + quad_length_on_level,
963 copy_vertex(vertex_index);
970 forest->connectivity,
973 quadrant->y + quad_length_on_level,
977 copy_vertex(vertex_index);
984 forest->connectivity,
986 quadrant->x + quad_length_on_level,
987 quadrant->y + quad_length_on_level,
991 copy_vertex(vertex_index);
993 are_vertices_initialized =
true;
1000 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1005 quad_length_on_level)
1007 constexpr
unsigned int dim = 3;
1009 double corner_point[dim] = {0};
1012 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1014 for (
unsigned int d = 0;
d < dim; ++
d)
1016 cell_vertices[vertex_index](
d) = corner_point[
d];
1018 corner_point[
d] = 0;
1026 unsigned int vertex_index = 0;
1028 forest->connectivity,
1036 copy_vertex(vertex_index);
1044 forest->connectivity,
1046 quadrant->x + quad_length_on_level,
1052 copy_vertex(vertex_index);
1059 forest->connectivity,
1062 quadrant->y + quad_length_on_level,
1067 copy_vertex(vertex_index);
1074 forest->connectivity,
1076 quadrant->x + quad_length_on_level,
1077 quadrant->y + quad_length_on_level,
1082 copy_vertex(vertex_index);
1089 forest->connectivity,
1093 quadrant->z + quad_length_on_level,
1097 copy_vertex(vertex_index);
1104 forest->connectivity,
1106 quadrant->x + quad_length_on_level,
1108 quadrant->z + quad_length_on_level,
1112 copy_vertex(vertex_index);
1119 forest->connectivity,
1122 quadrant->y + quad_length_on_level,
1123 quadrant->z + quad_length_on_level,
1127 copy_vertex(vertex_index);
1134 forest->connectivity,
1136 quadrant->x + quad_length_on_level,
1137 quadrant->y + quad_length_on_level,
1138 quadrant->z + quad_length_on_level,
1142 copy_vertex(vertex_index);
1145 are_vertices_initialized =
true;
1155 template <
int dim,
int spacedim>
1156 class RefineAndCoarsenList
1160 const std::vector<types::global_dof_index>
1161 &p4est_tree_to_coarse_cell_permutation,
1189 pointers_are_at_end()
const;
1192 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1193 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1194 const_iterator current_refine_pointer;
1196 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1197 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1198 const_iterator current_coarsen_pointer;
1209 template <
int dim,
int spacedim>
1211 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1213 return ((current_refine_pointer == refine_list.end()) &&
1214 (current_coarsen_pointer == coarsen_list.end()));
1219 template <
int dim,
int spacedim>
1220 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1222 const std::vector<types::global_dof_index>
1223 & p4est_tree_to_coarse_cell_permutation,
1227 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1228 for (
const auto &cell :
triangulation.active_cell_iterators())
1231 if (cell->subdomain_id() != my_subdomain)
1234 if (cell->refine_flag_set())
1236 else if (cell->coarsen_flag_set())
1240 refine_list.reserve(n_refine_flags);
1241 coarsen_list.reserve(n_coarsen_flags);
1253 unsigned int coarse_cell_index =
1254 p4est_tree_to_coarse_cell_permutation[c];
1263 p4est_cell.p.which_tree = c;
1264 build_lists(cell, p4est_cell, my_subdomain);
1272 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1273 Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1275 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1276 Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1279 current_refine_pointer = refine_list.begin();
1280 current_coarsen_pointer = coarsen_list.begin();
1285 template <
int dim,
int spacedim>
1287 RefineAndCoarsenList<dim, spacedim>::build_lists(
1292 if (cell->is_active())
1294 if (cell->subdomain_id() == my_subdomain)
1296 if (cell->refine_flag_set())
1297 refine_list.push_back(p4est_cell);
1298 else if (cell->coarsen_flag_set())
1299 coarsen_list.push_back(p4est_cell);
1306 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1311 P4EST_QUADRANT_INIT(&p4est_child[c]);
1314 P8EST_QUADRANT_INIT(&p4est_child[c]);
1321 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1324 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1325 build_lists(cell->child(c), p4est_child[c], my_subdomain);
1331 template <
int dim,
int spacedim>
1333 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1338 RefineAndCoarsenList<dim, spacedim> *this_object =
1339 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1340 forest->user_pointer);
1344 if (this_object->current_refine_pointer == this_object->refine_list.end())
1347 Assert(coarse_cell_index <=
1348 this_object->current_refine_pointer->p.which_tree,
1353 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1357 Assert(coarse_cell_index <=
1358 this_object->current_refine_pointer->p.which_tree,
1364 quadrant, &*this_object->current_refine_pointer) <= 0,
1369 quadrant, &*this_object->current_refine_pointer))
1371 ++this_object->current_refine_pointer;
1381 template <
int dim,
int spacedim>
1383 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1388 RefineAndCoarsenList<dim, spacedim> *this_object =
1389 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1390 forest->user_pointer);
1394 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1397 Assert(coarse_cell_index <=
1398 this_object->current_coarsen_pointer->p.which_tree,
1403 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1407 Assert(coarse_cell_index <=
1408 this_object->current_coarsen_pointer->p.which_tree,
1414 children[0], &*this_object->current_coarsen_pointer) <= 0,
1420 children[0], &*this_object->current_coarsen_pointer))
1423 ++this_object->current_coarsen_pointer;
1427 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1431 children[c], &*this_object->current_coarsen_pointer),
1433 ++this_object->current_coarsen_pointer;
1451 template <
int dim,
int spacedim>
1452 class PartitionWeights
1460 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1475 std::vector<unsigned int> cell_weights_list;
1476 std::vector<unsigned int>::const_iterator current_pointer;
1480 template <
int dim,
int spacedim>
1481 PartitionWeights<dim, spacedim>::PartitionWeights(
1482 const std::vector<unsigned int> &cell_weights)
1483 : cell_weights_list(cell_weights)
1487 current_pointer = cell_weights_list.begin();
1491 template <
int dim,
int spacedim>
1493 PartitionWeights<dim, spacedim>::cell_weight(
1502 PartitionWeights<dim, spacedim> *this_object =
1503 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1505 Assert(this_object->current_pointer >=
1506 this_object->cell_weights_list.begin(),
1508 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1514 const unsigned int weight = *this_object->current_pointer;
1515 ++this_object->current_pointer;
1518 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1519 "weights for cells. The weight provided here exceeds "
1520 "the maximum value represented as a 'signed int'."));
1521 return static_cast<int>(weight);
1524 template <
int dim,
int spacedim>
1525 using cell_relation_t =
typename std::pair<
1526 typename ::Triangulation<dim, spacedim>::cell_iterator,
1527 typename ::Triangulation<dim, spacedim>::CellStatus>;
1538 template <
int dim,
int spacedim>
1540 add_single_cell_relation(
1541 std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
1542 const typename ::internal::p4est::types<dim>::tree & tree,
1543 const unsigned int idx,
1547 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1553 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1567 template <
int dim,
int spacedim>
1569 update_cell_relations_recursively(
1570 std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
1571 const typename ::internal::p4est::types<dim>::tree & tree,
1573 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1576 const int idx = sc_array_bsearch(
1577 const_cast<sc_array_t *
>(&tree.quadrants),
1582 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1584 &p4est_cell) ==
false))
1589 const bool p4est_has_children = (idx == -1);
1590 if (p4est_has_children && dealii_cell->has_children())
1593 typename ::internal::p4est::types<dim>::quadrant
1596 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1601 P4EST_QUADRANT_INIT(&p4est_child[c]);
1604 P8EST_QUADRANT_INIT(&p4est_child[c]);
1611 &p4est_cell, p4est_child);
1613 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1616 update_cell_relations_recursively<dim, spacedim>(
1617 cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1620 else if (!p4est_has_children && !dealii_cell->has_children())
1624 add_single_cell_relation<dim, spacedim>(
1631 else if (p4est_has_children)
1639 typename ::internal::p4est::types<dim>::quadrant
1641 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1646 P4EST_QUADRANT_INIT(&p4est_child[c]);
1649 P8EST_QUADRANT_INIT(&p4est_child[c]);
1656 &p4est_cell, p4est_child);
1663 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1666 child_idx = sc_array_bsearch(
1667 const_cast<sc_array_t *
>(&tree.quadrants),
1674 add_single_cell_relation<dim, spacedim>(
1675 cell_rel, tree, child_idx, dealii_cell, cell_status);
1683 add_single_cell_relation<dim, spacedim>(
1697 namespace distributed
1700 template <
int dim,
int spacedim>
1702 const MPI_Comm &mpi_communicator,
1703 const typename ::Triangulation<dim, spacedim>::MeshSmoothing
1715 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1719 , triangulation_has_content(false)
1720 , connectivity(nullptr)
1721 , parallel_forest(nullptr)
1728 template <
int dim,
int spacedim>
1748 template <
int dim,
int spacedim>
1761 const typename ::Triangulation<dim, spacedim>::DistortedCellList
1770 this->all_reference_cells_are_hyper_cube(),
1772 "The class parallel::distributed::Triangulation only supports meshes "
1773 "consisting only of hypercube-like cells."));
1778 triangulation_has_content =
true;
1780 setup_coarse_cell_to_p4est_tree_permutation();
1782 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1786 copy_local_forest_to_triangulation();
1795 this->update_periodic_face_map();
1796 this->update_number_cache();
1801 template <
int dim,
int spacedim>
1807 (void)construction_data;
1814 template <
int dim,
int spacedim>
1818 triangulation_has_content =
false;
1820 if (parallel_ghost !=
nullptr)
1824 parallel_ghost =
nullptr;
1827 if (parallel_forest !=
nullptr)
1830 parallel_forest =
nullptr;
1833 if (connectivity !=
nullptr)
1837 connectivity =
nullptr;
1840 coarse_cell_to_p4est_tree_permutation.resize(0);
1841 p4est_tree_to_coarse_cell_permutation.resize(0);
1845 this->update_number_cache();
1850 template <
int dim,
int spacedim>
1860 template <
int dim,
int spacedim>
1870 template <
int dim,
int spacedim>
1873 const typename ::internal::p4est::types<dim>::forest
1875 const typename ::internal::p4est::types<dim>::gloidx
1876 *previous_global_first_quadrant)
1878 Assert(this->data_transfer.sizes_fixed_cumulative.size() > 0,
1882 this->data_transfer.dest_data_fixed.resize(
1883 parallel_forest->local_num_quadrants *
1884 this->data_transfer.sizes_fixed_cumulative.back());
1887 typename ::internal::p4est::types<dim>::transfer_context
1891 parallel_forest->global_first_quadrant,
1892 previous_global_first_quadrant,
1893 parallel_forest->mpicomm,
1895 this->data_transfer.dest_data_fixed.data(),
1896 this->data_transfer.src_data_fixed.data(),
1897 this->data_transfer.sizes_fixed_cumulative.back());
1899 if (this->data_transfer.variable_size_data_stored)
1902 this->data_transfer.dest_sizes_variable.resize(
1903 parallel_forest->local_num_quadrants);
1908 parallel_forest->global_first_quadrant,
1909 previous_global_first_quadrant,
1910 parallel_forest->mpicomm,
1912 this->data_transfer.dest_sizes_variable.data(),
1913 this->data_transfer.src_sizes_variable.data(),
1914 sizeof(
unsigned int));
1920 this->data_transfer.src_data_fixed.clear();
1921 this->data_transfer.src_data_fixed.shrink_to_fit();
1923 if (this->data_transfer.variable_size_data_stored)
1926 this->data_transfer.dest_data_variable.resize(
1927 std::accumulate(this->data_transfer.dest_sizes_variable.begin(),
1928 this->data_transfer.dest_sizes_variable.end(),
1931 # if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1938 if (this->data_transfer.src_sizes_variable.size() == 0)
1939 this->data_transfer.src_sizes_variable.resize(1);
1940 if (this->data_transfer.dest_sizes_variable.size() == 0)
1941 this->data_transfer.dest_sizes_variable.resize(1);
1946 parallel_forest->global_first_quadrant,
1947 previous_global_first_quadrant,
1948 parallel_forest->mpicomm,
1950 this->data_transfer.dest_data_variable.data(),
1951 this->data_transfer.dest_sizes_variable.data(),
1952 this->data_transfer.src_data_variable.data(),
1953 this->data_transfer.src_sizes_variable.data());
1956 this->data_transfer.src_sizes_variable.clear();
1957 this->data_transfer.src_sizes_variable.shrink_to_fit();
1958 this->data_transfer.src_data_variable.clear();
1959 this->data_transfer.src_data_variable.shrink_to_fit();
1965 template <
int dim,
int spacedim>
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.c_str());
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)
2042 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2043 parallel_forest->global_num_quadrants,
2051 this->signals.post_distributed_save();
2056 template <
int dim,
int spacedim>
2063 "load() only works if the Triangulation already contains a coarse mesh!"));
2065 this->n_levels() == 1,
2067 "Triangulation may only contain coarse cells when calling load()."));
2073 this->signals.pre_distributed_load();
2075 if (parallel_ghost !=
nullptr)
2079 parallel_ghost =
nullptr;
2082 parallel_forest =
nullptr;
2085 connectivity =
nullptr;
2087 unsigned int version, numcpus, attached_count_fixed,
2088 attached_count_variable, n_coarse_cells;
2090 std::string fname = std::string(filename) +
".info";
2091 std::ifstream f(fname.c_str());
2093 std::string firstline;
2094 getline(f, firstline);
2095 f >> version >> numcpus >> attached_count_fixed >>
2096 attached_count_variable >> n_coarse_cells;
2100 ExcMessage(
"Incompatible version found in .info file."));
2102 ExcMessage(
"Number of coarse cells differ!"));
2106 this->cell_attached_data.n_attached_data_sets = 0;
2107 this->cell_attached_data.n_attached_deserialize =
2108 attached_count_fixed + attached_count_variable;
2112 this->mpi_communicator,
2130 copy_local_forest_to_triangulation();
2140 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2141 parallel_forest->global_num_quadrants,
2142 parallel_forest->local_num_quadrants,
2144 attached_count_fixed,
2145 attached_count_variable);
2148 this->signals.post_distributed_load();
2150 this->update_periodic_face_map();
2151 this->update_number_cache();
2156 template <
int dim,
int spacedim>
2159 const bool autopartition)
2161 (void)autopartition;
2167 template <
int dim,
int spacedim>
2170 const typename ::internal::p4est::types<dim>::forest *forest)
2174 "load() only works if the Triangulation already contains "
2178 "Coarse mesh of the Triangulation does not match the one "
2179 "of the provided forest!"));
2182 if (parallel_ghost !=
nullptr)
2186 parallel_ghost =
nullptr;
2189 parallel_forest =
nullptr;
2195 typename ::internal::p4est::types<dim>::forest *temp =
2196 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2200 parallel_forest->connectivity = connectivity;
2201 parallel_forest->user_pointer =
this;
2205 copy_local_forest_to_triangulation();
2214 this->update_periodic_face_map();
2215 this->update_number_cache();
2220 template <
int dim,
int spacedim>
2224 Assert(parallel_forest !=
nullptr,
2226 "Can't produce a check sum when no forest is created yet."));
2227 return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2232 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>
2244 typename ::internal::p4est::types<dim>::tree *
2246 const int dealii_coarse_cell_index)
const
2248 const unsigned int tree_index =
2249 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2250 typename ::internal::p4est::types<dim>::tree *tree =
2251 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2252 sc_array_index(parallel_forest->trees, tree_index));
2265 std::integral_constant<int, 2>)
2267 const unsigned int dim = 2, spacedim = 2;
2275 std::vector<unsigned int> vertex_touch_count;
2277 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2280 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2281 const ::internal::p4est::types<2>::locidx num_vtt =
2282 std::accumulate(vertex_touch_count.begin(),
2283 vertex_touch_count.end(),
2289 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2292 (set_vertex_info ==
true ? this->n_vertices() : 0),
2297 set_vertex_and_cell_info(*
this,
2300 coarse_cell_to_p4est_tree_permutation,
2304 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2309 this->mpi_communicator,
2325 std::integral_constant<int, 2>)
2327 const unsigned int dim = 2, spacedim = 3;
2335 std::vector<unsigned int> vertex_touch_count;
2337 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2340 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2341 const ::internal::p4est::types<2>::locidx num_vtt =
2342 std::accumulate(vertex_touch_count.begin(),
2343 vertex_touch_count.end(),
2349 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2352 (set_vertex_info ==
true ? this->n_vertices() : 0),
2357 set_vertex_and_cell_info(*
this,
2360 coarse_cell_to_p4est_tree_permutation,
2364 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2369 this->mpi_communicator,
2383 std::integral_constant<int, 3>)
2385 const int dim = 3, spacedim = 3;
2393 std::vector<unsigned int> vertex_touch_count;
2394 std::vector<std::list<
2395 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2397 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2398 const ::internal::p4est::types<2>::locidx num_vtt =
2399 std::accumulate(vertex_touch_count.begin(),
2400 vertex_touch_count.end(),
2403 std::vector<unsigned int> edge_touch_count;
2404 std::vector<std::list<
2405 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2407 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2408 const ::internal::p4est::types<2>::locidx num_ett =
2409 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2412 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2415 (set_vertex_info ==
true ? this->n_vertices() : 0),
2417 this->n_active_lines(),
2422 set_vertex_and_cell_info(*
this,
2425 coarse_cell_to_p4est_tree_permutation,
2454 const unsigned int deal_to_p4est_line_index[12] = {
2455 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2458 this->begin_active();
2459 cell != this->
end();
2462 const unsigned int index =
2463 coarse_cell_to_p4est_tree_permutation[cell->index()];
2464 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2466 deal_to_p4est_line_index[
e]] =
2467 cell->line(
e)->index();
2472 connectivity->ett_offset[0] = 0;
2473 std::partial_sum(edge_touch_count.begin(),
2474 edge_touch_count.end(),
2475 &connectivity->ett_offset[1]);
2477 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2480 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2482 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2486 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2487 unsigned int>>::const_iterator p =
2488 edge_to_cell[v].begin();
2489 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2491 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2492 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2493 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2494 deal_to_p4est_line_index[p->second];
2498 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2503 this->mpi_communicator,
2520 template <
int dim,
int spacedim>
2522 enforce_mesh_balance_over_periodic_boundaries(
2528 std::vector<bool> flags_before[2];
2532 std::vector<unsigned int> topological_vertex_numbering(
2534 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2535 topological_vertex_numbering[i] = i;
2553 using cell_iterator =
2555 typename std::map<std::pair<cell_iterator, unsigned int>,
2556 std::pair<std::pair<cell_iterator, unsigned int>,
2557 std::bitset<3>>>::const_iterator it;
2562 const cell_iterator &cell_1 = it->first.first;
2563 const unsigned int face_no_1 = it->first.second;
2564 const cell_iterator &cell_2 = it->second.first.first;
2565 const unsigned int face_no_2 = it->second.first.second;
2566 const std::bitset<3> face_orientation = it->second.second;
2568 if (cell_1->level() == cell_2->level())
2570 for (
unsigned int v = 0;
2576 const unsigned int vface0 =
2579 face_orientation[0],
2580 face_orientation[1],
2581 face_orientation[2]);
2582 const unsigned int vi0 =
2583 topological_vertex_numbering[cell_1->face(face_no_1)
2584 ->vertex_index(vface0)];
2585 const unsigned int vi1 =
2586 topological_vertex_numbering[cell_2->face(face_no_2)
2588 const unsigned int min_index =
std::min(vi0, vi1);
2589 topological_vertex_numbering[cell_1->face(face_no_1)
2590 ->vertex_index(vface0)] =
2591 topological_vertex_numbering[cell_2->face(face_no_2)
2592 ->vertex_index(v)] =
2599 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2601 const unsigned int j = topological_vertex_numbering[i];
2609 bool continue_iterating =
true;
2611 while (continue_iterating)
2615 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2619 for (; cell != endc; ++cell)
2621 if (cell->refine_flag_set())
2622 for (
const unsigned int vertex :
2624 vertex_level[topological_vertex_numbering
2625 [cell->vertex_index(vertex)]] =
2626 std::max(vertex_level[topological_vertex_numbering
2627 [cell->vertex_index(vertex)]],
2629 else if (!cell->coarsen_flag_set())
2630 for (
const unsigned int vertex :
2632 vertex_level[topological_vertex_numbering
2633 [cell->vertex_index(vertex)]] =
2634 std::max(vertex_level[topological_vertex_numbering
2635 [cell->vertex_index(vertex)]],
2646 for (
const unsigned int vertex :
2648 vertex_level[topological_vertex_numbering
2649 [cell->vertex_index(vertex)]] =
2650 std::max(vertex_level[topological_vertex_numbering
2651 [cell->vertex_index(vertex)]],
2656 continue_iterating =
false;
2668 if (cell->refine_flag_set() ==
false)
2670 for (
const unsigned int vertex :
2672 if (vertex_level[topological_vertex_numbering
2673 [cell->vertex_index(vertex)]] >=
2677 cell->clear_coarsen_flag();
2682 if (vertex_level[topological_vertex_numbering
2683 [cell->vertex_index(vertex)]] >
2686 cell->set_refine_flag();
2687 continue_iterating =
true;
2689 for (
const unsigned int v :
2691 vertex_level[topological_vertex_numbering
2692 [cell->vertex_index(v)]] =
2694 vertex_level[topological_vertex_numbering
2695 [cell->vertex_index(v)]],
2709 if (cell->is_active())
2712 const unsigned int n_children = cell->n_children();
2713 unsigned int flagged_children = 0;
2714 for (
unsigned int child = 0; child < n_children; ++child)
2715 if (cell->child(child)->is_active() &&
2716 cell->child(child)->coarsen_flag_set())
2721 if (flagged_children < n_children)
2722 for (
unsigned int child = 0; child < n_children; ++child)
2723 if (cell->child(child)->is_active())
2724 cell->child(child)->clear_coarsen_flag();
2727 std::vector<bool> flags_after[2];
2730 return ((flags_before[0] != flags_after[0]) ||
2731 (flags_before[1] != flags_after[1]));
2737 template <
int dim,
int spacedim>
2741 bool mesh_changed =
false;
2742 unsigned int loop_counter = 0;
2743 unsigned int n_changes = 0;
2748 this->update_periodic_face_map();
2750 mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*
this);
2751 n_changes += mesh_changed;
2762 "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2763 "for periodic boundaries detected. Aborting."));
2765 while (mesh_changed);
2768 return n_changes > 0;
2773 template <
int dim,
int spacedim>
2797 spacedim>::limit_level_difference_at_vertices;
2801 bool mesh_changed =
false;
2809 if (
settings & mesh_reconstruction_after_repartitioning)
2810 while (this->n_levels() > 1)
2817 for (
const auto &cell :
2818 this->active_cell_iterators_on_level(this->n_levels() - 1))
2820 cell->set_coarsen_flag();
2838 if (parallel_ghost !=
nullptr)
2842 parallel_ghost =
nullptr;
2846 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2847 P4EST_CONNECT_CORNER) :
2848 typename ::internal::p4est::types<dim>::balance_type(
2849 P8EST_CONNECT_CORNER)));
2856 for (
const auto &cell : this->cell_iterators_on_level(0))
2861 for (
const auto &cell : this->cell_iterators_on_level(0))
2866 if (tree_exists_locally<dim, spacedim>(
2868 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2871 delete_all_children<dim, spacedim>(cell);
2872 if (cell->is_active())
2881 typename ::internal::p4est::types<dim>::quadrant
2883 typename ::internal::p4est::types<dim>::tree *tree =
2884 init_tree(cell->index());
2886 ::internal::p4est::init_coarse_quadrant<dim>(
2889 match_tree_recursively<dim, spacedim>(*tree,
2893 this->my_subdomain);
2900 typename ::internal::p4est::types<dim>::quadrant *quadr;
2902 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2904 for (
unsigned int g_idx = 0;
2905 g_idx < parallel_ghost->ghosts.elem_count;
2908 while (g_idx >=
static_cast<unsigned int>(
2909 parallel_ghost->proc_offsets[ghost_owner + 1]))
2911 while (g_idx >=
static_cast<unsigned int>(
2912 parallel_ghost->tree_offsets[ghost_tree + 1]))
2915 quadr =
static_cast<
2916 typename ::internal::p4est::types<dim>::quadrant *
>(
2917 sc_array_index(¶llel_ghost->ghosts, g_idx));
2919 unsigned int coarse_cell_index =
2920 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2922 match_quadrant<dim, spacedim>(
this,
2929 this->prepare_coarsening_and_refinement();
2933 std::any_of(this->begin_active(),
2936 return cell.refine_flag_set() ||
2937 cell.coarsen_flag_set();
2955 while (mesh_changed);
2959 unsigned int num_ghosts = 0;
2961 for (
const auto &cell : this->active_cell_iterators())
2963 if (cell->subdomain_id() != this->my_subdomain &&
2968 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2985 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
2988 for (
const auto &cell : this->cell_iterators_on_level(lvl))
2990 if ((cell->is_active() &&
2991 cell->subdomain_id() ==
2993 (cell->has_children() &&
2994 cell->child(0)->level_subdomain_id() ==
2996 cell->set_level_subdomain_id(
2997 this->locally_owned_subdomain());
3001 cell->set_level_subdomain_id(
3009 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3010 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3011 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3013 for (
const auto &cell : this->cell_iterators_on_level(0))
3015 typename ::internal::p4est::types<dim>::quadrant
3017 const unsigned int tree_index =
3018 coarse_cell_to_p4est_tree_permutation[cell->index()];
3019 typename ::internal::p4est::types<dim>::tree *tree =
3020 init_tree(cell->index());
3022 ::internal::p4est::init_coarse_quadrant<dim>(
3025 determine_level_subdomain_id_recursively<dim, spacedim>(
3036 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3039 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3041 if (cell->has_children())
3042 for (
unsigned int c = 0;
3043 c < GeometryInfo<dim>::max_children_per_cell;
3046 if (cell->child(c)->level_subdomain_id() ==
3052 cell->child(0)->level_subdomain_id();
3056 cell->set_level_subdomain_id(mark);
3073 (void)total_local_cells;
3077 Assert(
static_cast<unsigned int>(
3078 parallel_forest->local_num_quadrants) == total_local_cells,
3083 Assert(
static_cast<unsigned int>(
3084 parallel_forest->local_num_quadrants) <= total_local_cells,
3089 unsigned int n_owned = 0;
3090 for (
const auto &cell : this->active_cell_iterators())
3092 if (cell->subdomain_id() == this->my_subdomain)
3096 Assert(
static_cast<unsigned int>(
3097 parallel_forest->local_num_quadrants) == n_owned,
3101 this->smooth_grid = save_smooth;
3107 update_cell_relations();
3112 template <
int dim,
int spacedim>
3117 std::vector<Point<dim>>
point{p};
3118 std::vector<types::subdomain_id> owner = find_point_owner_rank(
point);
3125 template <
int dim,
int spacedim>
3126 std::vector<types::subdomain_id>
3130 # ifndef P4EST_SEARCH_LOCAL
3135 "This function is only available if p4est is version 2.2 and higher."));
3137 return std::vector<unsigned int>(1,
3141 AssertThrow(this->are_vertices_communicated_to_p4est(),
3143 "Vertices need to be communicated to p4est to use this "
3144 "function. This must explicitly be turned on in the "
3145 "settings of the triangulation's constructor."));
3148 for (
const auto &
manifold_id : this->get_manifold_ids())
3153 "This function can only be used if the triangulation "
3154 "has no other manifold than a Cartesian (flat) manifold attached."));
3158 PartitionSearch<dim> partition_search;
3164 parallel_forest->user_pointer = &partition_search;
3170 sc_array_t *point_sc_array;
3174 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3177 for (
size_t i = 0; i < points.size(); ++i)
3182 double *this_sc_point =
3183 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3185 for (
unsigned int d = 0;
d < dim; ++
d)
3187 this_sc_point[
d] = p(
d);
3189 this_sc_point[dim] = -1.0;
3195 static_cast<int>(
false),
3196 &PartitionSearch<dim>::local_quadrant_fn,
3197 &PartitionSearch<dim>::local_point_fn,
3201 std::vector<types::subdomain_id> owner_rank(
3205 for (
size_t i = 0; i < points.size(); ++i)
3208 double *this_sc_point =
3209 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3214 parallel_forest->user_pointer =
this;
3217 sc_array_destroy_null(&point_sc_array);
3225 template <
int dim,
int spacedim>
3231 for (
const auto &cell : this->active_cell_iterators())
3232 if (cell->is_locally_owned() && cell->refine_flag_set())
3233 Assert(cell->refine_flag_set() ==
3236 "This class does not support anisotropic refinement"));
3241 if (this->n_levels() ==
3245 cell = this->begin_active(
3253 !(cell->refine_flag_set()),
3255 "Fatal Error: maximum refinement level of p4est reached."));
3259 this->prepare_coarsening_and_refinement();
3262 this->signals.pre_distributed_refinement();
3267 for (
const auto &cell : this->active_cell_iterators())
3268 if (cell->is_ghost() || cell->is_artificial())
3270 cell->clear_refine_flag();
3271 cell->clear_coarsen_flag();
3277 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3278 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3285 parallel_forest->user_pointer = &refine_and_coarsen_list;
3287 if (parallel_ghost !=
nullptr)
3291 parallel_ghost =
nullptr;
3296 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3301 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3308 parallel_forest->user_pointer =
this;
3314 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3315 P4EST_CONNECT_FULL) :
3316 typename ::internal::p4est::types<dim>::balance_type(
3317 P8EST_CONNECT_FULL)),
3322 update_cell_relations();
3326 this->signals.post_p4est_refinement();
3330 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3331 previous_global_first_quadrant;
3333 if (this->cell_attached_data.n_attached_data_sets > 0)
3335 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3336 std::memcpy(previous_global_first_quadrant.data(),
3337 parallel_forest->global_first_quadrant,
3339 typename ::internal::p4est::types<dim>::gloidx) *
3340 (parallel_forest->mpisize + 1));
3343 if (!(
settings & no_automatic_repartitioning))
3347 if (this->signals.weight.empty())
3355 const std::vector<unsigned int> cell_weights = get_cell_weights();
3361 this->mpi_communicator) > 0,
3363 "The global sum of weights over all active cells "
3364 "is zero. Please verify how you generate weights."));
3366 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3371 parallel_forest->user_pointer = &partition_weights;
3377 &PartitionWeights<dim, spacedim>::cell_weight);
3381 parallel_forest, 0,
nullptr,
nullptr);
3383 parallel_forest->user_pointer =
this;
3388 if (this->cell_attached_data.n_attached_data_sets > 0)
3390 this->data_transfer.pack_data(
3391 this->local_cell_relations,
3392 this->cell_attached_data.pack_callbacks_fixed,
3393 this->cell_attached_data.pack_callbacks_variable);
3399 for (
const auto &cell : this->active_cell_iterators())
3401 cell->clear_refine_flag();
3402 cell->clear_coarsen_flag();
3407 copy_local_forest_to_triangulation();
3417 if (this->cell_attached_data.n_attached_data_sets > 0)
3419 this->execute_transfer(parallel_forest,
3420 previous_global_first_quadrant.data());
3423 this->data_transfer.unpack_cell_status(this->local_cell_relations);
3444 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3446 std::vector<bool> active_verts =
3447 this->mark_locally_active_vertices_on_level(lvl);
3449 const unsigned int maybe_coarser_lvl =
3450 (lvl > 0) ? (lvl - 1) : lvl;
3452 cell = this->
begin(maybe_coarser_lvl),
3453 endc = this->
end(lvl);
3454 for (; cell != endc; ++cell)
3455 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3457 const bool is_level_artificial =
3458 (cell->level_subdomain_id() ==
3460 bool need_to_know =
false;
3461 for (
const unsigned int vertex :
3463 if (active_verts[cell->vertex_index(vertex)])
3465 need_to_know =
true;
3470 !need_to_know || !is_level_artificial,
3472 "Internal error: the owner of cell" +
3473 cell->id().to_string() +
3474 " is unknown even though it is needed for geometric multigrid."));
3480 this->update_periodic_face_map();
3481 this->update_number_cache();
3484 this->signals.post_distributed_refinement();
3489 template <
int dim,
int spacedim>
3494 for (
const auto &cell : this->active_cell_iterators())
3495 if (cell->is_locally_owned())
3497 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3499 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3503 this->signals.pre_distributed_repartition();
3507 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3508 previous_global_first_quadrant;
3510 if (this->cell_attached_data.n_attached_data_sets > 0)
3512 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3513 std::memcpy(previous_global_first_quadrant.data(),
3514 parallel_forest->global_first_quadrant,
3516 typename ::internal::p4est::types<dim>::gloidx) *
3517 (parallel_forest->mpisize + 1));
3520 if (this->signals.weight.empty())
3532 const std::vector<unsigned int> cell_weights = get_cell_weights();
3538 this->mpi_communicator) > 0,
3540 "The global sum of weights over all active cells "
3541 "is zero. Please verify how you generate weights."));
3543 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3548 parallel_forest->user_pointer = &partition_weights;
3554 &PartitionWeights<dim, spacedim>::cell_weight);
3557 parallel_forest->user_pointer =
this;
3561 if (this->cell_attached_data.n_attached_data_sets > 0)
3563 this->data_transfer.pack_data(
3564 this->local_cell_relations,
3565 this->cell_attached_data.pack_callbacks_fixed,
3566 this->cell_attached_data.pack_callbacks_variable);
3571 copy_local_forest_to_triangulation();
3581 if (this->cell_attached_data.n_attached_data_sets > 0)
3583 this->execute_transfer(parallel_forest,
3584 previous_global_first_quadrant.data());
3587 this->update_periodic_face_map();
3590 this->update_number_cache();
3593 this->signals.post_distributed_repartition();
3598 template <
int dim,
int spacedim>
3599 const std::vector<types::global_dof_index> &
3603 return p4est_tree_to_coarse_cell_permutation;
3608 template <
int dim,
int spacedim>
3609 const std::vector<types::global_dof_index> &
3613 return coarse_cell_to_p4est_tree_permutation;
3618 template <
int dim,
int spacedim>
3621 const int level)
const
3625 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3626 for (
const auto &cell : this->cell_iterators_on_level(
level))
3629 marked_vertices[cell->vertex_index(v)] =
true;
3646 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3647 for (
const auto &it : this->get_periodic_face_map())
3650 const unsigned int face_no_1 = it.first.second;
3652 const unsigned int face_no_2 = it.second.first.second;
3653 const std::bitset<3> &face_orientation = it.second.second;
3655 if (cell_1->level() ==
level && cell_2->level() ==
level)
3657 for (
unsigned int v = 0;
3663 const unsigned int vface0 =
3666 face_orientation[0],
3667 face_orientation[1],
3668 face_orientation[2]);
3669 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3671 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3673 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3675 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3681 return marked_vertices;
3686 template <
int dim,
int spacedim>
3696 template <
int dim,
int spacedim>
3699 const unsigned int coarse_cell_index)
const
3701 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3706 template <
int dim,
int spacedim>
3710 &periodicity_vector)
3712 Assert(triangulation_has_content ==
true,
3714 Assert(this->n_levels() == 1,
3715 ExcMessage(
"The triangulation is refined!"));
3722 for (
const auto &face_pair : periodicity_vector)
3726 const unsigned int face_left = face_pair.face_idx[0];
3727 const unsigned int face_right = face_pair.face_idx[1];
3730 const unsigned int tree_left =
3731 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3732 const unsigned int tree_right =
3733 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3742 unsigned int p4est_orientation = 0;
3744 p4est_orientation = face_pair.orientation[1];
3747 const unsigned int face_idx_list[] = {face_left, face_right};
3748 const cell_iterator cell_list[] = {first_cell, second_cell};
3749 unsigned int lower_idx, higher_idx;
3750 if (face_left <= face_right)
3763 unsigned int first_p4est_idx_on_cell =
3764 p8est_face_corners[face_idx_list[lower_idx]][0];
3765 unsigned int first_dealii_idx_on_face =
3767 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3770 const unsigned int first_dealii_idx_on_cell =
3772 face_idx_list[lower_idx],
3774 cell_list[lower_idx]->face_orientation(
3775 face_idx_list[lower_idx]),
3776 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3777 cell_list[lower_idx]->face_rotation(
3778 face_idx_list[lower_idx]));
3779 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3781 first_dealii_idx_on_face = i;
3788 constexpr
unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
3796 constexpr
unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
3804 const unsigned int second_dealii_idx_on_face =
3805 lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
3806 [first_dealii_idx_on_face] :
3807 right_to_left[face_pair.orientation.to_ulong()]
3808 [first_dealii_idx_on_face];
3809 const unsigned int second_dealii_idx_on_cell =
3811 face_idx_list[higher_idx],
3812 second_dealii_idx_on_face,
3813 cell_list[higher_idx]->face_orientation(
3814 face_idx_list[higher_idx]),
3815 cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
3816 cell_list[higher_idx]->face_rotation(
3817 face_idx_list[higher_idx]));
3819 const unsigned int second_p4est_idx_on_face =
3820 p8est_corner_face_corners[second_dealii_idx_on_cell]
3821 [face_idx_list[higher_idx]];
3822 p4est_orientation = second_p4est_idx_on_face;
3842 this->mpi_communicator,
3853 copy_local_forest_to_triangulation();
3864 this->update_number_cache();
3869 template <
int dim,
int spacedim>
3880 this->cell_attached_data.n_attached_data_sets) +
3887 coarse_cell_to_p4est_tree_permutation) +
3889 p4est_tree_to_coarse_cell_permutation) +
3890 memory_consumption_p4est();
3897 template <
int dim,
int spacedim>
3901 return ::internal::p4est::functions<dim>::forest_memory_used(
3909 template <
int dim,
int spacedim>
3912 const ::Triangulation<dim, spacedim> &other_tria)
3916 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
3918 (other_tria.n_global_levels() == 1),
3929 const typename ::Triangulation<dim, spacedim>::DistortedCellList
3937 if (const ::parallel::distributed::Triangulation<dim, spacedim>
3938 *other_distributed =
3939 dynamic_cast<const ::parallel::distributed::
3940 Triangulation<dim, spacedim> *
>(&other_tria))
3943 settings = other_distributed->settings;
3944 triangulation_has_content =
3945 other_distributed->triangulation_has_content;
3946 coarse_cell_to_p4est_tree_permutation =
3947 other_distributed->coarse_cell_to_p4est_tree_permutation;
3948 p4est_tree_to_coarse_cell_permutation =
3949 other_distributed->p4est_tree_to_coarse_cell_permutation;
3952 typename ::internal::p4est::types<dim>::connectivity
3953 *temp_connectivity =
const_cast<
3954 typename ::internal::p4est::types<dim>::connectivity *
>(
3955 other_distributed->connectivity);
3957 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
3960 typename ::internal::p4est::types<dim>::forest *temp_forest =
3961 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
3962 other_distributed->parallel_forest);
3966 parallel_forest->connectivity = connectivity;
3967 parallel_forest->user_pointer =
this;
3971 triangulation_has_content =
true;
3972 setup_coarse_cell_to_p4est_tree_permutation();
3973 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
3978 copy_local_forest_to_triangulation();
3987 this->update_periodic_face_map();
3988 this->update_number_cache();
3993 template <
int dim,
int spacedim>
3998 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
3999 this->local_cell_relations.shrink_to_fit();
4002 for (
const auto &cell : this->cell_iterators_on_level(0))
4005 if (tree_exists_locally<dim, spacedim>(
4007 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4011 typename ::internal::p4est::types<dim>::quadrant
4013 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4016 typename ::internal::p4est::types<dim>::tree *tree =
4017 init_tree(cell->index());
4019 update_cell_relations_recursively<dim, spacedim>(
4020 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4026 template <
int dim,
int spacedim>
4027 std::vector<unsigned int>
4032 Assert(this->local_cell_relations.size() ==
4033 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4038 std::vector<unsigned int> weights;
4039 weights.reserve(this->local_cell_relations.size());
4047 for (
const auto &cell_rel : this->local_cell_relations)
4049 const auto &cell_it = cell_rel.first;
4050 const auto &cell_status = cell_rel.second;
4052 weights.push_back(this->signals.weight(cell_it, cell_status));
4060 template <
int spacedim>
4062 const MPI_Comm &mpi_communicator,
4063 const typename ::Triangulation<1, spacedim>::MeshSmoothing
4075 template <
int spacedim>
4083 template <
int spacedim>
4084 const std::vector<types::global_dof_index> &
4088 static std::vector<types::global_dof_index> a;
4094 template <
int spacedim>
4095 std::map<unsigned int, std::set<::types::subdomain_id>>
4097 const unsigned int )
const
4101 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4106 template <
int spacedim>
4109 const unsigned int)
const
4112 return std::vector<bool>();
4117 template <
int spacedim>
4128 template <
int spacedim>
4131 const unsigned int)
const
4139 template <
int spacedim>
4148 template <
int spacedim>
4157 template <
int spacedim>
4166 template <
int spacedim>
4176 template <
int spacedim>
4186 template <
int spacedim>
4203 namespace distributed
4205 template <
int dim,
int spacedim>
4213 #ifdef DEAL_II_WITH_P4EST
4223 const auto &cell = pair.first;
4224 const auto &status = pair.second;
4228 case ::Triangulation<dim, spacedim>::CELL_PERSIST:
4230 cell->clear_refine_flag();
4231 cell->clear_coarsen_flag();
4234 case ::Triangulation<dim, spacedim>::CELL_REFINE:
4236 cell->clear_coarsen_flag();
4237 cell->set_refine_flag();
4240 case ::Triangulation<dim, spacedim>::CELL_COARSEN:
4242 for (
const auto &child : cell->child_iterators())
4244 child->clear_refine_flag();
4245 child->set_coarsen_flag();
4249 case ::Triangulation<dim, spacedim>::CELL_INVALID:
4264 template <
int dim,
int spacedim>
4267 #ifdef DEAL_II_WITH_P4EST
4268 if (distributed_tria)
4271 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4272 distributed_tria->load_refine_flags(saved_refine_flags);
4276 (void)distributed_tria;
4285 #include "tria.inst"
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 void execute_coarsening_and_refinement() override
void setup_coarse_cell_to_p4est_tree_permutation()
bool are_vertices_communicated_to_p4est() const
virtual types::coarse_cell_id coarse_cell_index_to_coarse_cell_id(const unsigned int coarse_cell_index) const override
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &other_tria) override
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
typename ::internal::p4est::types< dim >::ghost * parallel_ghost
Triangulation(const MPI_Comm &mpi_communicator, const typename ::Triangulation< dim, spacedim >::MeshSmoothing smooth_grid=(::Triangulation< dim, spacedim >::none), const Settings settings=default_setting)
std::vector< unsigned int > get_cell_weights() const
unsigned int get_checksum() const
virtual ~Triangulation() override
virtual void update_cell_relations() override
void execute_transfer(const typename ::internal::p4est::types< dim >::forest *parallel_forest, const typename ::internal::p4est::types< dim >::gloidx *previous_global_first_quadrant)
typename ::internal::p4est::types< dim >::tree * init_tree(const int dealii_coarse_cell_index) const
types::subdomain_id find_point_owner_rank(const Point< dim > &p)
virtual bool prepare_coarsening_and_refinement() override
std::vector< bool > mark_locally_active_vertices_on_level(const int level) const
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator >> &) override
const ::internal::p4est::types< dim >::forest * get_p4est() const
const std::vector< types::global_dof_index > & get_coarse_cell_to_p4est_tree_permutation() const
void write_mesh_vtk(const std::string &file_basename) const
void copy_local_forest_to_triangulation()
virtual unsigned int coarse_cell_id_to_coarse_cell_index(const types::coarse_cell_id coarse_cell_id) const override
void copy_new_triangulation_to_p4est(std::integral_constant< int, 2 >)
virtual void save(const std::string &filename) const override
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata) override
virtual std::size_t memory_consumption_p4est() const
virtual void clear() override
bool is_multilevel_hierarchy_constructed() const override
virtual void load(const std::string &filename) override
virtual std::size_t memory_consumption() const override
#define DEAL_II_NAMESPACE_OPEN
#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 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< T >::value, 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)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(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