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;
1517 Assert(weight <
static_cast<unsigned int>(std::numeric_limits<int>::max()),
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>
1703 const typename ::Triangulation<dim, spacedim>::MeshSmoothing
1711 (
settings & construct_multigrid_hierarchy) ?
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(),
1929 std::vector<int>::size_type(0)));
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()."));
2009 Assert(this->n_cells() > 0,
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>
2061 this->n_cells() > 0,
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."));
2101 Assert(this->n_cells(0) == n_coarse_cells,
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)
2172 Assert(this->n_cells() > 0,
2174 "load() only works if the Triangulation already contains "
2176 Assert(this->n_cells() == forest->trees->elem_count,
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 std::vector<bool> flags_before[2];
2742 this->save_coarsen_flags(flags_before[0]);
2743 this->save_refine_flags(flags_before[1]);
2745 bool mesh_changed =
false;
2746 unsigned int loop_counter = 0;
2751 this->update_periodic_face_map();
2753 mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*
this);
2764 "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2765 "for periodic boundaries detected. Aborting."));
2767 while (mesh_changed);
2771 std::vector<bool> flags_after[2];
2772 this->save_coarsen_flags(flags_after[0]);
2773 this->save_refine_flags(flags_after[1]);
2774 return ((flags_before[0] != flags_after[0]) ||
2775 (flags_before[1] != flags_after[1]));
2780 template <
int dim,
int spacedim>
2801 if (
settings & construct_multigrid_hierarchy)
2804 spacedim>::limit_level_difference_at_vertices;
2808 bool mesh_changed =
false;
2816 if (
settings & mesh_reconstruction_after_repartitioning)
2817 while (this->n_levels() > 1)
2824 for (
const auto &cell :
2825 this->active_cell_iterators_on_level(this->n_levels() - 1))
2827 cell->set_coarsen_flag();
2845 if (parallel_ghost !=
nullptr)
2849 parallel_ghost =
nullptr;
2853 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2854 P4EST_CONNECT_CORNER) :
2855 typename ::internal::p4est::types<dim>::balance_type(
2856 P8EST_CONNECT_CORNER)));
2863 for (
const auto &cell : this->cell_iterators_on_level(0))
2868 for (
const auto &cell : this->cell_iterators_on_level(0))
2873 if (tree_exists_locally<dim, spacedim>(
2875 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2878 delete_all_children<dim, spacedim>(cell);
2879 if (cell->is_active())
2888 typename ::internal::p4est::types<dim>::quadrant
2890 typename ::internal::p4est::types<dim>::tree *tree =
2891 init_tree(cell->index());
2893 ::internal::p4est::init_coarse_quadrant<dim>(
2896 match_tree_recursively<dim, spacedim>(*tree,
2900 this->my_subdomain);
2907 typename ::internal::p4est::types<dim>::quadrant *quadr;
2909 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2911 for (
unsigned int g_idx = 0;
2912 g_idx < parallel_ghost->ghosts.elem_count;
2915 while (g_idx >=
static_cast<unsigned int>(
2916 parallel_ghost->proc_offsets[ghost_owner + 1]))
2918 while (g_idx >=
static_cast<unsigned int>(
2919 parallel_ghost->tree_offsets[ghost_tree + 1]))
2922 quadr =
static_cast<
2923 typename ::internal::p4est::types<dim>::quadrant *
>(
2924 sc_array_index(¶llel_ghost->ghosts, g_idx));
2926 unsigned int coarse_cell_index =
2927 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2929 match_quadrant<dim, spacedim>(
this,
2936 this->prepare_coarsening_and_refinement();
2940 std::any_of(this->begin_active(),
2943 return cell.refine_flag_set() ||
2944 cell.coarsen_flag_set();
2962 while (mesh_changed);
2966 unsigned int num_ghosts = 0;
2968 for (
const auto &cell : this->active_cell_iterators())
2970 if (cell->subdomain_id() != this->my_subdomain &&
2975 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2986 if (
settings & construct_multigrid_hierarchy)
2992 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
2995 for (
const auto &cell : this->cell_iterators_on_level(lvl))
2997 if ((cell->is_active() &&
2998 cell->subdomain_id() ==
2999 this->locally_owned_subdomain()) ||
3000 (cell->has_children() &&
3001 cell->child(0)->level_subdomain_id() ==
3002 this->locally_owned_subdomain()))
3003 cell->set_level_subdomain_id(
3004 this->locally_owned_subdomain());
3008 cell->set_level_subdomain_id(
3016 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3017 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3018 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3020 for (
const auto &cell : this->cell_iterators_on_level(0))
3022 typename ::internal::p4est::types<dim>::quadrant
3024 const unsigned int tree_index =
3025 coarse_cell_to_p4est_tree_permutation[cell->index()];
3026 typename ::internal::p4est::types<dim>::tree *tree =
3027 init_tree(cell->index());
3029 ::internal::p4est::init_coarse_quadrant<dim>(
3032 determine_level_subdomain_id_recursively<dim, spacedim>(
3043 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3046 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3048 if (cell->has_children())
3049 for (
unsigned int c = 0;
3050 c < GeometryInfo<dim>::max_children_per_cell;
3053 if (cell->child(c)->level_subdomain_id() ==
3054 this->locally_owned_subdomain())
3059 cell->child(0)->level_subdomain_id();
3063 cell->set_level_subdomain_id(mark);
3079 const unsigned int total_local_cells = this->n_active_cells();
3080 (void)total_local_cells;
3084 Assert(
static_cast<unsigned int>(
3085 parallel_forest->local_num_quadrants) == total_local_cells,
3090 Assert(
static_cast<unsigned int>(
3091 parallel_forest->local_num_quadrants) <= total_local_cells,
3096 unsigned int n_owned = 0;
3097 for (
const auto &cell : this->active_cell_iterators())
3099 if (cell->subdomain_id() == this->my_subdomain)
3103 Assert(
static_cast<unsigned int>(
3104 parallel_forest->local_num_quadrants) == n_owned,
3108 this->smooth_grid = save_smooth;
3114 update_cell_relations();
3119 template <
int dim,
int spacedim>
3124 std::vector<Point<dim>> point{p};
3125 std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
3132 template <
int dim,
int spacedim>
3133 std::vector<types::subdomain_id>
3137# ifndef P4EST_SEARCH_LOCAL
3142 "This function is only available if p4est is version 2.2 and higher."));
3144 return std::vector<unsigned int>(1,
3148 AssertThrow(this->are_vertices_communicated_to_p4est(),
3150 "Vertices need to be communicated to p4est to use this "
3151 "function. This must explicitly be turned on in the "
3152 "settings of the triangulation's constructor."));
3155 for (
const auto &manifold_id : this->get_manifold_ids())
3160 "This function can only be used if the triangulation "
3161 "has no other manifold than a Cartesian (flat) manifold attached."));
3165 PartitionSearch<dim> partition_search;
3171 parallel_forest->user_pointer = &partition_search;
3177 sc_array_t *point_sc_array;
3181 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3184 for (
size_t i = 0; i < points.size(); ++i)
3189 double *this_sc_point =
3190 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3192 for (
unsigned int d = 0; d < dim; ++d)
3194 this_sc_point[d] = p(d);
3196 this_sc_point[dim] = -1.0;
3202 static_cast<int>(
false),
3203 &PartitionSearch<dim>::local_quadrant_fn,
3204 &PartitionSearch<dim>::local_point_fn,
3208 std::vector<types::subdomain_id> owner_rank(
3212 for (
size_t i = 0; i < points.size(); ++i)
3215 double *this_sc_point =
3216 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3221 parallel_forest->user_pointer =
this;
3224 sc_array_destroy_null(&point_sc_array);
3232 template <
int dim,
int spacedim>
3238 for (
const auto &cell : this->active_cell_iterators())
3239 if (cell->is_locally_owned() && cell->refine_flag_set())
3240 Assert(cell->refine_flag_set() ==
3243 "This class does not support anisotropic refinement"));
3248 if (this->n_levels() ==
3252 cell = this->begin_active(
3260 !(cell->refine_flag_set()),
3262 "Fatal Error: maximum refinement level of p4est reached."));
3266 this->prepare_coarsening_and_refinement();
3269 this->signals.pre_distributed_refinement();
3274 for (
const auto &cell : this->active_cell_iterators())
3275 if (cell->is_ghost() || cell->is_artificial())
3277 cell->clear_refine_flag();
3278 cell->clear_coarsen_flag();
3284 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3285 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3292 parallel_forest->user_pointer = &refine_and_coarsen_list;
3294 if (parallel_ghost !=
nullptr)
3298 parallel_ghost =
nullptr;
3303 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3308 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3315 parallel_forest->user_pointer =
this;
3321 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3322 P4EST_CONNECT_FULL) :
3323 typename ::internal::p4est::types<dim>::balance_type(
3324 P8EST_CONNECT_FULL)),
3329 update_cell_relations();
3333 this->signals.post_p4est_refinement();
3337 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3338 previous_global_first_quadrant;
3340 if (this->cell_attached_data.n_attached_data_sets > 0)
3342 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3343 std::memcpy(previous_global_first_quadrant.data(),
3344 parallel_forest->global_first_quadrant,
3346 typename ::internal::p4est::types<dim>::gloidx) *
3347 (parallel_forest->mpisize + 1));
3350 if (!(
settings & no_automatic_repartitioning))
3354 if (this->signals.weight.empty())
3362 const std::vector<unsigned int> cell_weights = get_cell_weights();
3368 this->mpi_communicator) > 0,
3370 "The global sum of weights over all active cells "
3371 "is zero. Please verify how you generate weights."));
3373 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3378 parallel_forest->user_pointer = &partition_weights;
3384 &PartitionWeights<dim, spacedim>::cell_weight);
3388 parallel_forest, 0,
nullptr,
nullptr);
3390 parallel_forest->user_pointer =
this;
3395 if (this->cell_attached_data.n_attached_data_sets > 0)
3397 this->data_transfer.pack_data(
3398 this->local_cell_relations,
3399 this->cell_attached_data.pack_callbacks_fixed,
3400 this->cell_attached_data.pack_callbacks_variable);
3406 for (
const auto &cell : this->active_cell_iterators())
3408 cell->clear_refine_flag();
3409 cell->clear_coarsen_flag();
3414 copy_local_forest_to_triangulation();
3424 if (this->cell_attached_data.n_attached_data_sets > 0)
3426 this->execute_transfer(parallel_forest,
3427 previous_global_first_quadrant.data());
3430 this->data_transfer.unpack_cell_status(this->local_cell_relations);
3449 if (
settings & construct_multigrid_hierarchy)
3451 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3453 std::vector<bool> active_verts =
3454 this->mark_locally_active_vertices_on_level(lvl);
3456 const unsigned int maybe_coarser_lvl =
3457 (lvl > 0) ? (lvl - 1) : lvl;
3459 cell = this->begin(maybe_coarser_lvl),
3460 endc = this->end(lvl);
3461 for (; cell != endc; ++cell)
3462 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3464 const bool is_level_artificial =
3465 (cell->level_subdomain_id() ==
3467 bool need_to_know =
false;
3468 for (
const unsigned int vertex :
3470 if (active_verts[cell->vertex_index(vertex)])
3472 need_to_know =
true;
3477 !need_to_know || !is_level_artificial,
3479 "Internal error: the owner of cell" +
3480 cell->id().to_string() +
3481 " is unknown even though it is needed for geometric multigrid."));
3487 this->update_periodic_face_map();
3488 this->update_number_cache();
3491 this->signals.post_distributed_refinement();
3496 template <
int dim,
int spacedim>
3501 for (
const auto &cell : this->active_cell_iterators())
3502 if (cell->is_locally_owned())
3504 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3506 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3510 this->signals.pre_distributed_repartition();
3514 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3515 previous_global_first_quadrant;
3517 if (this->cell_attached_data.n_attached_data_sets > 0)
3519 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3520 std::memcpy(previous_global_first_quadrant.data(),
3521 parallel_forest->global_first_quadrant,
3523 typename ::internal::p4est::types<dim>::gloidx) *
3524 (parallel_forest->mpisize + 1));
3527 if (this->signals.weight.empty())
3539 const std::vector<unsigned int> cell_weights = get_cell_weights();
3545 this->mpi_communicator) > 0,
3547 "The global sum of weights over all active cells "
3548 "is zero. Please verify how you generate weights."));
3550 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3555 parallel_forest->user_pointer = &partition_weights;
3561 &PartitionWeights<dim, spacedim>::cell_weight);
3564 parallel_forest->user_pointer =
this;
3568 if (this->cell_attached_data.n_attached_data_sets > 0)
3570 this->data_transfer.pack_data(
3571 this->local_cell_relations,
3572 this->cell_attached_data.pack_callbacks_fixed,
3573 this->cell_attached_data.pack_callbacks_variable);
3578 copy_local_forest_to_triangulation();
3588 if (this->cell_attached_data.n_attached_data_sets > 0)
3590 this->execute_transfer(parallel_forest,
3591 previous_global_first_quadrant.data());
3594 this->update_periodic_face_map();
3597 this->update_number_cache();
3600 this->signals.post_distributed_repartition();
3605 template <
int dim,
int spacedim>
3606 const std::vector<types::global_dof_index> &
3610 return p4est_tree_to_coarse_cell_permutation;
3615 template <
int dim,
int spacedim>
3616 const std::vector<types::global_dof_index> &
3620 return coarse_cell_to_p4est_tree_permutation;
3625 template <
int dim,
int spacedim>
3628 const int level)
const
3632 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3633 for (
const auto &cell : this->cell_iterators_on_level(
level))
3634 if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3636 marked_vertices[cell->vertex_index(v)] =
true;
3653 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3654 for (
const auto &it : this->get_periodic_face_map())
3657 const unsigned int face_no_1 = it.first.second;
3659 const unsigned int face_no_2 = it.second.first.second;
3660 const std::bitset<3> &face_orientation = it.second.second;
3662 if (cell_1->level() ==
level && cell_2->level() ==
level)
3664 for (
unsigned int v = 0;
3670 const unsigned int vface0 =
3673 face_orientation[0],
3674 face_orientation[1],
3675 face_orientation[2]);
3676 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3678 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3680 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3682 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3688 return marked_vertices;
3693 template <
int dim,
int spacedim>
3698 return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3703 template <
int dim,
int spacedim>
3706 const unsigned int coarse_cell_index)
const
3708 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];
3713 template <
int dim,
int spacedim>
3717 &periodicity_vector)
3719 Assert(triangulation_has_content ==
true,
3721 Assert(this->n_levels() == 1,
3722 ExcMessage(
"The triangulation is refined!"));
3729 for (
const auto &face_pair : periodicity_vector)
3733 const unsigned int face_left = face_pair.face_idx[0];
3734 const unsigned int face_right = face_pair.face_idx[1];
3737 const unsigned int tree_left =
3738 coarse_cell_to_p4est_tree_permutation[first_cell->index()];
3739 const unsigned int tree_right =
3740 coarse_cell_to_p4est_tree_permutation[second_cell->index()];
3749 unsigned int p4est_orientation = 0;
3751 p4est_orientation = face_pair.orientation[1];
3754 const unsigned int face_idx_list[] = {face_left, face_right};
3755 const cell_iterator cell_list[] = {first_cell, second_cell};
3756 unsigned int lower_idx, higher_idx;
3757 if (face_left <= face_right)
3770 unsigned int first_p4est_idx_on_cell =
3771 p8est_face_corners[face_idx_list[lower_idx]][0];
3772 unsigned int first_dealii_idx_on_face =
3774 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
3777 const unsigned int first_dealii_idx_on_cell =
3779 face_idx_list[lower_idx],
3781 cell_list[lower_idx]->face_orientation(
3782 face_idx_list[lower_idx]),
3783 cell_list[lower_idx]->face_flip(face_idx_list[lower_idx]),
3784 cell_list[lower_idx]->face_rotation(
3785 face_idx_list[lower_idx]));
3786 if (first_p4est_idx_on_cell == first_dealii_idx_on_cell)
3788 first_dealii_idx_on_face = i;
3795 constexpr unsigned int left_to_right[8][4] = {{0, 2, 1, 3},
3803 constexpr unsigned int right_to_left[8][4] = {{0, 2, 1, 3},
3811 const unsigned int second_dealii_idx_on_face =
3812 lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()]
3813 [first_dealii_idx_on_face] :
3814 right_to_left[face_pair.orientation.to_ulong()]
3815 [first_dealii_idx_on_face];
3816 const unsigned int second_dealii_idx_on_cell =
3818 face_idx_list[higher_idx],
3819 second_dealii_idx_on_face,
3820 cell_list[higher_idx]->face_orientation(
3821 face_idx_list[higher_idx]),
3822 cell_list[higher_idx]->face_flip(face_idx_list[higher_idx]),
3823 cell_list[higher_idx]->face_rotation(
3824 face_idx_list[higher_idx]));
3826 const unsigned int second_p4est_idx_on_face =
3827 p8est_corner_face_corners[second_dealii_idx_on_cell]
3828 [face_idx_list[higher_idx]];
3829 p4est_orientation = second_p4est_idx_on_face;
3849 this->mpi_communicator,
3860 copy_local_forest_to_triangulation();
3871 this->update_number_cache();
3876 template <
int dim,
int spacedim>
3887 this->cell_attached_data.n_attached_data_sets) +
3894 coarse_cell_to_p4est_tree_permutation) +
3896 p4est_tree_to_coarse_cell_permutation) +
3897 memory_consumption_p4est();
3904 template <
int dim,
int spacedim>
3908 return ::internal::p4est::functions<dim>::forest_memory_used(
3916 template <
int dim,
int spacedim>
3919 const ::Triangulation<dim, spacedim> &other_tria)
3923 const ::parallel::distributed::Triangulation<dim, spacedim> *
>(
3925 (other_tria.n_global_levels() == 1),
3936 const typename ::Triangulation<dim, spacedim>::DistortedCellList
3944 if (const ::parallel::distributed::Triangulation<dim, spacedim>
3945 *other_distributed =
3946 dynamic_cast<const ::parallel::distributed::
3947 Triangulation<dim, spacedim> *
>(&other_tria))
3950 settings = other_distributed->settings;
3951 triangulation_has_content =
3952 other_distributed->triangulation_has_content;
3953 coarse_cell_to_p4est_tree_permutation =
3954 other_distributed->coarse_cell_to_p4est_tree_permutation;
3955 p4est_tree_to_coarse_cell_permutation =
3956 other_distributed->p4est_tree_to_coarse_cell_permutation;
3959 typename ::internal::p4est::types<dim>::connectivity
3960 *temp_connectivity =
const_cast<
3961 typename ::internal::p4est::types<dim>::connectivity *
>(
3962 other_distributed->connectivity);
3964 ::internal::p4est::copy_connectivity<dim>(temp_connectivity);
3967 typename ::internal::p4est::types<dim>::forest *temp_forest =
3968 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
3969 other_distributed->parallel_forest);
3973 parallel_forest->connectivity = connectivity;
3974 parallel_forest->user_pointer =
this;
3978 triangulation_has_content =
true;
3979 setup_coarse_cell_to_p4est_tree_permutation();
3980 copy_new_triangulation_to_p4est(std::integral_constant<int, dim>());
3985 copy_local_forest_to_triangulation();
3994 this->update_periodic_face_map();
3995 this->update_number_cache();
4000 template <
int dim,
int spacedim>
4005 this->local_cell_relations.resize(parallel_forest->local_num_quadrants);
4006 this->local_cell_relations.shrink_to_fit();
4009 for (
const auto &cell : this->cell_iterators_on_level(0))
4012 if (tree_exists_locally<dim, spacedim>(
4014 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
false)
4018 typename ::internal::p4est::types<dim>::quadrant
4020 ::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4023 typename ::internal::p4est::types<dim>::tree *tree =
4024 init_tree(cell->index());
4026 update_cell_relations_recursively<dim, spacedim>(
4027 this->local_cell_relations, *tree, cell, p4est_coarse_cell);
4033 template <
int dim,
int spacedim>
4034 std::vector<unsigned int>
4039 Assert(this->local_cell_relations.size() ==
4040 static_cast<unsigned int>(parallel_forest->local_num_quadrants),
4045 std::vector<unsigned int> weights;
4046 weights.reserve(this->local_cell_relations.size());
4054 for (
const auto &cell_rel : this->local_cell_relations)
4056 const auto &cell_it = cell_rel.first;
4057 const auto &cell_status = cell_rel.second;
4059 weights.push_back(this->signals.weight(cell_it, cell_status));
4067 template <
int spacedim>
4070 const typename ::Triangulation<1, spacedim>::MeshSmoothing
4082 template <
int spacedim>
4090 template <
int spacedim>
4091 const std::vector<types::global_dof_index> &
4095 static std::vector<types::global_dof_index> a;
4101 template <
int spacedim>
4102 std::map<unsigned int, std::set<::types::subdomain_id>>
4104 const unsigned int )
const
4108 return std::map<unsigned int, std::set<::types::subdomain_id>>();
4113 template <
int spacedim>
4116 const unsigned int)
const
4119 return std::vector<bool>();
4124 template <
int spacedim>
4135 template <
int spacedim>
4138 const unsigned int)
const
4146 template <
int spacedim>
4155 template <
int spacedim>
4164 template <
int spacedim>
4173 template <
int spacedim>
4183 template <
int spacedim>
4193 template <
int spacedim>
4210 namespace distributed
4212 template <
int dim,
int spacedim>
4220#ifdef DEAL_II_WITH_P4EST
4230 const auto &cell = pair.first;
4231 const auto &status = pair.second;
4235 case ::Triangulation<dim, spacedim>::CELL_PERSIST:
4237 cell->clear_refine_flag();
4238 cell->clear_coarsen_flag();
4241 case ::Triangulation<dim, spacedim>::CELL_REFINE:
4243 cell->clear_coarsen_flag();
4244 cell->set_refine_flag();
4247 case ::Triangulation<dim, spacedim>::CELL_COARSEN:
4249 for (
const auto &child : cell->child_iterators())
4251 child->clear_refine_flag();
4252 child->set_coarsen_flag();
4256 case ::Triangulation<dim, spacedim>::CELL_INVALID:
4271 template <
int dim,
int spacedim>
4274#ifdef DEAL_II_WITH_P4EST
4275 if (distributed_tria)
4278 distributed_tria->load_coarsen_flags(saved_coarsen_flags);
4279 distributed_tria->load_refine_flags(saved_refine_flags);
4283 (void)distributed_tria;
active_cell_iterator last_active() const
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
cell_iterator end() const
void save_refine_flags(std::ostream &out) const
unsigned int n_vertices() const
void save_coarsen_flags(std::ostream &out) const
active_cell_iterator begin_active(const unsigned int level=0) const
virtual void clear() override
virtual std::size_t memory_consumption() const override
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &old_tria) override
TemporarilyMatchRefineFlags(::Triangulation< dim, spacedim > &tria)
~TemporarilyMatchRefineFlags()
const SmartPointer< ::parallel::distributed::Triangulation< dim, spacedim > > distributed_tria
std::vector< bool > saved_refine_flags
std::vector< bool > saved_coarsen_flags
virtual void execute_coarsening_and_refinement() override
void setup_coarse_cell_to_p4est_tree_permutation()
virtual void add_periodicity(const std::vector<::GridTools::PeriodicFacePair< cell_iterator > > &) override
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
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
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 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 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() const override
typename::internal::p4est::types< dim >::ghost * parallel_ghost
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ::Triangulation< dim, spacedim >::cell_iterator cell_iterator
typename ::Triangulation< dim, spacedim >::active_cell_iterator active_cell_iterator
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
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
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static unsigned int standard_to_real_face_vertex(const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static bool is_inside_unit_cell(const Point< dim > &p)
static Point< dim > unit_cell_vertex(const unsigned int vertex)
const TriangulationDescription::Settings settings
const ::Triangulation< dim, spacedim > & tria