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,
1525 typename ::Triangulation<dim, spacedim>::CellStatus>;
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>(
1629 else if (p4est_has_children)
1637 typename ::internal::p4est::types<dim>::quadrant
1639 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1644 P4EST_QUADRANT_INIT(&p4est_child[c]);
1647 P8EST_QUADRANT_INIT(&p4est_child[c]);
1654 &p4est_cell, p4est_child);
1661 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1664 child_idx = sc_array_bsearch(
1665 const_cast<sc_array_t *
>(&tree.quadrants),
1672 add_single_cell_relation<dim, spacedim>(
1673 cell_rel, tree, child_idx, dealii_cell, cell_status);
1681 add_single_cell_relation<dim, spacedim>(
1695 namespace distributed
1698 template <
int dim,
int spacedim>
1700 const MPI_Comm &mpi_communicator,
1701 const typename ::Triangulation<dim, spacedim>::MeshSmoothing
1713 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1717 , triangulation_has_content(false)
1718 , connectivity(nullptr)
1719 , parallel_forest(nullptr)
1726 template <
int dim,
int spacedim>
1746 template <
int dim,
int spacedim>
1759 const typename ::Triangulation<dim, spacedim>::DistortedCellList
1768 this->all_reference_cells_are_hyper_cube(),
1770 "The class parallel::distributed::Triangulation only supports meshes "
1771 "consisting only of hypercube-like cells."));
1776 triangulation_has_content =
true;
1778 setup_coarse_cell_to_p4est_tree_permutation();
1780 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
1784 copy_local_forest_to_triangulation();
1793 this->update_periodic_face_map();
1794 this->update_number_cache();
1799 template <
int dim,
int spacedim>
1805 (void)construction_data;
1812 template <
int dim,
int spacedim>
1816 triangulation_has_content =
false;
1818 if (parallel_ghost !=
nullptr)
1822 parallel_ghost =
nullptr;
1825 if (parallel_forest !=
nullptr)
1828 parallel_forest =
nullptr;
1831 if (connectivity !=
nullptr)
1835 connectivity =
nullptr;
1838 coarse_cell_to_p4est_tree_permutation.resize(0);
1839 p4est_tree_to_coarse_cell_permutation.resize(0);
1843 this->update_number_cache();
1848 template <
int dim,
int spacedim>
1858 template <
int dim,
int spacedim>
1868 template <
int dim,
int spacedim>
1871 const typename ::internal::p4est::types<dim>::forest
1873 const typename ::internal::p4est::types<dim>::gloidx
1874 *previous_global_first_quadrant)
1876 Assert(this->data_transfer.sizes_fixed_cumulative.size() > 0,
1880 this->data_transfer.dest_data_fixed.resize(
1881 parallel_forest->local_num_quadrants *
1882 this->data_transfer.sizes_fixed_cumulative.back());
1885 typename ::internal::p4est::types<dim>::transfer_context
1889 parallel_forest->global_first_quadrant,
1890 previous_global_first_quadrant,
1891 parallel_forest->mpicomm,
1893 this->data_transfer.dest_data_fixed.data(),
1894 this->data_transfer.src_data_fixed.data(),
1895 this->data_transfer.sizes_fixed_cumulative.back());
1897 if (this->data_transfer.variable_size_data_stored)
1900 this->data_transfer.dest_sizes_variable.resize(
1901 parallel_forest->local_num_quadrants);
1906 parallel_forest->global_first_quadrant,
1907 previous_global_first_quadrant,
1908 parallel_forest->mpicomm,
1910 this->data_transfer.dest_sizes_variable.data(),
1911 this->data_transfer.src_sizes_variable.data(),
1912 sizeof(
unsigned int));
1918 this->data_transfer.src_data_fixed.clear();
1919 this->data_transfer.src_data_fixed.shrink_to_fit();
1921 if (this->data_transfer.variable_size_data_stored)
1924 this->data_transfer.dest_data_variable.resize(
1925 std::accumulate(this->data_transfer.dest_sizes_variable.begin(),
1926 this->data_transfer.dest_sizes_variable.end(),
1929 # if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1936 if (this->data_transfer.src_sizes_variable.size() == 0)
1937 this->data_transfer.src_sizes_variable.resize(1);
1938 if (this->data_transfer.dest_sizes_variable.size() == 0)
1939 this->data_transfer.dest_sizes_variable.resize(1);
1944 parallel_forest->global_first_quadrant,
1945 previous_global_first_quadrant,
1946 parallel_forest->mpicomm,
1948 this->data_transfer.dest_data_variable.data(),
1949 this->data_transfer.dest_sizes_variable.data(),
1950 this->data_transfer.src_data_variable.data(),
1951 this->data_transfer.src_sizes_variable.data());
1954 this->data_transfer.src_sizes_variable.clear();
1955 this->data_transfer.src_sizes_variable.shrink_to_fit();
1956 this->data_transfer.src_data_variable.clear();
1957 this->data_transfer.src_data_variable.shrink_to_fit();
1963 template <
int dim,
int spacedim>
1970 coarse_cell_to_p4est_tree_permutation.resize(this->
n_cells(0));
1972 cell_connectivity, coarse_cell_to_p4est_tree_permutation);
1974 p4est_tree_to_coarse_cell_permutation =
1980 template <
int dim,
int spacedim>
1983 const std::string &file_basename)
const
1985 Assert(parallel_forest !=
nullptr,
1986 ExcMessage(
"Can't produce output when no forest is created yet."));
1990 "To use this function the triangulation's flag "
1991 "Settings::communicate_vertices_to_p4est must be set."));
1994 parallel_forest,
nullptr, file_basename.c_str());
1999 template <
int dim,
int spacedim>
2004 this->cell_attached_data.n_attached_deserialize == 0,
2006 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2008 ExcMessage(
"Can not save() an empty Triangulation."));
2014 this->signals.pre_distributed_save();
2016 if (this->my_subdomain == 0)
2018 std::string fname = std::string(filename) +
".info";
2019 std::ofstream f(fname.c_str());
2020 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2024 << this->cell_attached_data.pack_callbacks_fixed.size() <<
" "
2025 << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
2026 << this->
n_cells(0) << std::endl;
2030 for (
const auto &cell_rel : this->local_cell_relations)
2040 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2041 parallel_forest->global_num_quadrants,
2049 this->signals.post_distributed_save();
2054 template <
int dim,
int spacedim>
2061 "load() only works if the Triangulation already contains a coarse mesh!"));
2063 this->n_levels() == 1,
2065 "Triangulation may only contain coarse cells when calling load()."));
2071 this->signals.pre_distributed_load();
2073 if (parallel_ghost !=
nullptr)
2077 parallel_ghost =
nullptr;
2080 parallel_forest =
nullptr;
2083 connectivity =
nullptr;
2085 unsigned int version, numcpus, attached_count_fixed,
2086 attached_count_variable, n_coarse_cells;
2088 std::string fname = std::string(filename) +
".info";
2089 std::ifstream f(fname.c_str());
2091 std::string firstline;
2092 getline(f, firstline);
2093 f >> version >> numcpus >> attached_count_fixed >>
2094 attached_count_variable >> n_coarse_cells;
2098 ExcMessage(
"Incompatible version found in .info file."));
2100 ExcMessage(
"Number of coarse cells differ!"));
2104 this->cell_attached_data.n_attached_data_sets = 0;
2105 this->cell_attached_data.n_attached_deserialize =
2106 attached_count_fixed + attached_count_variable;
2110 this->mpi_communicator,
2128 copy_local_forest_to_triangulation();
2138 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2139 parallel_forest->global_num_quadrants,
2140 parallel_forest->local_num_quadrants,
2142 attached_count_fixed,
2143 attached_count_variable);
2146 this->signals.post_distributed_load();
2148 this->update_periodic_face_map();
2149 this->update_number_cache();
2154 template <
int dim,
int spacedim>
2157 const bool autopartition)
2159 (void)autopartition;
2165 template <
int dim,
int spacedim>
2168 const typename ::internal::p4est::types<dim>::forest *forest)
2172 "load() only works if the Triangulation already contains "
2176 "Coarse mesh of the Triangulation does not match the one "
2177 "of the provided forest!"));
2180 if (parallel_ghost !=
nullptr)
2184 parallel_ghost =
nullptr;
2187 parallel_forest =
nullptr;
2193 typename ::internal::p4est::types<dim>::forest *temp =
2194 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2198 parallel_forest->connectivity = connectivity;
2199 parallel_forest->user_pointer =
this;
2203 copy_local_forest_to_triangulation();
2212 this->update_periodic_face_map();
2213 this->update_number_cache();
2218 template <
int dim,
int spacedim>
2222 Assert(parallel_forest !=
nullptr,
2224 "Can't produce a check sum when no forest is created yet."));
2225 return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2230 template <
int dim,
int spacedim>
2231 const typename ::internal::p4est::types<dim>::forest *
2234 Assert(parallel_forest !=
nullptr,
2235 ExcMessage(
"The forest has not been allocated yet."));
2236 return parallel_forest;
2241 template <
int dim,
int spacedim>
2242 typename ::internal::p4est::types<dim>::tree *
2244 const int dealii_coarse_cell_index)
const
2246 const unsigned int tree_index =
2247 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2248 typename ::internal::p4est::types<dim>::tree *tree =
2249 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2250 sc_array_index(parallel_forest->trees, tree_index));
2263 std::integral_constant<int, 2>)
2265 const unsigned int dim = 2, spacedim = 2;
2273 std::vector<unsigned int> vertex_touch_count;
2275 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2278 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2279 const ::internal::p4est::types<2>::locidx num_vtt =
2280 std::accumulate(vertex_touch_count.begin(),
2281 vertex_touch_count.end(),
2287 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2290 (set_vertex_info ==
true ? this->n_vertices() : 0),
2295 set_vertex_and_cell_info(*
this,
2298 coarse_cell_to_p4est_tree_permutation,
2302 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2307 this->mpi_communicator,
2323 std::integral_constant<int, 2>)
2325 const unsigned int dim = 2, spacedim = 3;
2333 std::vector<unsigned int> vertex_touch_count;
2335 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2338 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2339 const ::internal::p4est::types<2>::locidx num_vtt =
2340 std::accumulate(vertex_touch_count.begin(),
2341 vertex_touch_count.end(),
2347 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2350 (set_vertex_info ==
true ? this->n_vertices() : 0),
2355 set_vertex_and_cell_info(*
this,
2358 coarse_cell_to_p4est_tree_permutation,
2362 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2367 this->mpi_communicator,
2381 std::integral_constant<int, 3>)
2383 const int dim = 3, spacedim = 3;
2391 std::vector<unsigned int> vertex_touch_count;
2392 std::vector<std::list<
2393 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2395 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2396 const ::internal::p4est::types<2>::locidx num_vtt =
2397 std::accumulate(vertex_touch_count.begin(),
2398 vertex_touch_count.end(),
2401 std::vector<unsigned int> edge_touch_count;
2402 std::vector<std::list<
2403 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2405 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2406 const ::internal::p4est::types<2>::locidx num_ett =
2407 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2410 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2413 (set_vertex_info ==
true ? this->n_vertices() : 0),
2415 this->n_active_lines(),
2420 set_vertex_and_cell_info(*
this,
2423 coarse_cell_to_p4est_tree_permutation,
2452 const unsigned int deal_to_p4est_line_index[12] = {
2453 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2456 this->begin_active();
2457 cell != this->
end();
2460 const unsigned int index =
2461 coarse_cell_to_p4est_tree_permutation[cell->index()];
2462 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2464 deal_to_p4est_line_index[
e]] =
2465 cell->line(
e)->index();
2470 connectivity->ett_offset[0] = 0;
2471 std::partial_sum(edge_touch_count.begin(),
2472 edge_touch_count.end(),
2473 &connectivity->ett_offset[1]);
2475 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2478 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2480 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2484 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2485 unsigned int>>::const_iterator p =
2486 edge_to_cell[v].begin();
2487 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2489 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2490 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2491 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2492 deal_to_p4est_line_index[p->second];
2496 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2501 this->mpi_communicator,
2518 template <
int dim,
int spacedim>
2520 enforce_mesh_balance_over_periodic_boundaries(
2526 std::vector<bool> flags_before[2];
2530 std::vector<unsigned int> topological_vertex_numbering(
2532 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2533 topological_vertex_numbering[i] = i;
2551 using cell_iterator =
2553 typename std::map<std::pair<cell_iterator, unsigned int>,
2554 std::pair<std::pair<cell_iterator, unsigned int>,
2555 std::bitset<3>>>::const_iterator it;
2560 const cell_iterator &cell_1 = it->first.first;
2561 const unsigned int face_no_1 = it->first.second;
2562 const cell_iterator &cell_2 = it->second.first.first;
2563 const unsigned int face_no_2 = it->second.first.second;
2564 const std::bitset<3> face_orientation = it->second.second;
2566 if (cell_1->level() == cell_2->level())
2568 for (
unsigned int v = 0;
2574 const unsigned int vface0 =
2577 face_orientation[0],
2578 face_orientation[1],
2579 face_orientation[2]);
2580 const unsigned int vi0 =
2581 topological_vertex_numbering[cell_1->face(face_no_1)
2582 ->vertex_index(vface0)];
2583 const unsigned int vi1 =
2584 topological_vertex_numbering[cell_2->face(face_no_2)
2586 const unsigned int min_index =
std::min(vi0, vi1);
2587 topological_vertex_numbering[cell_1->face(face_no_1)
2588 ->vertex_index(vface0)] =
2589 topological_vertex_numbering[cell_2->face(face_no_2)
2590 ->vertex_index(v)] =
2597 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2599 const unsigned int j = topological_vertex_numbering[i];
2607 bool continue_iterating =
true;
2609 while (continue_iterating)
2613 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2617 for (; cell != endc; ++cell)
2619 if (cell->refine_flag_set())
2620 for (
const unsigned int vertex :
2622 vertex_level[topological_vertex_numbering
2623 [cell->vertex_index(vertex)]] =
2624 std::max(vertex_level[topological_vertex_numbering
2625 [cell->vertex_index(vertex)]],
2627 else if (!cell->coarsen_flag_set())
2628 for (
const unsigned int vertex :
2630 vertex_level[topological_vertex_numbering
2631 [cell->vertex_index(vertex)]] =
2632 std::max(vertex_level[topological_vertex_numbering
2633 [cell->vertex_index(vertex)]],
2644 for (
const unsigned int vertex :
2646 vertex_level[topological_vertex_numbering
2647 [cell->vertex_index(vertex)]] =
2648 std::max(vertex_level[topological_vertex_numbering
2649 [cell->vertex_index(vertex)]],
2654 continue_iterating =
false;
2666 if (cell->refine_flag_set() ==
false)
2668 for (
const unsigned int vertex :
2670 if (vertex_level[topological_vertex_numbering
2671 [cell->vertex_index(vertex)]] >=
2675 cell->clear_coarsen_flag();
2680 if (vertex_level[topological_vertex_numbering
2681 [cell->vertex_index(vertex)]] >
2684 cell->set_refine_flag();
2685 continue_iterating =
true;
2687 for (
const unsigned int v :
2689 vertex_level[topological_vertex_numbering
2690 [cell->vertex_index(v)]] =
2692 vertex_level[topological_vertex_numbering
2693 [cell->vertex_index(v)]],
2707 if (cell->is_active())
2710 const unsigned int n_children = cell->n_children();
2711 unsigned int flagged_children = 0;
2712 for (
unsigned int child = 0; child < n_children; ++child)
2713 if (cell->child(child)->is_active() &&
2714 cell->child(child)->coarsen_flag_set())
2719 if (flagged_children < n_children)
2720 for (
unsigned int child = 0; child < n_children; ++child)
2721 if (cell->child(child)->is_active())
2722 cell->child(child)->clear_coarsen_flag();
2725 std::vector<bool> flags_after[2];
2728 return ((flags_before[0] != flags_after[0]) ||
2729 (flags_before[1] != flags_after[1]));
2735 template <
int dim,
int spacedim>
2739 bool mesh_changed =
false;
2740 unsigned int loop_counter = 0;
2741 unsigned int n_changes = 0;
2746 this->update_periodic_face_map();
2748 mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*
this);
2749 n_changes += mesh_changed;
2760 "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2761 "for periodic boundaries detected. Aborting."));
2763 while (mesh_changed);
2766 return n_changes > 0;
2771 template <
int dim,
int spacedim>
2795 spacedim>::limit_level_difference_at_vertices;
2799 bool mesh_changed =
false;
2807 if (
settings & mesh_reconstruction_after_repartitioning)
2808 while (this->n_levels() > 1)
2815 for (
const auto &cell :
2816 this->active_cell_iterators_on_level(this->n_levels() - 1))
2818 cell->set_coarsen_flag();
2836 if (parallel_ghost !=
nullptr)
2840 parallel_ghost =
nullptr;
2844 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2845 P4EST_CONNECT_CORNER) :
2846 typename ::internal::p4est::types<dim>::balance_type(
2847 P8EST_CONNECT_CORNER)));
2854 for (
const auto &cell : this->cell_iterators_on_level(0))
2859 for (
const auto &cell : this->cell_iterators_on_level(0))
2864 if (tree_exists_locally<dim, spacedim>(
2866 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2869 delete_all_children<dim, spacedim>(cell);
2870 if (cell->is_active())
2879 typename ::internal::p4est::types<dim>::quadrant
2881 typename ::internal::p4est::types<dim>::tree *tree =
2882 init_tree(cell->index());
2884 ::internal::p4est::init_coarse_quadrant<dim>(
2887 match_tree_recursively<dim, spacedim>(*tree,
2891 this->my_subdomain);
2898 typename ::internal::p4est::types<dim>::quadrant *quadr;
2900 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2902 for (
unsigned int g_idx = 0;
2903 g_idx < parallel_ghost->ghosts.elem_count;
2906 while (g_idx >=
static_cast<unsigned int>(
2907 parallel_ghost->proc_offsets[ghost_owner + 1]))
2909 while (g_idx >=
static_cast<unsigned int>(
2910 parallel_ghost->tree_offsets[ghost_tree + 1]))
2913 quadr =
static_cast<
2914 typename ::internal::p4est::types<dim>::quadrant *
>(
2915 sc_array_index(¶llel_ghost->ghosts, g_idx));
2917 unsigned int coarse_cell_index =
2918 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2920 match_quadrant<dim, spacedim>(
this,
2927 this->prepare_coarsening_and_refinement();
2931 std::any_of(this->begin_active(),
2934 return cell.refine_flag_set() ||
2935 cell.coarsen_flag_set();
2953 while (mesh_changed);
2957 unsigned int num_ghosts = 0;
2959 for (
const auto &cell : this->active_cell_iterators())
2961 if (cell->subdomain_id() != this->my_subdomain &&
2966 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2983 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
2986 for (
const auto &cell : this->cell_iterators_on_level(lvl))
2988 if ((cell->is_active() &&
2989 cell->subdomain_id() ==
2991 (cell->has_children() &&
2992 cell->child(0)->level_subdomain_id() ==
2994 cell->set_level_subdomain_id(
2995 this->locally_owned_subdomain());
2999 cell->set_level_subdomain_id(
3007 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3008 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3009 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3011 for (
const auto &cell : this->cell_iterators_on_level(0))
3013 typename ::internal::p4est::types<dim>::quadrant
3015 const unsigned int tree_index =
3016 coarse_cell_to_p4est_tree_permutation[cell->index()];
3017 typename ::internal::p4est::types<dim>::tree *tree =
3018 init_tree(cell->index());
3020 ::internal::p4est::init_coarse_quadrant<dim>(
3023 determine_level_subdomain_id_recursively<dim, spacedim>(
3034 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3037 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3039 if (cell->has_children())
3040 for (
unsigned int c = 0;
3041 c < GeometryInfo<dim>::max_children_per_cell;
3044 if (cell->child(c)->level_subdomain_id() ==
3050 cell->child(0)->level_subdomain_id();
3054 cell->set_level_subdomain_id(mark);
3071 (void)total_local_cells;
3075 Assert(
static_cast<unsigned int>(
3076 parallel_forest->local_num_quadrants) == total_local_cells,
3081 Assert(
static_cast<unsigned int>(
3082 parallel_forest->local_num_quadrants) <= total_local_cells,
3087 unsigned int n_owned = 0;
3088 for (
const auto &cell : this->active_cell_iterators())
3090 if (cell->subdomain_id() == this->my_subdomain)
3094 Assert(
static_cast<unsigned int>(
3095 parallel_forest->local_num_quadrants) == n_owned,
3099 this->smooth_grid = save_smooth;
3105 update_cell_relations();
3110 template <
int dim,
int spacedim>
3115 std::vector<Point<dim>>
point{p};
3116 std::vector<types::subdomain_id> owner = find_point_owner_rank(
point);
3123 template <
int dim,
int spacedim>
3124 std::vector<types::subdomain_id>
3128 # ifndef P4EST_SEARCH_LOCAL
3133 "This function is only available if p4est is version 2.2 and higher."));
3135 return std::vector<unsigned int>(1,
3139 AssertThrow(this->are_vertices_communicated_to_p4est(),
3141 "Vertices need to be communicated to p4est to use this "
3142 "function. This must explicitly be turned on in the "
3143 "settings of the triangulation's constructor."));
3146 for (
const auto &
manifold_id : this->get_manifold_ids())
3151 "This function can only be used if the triangulation "
3152 "has no other manifold than a Cartesian (flat) manifold attached."));
3156 PartitionSearch<dim> partition_search;
3162 parallel_forest->user_pointer = &partition_search;
3168 sc_array_t *point_sc_array;
3172 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3175 for (
size_t i = 0; i < points.size(); ++i)
3180 double *this_sc_point =
3181 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3183 for (
unsigned int d = 0;
d < dim; ++
d)
3185 this_sc_point[
d] = p(
d);
3187 this_sc_point[dim] = -1.0;
3193 static_cast<int>(
false),
3194 &PartitionSearch<dim>::local_quadrant_fn,
3195 &PartitionSearch<dim>::local_point_fn,
3199 std::vector<types::subdomain_id> owner_rank(
3203 for (
size_t i = 0; i < points.size(); ++i)
3206 double *this_sc_point =
3207 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3212 parallel_forest->user_pointer =
this;
3215 sc_array_destroy_null(&point_sc_array);
3223 template <
int dim,
int spacedim>
3229 for (
const auto &cell : this->active_cell_iterators())
3230 if (cell->is_locally_owned() && cell->refine_flag_set())
3231 Assert(cell->refine_flag_set() ==
3234 "This class does not support anisotropic refinement"));
3239 if (this->n_levels() ==
3243 cell = this->begin_active(
3251 !(cell->refine_flag_set()),
3253 "Fatal Error: maximum refinement level of p4est reached."));
3257 this->prepare_coarsening_and_refinement();
3260 this->signals.pre_distributed_refinement();
3265 for (
const auto &cell : this->active_cell_iterators())
3266 if (cell->is_ghost() || cell->is_artificial())
3268 cell->clear_refine_flag();
3269 cell->clear_coarsen_flag();
3275 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3276 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3283 parallel_forest->user_pointer = &refine_and_coarsen_list;
3285 if (parallel_ghost !=
nullptr)
3289 parallel_ghost =
nullptr;
3294 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3299 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3306 parallel_forest->user_pointer =
this;
3312 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3313 P4EST_CONNECT_FULL) :
3314 typename ::internal::p4est::types<dim>::balance_type(
3315 P8EST_CONNECT_FULL)),
3320 update_cell_relations();
3324 this->signals.post_p4est_refinement();
3328 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3329 previous_global_first_quadrant;
3331 if (this->cell_attached_data.n_attached_data_sets > 0)
3333 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3334 std::memcpy(previous_global_first_quadrant.data(),
3335 parallel_forest->global_first_quadrant,
3337 typename ::internal::p4est::types<dim>::gloidx) *
3338 (parallel_forest->mpisize + 1));
3341 if (!(
settings & no_automatic_repartitioning))
3345 if (this->signals.weight.empty())
3353 const std::vector<unsigned int> cell_weights = get_cell_weights();
3359 this->mpi_communicator) > 0,
3361 "The global sum of weights over all active cells "
3362 "is zero. Please verify how you generate weights."));
3364 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3369 parallel_forest->user_pointer = &partition_weights;
3375 &PartitionWeights<dim, spacedim>::cell_weight);
3379 parallel_forest, 0,
nullptr,
nullptr);
3381 parallel_forest->user_pointer =
this;
3386 if (this->cell_attached_data.n_attached_data_sets > 0)
3388 this->data_transfer.pack_data(
3389 this->local_cell_relations,
3390 this->cell_attached_data.pack_callbacks_fixed,
3391 this->cell_attached_data.pack_callbacks_variable);
3397 for (
const auto &cell : this->active_cell_iterators())
3399 cell->clear_refine_flag();
3400 cell->clear_coarsen_flag();
3405 copy_local_forest_to_triangulation();
3415 if (this->cell_attached_data.n_attached_data_sets > 0)
3417 this->execute_transfer(parallel_forest,
3418 previous_global_first_quadrant.data());
3421 this->data_transfer.unpack_cell_status(this->local_cell_relations);
3442 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3444 std::vector<bool> active_verts =
3445 this->mark_locally_active_vertices_on_level(lvl);
3447 const unsigned int maybe_coarser_lvl =
3448 (lvl > 0) ? (lvl - 1) : lvl;
3450 cell = this->
begin(maybe_coarser_lvl),
3451 endc = this->
end(lvl);
3452 for (; cell != endc; ++cell)
3453 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3455 const bool is_level_artificial =
3456 (cell->level_subdomain_id() ==
3458 bool need_to_know =
false;
3459 for (
const unsigned int vertex :
3461 if (active_verts[cell->vertex_index(vertex)])
3463 need_to_know =
true;
3468 !need_to_know || !is_level_artificial,
3470 "Internal error: the owner of cell" +
3471 cell->id().to_string() +
3472 " is unknown even though it is needed for geometric multigrid."));
3478 this->update_periodic_face_map();
3479 this->update_number_cache();
3482 this->signals.post_distributed_refinement();
3487 template <
int dim,
int spacedim>
3492 for (
const auto &cell : this->active_cell_iterators())
3493 if (cell->is_locally_owned())
3495 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3497 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3501 this->signals.pre_distributed_repartition();
3505 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3506 previous_global_first_quadrant;
3508 if (this->cell_attached_data.n_attached_data_sets > 0)
3510 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3511 std::memcpy(previous_global_first_quadrant.data(),
3512 parallel_forest->global_first_quadrant,
3514 typename ::internal::p4est::types<dim>::gloidx) *
3515 (parallel_forest->mpisize + 1));
3518 if (this->signals.weight.empty())
3530 const std::vector<unsigned int> cell_weights = get_cell_weights();
3536 this->mpi_communicator) > 0,
3538 "The global sum of weights over all active cells "
3539 "is zero. Please verify how you generate weights."));
3541 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3546 parallel_forest->user_pointer = &partition_weights;
3552 &PartitionWeights<dim, spacedim>::cell_weight);
3555 parallel_forest->user_pointer =
this;
3559 if (this->cell_attached_data.n_attached_data_sets > 0)
3561 this->data_transfer.pack_data(
3562 this->local_cell_relations,
3563 this->cell_attached_data.pack_callbacks_fixed,
3564 this->cell_attached_data.pack_callbacks_variable);
3569 copy_local_forest_to_triangulation();
3579 if (this->cell_attached_data.n_attached_data_sets > 0)
3581 this->execute_transfer(parallel_forest,
3582 previous_global_first_quadrant.data());
3585 this->update_periodic_face_map();
3588 this->update_number_cache();
3591 this->signals.post_distributed_repartition();
3596 template <
int dim,
int spacedim>
3597 const std::vector<types::global_dof_index> &
3601 return p4est_tree_to_coarse_cell_permutation;
3606 template <
int dim,
int spacedim>
3607 const std::vector<types::global_dof_index> &
3611 return coarse_cell_to_p4est_tree_permutation;
3616 template <
int dim,
int spacedim>
3619 const int level)
const
3623 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3624 for (
const auto &cell : this->cell_iterators_on_level(
level))
3627 marked_vertices[cell->vertex_index(v)] =
true;
3644 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3645 for (
const auto &it : this->get_periodic_face_map())
3648 const unsigned int face_no_1 = it.first.second;
3650 const unsigned int face_no_2 = it.second.first.second;
3651 const std::bitset<3> &face_orientation = it.second.second;
3653 if (cell_1->level() ==
level && cell_2->level() ==
level)
3655 for (
unsigned int v = 0;
3661 const unsigned int vface0 =
3664 face_orientation[0],
3665 face_orientation[1],
3666 face_orientation[2]);
3667 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3669 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3671 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3673 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3679 return marked_vertices;
3684 template <
int dim,
int spacedim>
3694 template <
int dim,
int spacedim>
3697 const unsigned int coarse_cell_index)
const
3699 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3704 template <
int dim,
int spacedim>
3708 &periodicity_vector)
3710 Assert(triangulation_has_content ==
true,
3712 Assert(this->n_levels() == 1,
3713 ExcMessage(
"The triangulation is refined!"));
3720 for (
const auto &face_pair : periodicity_vector)
3724 const unsigned int face_left = face_pair.face_idx[0];
3725 const unsigned int face_right = face_pair.face_idx[1];
3728 const unsigned int tree_left =
3729 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3730 const unsigned int tree_right =
3731 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3740 unsigned int p4est_orientation = 0;
3742 p4est_orientation = face_pair.orientation[1];
3745 const unsigned int face_idx_list[] = {face_left, face_right};
3746 const cell_iterator cell_list[] = {first_cell, second_cell};
3747 unsigned int lower_idx, higher_idx;
3748 if (face_left <= face_right)
3761 unsigned int first_p4est_idx_on_cell =
3762 p8est_face_corners[face_idx_list[lower_idx]][0];
3763 unsigned int first_dealii_idx_on_face =
3765 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3768 const unsigned int first_dealii_idx_on_cell =
3770 face_idx_list[lower_idx],
3772 cell_list[lower_idx]->face_orientation(
3773 face_idx_list[lower_idx]),
3774 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3775 cell_list[lower_idx]->face_rotation(
3776 face_idx_list[lower_idx]));
3777 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3779 first_dealii_idx_on_face = i;
3786 constexpr
unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
3794 constexpr
unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
3802 const unsigned int second_dealii_idx_on_face =
3803 lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
3804 [first_dealii_idx_on_face] :
3805 right_to_left[face_pair.orientation.to_ulong()]
3806 [first_dealii_idx_on_face];
3807 const unsigned int second_dealii_idx_on_cell =
3809 face_idx_list[higher_idx],
3810 second_dealii_idx_on_face,
3811 cell_list[higher_idx]->face_orientation(
3812 face_idx_list[higher_idx]),
3813 cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
3814 cell_list[higher_idx]->face_rotation(
3815 face_idx_list[higher_idx]));
3817 const unsigned int second_p4est_idx_on_face =
3818 p8est_corner_face_corners[second_dealii_idx_on_cell]
3819 [face_idx_list[higher_idx]];
3820 p4est_orientation = second_p4est_idx_on_face;
3840 this->mpi_communicator,
3851 copy_local_forest_to_triangulation();
3862 this->update_number_cache();
3867 template <
int dim,
int spacedim>
3878 this->cell_attached_data.n_attached_data_sets) +
3885 coarse_cell_to_p4est_tree_permutation) +
3887 p4est_tree_to_coarse_cell_permutation) +
3888 memory_consumption_p4est();
3895 template <
int dim,
int spacedim>
3899 return ::internal::p4est::functions<dim>::forest_memory_used(
3907 template <
int dim,
int spacedim>
3910 const ::Triangulation<dim, spacedim> &other_tria)
3914 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
3916 (other_tria.n_global_levels() == 1),
3927 const typename ::Triangulation<dim, spacedim>::DistortedCellList
3935 if (const ::parallel::distributed::Triangulation<dim, spacedim>
3936 *other_distributed =
3937 dynamic_cast<const ::parallel::distributed::
3938 Triangulation<dim, spacedim> *
>(&other_tria))
3941 settings = other_distributed->settings;
3942 triangulation_has_content =
3943 other_distributed->triangulation_has_content;
3944 coarse_cell_to_p4est_tree_permutation =
3945 other_distributed->coarse_cell_to_p4est_tree_permutation;
3946 p4est_tree_to_coarse_cell_permutation =
3947 other_distributed->p4est_tree_to_coarse_cell_permutation;
3950 typename ::internal::p4est::types<dim>::connectivity
3951 *temp_connectivity =
const_cast<
3952 typename ::internal::p4est::types<dim>::connectivity *
>(
3953 other_distributed->connectivity);
3955 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
3958 typename ::internal::p4est::types<dim>::forest *temp_forest =
3959 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
3960 other_distributed->parallel_forest);
3964 parallel_forest->connectivity = connectivity;
3965 parallel_forest->user_pointer =
this;
3969 triangulation_has_content =
true;
3970 setup_coarse_cell_to_p4est_tree_permutation();
3971 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
3976 copy_local_forest_to_triangulation();
3985 this->update_periodic_face_map();
3986 this->update_number_cache();
3991 template <
int dim,
int spacedim>
3996 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
3997 this->local_cell_relations.shrink_to_fit();
4000 for (
const auto &cell : this->cell_iterators_on_level(0))
4003 if (tree_exists_locally<dim, spacedim>(
4005 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4009 typename ::internal::p4est::types<dim>::quadrant
4011 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4014 typename ::internal::p4est::types<dim>::tree *tree =
4015 init_tree(cell->index());
4017 update_cell_relations_recursively<dim, spacedim>(
4018 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4024 template <
int dim,
int spacedim>
4025 std::vector<unsigned int>
4030 Assert(this->local_cell_relations.size() ==
4031 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4036 std::vector<unsigned int> weights;
4037 weights.reserve(this->local_cell_relations.size());
4045 for (
const auto &cell_rel : this->local_cell_relations)
4047 const auto &cell_it = cell_rel.first;
4048 const auto &cell_status = cell_rel.second;
4050 weights.push_back(this->signals.weight(cell_it, cell_status));
4058 template <
int spacedim>
4060 const MPI_Comm &mpi_communicator,
4061 const typename ::Triangulation<1, spacedim>::MeshSmoothing
4073 template <
int spacedim>
4081 template <
int spacedim>
4082 const std::vector<types::global_dof_index> &
4086 static std::vector<types::global_dof_index> a;
4092 template <
int spacedim>
4093 std::map<unsigned int, std::set<::types::subdomain_id>>
4095 const unsigned int )
const
4099 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4104 template <
int spacedim>
4107 const unsigned int)
const
4110 return std::vector<bool>();
4115 template <
int spacedim>
4126 template <
int spacedim>
4129 const unsigned int)
const
4137 template <
int spacedim>
4146 template <
int spacedim>
4155 template <
int spacedim>
4164 template <
int spacedim>
4174 template <
int spacedim>
4184 template <
int spacedim>
4201 namespace distributed
4203 template <
int dim,
int spacedim>
4211 #ifdef DEAL_II_WITH_P4EST
4221 const auto &cell = pair.first;
4222 const auto &status = pair.second;
4226 case ::Triangulation<dim, spacedim>::CELL_PERSIST:
4228 cell->clear_refine_flag();
4229 cell->clear_coarsen_flag();
4232 case ::Triangulation<dim, spacedim>::CELL_REFINE:
4234 cell->clear_coarsen_flag();
4235 cell->set_refine_flag();
4238 case ::Triangulation<dim, spacedim>::CELL_COARSEN:
4240 for (
const auto &child : cell->child_iterators())
4242 child->clear_refine_flag();
4243 child->set_coarsen_flag();
4247 case ::Triangulation<dim, spacedim>::CELL_INVALID:
4262 template <
int dim,
int spacedim>
4265 #ifdef DEAL_II_WITH_P4EST
4266 if (distributed_tria)
4269 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4270 distributed_tria->load_refine_flags(saved_refine_flags);
4274 (void)distributed_tria;
4283 #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