870 p_out += (quadrant_mapping_matrix(alpha, 0) +
871 quadrant_mapping_matrix(alpha, 1) * p(0) +
872 quadrant_mapping_matrix(alpha, 2) * p(1) +
873 quadrant_mapping_matrix(alpha, 3) * p(2) +
874 quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
875 quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
876 quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
877 quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
887 PartitionSearch<dim>::QuadrantData::QuadrantData()
891 , are_vertices_initialized(false)
892 , is_reference_mapping_initialized(false)
899 PartitionSearch<dim>::QuadrantData::initialize_mapping()
902 are_vertices_initialized,
904 "Cell vertices must be initialized before the cell mapping can be filled."));
911 for (
unsigned int alpha = 0;
912 alpha < GeometryInfo<dim>::vertices_per_cell;
916 point_matrix(0, alpha) = 1;
917 point_matrix(1, alpha) = cell_vertices[alpha](0);
918 point_matrix(2, alpha) = cell_vertices[alpha](1);
919 point_matrix(3, alpha) =
920 cell_vertices[alpha](0) * cell_vertices[alpha](1);
927 quadrant_mapping_matrix.invert(point_matrix);
931 for (
unsigned int alpha = 0;
932 alpha < GeometryInfo<dim>::vertices_per_cell;
936 point_matrix(0, alpha) = 1;
937 point_matrix(1, alpha) = cell_vertices[alpha](0);
938 point_matrix(2, alpha) = cell_vertices[alpha](1);
939 point_matrix(3, alpha) = cell_vertices[alpha](2);
940 point_matrix(4, alpha) =
941 cell_vertices[alpha](0) * cell_vertices[alpha](1);
942 point_matrix(5, alpha) =
943 cell_vertices[alpha](1) * cell_vertices[alpha](2);
944 point_matrix(6, alpha) =
945 cell_vertices[alpha](0) * cell_vertices[alpha](2);
946 point_matrix(7, alpha) = cell_vertices[alpha](0) *
947 cell_vertices[alpha](1) *
948 cell_vertices[alpha](2);
955 quadrant_mapping_matrix.invert(point_matrix);
958 is_reference_mapping_initialized =
true;
965 PartitionSearch<2>::QuadrantData::set_cell_vertices(
966 typename internal::p4est::types<2>::forest *forest,
967 typename internal::p4est::types<2>::topidx which_tree,
968 typename internal::p4est::types<2>::quadrant *quadrant,
969 const typename internal::p4est::types<2>::quadrant_coord
970 quad_length_on_level)
972 constexpr unsigned int dim = 2;
976 double corner_point[dim + 1] = {0};
979 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
981 for (
unsigned int d = 0;
d < dim; ++
d)
983 cell_vertices[vertex_index](
d) = corner_point[d];
993 unsigned int vertex_index = 0;
994 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
995 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
998 copy_vertex(vertex_index);
1004 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1005 forest->connectivity,
1007 quadrant->x + quad_length_on_level,
1012 copy_vertex(vertex_index);
1018 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1019 forest->connectivity,
1022 quadrant->y + quad_length_on_level,
1026 copy_vertex(vertex_index);
1032 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1033 forest->connectivity,
1035 quadrant->x + quad_length_on_level,
1036 quadrant->y + quad_length_on_level,
1040 copy_vertex(vertex_index);
1042 are_vertices_initialized =
true;
1049 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1050 typename internal::p4est::types<3>::forest *forest,
1051 typename internal::p4est::types<3>::topidx which_tree,
1052 typename internal::p4est::types<3>::quadrant *quadrant,
1053 const typename internal::p4est::types<3>::quadrant_coord
1054 quad_length_on_level)
1056 constexpr unsigned int dim = 3;
1058 double corner_point[dim] = {0};
1061 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1063 for (
unsigned int d = 0;
d < dim; ++
d)
1065 cell_vertices[vertex_index](
d) = corner_point[d];
1067 corner_point[
d] = 0;
1075 unsigned int vertex_index = 0;
1076 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1077 forest->connectivity,
1085 copy_vertex(vertex_index);
1092 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1093 forest->connectivity,
1095 quadrant->x + quad_length_on_level,
1101 copy_vertex(vertex_index);
1107 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1108 forest->connectivity,
1111 quadrant->y + quad_length_on_level,
1116 copy_vertex(vertex_index);
1122 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1123 forest->connectivity,
1125 quadrant->x + quad_length_on_level,
1126 quadrant->y + quad_length_on_level,
1131 copy_vertex(vertex_index);
1137 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1138 forest->connectivity,
1142 quadrant->z + quad_length_on_level,
1146 copy_vertex(vertex_index);
1152 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1153 forest->connectivity,
1155 quadrant->x + quad_length_on_level,
1157 quadrant->z + quad_length_on_level,
1161 copy_vertex(vertex_index);
1167 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1168 forest->connectivity,
1171 quadrant->y + quad_length_on_level,
1172 quadrant->z + quad_length_on_level,
1176 copy_vertex(vertex_index);
1182 internal::p4est::functions<dim>::quadrant_coord_to_vertex(
1183 forest->connectivity,
1185 quadrant->x + quad_length_on_level,
1186 quadrant->y + quad_length_on_level,
1187 quadrant->z + quad_length_on_level,
1191 copy_vertex(vertex_index);
1194 are_vertices_initialized =
true;
1204 template <
int dim,
int spacedim>
1205 class RefineAndCoarsenList
1208 RefineAndCoarsenList(
const Triangulation<dim, spacedim> &triangulation,
1209 const std::vector<types::global_dof_index>
1210 &p4est_tree_to_coarse_cell_permutation,
1223 typename internal::p4est::types<dim>::forest *forest,
1224 typename internal::p4est::types<dim>::topidx coarse_cell_index,
1225 typename internal::p4est::types<dim>::quadrant *quadrant);
1233 typename internal::p4est::types<dim>::forest *forest,
1234 typename internal::p4est::types<dim>::topidx coarse_cell_index,
1235 typename internal::p4est::types<dim>::quadrant *children[]);
1238 pointers_are_at_end()
const;
1241 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1242 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1243 const_iterator current_refine_pointer;
1245 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1246 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1247 const_iterator current_coarsen_pointer;
1252 const typename internal::p4est::types<dim>::quadrant &p4est_cell,
1258 template <
int dim,
int spacedim>
1260 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1262 return ((current_refine_pointer == refine_list.end()) &&
1263 (current_coarsen_pointer == coarsen_list.end()));
1268 template <
int dim,
int spacedim>
1269 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1271 const std::vector<types::global_dof_index>
1272 &p4est_tree_to_coarse_cell_permutation,
1276 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1289 refine_list.reserve(n_refine_flags);
1290 coarsen_list.reserve(n_coarsen_flags);
1300 for (
unsigned int c = 0; c < triangulation.
n_cells(0); ++c)
1302 unsigned int coarse_cell_index =
1303 p4est_tree_to_coarse_cell_permutation[c];
1306 &triangulation, 0, coarse_cell_index);
1308 typename internal::p4est::types<dim>::quadrant p4est_cell;
1309 internal::p4est::functions<dim>::quadrant_set_morton(&p4est_cell,
1312 p4est_cell.p.which_tree = c;
1313 build_lists(cell, p4est_cell, my_subdomain);
1321 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1322 Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1324 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1325 Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1328 current_refine_pointer = refine_list.begin();
1329 current_coarsen_pointer = coarsen_list.begin();
1334 template <
int dim,
int spacedim>
1336 RefineAndCoarsenList<dim, spacedim>::build_lists(
1338 const typename internal::p4est::types<dim>::quadrant &p4est_cell,
1346 refine_list.push_back(p4est_cell);
1348 coarsen_list.push_back(p4est_cell);
1353 typename internal::p4est::types<dim>::quadrant
1355 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1360 P4EST_QUADRANT_INIT(&p4est_child[c]);
1363 P8EST_QUADRANT_INIT(&p4est_child[c]);
1368 internal::p4est::functions<dim>::quadrant_childrenv(&p4est_cell,
1370 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1373 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1374 build_lists(cell->
child(c), p4est_child[c], my_subdomain);
1380 template <
int dim,
int spacedim>
1382 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1383 typename internal::p4est::types<dim>::forest *forest,
1384 typename internal::p4est::types<dim>::topidx coarse_cell_index,
1385 typename internal::p4est::types<dim>::quadrant *quadrant)
1387 RefineAndCoarsenList<dim, spacedim> *this_object =
1388 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1389 forest->user_pointer);
1393 if (this_object->current_refine_pointer == this_object->refine_list.end())
1396 Assert(coarse_cell_index <=
1397 this_object->current_refine_pointer->p.which_tree,
1402 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1406 Assert(coarse_cell_index <=
1407 this_object->current_refine_pointer->p.which_tree,
1412 Assert(internal::p4est::functions<dim>::quadrant_compare(
1413 quadrant, &*this_object->current_refine_pointer) <= 0,
1417 if (internal::p4est::functions<dim>::quadrant_is_equal(
1418 quadrant, &*this_object->current_refine_pointer))
1420 ++this_object->current_refine_pointer;
1430 template <
int dim,
int spacedim>
1432 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1433 typename internal::p4est::types<dim>::forest *forest,
1434 typename internal::p4est::types<dim>::topidx coarse_cell_index,
1435 typename internal::p4est::types<dim>::quadrant *children[])
1437 RefineAndCoarsenList<dim, spacedim> *this_object =
1438 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1439 forest->user_pointer);
1443 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1446 Assert(coarse_cell_index <=
1447 this_object->current_coarsen_pointer->p.which_tree,
1452 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1456 Assert(coarse_cell_index <=
1457 this_object->current_coarsen_pointer->p.which_tree,
1462 Assert(internal::p4est::functions<dim>::quadrant_compare(
1463 children[0], &*this_object->current_coarsen_pointer) <= 0,
1468 if (internal::p4est::functions<dim>::quadrant_is_equal(
1469 children[0], &*this_object->current_coarsen_pointer))
1472 ++this_object->current_coarsen_pointer;
1476 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1479 Assert(internal::p4est::functions<dim>::quadrant_is_equal(
1480 children[c], &*this_object->current_coarsen_pointer),
1482 ++this_object->current_coarsen_pointer;
1500 template <
int dim,
int spacedim>
1501 class PartitionWeights
1509 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1519 cell_weight(
typename internal::p4est::types<dim>::forest *forest,
1520 typename internal::p4est::types<dim>::topidx coarse_cell_index,
1521 typename internal::p4est::types<dim>::quadrant *quadrant);
1524 std::vector<unsigned int> cell_weights_list;
1525 std::vector<unsigned int>::const_iterator current_pointer;
1529 template <
int dim,
int spacedim>
1530 PartitionWeights<dim, spacedim>::PartitionWeights(
1531 const std::vector<unsigned int> &cell_weights)
1532 : cell_weights_list(cell_weights)
1536 current_pointer = cell_weights_list.begin();
1540 template <
int dim,
int spacedim>
1542 PartitionWeights<dim, spacedim>::cell_weight(
1543 typename internal::p4est::types<dim>::forest *forest,
1544 typename internal::p4est::types<dim>::topidx,
1545 typename internal::p4est::types<dim>::quadrant *)
1551 PartitionWeights<dim, spacedim> *this_object =
1552 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1554 Assert(this_object->current_pointer >=
1555 this_object->cell_weights_list.begin(),
1557 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1563 const unsigned int weight = *this_object->current_pointer;
1564 ++this_object->current_pointer;
1566 Assert(weight <
static_cast<unsigned int>(std::numeric_limits<int>::max()),
1567 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1568 "weights for cells. The weight provided here exceeds "
1569 "the maximum value represented as a 'signed int'."));
1570 return static_cast<int>(weight);
1573 template <
int dim,
int spacedim>
1574 using cell_relation_t =
typename std::pair<
1575 typename ::Triangulation<dim, spacedim>::cell_iterator,
1587 template <
int dim,
int spacedim>
1589 add_single_cell_relation(
1590 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1591 const typename ::internal::p4est::types<dim>::tree &tree,
1592 const unsigned int idx,
1596 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1602 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1616 template <
int dim,
int spacedim>
1618 update_cell_relations_recursively(
1619 std::vector<cell_relation_t<dim, spacedim>> &cell_rel,
1620 const typename ::internal::p4est::types<dim>::tree &tree,
1622 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1625 const int idx = sc_array_bsearch(
1626 const_cast<sc_array_t *
>(&tree.quadrants),
1628 ::internal::p4est::functions<dim>::quadrant_compare);
1630 (::internal::p4est::functions<dim>::quadrant_overlaps_tree(
1631 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1633 &p4est_cell) ==
false))
1638 const bool p4est_has_children = (idx == -1);
1642 typename ::internal::p4est::types<dim>::quadrant
1645 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1650 P4EST_QUADRANT_INIT(&p4est_child[c]);
1653 P8EST_QUADRANT_INIT(&p4est_child[c]);
1659 ::internal::p4est::functions<dim>::quadrant_childrenv(
1660 &p4est_cell, p4est_child);
1662 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1665 update_cell_relations_recursively<dim, spacedim>(
1666 cell_rel, tree, dealii_cell->
child(c), p4est_child[c]);
1669 else if (!p4est_has_children && !dealii_cell->
has_children())
1673 add_single_cell_relation<dim, spacedim>(
1676 else if (p4est_has_children)
1684 typename ::internal::p4est::types<dim>::quadrant
1686 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1691 P4EST_QUADRANT_INIT(&p4est_child[c]);
1694 P8EST_QUADRANT_INIT(&p4est_child[c]);
1700 ::internal::p4est::functions<dim>::quadrant_childrenv(
1701 &p4est_cell, p4est_child);
1709 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1712 child_idx = sc_array_bsearch(
1713 const_cast<sc_array_t *
>(&tree.quadrants),
1715 ::internal::p4est::functions<dim>::quadrant_compare);
1720 add_single_cell_relation<dim, spacedim>(
1721 cell_rel, tree, child_idx, dealii_cell, cell_status);
1729 add_single_cell_relation<dim, spacedim>(
1746 template <
int dim,
int spacedim>
1775 template <
int dim,
int spacedim>
1796 template <
int dim,
int spacedim>
1809 const typename ::Triangulation<dim, spacedim>::DistortedCellList
1820 "The class parallel::distributed::Triangulation only supports meshes "
1821 "consisting only of hypercube-like cells."));
1836 catch (
const typename Triangulation<dim>::DistortedCellList &)
1849 template <
int dim,
int spacedim>
1855 (void)construction_data;
1862 template <
int dim,
int spacedim>
1898 template <
int dim,
int spacedim>
1909 template <
int dim,
int spacedim>
1920 template <
int dim,
int spacedim>
1926 *previous_global_first_quadrant)
1934 this->data_serializer.sizes_fixed_cumulative.back());
1937 typename ::internal::p4est::types<dim>::transfer_context
1942 previous_global_first_quadrant,
1945 this->data_serializer.dest_data_fixed.data(),
1946 this->data_serializer.src_data_fixed.data(),
1947 this->data_serializer.sizes_fixed_cumulative.back());
1959 previous_global_first_quadrant,
1962 this->data_serializer.dest_sizes_variable.data(),
1963 this->data_serializer.src_sizes_variable.data(),
1964 sizeof(
unsigned int));
1978 this->data_serializer.dest_sizes_variable.end(),
1979 std::vector<int>::size_type(0)));
1981# if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1997 previous_global_first_quadrant,
2000 this->data_serializer.dest_data_variable.data(),
2001 this->data_serializer.dest_sizes_variable.data(),
2002 this->data_serializer.src_data_variable.data(),
2003 this->data_serializer.src_sizes_variable.data());
2015 template <
int dim,
int spacedim>
2033 template <
int dim,
int spacedim>
2036 const
std::
string &file_basename)
const
2039 ExcMessage(
"Can't produce output when no forest is created yet."));
2043 "To use this function the triangulation's flag "
2044 "Settings::communicate_vertices_to_p4est must be set."));
2052 template <
int dim,
int spacedim>
2055 const
std::
string &file_basename)
const
2060 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2062 ExcMessage(
"Can not save() an empty Triangulation."));
2068 this->
signals.pre_distributed_save();
2070 if (this->my_subdomain == 0)
2072 std::string fname = file_basename +
".info";
2073 std::ofstream f(fname);
2074 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2080 << this->
n_cells(0) << std::endl;
2087 Assert((cell_rel.second ==
2102 this->
signals.post_distributed_save();
2107 template <
int dim,
int spacedim>
2114 "load() only works if the Triangulation already contains a coarse mesh!"));
2118 "Triangulation may only contain coarse cells when calling load()."));
2124 this->
signals.pre_distributed_load();
2138 unsigned int version, numcpus, attached_count_fixed,
2139 attached_count_variable, n_coarse_cells;
2141 std::string fname = std::string(file_basename) +
".info";
2142 std::ifstream f(fname);
2144 std::string firstline;
2145 getline(f, firstline);
2146 f >> version >> numcpus >> attached_count_fixed >>
2147 attached_count_variable >> n_coarse_cells;
2151 ExcMessage(
"Incompatible version found in .info file."));
2153 ExcMessage(
"Number of coarse cells differ!"));
2159 attached_count_fixed + attached_count_variable;
2162 file_basename.c_str(),
2163 this->mpi_communicator,
2183 catch (
const typename Triangulation<dim>::DistortedCellList &)
2195 attached_count_fixed,
2196 attached_count_variable);
2199 this->
signals.post_distributed_load();
2207 template <
int dim,
int spacedim>
2210 const
bool autopartition)
2212 (void)autopartition;
2218 template <
int dim,
int spacedim>
2221 const typename ::
internal::p4est::
types<dim>::forest *forest)
2225 "load() only works if the Triangulation already contains "
2229 "Coarse mesh of the Triangulation does not match the one "
2230 "of the provided forest!"));
2246 typename ::internal::p4est::types<dim>::forest *temp =
2247 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2258 catch (
const typename Triangulation<dim>::DistortedCellList &)
2271 template <
int dim,
int spacedim>
2277 "Can't produce a check sum when no forest is created yet."));
2282# if !DEAL_II_P4EST_VERSION_GTE(2, 8, 6, 0)
2298 template <
int dim,
int spacedim>
2300 const typename ::internal::p4est::types<dim>::forest
2304 ExcMessage(
"The forest has not been allocated yet."));
2310 template <
int dim,
int spacedim>
2312 typename ::internal::p4est::types<dim>::tree
2314 const int dealii_coarse_cell_index)
const
2316 const unsigned int tree_index =
2318 typename ::internal::p4est::types<dim>::tree *tree =
2319 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2334 std::integral_constant<int, 2>)
2336 const unsigned int dim = 2, spacedim = 2;
2344 std::vector<unsigned int> vertex_touch_count;
2346 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2349 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2350 const ::internal::p4est::types<2>::locidx num_vtt =
2351 std::accumulate(vertex_touch_count.begin(),
2352 vertex_touch_count.end(),
2358 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2361 (set_vertex_info ==
true ? this->n_vertices() : 0),
2366 set_vertex_and_cell_info(*
this,
2369 coarse_cell_to_p4est_tree_permutation,
2373 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2378 this->mpi_communicator,
2395 std::integral_constant<int, 2>)
2397 const unsigned int dim = 2, spacedim = 3;
2405 std::vector<unsigned int> vertex_touch_count;
2407 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2410 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2411 const ::internal::p4est::types<2>::locidx num_vtt =
2412 std::accumulate(vertex_touch_count.begin(),
2413 vertex_touch_count.end(),
2419 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2422 (set_vertex_info ==
true ? this->n_vertices() : 0),
2427 set_vertex_and_cell_info(*
this,
2430 coarse_cell_to_p4est_tree_permutation,
2434 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2439 this->mpi_communicator,
2454 std::integral_constant<int, 3>)
2456 const int dim = 3, spacedim = 3;
2464 std::vector<unsigned int> vertex_touch_count;
2465 std::vector<std::list<
2466 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2468 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2469 const ::internal::p4est::types<2>::locidx num_vtt =
2470 std::accumulate(vertex_touch_count.begin(),
2471 vertex_touch_count.end(),
2474 std::vector<unsigned int> edge_touch_count;
2475 std::vector<std::list<
2476 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2478 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2479 const ::internal::p4est::types<2>::locidx num_ett =
2480 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2483 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2486 (set_vertex_info ==
true ? this->n_vertices() : 0),
2488 this->n_active_lines(),
2493 set_vertex_and_cell_info(*
this,
2496 coarse_cell_to_p4est_tree_permutation,
2525 const unsigned int deal_to_p4est_line_index[12] = {
2526 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2529 this->begin_active();
2530 cell != this->
end();
2533 const unsigned int index =
2534 coarse_cell_to_p4est_tree_permutation[cell->
index()];
2535 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2537 deal_to_p4est_line_index[e]] =
2538 cell->
line(e)->index();
2543 connectivity->ett_offset[0] = 0;
2544 std::partial_sum(edge_touch_count.begin(),
2545 edge_touch_count.end(),
2546 &connectivity->ett_offset[1]);
2548 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2551 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2553 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2557 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2558 unsigned int>>::const_iterator p =
2559 edge_to_cell[v].begin();
2560 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2562 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2563 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2564 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2565 deal_to_p4est_line_index[p->second];
2569 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2573 parallel_forest = ::internal::p4est::functions<3>::new_forest(
2574 this->mpi_communicator,
2591 template <
int dim,
int spacedim>
2593 enforce_mesh_balance_over_periodic_boundaries(
2599 std::vector<bool> flags_before[2];
2603 std::vector<unsigned int> topological_vertex_numbering(
2605 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2606 topological_vertex_numbering[i] = i;
2624 using cell_iterator =
2628 const cell_iterator &cell_1 = it.first.first;
2629 const unsigned int face_no_1 = it.first.second;
2630 const cell_iterator &cell_2 = it.second.first.first;
2631 const unsigned int face_no_2 = it.second.first.second;
2632 const unsigned char combined_orientation = it.second.second;
2634 if (cell_1->level() == cell_2->level())
2636 for (
const unsigned int v :
2637 cell_1->face(face_no_1)->vertex_indices())
2641 const unsigned int vface1 =
2642 cell_1->reference_cell().standard_to_real_face_vertex(
2643 v, face_no_1, combined_orientation);
2644 const unsigned int vi1 =
2645 topological_vertex_numbering[cell_1->face(face_no_1)
2646 ->vertex_index(vface1)];
2647 const unsigned int vi2 =
2648 topological_vertex_numbering[cell_2->face(face_no_2)
2650 const unsigned int min_index =
std::min(vi1, vi2);
2651 topological_vertex_numbering[cell_1->face(face_no_1)
2652 ->vertex_index(vface1)] =
2653 topological_vertex_numbering[cell_2->face(face_no_2)
2654 ->vertex_index(v)] =
2661 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2663 const unsigned int j = topological_vertex_numbering[i];
2664 Assert(j == i || topological_vertex_numbering[j] == j,
2665 ExcMessage(
"Got inconclusive constraints with chain: " +
2666 std::to_string(i) +
" vs " + std::to_string(j) +
2667 " which should be equal to " +
2668 std::to_string(topological_vertex_numbering[j])));
2675 bool continue_iterating =
true;
2676 std::vector<int> vertex_level(tria.
n_vertices());
2677 while (continue_iterating)
2681 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2685 for (; cell != endc; ++cell)
2687 if (cell->refine_flag_set())
2688 for (
const unsigned int vertex :
2690 vertex_level[topological_vertex_numbering
2691 [cell->vertex_index(vertex)]] =
2692 std::max(vertex_level[topological_vertex_numbering
2693 [cell->vertex_index(vertex)]],
2695 else if (!cell->coarsen_flag_set())
2696 for (
const unsigned int vertex :
2698 vertex_level[topological_vertex_numbering
2699 [cell->vertex_index(vertex)]] =
2700 std::max(vertex_level[topological_vertex_numbering
2701 [cell->vertex_index(vertex)]],
2712 for (
const unsigned int vertex :
2714 vertex_level[topological_vertex_numbering
2715 [cell->vertex_index(vertex)]] =
2716 std::max(vertex_level[topological_vertex_numbering
2717 [cell->vertex_index(vertex)]],
2722 continue_iterating =
false;
2733 for (cell = tria.
last_active(); cell != endc; --cell)
2734 if (cell->refine_flag_set() ==
false)
2736 for (
const unsigned int vertex :
2738 if (vertex_level[topological_vertex_numbering
2739 [cell->vertex_index(vertex)]] >=
2743 cell->clear_coarsen_flag();
2748 if (vertex_level[topological_vertex_numbering
2749 [cell->vertex_index(vertex)]] >
2752 cell->set_refine_flag();
2753 continue_iterating =
true;
2755 for (
const unsigned int v :
2757 vertex_level[topological_vertex_numbering
2758 [cell->vertex_index(v)]] =
2760 vertex_level[topological_vertex_numbering
2761 [cell->vertex_index(v)]],
2775 if (cell->is_active())
2778 const unsigned int n_children = cell->n_children();
2779 unsigned int flagged_children = 0;
2780 for (
unsigned int child = 0; child < n_children; ++child)
2781 if (cell->child(child)->is_active() &&
2782 cell->child(child)->coarsen_flag_set())
2787 if (flagged_children < n_children)
2788 for (
unsigned int child = 0; child < n_children; ++child)
2789 if (cell->child(child)->is_active())
2790 cell->child(child)->clear_coarsen_flag();
2793 std::vector<bool> flags_after[2];
2796 return ((flags_before[0] != flags_after[0]) ||
2797 (flags_before[1] != flags_after[1]));
2803 template <
int dim,
int spacedim>
2822 template <
int dim,
int spacedim>
2850 bool mesh_changed =
false;
2866 for (
const auto &cell :
2869 cell->set_coarsen_flag();
2877 const typename Triangulation<dim, spacedim>::DistortedCellList &)
2895 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2896 P4EST_CONNECT_CORNER) :
2897 typename ::internal::p4est::types<dim>::balance_type(
2898 P8EST_CONNECT_CORNER)));
2915 if (tree_exists_locally<dim, spacedim>(
2920 delete_all_children<dim, spacedim>(cell);
2921 if (cell->is_active())
2930 typename ::internal::p4est::types<dim>::quadrant
2932 typename ::internal::p4est::types<dim>::tree *tree =
2938 match_tree_recursively<dim, spacedim>(*tree,
2942 this->my_subdomain);
2949 typename ::internal::p4est::types<dim>::quadrant *quadr;
2951 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2953 for (
unsigned int g_idx = 0;
2957 while (g_idx >=
static_cast<unsigned int>(
2960 while (g_idx >=
static_cast<unsigned int>(
2964 quadr =
static_cast<
2965 typename ::internal::p4est::types<dim>::quadrant *
>(
2968 unsigned int coarse_cell_index =
2971 match_quadrant<dim, spacedim>(
this,
2986 bool mesh_changed =
true;
2987 unsigned int loop_counter = 0;
2997 enforce_mesh_balance_over_periodic_boundaries(*
this);
3009 "parallel::distributed::Triangulation::copy_local_forest_to_triangulation() "
3010 "for periodic boundaries detected. Aborting."));
3012 while (mesh_changed);
3020 return cell.refine_flag_set() ||
3021 cell.coarsen_flag_set();
3032 const typename Triangulation<dim, spacedim>::DistortedCellList &)
3039 while (mesh_changed);
3043 unsigned int num_ghosts = 0;
3047 if (cell->subdomain_id() != this->my_subdomain &&
3069 for (
unsigned int lvl = this->
n_levels(); lvl > 0;)
3074 if ((cell->is_active() &&
3075 cell->subdomain_id() ==
3076 this->locally_owned_subdomain()) ||
3077 (cell->has_children() &&
3078 cell->child(0)->level_subdomain_id() ==
3079 this->locally_owned_subdomain()))
3080 cell->set_level_subdomain_id(
3085 cell->set_level_subdomain_id(
3093 std::vector<std::vector<bool>> marked_vertices(this->
n_levels());
3094 for (
unsigned int lvl = 0; lvl < this->
n_levels(); ++lvl)
3099 typename ::internal::p4est::types<dim>::quadrant
3101 const unsigned int tree_index =
3103 typename ::internal::p4est::types<dim>::tree *tree =
3109 determine_level_subdomain_id_recursively<dim, spacedim>(
3120 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3125 if (cell->has_children())
3126 for (
unsigned int c = 0;
3127 c < GeometryInfo<dim>::max_children_per_cell;
3130 if (cell->child(c)->level_subdomain_id() ==
3131 this->locally_owned_subdomain())
3136 cell->child(0)->level_subdomain_id();
3140 cell->set_level_subdomain_id(mark);
3157 (void)total_local_cells;
3161 Assert(
static_cast<unsigned int>(
3167 Assert(
static_cast<unsigned int>(
3174 unsigned int n_owned = 0;
3177 if (cell->subdomain_id() == this->my_subdomain)
3181 Assert(
static_cast<unsigned int>(
3198 template <
int dim,
int spacedim>
3204 std::vector<Point<dim>> point{p};
3212 template <
int dim,
int spacedim>
3217# ifndef P4EST_SEARCH_LOCAL
3222 "This function is only available if p4est is version 2.2 and higher."));
3224 return std::vector<unsigned int>(1,
3230 "Vertices need to be communicated to p4est to use this "
3231 "function. This must explicitly be turned on in the "
3232 "settings of the triangulation's constructor."));
3240 "This function can only be used if the triangulation "
3241 "has no other manifold than a Cartesian (flat) manifold attached."));
3245 PartitionSearch<dim> partition_search;
3257 sc_array_t *point_sc_array;
3261 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3264 for (
size_t i = 0; i < points.size(); ++i)
3269 double *this_sc_point =
3270 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3272 for (
unsigned int d = 0; d < dim; ++d)
3274 this_sc_point[d] = p(d);
3276 this_sc_point[dim] = -1.0;
3282 static_cast<int>(
false),
3283 &PartitionSearch<dim>::local_quadrant_fn,
3284 &PartitionSearch<dim>::local_point_fn,
3288 std::vector<types::subdomain_id> owner_rank(
3292 for (
size_t i = 0; i < points.size(); ++i)
3295 double *this_sc_point =
3296 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3297 Assert(this_sc_point[dim] >= 0. || this_sc_point[dim] == -1.,
3299 if (this_sc_point[dim] < 0.)
3310 sc_array_destroy_null(&point_sc_array);
3318 template <
int dim,
int spacedim>
3325 if (cell->is_locally_owned() && cell->refine_flag_set())
3326 Assert(cell->refine_flag_set() ==
3329 "This class does not support anisotropic refinement"));
3346 !(cell->refine_flag_set()),
3348 "Fatal Error: maximum refinement level of p4est reached."));
3355 this->
signals.pre_distributed_refinement();
3361 if (cell->is_ghost() || cell->is_artificial())
3363 cell->clear_refine_flag();
3364 cell->clear_coarsen_flag();
3370 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3389 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3394 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3407 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3408 P4EST_CONNECT_FULL) :
3409 typename ::internal::p4est::types<dim>::balance_type(
3410 P8EST_CONNECT_FULL)),
3419 this->
signals.post_p4est_refinement();
3423 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3424 previous_global_first_quadrant;
3429 std::memcpy(previous_global_first_quadrant.data(),
3432 typename ::internal::p4est::types<dim>::gloidx) *
3440 if (this->
signals.weight.empty())
3454 this->mpi_communicator) > 0,
3456 "The global sum of weights over all active cells "
3457 "is zero. Please verify how you generate weights."));
3459 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3470 &PartitionWeights<dim, spacedim>::cell_weight);
3486 this->cell_attached_data.pack_callbacks_variable,
3487 this->get_communicator());
3495 cell->clear_refine_flag();
3496 cell->clear_coarsen_flag();
3503 catch (
const typename Triangulation<dim>::DistortedCellList &)
3514 previous_global_first_quadrant.data());
3540 std::vector<bool> active_verts =
3543 const unsigned int maybe_coarser_lvl =
3544 (lvl > 0) ? (lvl - 1) : lvl;
3546 cell = this->
begin(maybe_coarser_lvl),
3547 endc = this->
end(lvl);
3548 for (; cell != endc; ++cell)
3549 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3551 const bool is_level_artificial =
3552 (cell->level_subdomain_id() ==
3554 bool need_to_know =
false;
3555 for (
const unsigned int vertex :
3557 if (active_verts[cell->vertex_index(vertex)])
3559 need_to_know =
true;
3564 !need_to_know || !is_level_artificial,
3566 "Internal error: the owner of cell" +
3567 cell->id().to_string() +
3568 " is unknown even though it is needed for geometric multigrid."));
3578 this->
signals.post_distributed_refinement();
3583 template <
int dim,
int spacedim>
3589 if (cell->is_locally_owned())
3591 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3593 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3597 this->
signals.pre_distributed_repartition();
3601 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3602 previous_global_first_quadrant;
3606 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3607 std::memcpy(previous_global_first_quadrant.data(),
3608 parallel_forest->global_first_quadrant,
3610 typename ::internal::p4est::types<dim>::gloidx) *
3611 (parallel_forest->mpisize + 1));
3614 if (this->
signals.weight.empty())
3618 ::internal::p4est::functions<dim>::partition(
3626 const std::vector<unsigned int> cell_weights = get_cell_weights();
3629 Assert(Utilities::MPI::sum(std::accumulate(cell_weights.begin(),
3632 this->mpi_communicator) > 0,
3634 "The global sum of weights over all active cells "
3635 "is zero. Please verify how you generate weights."));
3637 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3641 Assert(parallel_forest->user_pointer == this, ExcInternalError());
3642 parallel_forest->user_pointer = &partition_weights;
3644 ::internal::p4est::functions<dim>::partition(
3648 &PartitionWeights<dim, spacedim>::cell_weight);
3651 parallel_forest->user_pointer = this;
3657 this->data_serializer.pack_data(
3658 this->local_cell_relations,
3659 this->cell_attached_data.pack_callbacks_fixed,
3660 this->cell_attached_data.pack_callbacks_variable,
3661 this->get_communicator());
3666 copy_local_forest_to_triangulation();
3668 catch (
const typename Triangulation<dim>::DistortedCellList &)
3678 this->execute_transfer(parallel_forest,
3679 previous_global_first_quadrant.data());
3682 this->update_periodic_face_map();
3685 this->update_number_cache();
3688 this->signals.post_distributed_repartition();
3693 template <
int dim,
int spacedim>
3695 const std::vector<types::global_dof_index>
3704 template <
int dim,
int spacedim>
3706 const std::vector<types::global_dof_index>
3715 template <
int dim,
int spacedim>
3722 std::vector<bool> marked_vertices(this->
n_vertices(),
false);
3724 if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3726 marked_vertices[cell->vertex_index(v)] =
true;
3743 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3747 const unsigned int face_no_1 = it.first.second;
3749 const unsigned int face_no_2 = it.second.first.second;
3750 const unsigned char combined_orientation = it.second.second;
3751 const auto [orientation, rotation, flip] =
3754 if (cell_1->level() == level && cell_2->level() == level)
3756 for (
unsigned int v = 0;
3762 const unsigned int vface0 =
3764 v, orientation, flip, rotation);
3765 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3767 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3769 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3771 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3777 return marked_vertices;
3782 template <
int dim,
int spacedim>
3786 const
types::coarse_cell_id coarse_cell_id)
const
3793 template <
int dim,
int spacedim>
3797 const unsigned int coarse_cell_index)
const