841 , is_reference_mapping_initialized(false)
848 PartitionSearch<dim>::QuadrantData::initialize_mapping()
851 are_vertices_initialized,
853 "Cell vertices must be initialized before the cell mapping can be filled."));
860 for (
unsigned int alpha = 0;
861 alpha < GeometryInfo<dim>::vertices_per_cell;
865 point_matrix(0, alpha) = 1;
866 point_matrix(1, alpha) = cell_vertices[alpha](0);
867 point_matrix(2, alpha) = cell_vertices[alpha](1);
868 point_matrix(3, alpha) =
869 cell_vertices[alpha](0) * cell_vertices[alpha](1);
876 quadrant_mapping_matrix.invert(point_matrix);
880 for (
unsigned int alpha = 0;
881 alpha < GeometryInfo<dim>::vertices_per_cell;
885 point_matrix(0, alpha) = 1;
886 point_matrix(1, alpha) = cell_vertices[alpha](0);
887 point_matrix(2, alpha) = cell_vertices[alpha](1);
888 point_matrix(3, alpha) = cell_vertices[alpha](2);
889 point_matrix(4, alpha) =
890 cell_vertices[alpha](0) * cell_vertices[alpha](1);
891 point_matrix(5, alpha) =
892 cell_vertices[alpha](1) * cell_vertices[alpha](2);
893 point_matrix(6, alpha) =
894 cell_vertices[alpha](0) * cell_vertices[alpha](2);
895 point_matrix(7, alpha) = cell_vertices[alpha](0) *
896 cell_vertices[alpha](1) *
897 cell_vertices[alpha](2);
904 quadrant_mapping_matrix.invert(point_matrix);
907 is_reference_mapping_initialized =
true;
914 PartitionSearch<2>::QuadrantData::set_cell_vertices(
919 quad_length_on_level)
921 constexpr unsigned int dim = 2;
925 double corner_point[dim + 1] = {0};
928 const auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
930 for (
unsigned int d = 0;
d < dim; ++
d)
932 cell_vertices[vertex_index](
d) = corner_point[d];
942 unsigned int vertex_index = 0;
944 forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
947 copy_vertex(vertex_index);
954 forest->connectivity,
956 quadrant->x + quad_length_on_level,
961 copy_vertex(vertex_index);
968 forest->connectivity,
971 quadrant->y + quad_length_on_level,
975 copy_vertex(vertex_index);
982 forest->connectivity,
984 quadrant->x + quad_length_on_level,
985 quadrant->y + quad_length_on_level,
989 copy_vertex(vertex_index);
991 are_vertices_initialized =
true;
998 PartitionSearch<3>::QuadrantData::set_cell_vertices(
1003 quad_length_on_level)
1005 constexpr unsigned int dim = 3;
1007 double corner_point[dim] = {0};
1010 auto copy_vertex = [&](
unsigned int vertex_index) ->
void {
1012 for (
unsigned int d = 0;
d < dim; ++
d)
1014 cell_vertices[vertex_index](
d) = corner_point[d];
1016 corner_point[
d] = 0;
1024 unsigned int vertex_index = 0;
1026 forest->connectivity,
1034 copy_vertex(vertex_index);
1042 forest->connectivity,
1044 quadrant->x + quad_length_on_level,
1050 copy_vertex(vertex_index);
1057 forest->connectivity,
1060 quadrant->y + quad_length_on_level,
1065 copy_vertex(vertex_index);
1072 forest->connectivity,
1074 quadrant->x + quad_length_on_level,
1075 quadrant->y + quad_length_on_level,
1080 copy_vertex(vertex_index);
1087 forest->connectivity,
1091 quadrant->z + quad_length_on_level,
1095 copy_vertex(vertex_index);
1102 forest->connectivity,
1104 quadrant->x + quad_length_on_level,
1106 quadrant->z + quad_length_on_level,
1110 copy_vertex(vertex_index);
1117 forest->connectivity,
1120 quadrant->y + quad_length_on_level,
1121 quadrant->z + quad_length_on_level,
1125 copy_vertex(vertex_index);
1132 forest->connectivity,
1134 quadrant->x + quad_length_on_level,
1135 quadrant->y + quad_length_on_level,
1136 quadrant->z + quad_length_on_level,
1140 copy_vertex(vertex_index);
1143 are_vertices_initialized =
true;
1153 template <
int dim,
int spacedim>
1154 class RefineAndCoarsenList
1158 const std::vector<types::global_dof_index>
1159 &p4est_tree_to_coarse_cell_permutation,
1187 pointers_are_at_end()
const;
1190 std::vector<typename internal::p4est::types<dim>::quadrant> refine_list;
1191 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1192 const_iterator current_refine_pointer;
1194 std::vector<typename internal::p4est::types<dim>::quadrant> coarsen_list;
1195 typename std::vector<typename internal::p4est::types<dim>::quadrant>::
1196 const_iterator current_coarsen_pointer;
1207 template <
int dim,
int spacedim>
1209 RefineAndCoarsenList<dim, spacedim>::pointers_are_at_end()
const
1211 return ((current_refine_pointer == refine_list.end()) &&
1212 (current_coarsen_pointer == coarsen_list.end()));
1217 template <
int dim,
int spacedim>
1218 RefineAndCoarsenList<dim, spacedim>::RefineAndCoarsenList(
1220 const std::vector<types::global_dof_index>
1221 & p4est_tree_to_coarse_cell_permutation,
1225 unsigned int n_refine_flags = 0, n_coarsen_flags = 0;
1226 for (
const auto &cell :
triangulation.active_cell_iterators())
1229 if (cell->subdomain_id() != my_subdomain)
1232 if (cell->refine_flag_set())
1234 else if (cell->coarsen_flag_set())
1238 refine_list.reserve(n_refine_flags);
1239 coarsen_list.reserve(n_coarsen_flags);
1251 unsigned int coarse_cell_index =
1252 p4est_tree_to_coarse_cell_permutation[c];
1261 p4est_cell.p.which_tree = c;
1262 build_lists(cell, p4est_cell, my_subdomain);
1270 for (
unsigned int i = 1; i < refine_list.size(); ++i)
1271 Assert(refine_list[i].p.which_tree >= refine_list[i - 1].p.which_tree,
1273 for (
unsigned int i = 1; i < coarsen_list.size(); ++i)
1274 Assert(coarsen_list[i].p.which_tree >= coarsen_list[i - 1].p.which_tree,
1277 current_refine_pointer = refine_list.begin();
1278 current_coarsen_pointer = coarsen_list.begin();
1283 template <
int dim,
int spacedim>
1285 RefineAndCoarsenList<dim, spacedim>::build_lists(
1290 if (cell->is_active())
1292 if (cell->subdomain_id() == my_subdomain)
1294 if (cell->refine_flag_set())
1295 refine_list.push_back(p4est_cell);
1296 else if (cell->coarsen_flag_set())
1297 coarsen_list.push_back(p4est_cell);
1304 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1309 P4EST_QUADRANT_INIT(&p4est_child[c]);
1312 P8EST_QUADRANT_INIT(&p4est_child[c]);
1319 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1322 p4est_child[c].p.which_tree = p4est_cell.p.which_tree;
1323 build_lists(cell->child(c), p4est_child[c], my_subdomain);
1329 template <
int dim,
int spacedim>
1331 RefineAndCoarsenList<dim, spacedim>::refine_callback(
1336 RefineAndCoarsenList<dim, spacedim> *this_object =
1337 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1338 forest->user_pointer);
1342 if (this_object->current_refine_pointer == this_object->refine_list.end())
1345 Assert(coarse_cell_index <=
1346 this_object->current_refine_pointer->p.which_tree,
1351 if (coarse_cell_index < this_object->current_refine_pointer->p.which_tree)
1355 Assert(coarse_cell_index <=
1356 this_object->current_refine_pointer->p.which_tree,
1362 quadrant, &*this_object->current_refine_pointer) <= 0,
1367 quadrant, &*this_object->current_refine_pointer))
1369 ++this_object->current_refine_pointer;
1379 template <
int dim,
int spacedim>
1381 RefineAndCoarsenList<dim, spacedim>::coarsen_callback(
1386 RefineAndCoarsenList<dim, spacedim> *this_object =
1387 reinterpret_cast<RefineAndCoarsenList<dim, spacedim> *
>(
1388 forest->user_pointer);
1392 if (this_object->current_coarsen_pointer == this_object->coarsen_list.end())
1395 Assert(coarse_cell_index <=
1396 this_object->current_coarsen_pointer->p.which_tree,
1401 if (coarse_cell_index < this_object->current_coarsen_pointer->p.which_tree)
1405 Assert(coarse_cell_index <=
1406 this_object->current_coarsen_pointer->p.which_tree,
1412 children[0], &*this_object->current_coarsen_pointer) <= 0,
1418 children[0], &*this_object->current_coarsen_pointer))
1421 ++this_object->current_coarsen_pointer;
1425 for (
unsigned int c = 1; c < GeometryInfo<dim>::max_children_per_cell;
1429 children[c], &*this_object->current_coarsen_pointer),
1431 ++this_object->current_coarsen_pointer;
1449 template <
int dim,
int spacedim>
1450 class PartitionWeights
1458 explicit PartitionWeights(
const std::vector<unsigned int> &cell_weights);
1473 std::vector<unsigned int> cell_weights_list;
1474 std::vector<unsigned int>::const_iterator current_pointer;
1478 template <
int dim,
int spacedim>
1479 PartitionWeights<dim, spacedim>::PartitionWeights(
1480 const std::vector<unsigned int> &cell_weights)
1481 : cell_weights_list(cell_weights)
1485 current_pointer = cell_weights_list.begin();
1489 template <
int dim,
int spacedim>
1491 PartitionWeights<dim, spacedim>::cell_weight(
1500 PartitionWeights<dim, spacedim> *this_object =
1501 reinterpret_cast<PartitionWeights<dim, spacedim> *
>(forest->user_pointer);
1503 Assert(this_object->current_pointer >=
1504 this_object->cell_weights_list.begin(),
1506 Assert(this_object->current_pointer < this_object->cell_weights_list.end(),
1512 const unsigned int weight = *this_object->current_pointer;
1513 ++this_object->current_pointer;
1515 Assert(weight <
static_cast<unsigned int>(std::numeric_limits<int>::max()),
1516 ExcMessage(
"p4est uses 'signed int' to represent the partition "
1517 "weights for cells. The weight provided here exceeds "
1518 "the maximum value represented as a 'signed int'."));
1519 return static_cast<int>(weight);
1522 template <
int dim,
int spacedim>
1523 using cell_relation_t =
typename std::pair<
1524 typename ::Triangulation<dim, spacedim>::cell_iterator,
1525 typename ::Triangulation<dim, spacedim>::CellStatus>;
1536 template <
int dim,
int spacedim>
1538 add_single_cell_relation(
1539 std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
1540 const typename ::internal::p4est::types<dim>::tree & tree,
1541 const unsigned int idx,
1545 const unsigned int local_quadrant_index = tree.quadrants_offset + idx;
1551 cell_rel[local_quadrant_index] = std::make_pair(dealii_cell, status);
1565 template <
int dim,
int spacedim>
1567 update_cell_relations_recursively(
1568 std::vector<cell_relation_t<dim, spacedim>> & cell_rel,
1569 const typename ::internal::p4est::types<dim>::tree & tree,
1571 const typename ::internal::p4est::types<dim>::quadrant &p4est_cell)
1574 const int idx = sc_array_bsearch(
1575 const_cast<sc_array_t *
>(&tree.quadrants),
1580 const_cast<typename ::internal::p4est::types<dim>::tree *
>(
1582 &p4est_cell) ==
false))
1587 const bool p4est_has_children = (idx == -1);
1588 if (p4est_has_children && dealii_cell->has_children())
1591 typename ::internal::p4est::types<dim>::quadrant
1594 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1599 P4EST_QUADRANT_INIT(&p4est_child[c]);
1602 P8EST_QUADRANT_INIT(&p4est_child[c]);
1609 &p4est_cell, p4est_child);
1611 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1614 update_cell_relations_recursively<dim, spacedim>(
1615 cell_rel, tree, dealii_cell->child(c), p4est_child[c]);
1618 else if (!p4est_has_children && !dealii_cell->has_children())
1622 add_single_cell_relation<dim, spacedim>(
1629 else if (p4est_has_children)
1637 typename ::internal::p4est::types<dim>::quadrant
1639 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1644 P4EST_QUADRANT_INIT(&p4est_child[c]);
1647 P8EST_QUADRANT_INIT(&p4est_child[c]);
1654 &p4est_cell, p4est_child);
1661 for (
unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_cell;
1664 child_idx = sc_array_bsearch(
1665 const_cast<sc_array_t *
>(&tree.quadrants),
1672 add_single_cell_relation<dim, spacedim>(
1673 cell_rel, tree, child_idx, dealii_cell, cell_status);
1681 add_single_cell_relation<dim, spacedim>(
1695 namespace distributed
1698 template <
int dim,
int spacedim>
1710 (
settings & construct_multigrid_hierarchy) ?
1714 Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
1718 , triangulation_has_content(false)
1719 , connectivity(
nullptr)
1720 , parallel_forest(
nullptr)
1722 parallel_ghost =
nullptr;
1727 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();
1788 catch (
const typename Triangulation<dim>::DistortedCellList &)
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>
1861 template <
int dim,
int spacedim>
1872 template <
int dim,
int spacedim>
1878 *previous_global_first_quadrant)
1880 Assert(this->data_transfer.sizes_fixed_cumulative.size() > 0,
1884 this->data_transfer.dest_data_fixed.resize(
1885 parallel_forest->local_num_quadrants *
1886 this->data_transfer.sizes_fixed_cumulative.back());
1889 typename ::internal::p4est::types<dim>::transfer_context
1893 parallel_forest->global_first_quadrant,
1894 previous_global_first_quadrant,
1895 parallel_forest->mpicomm,
1897 this->data_transfer.dest_data_fixed.data(),
1898 this->data_transfer.src_data_fixed.data(),
1899 this->data_transfer.sizes_fixed_cumulative.back());
1901 if (this->data_transfer.variable_size_data_stored)
1904 this->data_transfer.dest_sizes_variable.resize(
1905 parallel_forest->local_num_quadrants);
1910 parallel_forest->global_first_quadrant,
1911 previous_global_first_quadrant,
1912 parallel_forest->mpicomm,
1914 this->data_transfer.dest_sizes_variable.data(),
1915 this->data_transfer.src_sizes_variable.data(),
1916 sizeof(
unsigned int));
1922 this->data_transfer.src_data_fixed.clear();
1923 this->data_transfer.src_data_fixed.shrink_to_fit();
1925 if (this->data_transfer.variable_size_data_stored)
1928 this->data_transfer.dest_data_variable.resize(
1929 std::accumulate(this->data_transfer.dest_sizes_variable.begin(),
1930 this->data_transfer.dest_sizes_variable.end(),
1931 std::vector<int>::size_type(0)));
1933# if DEAL_II_P4EST_VERSION_GTE(2, 0, 65, 0)
1940 if (this->data_transfer.src_sizes_variable.size() == 0)
1941 this->data_transfer.src_sizes_variable.resize(1);
1942 if (this->data_transfer.dest_sizes_variable.size() == 0)
1943 this->data_transfer.dest_sizes_variable.resize(1);
1948 parallel_forest->global_first_quadrant,
1949 previous_global_first_quadrant,
1950 parallel_forest->mpicomm,
1952 this->data_transfer.dest_data_variable.data(),
1953 this->data_transfer.dest_sizes_variable.data(),
1954 this->data_transfer.src_data_variable.data(),
1955 this->data_transfer.src_sizes_variable.data());
1958 this->data_transfer.src_sizes_variable.clear();
1959 this->data_transfer.src_sizes_variable.shrink_to_fit();
1960 this->data_transfer.src_data_variable.clear();
1961 this->data_transfer.src_data_variable.shrink_to_fit();
1967 template <
int dim,
int spacedim>
1970 spacedim>::setup_coarse_cell_to_p4est_tree_permutation()
1975 coarse_cell_to_p4est_tree_permutation.resize(this->n_cells(0));
1977 cell_connectivity, coarse_cell_to_p4est_tree_permutation);
1979 p4est_tree_to_coarse_cell_permutation =
1985 template <
int dim,
int spacedim>
1988 const
std::
string &file_basename)
const
1990 Assert(parallel_forest !=
nullptr,
1991 ExcMessage(
"Can't produce output when no forest is created yet."));
1995 "To use this function the triangulation's flag "
1996 "Settings::communicate_vertices_to_p4est must be set."));
1999 parallel_forest,
nullptr, file_basename.c_str());
2004 template <
int dim,
int spacedim>
2009 this->cell_attached_data.n_attached_deserialize == 0,
2011 "Not all SolutionTransfer objects have been deserialized after the last call to load()."));
2012 Assert(this->n_cells() > 0,
2013 ExcMessage(
"Can not save() an empty Triangulation."));
2019 this->signals.pre_distributed_save();
2021 if (this->my_subdomain == 0)
2023 std::string fname = std::string(filename) +
".info";
2024 std::ofstream f(fname.c_str());
2025 f <<
"version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells"
2029 << this->cell_attached_data.pack_callbacks_fixed.size() <<
" "
2030 << this->cell_attached_data.pack_callbacks_variable.size() <<
" "
2031 << this->n_cells(0) << std::endl;
2035 for (
const auto &cell_rel : this->local_cell_relations)
2045 this->save_attached_data(parallel_forest->global_first_quadrant[myrank],
2046 parallel_forest->global_num_quadrants,
2054 this->signals.post_distributed_save();
2059 template <
int dim,
int spacedim>
2064 this->n_cells() > 0,
2066 "load() only works if the Triangulation already contains a coarse mesh!"));
2068 this->n_levels() == 1,
2070 "Triangulation may only contain coarse cells when calling load()."));
2076 this->signals.pre_distributed_load();
2078 if (parallel_ghost !=
nullptr)
2082 parallel_ghost =
nullptr;
2085 parallel_forest =
nullptr;
2088 connectivity =
nullptr;
2090 unsigned int version, numcpus, attached_count_fixed,
2091 attached_count_variable, n_coarse_cells;
2093 std::string fname = std::string(filename) +
".info";
2094 std::ifstream f(fname.c_str());
2096 std::string firstline;
2097 getline(f, firstline);
2098 f >> version >> numcpus >> attached_count_fixed >>
2099 attached_count_variable >> n_coarse_cells;
2103 ExcMessage(
"Incompatible version found in .info file."));
2104 Assert(this->n_cells(0) == n_coarse_cells,
2105 ExcMessage(
"Number of coarse cells differ!"));
2109 this->cell_attached_data.n_attached_data_sets = 0;
2110 this->cell_attached_data.n_attached_deserialize =
2111 attached_count_fixed + attached_count_variable;
2115 this->mpi_communicator,
2133 copy_local_forest_to_triangulation();
2135 catch (
const typename Triangulation<dim>::DistortedCellList &)
2143 this->load_attached_data(parallel_forest->global_first_quadrant[myrank],
2144 parallel_forest->global_num_quadrants,
2145 parallel_forest->local_num_quadrants,
2147 attached_count_fixed,
2148 attached_count_variable);
2151 this->signals.post_distributed_load();
2153 this->update_periodic_face_map();
2154 this->update_number_cache();
2159 template <
int dim,
int spacedim>
2162 const
bool autopartition)
2164 (void)autopartition;
2170 template <
int dim,
int spacedim>
2173 const typename ::
internal::p4est::
types<dim>::forest *forest)
2175 Assert(this->n_cells() > 0,
2177 "load() only works if the Triangulation already contains "
2179 Assert(this->n_cells() == forest->trees->elem_count,
2181 "Coarse mesh of the Triangulation does not match the one "
2182 "of the provided forest!"));
2185 if (parallel_ghost !=
nullptr)
2189 parallel_ghost =
nullptr;
2192 parallel_forest =
nullptr;
2198 typename ::internal::p4est::types<dim>::forest *temp =
2199 const_cast<typename ::internal::p4est::types<dim>::forest *
>(
2203 parallel_forest->connectivity = connectivity;
2204 parallel_forest->user_pointer =
this;
2208 copy_local_forest_to_triangulation();
2210 catch (
const typename Triangulation<dim>::DistortedCellList &)
2217 this->update_periodic_face_map();
2218 this->update_number_cache();
2223 template <
int dim,
int spacedim>
2227 Assert(parallel_forest !=
nullptr,
2229 "Can't produce a check sum when no forest is created yet."));
2230 return ::internal::p4est::functions<dim>::checksum(parallel_forest);
2235 template <
int dim,
int spacedim>
2237 const typename ::internal::p4est::types<dim>::forest
2240 Assert(parallel_forest !=
nullptr,
2241 ExcMessage(
"The forest has not been allocated yet."));
2242 return parallel_forest;
2247 template <
int dim,
int spacedim>
2249 typename ::internal::p4est::types<dim>::tree
2251 const int dealii_coarse_cell_index)
const
2253 const unsigned int tree_index =
2254 coarse_cell_to_p4est_tree_permutation[dealii_coarse_cell_index];
2255 typename ::internal::p4est::types<dim>::tree *tree =
2256 static_cast<typename ::internal::p4est::types<dim>::tree *
>(
2257 sc_array_index(parallel_forest->trees, tree_index));
2270 std::integral_constant<int, 2>)
2272 const unsigned int dim = 2, spacedim = 2;
2280 std::vector<unsigned int> vertex_touch_count;
2282 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2285 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2286 const ::internal::p4est::types<2>::locidx num_vtt =
2287 std::accumulate(vertex_touch_count.begin(),
2288 vertex_touch_count.end(),
2294 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2297 (set_vertex_info ==
true ? this->n_vertices() : 0),
2302 set_vertex_and_cell_info(*
this,
2305 coarse_cell_to_p4est_tree_permutation,
2309 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2314 this->mpi_communicator,
2330 std::integral_constant<int, 2>)
2332 const unsigned int dim = 2, spacedim = 3;
2340 std::vector<unsigned int> vertex_touch_count;
2342 std::list<std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2345 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2346 const ::internal::p4est::types<2>::locidx num_vtt =
2347 std::accumulate(vertex_touch_count.begin(),
2348 vertex_touch_count.end(),
2354 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2357 (set_vertex_info ==
true ? this->n_vertices() : 0),
2362 set_vertex_and_cell_info(*
this,
2365 coarse_cell_to_p4est_tree_permutation,
2369 Assert(p4est_connectivity_is_valid(connectivity) == 1,
2374 this->mpi_communicator,
2388 std::integral_constant<int, 3>)
2390 const int dim = 3, spacedim = 3;
2398 std::vector<unsigned int> vertex_touch_count;
2399 std::vector<std::list<
2400 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2402 get_vertex_to_cell_mappings(*
this, vertex_touch_count, vertex_to_cell);
2403 const ::internal::p4est::types<2>::locidx num_vtt =
2404 std::accumulate(vertex_touch_count.begin(),
2405 vertex_touch_count.end(),
2408 std::vector<unsigned int> edge_touch_count;
2409 std::vector<std::list<
2410 std::pair<Triangulation<3>::active_cell_iterator,
unsigned int>>>
2412 get_edge_to_cell_mappings(*
this, edge_touch_count, edge_to_cell);
2413 const ::internal::p4est::types<2>::locidx num_ett =
2414 std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
2417 const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
2420 (set_vertex_info ==
true ? this->n_vertices() : 0),
2422 this->n_active_lines(),
2427 set_vertex_and_cell_info(*
this,
2430 coarse_cell_to_p4est_tree_permutation,
2459 const unsigned int deal_to_p4est_line_index[12] = {
2460 4, 5, 0, 1, 6, 7, 2, 3, 8, 9, 10, 11};
2463 this->begin_active();
2464 cell != this->
end();
2467 const unsigned int index =
2468 coarse_cell_to_p4est_tree_permutation[cell->index()];
2469 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_cell; ++
e)
2471 deal_to_p4est_line_index[e]] =
2472 cell->line(e)->index();
2477 connectivity->ett_offset[0] = 0;
2478 std::partial_sum(edge_touch_count.begin(),
2479 edge_touch_count.end(),
2480 &connectivity->ett_offset[1]);
2482 Assert(connectivity->ett_offset[this->n_active_lines()] == num_ett,
2485 for (
unsigned int v = 0; v < this->n_active_lines(); ++v)
2487 Assert(edge_to_cell[v].size() == edge_touch_count[v],
2491 std::pair<Triangulation<dim, spacedim>::active_cell_iterator,
2492 unsigned int>>::const_iterator p =
2493 edge_to_cell[v].begin();
2494 for (
unsigned int c = 0; c < edge_touch_count[v]; ++c, ++p)
2496 connectivity->edge_to_tree[connectivity->ett_offset[v] + c] =
2497 coarse_cell_to_p4est_tree_permutation[p->first->index()];
2498 connectivity->edge_to_edge[connectivity->ett_offset[v] + c] =
2499 deal_to_p4est_line_index[p->second];
2503 Assert(p8est_connectivity_is_valid(connectivity) == 1,
2508 this->mpi_communicator,
2525 template <
int dim,
int spacedim>
2527 enforce_mesh_balance_over_periodic_boundaries(
2533 std::vector<bool> flags_before[2];
2537 std::vector<unsigned int> topological_vertex_numbering(
2539 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2540 topological_vertex_numbering[i] = i;
2558 using cell_iterator =
2560 typename std::map<std::pair<cell_iterator, unsigned int>,
2561 std::pair<std::pair<cell_iterator, unsigned int>,
2562 std::bitset<3>>>::const_iterator it;
2567 const cell_iterator &cell_1 = it->first.first;
2568 const unsigned int face_no_1 = it->first.second;
2569 const cell_iterator &cell_2 = it->second.first.first;
2570 const unsigned int face_no_2 = it->second.first.second;
2571 const std::bitset<3> face_orientation = it->second.second;
2573 if (cell_1->level() == cell_2->level())
2575 for (
unsigned int v = 0;
2581 const unsigned int vface0 =
2584 face_orientation[0],
2585 face_orientation[1],
2586 face_orientation[2]);
2587 const unsigned int vi0 =
2588 topological_vertex_numbering[cell_1->face(face_no_1)
2589 ->vertex_index(vface0)];
2590 const unsigned int vi1 =
2591 topological_vertex_numbering[cell_2->face(face_no_2)
2593 const unsigned int min_index =
std::min(vi0, vi1);
2594 topological_vertex_numbering[cell_1->face(face_no_1)
2595 ->vertex_index(vface0)] =
2596 topological_vertex_numbering[cell_2->face(face_no_2)
2597 ->vertex_index(v)] =
2604 for (
unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
2606 const unsigned int j = topological_vertex_numbering[i];
2614 bool continue_iterating =
true;
2616 while (continue_iterating)
2620 std::fill(vertex_level.begin(), vertex_level.end(), 0);
2624 for (; cell != endc; ++cell)
2626 if (cell->refine_flag_set())
2627 for (
const unsigned int vertex :
2629 vertex_level[topological_vertex_numbering
2630 [cell->vertex_index(vertex)]] =
2631 std::max(vertex_level[topological_vertex_numbering
2632 [cell->vertex_index(vertex)]],
2634 else if (!cell->coarsen_flag_set())
2635 for (
const unsigned int vertex :
2637 vertex_level[topological_vertex_numbering
2638 [cell->vertex_index(vertex)]] =
2639 std::max(vertex_level[topological_vertex_numbering
2640 [cell->vertex_index(vertex)]],
2651 for (
const unsigned int vertex :
2653 vertex_level[topological_vertex_numbering
2654 [cell->vertex_index(vertex)]] =
2655 std::max(vertex_level[topological_vertex_numbering
2656 [cell->vertex_index(vertex)]],
2661 continue_iterating =
false;
2673 if (cell->refine_flag_set() ==
false)
2675 for (
const unsigned int vertex :
2677 if (vertex_level[topological_vertex_numbering
2678 [cell->vertex_index(vertex)]] >=
2682 cell->clear_coarsen_flag();
2687 if (vertex_level[topological_vertex_numbering
2688 [cell->vertex_index(vertex)]] >
2691 cell->set_refine_flag();
2692 continue_iterating =
true;
2694 for (
const unsigned int v :
2696 vertex_level[topological_vertex_numbering
2697 [cell->vertex_index(v)]] =
2699 vertex_level[topological_vertex_numbering
2700 [cell->vertex_index(v)]],
2711 for (
const auto &cell :
tria.cell_iterators())
2714 if (cell->is_active())
2717 const unsigned int n_children = cell->n_children();
2718 unsigned int flagged_children = 0;
2719 for (
unsigned int child = 0; child < n_children; ++child)
2720 if (cell->child(child)->is_active() &&
2721 cell->child(child)->coarsen_flag_set())
2726 if (flagged_children < n_children)
2727 for (
unsigned int child = 0; child < n_children; ++child)
2728 if (cell->child(child)->is_active())
2729 cell->child(child)->clear_coarsen_flag();
2732 std::vector<bool> flags_after[2];
2735 return ((flags_before[0] != flags_after[0]) ||
2736 (flags_before[1] != flags_after[1]));
2742 template <
int dim,
int spacedim>
2746 bool mesh_changed =
false;
2747 unsigned int loop_counter = 0;
2748 unsigned int n_changes = 0;
2753 this->update_periodic_face_map();
2755 mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*
this);
2756 n_changes += mesh_changed;
2767 "parallel::distributed::Triangulation::prepare_coarsening_and_refinement() "
2768 "for periodic boundaries detected. Aborting."));
2770 while (mesh_changed);
2773 return n_changes > 0;
2778 template <
int dim,
int spacedim>
2799 if (
settings & construct_multigrid_hierarchy)
2802 spacedim>::limit_level_difference_at_vertices;
2806 bool mesh_changed =
false;
2814 if (
settings & mesh_reconstruction_after_repartitioning)
2815 while (this->n_levels() > 1)
2822 for (
const auto &cell :
2823 this->active_cell_iterators_on_level(this->n_levels() - 1))
2825 cell->set_coarsen_flag();
2833 const typename Triangulation<dim, spacedim>::DistortedCellList &)
2843 if (parallel_ghost !=
nullptr)
2847 parallel_ghost =
nullptr;
2851 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
2852 P4EST_CONNECT_CORNER) :
2853 typename ::internal::p4est::types<dim>::balance_type(
2854 P8EST_CONNECT_CORNER)));
2861 for (
const auto &cell : this->cell_iterators_on_level(0))
2866 for (
const auto &cell : this->cell_iterators_on_level(0))
2871 if (tree_exists_locally<dim, spacedim>(
2873 coarse_cell_to_p4est_tree_permutation[cell->index()]) ==
2876 delete_all_children<dim, spacedim>(cell);
2877 if (cell->is_active())
2886 typename ::internal::p4est::types<dim>::quadrant
2888 typename ::internal::p4est::types<dim>::tree *tree =
2889 init_tree(cell->index());
2894 match_tree_recursively<dim, spacedim>(*tree,
2898 this->my_subdomain);
2905 typename ::internal::p4est::types<dim>::quadrant *quadr;
2907 typename ::internal::p4est::types<dim>::topidx ghost_tree = 0;
2909 for (
unsigned int g_idx = 0;
2910 g_idx < parallel_ghost->ghosts.elem_count;
2913 while (g_idx >=
static_cast<unsigned int>(
2914 parallel_ghost->proc_offsets[ghost_owner + 1]))
2916 while (g_idx >=
static_cast<unsigned int>(
2917 parallel_ghost->tree_offsets[ghost_tree + 1]))
2920 quadr =
static_cast<
2921 typename ::internal::p4est::types<dim>::quadrant *
>(
2922 sc_array_index(¶llel_ghost->ghosts, g_idx));
2924 unsigned int coarse_cell_index =
2925 p4est_tree_to_coarse_cell_permutation[ghost_tree];
2927 match_quadrant<dim, spacedim>(
this,
2934 this->prepare_coarsening_and_refinement();
2938 std::any_of(this->begin_active(),
2941 return cell.refine_flag_set() ||
2942 cell.coarsen_flag_set();
2953 const typename Triangulation<dim, spacedim>::DistortedCellList &)
2960 while (mesh_changed);
2964 unsigned int num_ghosts = 0;
2966 for (
const auto &cell : this->active_cell_iterators())
2968 if (cell->subdomain_id() != this->my_subdomain &&
2973 Assert(num_ghosts == parallel_ghost->ghosts.elem_count,
2984 if (
settings & construct_multigrid_hierarchy)
2990 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
2993 for (
const auto &cell : this->cell_iterators_on_level(lvl))
2995 if ((cell->is_active() &&
2996 cell->subdomain_id() ==
2997 this->locally_owned_subdomain()) ||
2998 (cell->has_children() &&
2999 cell->child(0)->level_subdomain_id() ==
3000 this->locally_owned_subdomain()))
3001 cell->set_level_subdomain_id(
3002 this->locally_owned_subdomain());
3006 cell->set_level_subdomain_id(
3014 std::vector<std::vector<bool>> marked_vertices(this->n_levels());
3015 for (
unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
3016 marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
3018 for (
const auto &cell : this->cell_iterators_on_level(0))
3020 typename ::internal::p4est::types<dim>::quadrant
3022 const unsigned int tree_index =
3023 coarse_cell_to_p4est_tree_permutation[cell->index()];
3024 typename ::internal::p4est::types<dim>::tree *tree =
3025 init_tree(cell->index());
3030 determine_level_subdomain_id_recursively<dim, spacedim>(
3041 for (
unsigned int lvl = this->n_levels(); lvl > 0;)
3044 for (
const auto &cell : this->cell_iterators_on_level(lvl))
3046 if (cell->has_children())
3047 for (
unsigned int c = 0;
3048 c < GeometryInfo<dim>::max_children_per_cell;
3051 if (cell->child(c)->level_subdomain_id() ==
3052 this->locally_owned_subdomain())
3057 cell->child(0)->level_subdomain_id();
3061 cell->set_level_subdomain_id(mark);
3077 const unsigned int total_local_cells = this->n_active_cells();
3078 (void)total_local_cells;
3082 Assert(
static_cast<unsigned int>(
3083 parallel_forest->local_num_quadrants) == total_local_cells,
3088 Assert(
static_cast<unsigned int>(
3089 parallel_forest->local_num_quadrants) <= total_local_cells,
3095 unsigned int n_owned = 0;
3096 for (
const auto &cell : this->active_cell_iterators())
3098 if (cell->subdomain_id() == this->my_subdomain)
3102 Assert(
static_cast<unsigned int>(
3103 parallel_forest->local_num_quadrants) == n_owned,
3108 this->smooth_grid = save_smooth;
3114 update_cell_relations();
3119 template <
int dim,
int spacedim>
3125 std::vector<Point<dim>> point{p};
3126 std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
3133 template <
int dim,
int spacedim>
3136 find_point_owner_rank(const
std::vector<
Point<dim>> &points)
3138# ifndef P4EST_SEARCH_LOCAL
3143 "This function is only available if p4est is version 2.2 and higher."));
3145 return std::vector<unsigned int>(1,
3149 AssertThrow(this->are_vertices_communicated_to_p4est(),
3151 "Vertices need to be communicated to p4est to use this "
3152 "function. This must explicitly be turned on in the "
3153 "settings of the triangulation's constructor."));
3156 for (
const auto &manifold_id : this->get_manifold_ids())
3161 "This function can only be used if the triangulation "
3162 "has no other manifold than a Cartesian (flat) manifold attached."));
3166 PartitionSearch<dim> partition_search;
3172 parallel_forest->user_pointer = &partition_search;
3178 sc_array_t *point_sc_array;
3182 sc_array_new_count(
sizeof(
double[dim + 1]), points.size());
3185 for (
size_t i = 0; i < points.size(); ++i)
3190 double *this_sc_point =
3191 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3193 for (
unsigned int d = 0; d < dim; ++d)
3195 this_sc_point[d] = p(d);
3197 this_sc_point[dim] = -1.0;
3203 static_cast<int>(
false),
3204 &PartitionSearch<dim>::local_quadrant_fn,
3205 &PartitionSearch<dim>::local_point_fn,
3209 std::vector<types::subdomain_id> owner_rank(
3213 for (
size_t i = 0; i < points.size(); ++i)
3216 double *this_sc_point =
3217 static_cast<double *
>(sc_array_index_ssize_t(point_sc_array, i));
3222 parallel_forest->user_pointer =
this;
3225 sc_array_destroy_null(&point_sc_array);
3233 template <
int dim,
int spacedim>
3239 for (
const auto &cell : this->active_cell_iterators())
3240 if (cell->is_locally_owned() && cell->refine_flag_set())
3241 Assert(cell->refine_flag_set() ==
3244 "This class does not support anisotropic refinement"));
3249 if (this->n_levels() ==
3253 cell = this->begin_active(
3261 !(cell->refine_flag_set()),
3263 "Fatal Error: maximum refinement level of p4est reached."));
3267 this->prepare_coarsening_and_refinement();
3270 this->signals.pre_distributed_refinement();
3275 for (
const auto &cell : this->active_cell_iterators())
3276 if (cell->is_ghost() || cell->is_artificial())
3278 cell->clear_refine_flag();
3279 cell->clear_coarsen_flag();
3285 RefineAndCoarsenList<dim, spacedim> refine_and_coarsen_list(
3286 *
this, p4est_tree_to_coarse_cell_permutation, this->my_subdomain);
3293 parallel_forest->user_pointer = &refine_and_coarsen_list;
3295 if (parallel_ghost !=
nullptr)
3299 parallel_ghost =
nullptr;
3304 &RefineAndCoarsenList<dim, spacedim>::refine_callback,
3309 &RefineAndCoarsenList<dim, spacedim>::coarsen_callback,
3316 parallel_forest->user_pointer =
this;
3322 (dim == 2 ? typename ::internal::p4est::types<dim>::balance_type(
3323 P4EST_CONNECT_FULL) :
3324 typename ::internal::p4est::types<dim>::balance_type(
3325 P8EST_CONNECT_FULL)),
3330 update_cell_relations();
3334 this->signals.post_p4est_refinement();
3338 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3339 previous_global_first_quadrant;
3341 if (this->cell_attached_data.n_attached_data_sets > 0)
3343 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3344 std::memcpy(previous_global_first_quadrant.data(),
3345 parallel_forest->global_first_quadrant,
3347 typename ::internal::p4est::types<dim>::gloidx) *
3348 (parallel_forest->mpisize + 1));
3351 if (!(
settings & no_automatic_repartitioning))
3355 if (this->signals.weight.empty())
3363 const std::vector<unsigned int> cell_weights = get_cell_weights();
3369 this->mpi_communicator) > 0,
3371 "The global sum of weights over all active cells "
3372 "is zero. Please verify how you generate weights."));
3374 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3379 parallel_forest->user_pointer = &partition_weights;
3385 &PartitionWeights<dim, spacedim>::cell_weight);
3389 parallel_forest, 0,
nullptr,
nullptr);
3391 parallel_forest->user_pointer =
this;
3396 if (this->cell_attached_data.n_attached_data_sets > 0)
3398 this->data_transfer.pack_data(
3399 this->local_cell_relations,
3400 this->cell_attached_data.pack_callbacks_fixed,
3401 this->cell_attached_data.pack_callbacks_variable);
3407 for (
const auto &cell : this->active_cell_iterators())
3409 cell->clear_refine_flag();
3410 cell->clear_coarsen_flag();
3415 copy_local_forest_to_triangulation();
3417 catch (
const typename Triangulation<dim>::DistortedCellList &)
3425 if (this->cell_attached_data.n_attached_data_sets > 0)
3427 this->execute_transfer(parallel_forest,
3428 previous_global_first_quadrant.data());
3431 this->data_transfer.unpack_cell_status(this->local_cell_relations);
3450 if (
settings & construct_multigrid_hierarchy)
3452 for (
unsigned int lvl = 0; lvl < this->n_global_levels(); ++lvl)
3454 std::vector<bool> active_verts =
3455 this->mark_locally_active_vertices_on_level(lvl);
3457 const unsigned int maybe_coarser_lvl =
3458 (lvl > 0) ? (lvl - 1) : lvl;
3460 cell = this->begin(maybe_coarser_lvl),
3461 endc = this->end(lvl);
3462 for (; cell != endc; ++cell)
3463 if (cell->level() ==
static_cast<int>(lvl) || cell->is_active())
3465 const bool is_level_artificial =
3466 (cell->level_subdomain_id() ==
3468 bool need_to_know =
false;
3469 for (
const unsigned int vertex :
3471 if (active_verts[cell->vertex_index(vertex)])
3473 need_to_know =
true;
3478 !need_to_know || !is_level_artificial,
3480 "Internal error: the owner of cell" +
3481 cell->id().to_string() +
3482 " is unknown even though it is needed for geometric multigrid."));
3488 this->update_periodic_face_map();
3489 this->update_number_cache();
3492 this->signals.post_distributed_refinement();
3497 template <
int dim,
int spacedim>
3502 for (
const auto &cell : this->active_cell_iterators())
3503 if (cell->is_locally_owned())
3505 !cell->refine_flag_set() && !cell->coarsen_flag_set(),
3507 "Error: There shouldn't be any cells flagged for coarsening/refinement when calling repartition()."));
3511 this->signals.pre_distributed_repartition();
3515 std::vector<typename ::internal::p4est::types<dim>::gloidx>
3516 previous_global_first_quadrant;
3518 if (this->cell_attached_data.n_attached_data_sets > 0)
3520 previous_global_first_quadrant.resize(parallel_forest->mpisize + 1);
3521 std::memcpy(previous_global_first_quadrant.data(),
3522 parallel_forest->global_first_quadrant,
3524 typename ::internal::p4est::types<dim>::gloidx) *
3525 (parallel_forest->mpisize + 1));
3528 if (this->signals.weight.empty())
3540 const std::vector<unsigned int> cell_weights = get_cell_weights();
3546 this->mpi_communicator) > 0,
3548 "The global sum of weights over all active cells "
3549 "is zero. Please verify how you generate weights."));
3551 PartitionWeights<dim, spacedim> partition_weights(cell_weights);
3556 parallel_forest->user_pointer = &partition_weights;
3562 &PartitionWeights<dim, spacedim>::cell_weight);
3565 parallel_forest->user_pointer =
this;
3569 if (this->cell_attached_data.n_attached_data_sets > 0)
3571 this->data_transfer.pack_data(
3572 this->local_cell_relations,
3573 this->cell_attached_data.pack_callbacks_fixed,
3574 this->cell_attached_data.pack_callbacks_variable);
3579 copy_local_forest_to_triangulation();
3581 catch (
const typename Triangulation<dim>::DistortedCellList &)
3589 if (this->cell_attached_data.n_attached_data_sets > 0)
3591 this->execute_transfer(parallel_forest,
3592 previous_global_first_quadrant.data());
3595 this->update_periodic_face_map();
3598 this->update_number_cache();
3601 this->signals.post_distributed_repartition();
3606 template <
int dim,
int spacedim>
3608 const std::vector<types::global_dof_index>
3612 return p4est_tree_to_coarse_cell_permutation;
3617 template <
int dim,
int spacedim>
3619 const std::vector<types::global_dof_index>
3623 return coarse_cell_to_p4est_tree_permutation;
3628 template <
int dim,
int spacedim>
3631 mark_locally_active_vertices_on_level(const
int level)
const
3635 std::vector<bool> marked_vertices(this->n_vertices(),
false);
3636 for (
const auto &cell : this->cell_iterators_on_level(
level))
3637 if (cell->level_subdomain_id() == this->locally_owned_subdomain())
3639 marked_vertices[cell->vertex_index(v)] =
true;
3656 for (
unsigned int repetition = 0; repetition < dim; ++repetition)
3657 for (
const auto &it : this->get_periodic_face_map())
3660 const unsigned int face_no_1 = it.first.second;
3662 const unsigned int face_no_2 = it.second.first.second;
3663 const std::bitset<3> &face_orientation = it.second.second;
3665 if (cell_1->level() ==
level && cell_2->level() ==
level)
3667 for (
unsigned int v = 0;
3673 const unsigned int vface0 =
3676 face_orientation[0],
3677 face_orientation[1],
3678 face_orientation[2]);
3679 if (marked_vertices[cell_1->face(face_no_1)->vertex_index(
3681 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3683 marked_vertices[cell_1->face(face_no_1)->vertex_index(
3685 marked_vertices[cell_2->face(face_no_2)->vertex_index(
3691 return marked_vertices;
3696 template <
int dim,
int spacedim>
3699 coarse_cell_id_to_coarse_cell_index(
3700 const
types::coarse_cell_id coarse_cell_id)
const
3702 return p4est_tree_to_coarse_cell_permutation[coarse_cell_id];
3707 template <
int dim,
int spacedim>
3711 const unsigned int coarse_cell_index)
const
3713 return coarse_cell_to_p4est_tree_permutation[coarse_cell_index];