From 35424a773ea4ec8e34001ef74379a5995a5c4780 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Fri, 20 Apr 2018 20:34:41 -0600 Subject: [PATCH] Manage relations of p4est quadrants and dealii cells. Use them to rewrite some of the recursively crawling functions. --- include/deal.II/distributed/tria.h | 36 +- source/distributed/tria.cc | 821 +++++++++++++---------------- 2 files changed, 406 insertions(+), 451 deletions(-) diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 7bef03bfcb..a604bf3dad 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -911,6 +911,40 @@ namespace parallel CellAttachedData cell_attached_data; + /** + * This auxiliary data structure stores the relation between a p4est + * quadrant, a deal.II cell and its current CellStatus. For an extensive + * description of the latter, see the documentation for the member + * function register_data_attach. + */ + typedef typename std::tuple< + typename dealii::internal::p4est::types::quadrant *, + CellStatus, + cell_iterator> + quadrant_cell_relation_t; + + /** + * Vector of tuples, which each contain a p4est quadrant, a deal.II cell + * and their relation after refinement. To update its contents, use the + * compute_quadrant_cell_relations member function. + * + * The size of this vector is assumed to be equal to the number of locally + * owned quadrants in the parallel_forest object. + */ + std::vector local_quadrant_cell_relations; + + /** + * Go through all p4est trees and store the relations between locally + * owned quadrants and cells in the private member + * local_quadrant_cell_relations. + * + * The stored vector will be ordered by the occurrence of quadrants in + * the corresponding local sc_array of the parallel_forest. p4est requires + * this specific ordering for its transfer functions. + */ + void + update_quadrant_cell_relations(); + /** * Two arrays that store which p4est tree corresponds to which coarse * grid cell and vice versa. We need these arrays because p4est goes @@ -996,7 +1030,7 @@ namespace parallel * them. */ std::vector - get_cell_weights(); + get_cell_weights() const; /** * Override the implementation in parallel::Triangulation because diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 991f5a5b9c..3e3c22f5ca 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -574,391 +574,6 @@ namespace - template - void - attach_mesh_data_recursively( - const typename internal::p4est::types::tree & tree, - const typename Triangulation::cell_iterator &dealii_cell, - const typename internal::p4est::types::quadrant & p4est_cell, - const std::vector::cell_iterator, - typename parallel::distributed:: - Triangulation::CellStatus, - void *)>>> & pack_callbacks) - { - int idx = - sc_array_bsearch(const_cast(&tree.quadrants), - &p4est_cell, - internal::p4est::functions::quadrant_compare); - - if (idx == -1 && - (internal::p4est::functions::quadrant_overlaps_tree( - const_cast::tree *>(&tree), - &p4est_cell) == false)) - return; // this quadrant and none of its children belongs to us. - - bool p4est_has_children = (idx == -1); - - if (p4est_has_children && dealii_cell->has_children()) - { - // recurse further - typename internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - switch (dim) - { - case 2: - P4EST_QUADRANT_INIT(&p4est_child[c]); - break; - case 3: - P8EST_QUADRANT_INIT(&p4est_child[c]); - break; - default: - Assert(false, ExcNotImplemented()); - } - - internal::p4est::functions::quadrant_childrenv(&p4est_cell, - p4est_child); - - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - { - attach_mesh_data_recursively( - tree, dealii_cell->child(c), p4est_child[c], pack_callbacks); - } - } - else if (!p4est_has_children && !dealii_cell->has_children()) - { - // this active cell didn't change - typename internal::p4est::types::quadrant *q; - q = static_cast::quadrant *>( - sc_array_index(const_cast(&tree.quadrants), idx)); - *static_cast::CellStatus *>( - q->p.user_data) = - parallel::distributed::Triangulation::CELL_PERSIST; - - // double check the condition that we will only ever attach data - // to active cells when we get here - Assert(dealii_cell->active(), ExcInternalError()); - - // call all callbacks - for (const auto &it : pack_callbacks) - { - void *ptr = - static_cast(q->p.user_data) + it.first; // add offset - (it.second)( - dealii_cell, - parallel::distributed::Triangulation::CELL_PERSIST, - ptr); - } - } - else if (p4est_has_children) - { - // this cell got refined in p4est, but the dealii_cell has not yet been - // refined - - // attach to the first child, because we can only attach to active - // quadrants - typename internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - switch (dim) - { - case 2: - P4EST_QUADRANT_INIT(&p4est_child[c]); - break; - case 3: - P8EST_QUADRANT_INIT(&p4est_child[c]); - break; - default: - Assert(false, ExcNotImplemented()); - } - - internal::p4est::functions::quadrant_childrenv(&p4est_cell, - p4est_child); - int child0_idx = - sc_array_bsearch(const_cast(&tree.quadrants), - &p4est_child[0], - internal::p4est::functions::quadrant_compare); - Assert( - child0_idx != -1, - ExcMessage("the first child should exist as an active quadrant!")); - - typename internal::p4est::types::quadrant *q; - q = static_cast::quadrant *>( - sc_array_index(const_cast(&tree.quadrants), - child0_idx)); - *static_cast::CellStatus *>( - q->p.user_data) = - parallel::distributed::Triangulation::CELL_REFINE; - - // double check the condition that we will only ever attach data - // to active cells when we get here - Assert(dealii_cell->active(), ExcInternalError()); - - // call all callbacks - for (const auto &it : pack_callbacks) - { - // compute the address where this particular callback is allowed - // to write its data, using the correct offset - void *ptr = static_cast(q->p.user_data) + it.first; - - (it.second)( - dealii_cell, - parallel::distributed::Triangulation::CELL_REFINE, - ptr); - } - - // mark other children as invalid, so that unpack only happens once - for (unsigned int i = 1; i < GeometryInfo::max_children_per_cell; - ++i) - { - int child_idx = sc_array_bsearch( - const_cast(&tree.quadrants), - &p4est_child[i], - internal::p4est::functions::quadrant_compare); - q = static_cast::quadrant *>( - sc_array_index(const_cast(&tree.quadrants), - child_idx)); - *static_cast::CellStatus *>( - q->p.user_data) = - parallel::distributed::Triangulation::CELL_INVALID; - } - } - else - { - // its children got coarsened into this cell in p4est, but the - // dealii_cell still has its children - typename internal::p4est::types::quadrant *q; - q = static_cast::quadrant *>( - sc_array_index(const_cast(&tree.quadrants), idx)); - *static_cast::CellStatus *>( - q->p.user_data) = - parallel::distributed::Triangulation::CELL_COARSEN; - - // double check the condition that we will only ever attach data - // to cells with children when we get here. however, we can - // only tolerate one level of coarsening at a time, so check - // that the children are all active - Assert(dealii_cell->active() == false, ExcInternalError()); - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - Assert(dealii_cell->child(c)->active(), ExcInternalError()); - - // then call all callbacks - for (const auto &it : pack_callbacks) - { - void *ptr = - static_cast(q->p.user_data) + it.first; // add offset - (it.second)( - dealii_cell, - parallel::distributed::Triangulation::CELL_COARSEN, - ptr); - } - } - } - - template - void - get_cell_weights_recursively( - const typename internal::p4est::types::tree & tree, - const typename Triangulation::cell_iterator &dealii_cell, - const typename internal::p4est::types::quadrant & p4est_cell, - const typename Triangulation::Signals & signals, - std::vector & weight) - { - const int idx = - sc_array_bsearch(const_cast(&tree.quadrants), - &p4est_cell, - internal::p4est::functions::quadrant_compare); - - if (idx == -1 && - (internal::p4est::functions::quadrant_overlaps_tree( - const_cast::tree *>(&tree), - &p4est_cell) == false)) - return; // This quadrant and none of its children belongs to us. - - const bool p4est_has_children = (idx == -1); - - if (p4est_has_children && dealii_cell->has_children()) - { - // recurse further - typename internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - switch (dim) - { - case 2: - P4EST_QUADRANT_INIT(&p4est_child[c]); - break; - case 3: - P8EST_QUADRANT_INIT(&p4est_child[c]); - break; - default: - Assert(false, ExcNotImplemented()); - } - - internal::p4est::functions::quadrant_childrenv(&p4est_cell, - p4est_child); - - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - { - get_cell_weights_recursively( - tree, dealii_cell->child(c), p4est_child[c], signals, weight); - } - } - else if (!p4est_has_children && !dealii_cell->has_children()) - { - // This active cell didn't change - weight.push_back(1000); - weight.back() += signals.cell_weight( - dealii_cell, - parallel::distributed::Triangulation::CELL_PERSIST); - } - else if (p4est_has_children) - { - // This cell will be refined - unsigned int parent_weight(1000); - parent_weight += signals.cell_weight( - dealii_cell, - parallel::distributed::Triangulation::CELL_REFINE); - - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - { - // We assign the weight of the parent cell equally to all children - weight.push_back(parent_weight); - } - } - else - { - // This cell's children will be coarsened into this cell - weight.push_back(1000); - weight.back() += signals.cell_weight( - dealii_cell, - parallel::distributed::Triangulation::CELL_COARSEN); - } - } - - - template - void - post_mesh_data_recursively( - const typename internal::p4est::types::tree & tree, - const typename Triangulation::cell_iterator &dealii_cell, - const typename Triangulation::cell_iterator &parent_cell, - const typename internal::p4est::types::quadrant & p4est_cell, - const unsigned int offset, - const typename std::function::cell_iterator, - typename parallel::distributed::Triangulation::CellStatus, - void *)> & unpack_callback) - { - int idx = - sc_array_bsearch(const_cast(&tree.quadrants), - &p4est_cell, - internal::p4est::functions::quadrant_compare); - if (idx == -1 && - (internal::p4est::functions::quadrant_overlaps_tree( - const_cast::tree *>(&tree), - &p4est_cell) == false)) - // this quadrant and none of its children belong to us. - return; - - - const bool p4est_has_children = (idx == -1); - if (p4est_has_children) - { - Assert(dealii_cell->has_children(), ExcInternalError()); - - // recurse further - typename internal::p4est::types::quadrant - p4est_child[GeometryInfo::max_children_per_cell]; - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - switch (dim) - { - case 2: - P4EST_QUADRANT_INIT(&p4est_child[c]); - break; - case 3: - P8EST_QUADRANT_INIT(&p4est_child[c]); - break; - default: - Assert(false, ExcNotImplemented()); - } - - internal::p4est::functions::quadrant_childrenv(&p4est_cell, - p4est_child); - - for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; - ++c) - { - post_mesh_data_recursively(tree, - dealii_cell->child(c), - dealii_cell, - p4est_child[c], - offset, - unpack_callback); - } - } - else - { - Assert(!dealii_cell->has_children(), ExcInternalError()); - - typename internal::p4est::types::quadrant *q; - q = static_cast::quadrant *>( - sc_array_index(const_cast(&tree.quadrants), idx)); - - void *ptr = static_cast(q->p.user_data) + offset; - typename parallel::distributed::Triangulation::CellStatus - status = *static_cast::CellStatus *>( - q->p.user_data); - switch (status) - { - case parallel::distributed::Triangulation::CELL_PERSIST: - { - unpack_callback(dealii_cell, status, ptr); - break; - } - case parallel::distributed::Triangulation::CELL_REFINE: - { - unpack_callback(parent_cell, status, ptr); - break; - } - case parallel::distributed::Triangulation::CELL_COARSEN: - { - unpack_callback(dealii_cell, status, ptr); - break; - } - case parallel::distributed::Triangulation::CELL_INVALID: - { - break; - } - default: - Assert(false, ExcInternalError()); - } - } - } - - - /** * A data structure that we use to store which cells (indicated by * internal::p4est::types::quadrant objects) shall be refined and which @@ -1326,6 +941,184 @@ namespace // get the weight, increment the pointer, and return the weight return *this_object->current_pointer++; } + + + + template + using quadrant_cell_relation_t = typename std::tuple< + typename dealii::internal::p4est::types::quadrant *, + typename dealii::Triangulation::CellStatus, + typename dealii::Triangulation::cell_iterator>; + + + + /** + * Adds a tuple of a p4est quadrant, @p status and @p dealii_cell + * to the vector containing all relations @p quad_cell_rel. + * The tuple will be inserted in the position corresponding to the one + * of the p4est quadrant in the underlying p4est sc_array. The position + * will be determined from @p idx, which is the position of the quadrant + * in its corresponding @p tree. The p4est quadrant will be deduced from + * the @p tree by @p idx. + */ + template + inline void + add_single_quadrant_cell_relation( + std::vector> & quad_cell_rel, + const typename dealii::internal::p4est::types::tree & tree, + const unsigned int idx, + const typename Triangulation::cell_iterator &dealii_cell, + const typename Triangulation::CellStatus status) + { + const unsigned int local_quadrant_index = tree.quadrants_offset + idx; + + const auto q = + static_cast::quadrant *>( + sc_array_index(const_cast(&tree.quadrants), idx)); + + // check if we will be writing into valid memory + Assert(local_quadrant_index < quad_cell_rel.size(), ExcInternalError()); + + // store relation + quad_cell_rel[local_quadrant_index] = + std::make_tuple(q, status, dealii_cell); + } + + + + /** + * This is the recursive part of the member function + * update_quadrant_cell_relations(). + * + * Find the relation between the @p p4est_cell and the @p dealii_cell in the + * corresponding @p tree. Depending on the CellStatus relation between the two, + * a new entry will either be inserted in @p quad_cell_rel or the recursion + * will be continued. + */ + template + void + update_quadrant_cell_relations_recursively( + std::vector> & quad_cell_rel, + const typename dealii::internal::p4est::types::tree & tree, + const typename Triangulation::cell_iterator & dealii_cell, + const typename dealii::internal::p4est::types::quadrant &p4est_cell) + { + // find index of p4est_cell in the quadrants array of the corresponding tree + const int idx = sc_array_bsearch( + const_cast(&tree.quadrants), + &p4est_cell, + dealii::internal::p4est::functions::quadrant_compare); + if (idx == -1 && + (dealii::internal::p4est::functions::quadrant_overlaps_tree( + const_cast::tree *>( + &tree), + &p4est_cell) == false)) + // this quadrant and none of its children belong to us. + return; + + // recurse further if both p4est and dealii still have children + const bool p4est_has_children = (idx == -1); + if (p4est_has_children && dealii_cell->has_children()) + { + // recurse further + typename dealii::internal::p4est::types::quadrant + p4est_child[GeometryInfo::max_children_per_cell]; + + for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; + ++c) + switch (dim) + { + case 2: + P4EST_QUADRANT_INIT(&p4est_child[c]); + break; + case 3: + P8EST_QUADRANT_INIT(&p4est_child[c]); + break; + default: + Assert(false, ExcNotImplemented()); + } + + dealii::internal::p4est::functions::quadrant_childrenv( + &p4est_cell, p4est_child); + + for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; + ++c) + { + update_quadrant_cell_relations_recursively( + quad_cell_rel, tree, dealii_cell->child(c), p4est_child[c]); + } + } + else if (!p4est_has_children && !dealii_cell->has_children()) + { + // this active cell didn't change + // save tuple into corresponding position + add_single_quadrant_cell_relation( + quad_cell_rel, + tree, + idx, + dealii_cell, + Triangulation::CELL_PERSIST); + } + else if (p4est_has_children) // based on the conditions above, we know that + // dealii_cell has no children + { + // this cell got refined in p4est, but the dealii_cell has not yet been + // refined + + // this quadrant is not active + // generate its children, and store information in those + typename dealii::internal::p4est::types::quadrant + p4est_child[GeometryInfo::max_children_per_cell]; + for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; + ++c) + switch (dim) + { + case 2: + P4EST_QUADRANT_INIT(&p4est_child[c]); + break; + case 3: + P8EST_QUADRANT_INIT(&p4est_child[c]); + break; + default: + Assert(false, ExcNotImplemented()); + } + + dealii::internal::p4est::functions::quadrant_childrenv( + &p4est_cell, p4est_child); + + // mark first child with CELL_REFINE and the remaining children with + // CELL_INVALID, but associate them all with the parent cell unpack + // algorithm will be called only on CELL_REFINE flagged quadrant + int child_idx; + typename Triangulation::CellStatus cell_status; + for (unsigned int i = 0; i < GeometryInfo::max_children_per_cell; + ++i) + { + child_idx = sc_array_bsearch( + const_cast(&tree.quadrants), + &p4est_child[i], + dealii::internal::p4est::functions::quadrant_compare); + + cell_status = (i == 0) ? Triangulation::CELL_REFINE : + Triangulation::CELL_INVALID; + + add_single_quadrant_cell_relation( + quad_cell_rel, tree, child_idx, dealii_cell, cell_status); + } + } + else // based on the conditions above, we know that p4est_cell has no + // children, and the dealii_cell does + { + // its children got coarsened into this cell in p4est, + // but the dealii_cell still has its children + add_single_quadrant_cell_relation( + quad_cell_rel, + tree, + idx, + dealii_cell, + Triangulation::CELL_COARSEN); + } + } } // namespace @@ -2963,6 +2756,12 @@ namespace parallel } this->smooth_grid = save_smooth; + + // finally, after syncing the parallel_forest with the triangulation, + // also update the quadrant_cell_relations, which will be used for + // repartitioning, further refinement/coarsening, and unpacking + // of stored or transfered data. + update_quadrant_cell_relations(); } @@ -3068,6 +2867,10 @@ namespace parallel P8EST_CONNECT_FULL)), /*init_callback=*/nullptr); + // since refinement and/or coarsening on the parallel forest + // has happened, we need to update the quadrant cell relations + update_quadrant_cell_relations(); + // before repartitioning the mesh let others attach mesh related info // (such as SolutionTransfer data) to the p4est attach_mesh_data(); @@ -3484,35 +3287,53 @@ namespace parallel ExcMessage("The notify_ready_to_unpack() has already been called " "once for each registered callback.")); + // check if local_quadrant_cell_relations have been previously gathered + // correctly + Assert(local_quadrant_cell_relations.size() == + static_cast(parallel_forest->local_num_quadrants), + ExcInternalError()); + const unsigned int offset = std::get<0>(cell_attached_data.pack_callbacks[handle]); // check if unpack function has been previously called Assert(offset != dealii::numbers::invalid_unsigned_int, ExcInternalError()); - // Recurse over p4est and hand the caller the data back - for (typename Triangulation::cell_iterator cell = - this->begin(0); - cell != this->end(0); - ++cell) + for (const auto &it_rel : local_quadrant_cell_relations) { - // skip coarse cells that are not ours - if (tree_exists_locally( - parallel_forest, - coarse_cell_to_p4est_tree_permutation[cell->index()]) == false) - continue; + const auto &q = std::get<0>(it_rel); + const auto &cell_it = std::get<2>(it_rel); - typename dealii::internal::p4est::types::quadrant - p4est_coarse_cell; - typename dealii::internal::p4est::types::tree *tree = - init_tree(cell->index()); + const auto cell_status = + *static_cast::CellStatus *>( + q->p.user_data); - dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + void *ptr = static_cast(q->p.user_data) + offset; + + switch (cell_status) + { + case parallel::distributed::Triangulation::CELL_PERSIST: + case parallel::distributed::Triangulation::CELL_COARSEN: + unpack_callback(cell_it, cell_status, ptr); + break; + + case parallel::distributed::Triangulation::CELL_REFINE: + unpack_callback(cell_it->parent(), cell_status, ptr); + break; - // parent_cell is not correct here, but is only used in a refined - // cell - post_mesh_data_recursively( - *tree, cell, cell, p4est_coarse_cell, offset, unpack_callback); + case parallel::distributed::Triangulation::CELL_INVALID: + // do nothing on these cells + break; + + default: + Assert(false, ExcInternalError()); + break; + } } --cell_attached_data.n_attached_datas; @@ -3941,6 +3762,43 @@ namespace parallel + template + void + Triangulation::update_quadrant_cell_relations() + { + // reorganize memory for local_quadrant_cell_relations + local_quadrant_cell_relations.resize( + parallel_forest->local_num_quadrants); + local_quadrant_cell_relations.shrink_to_fit(); + + // recurse over p4est + for (typename Triangulation::cell_iterator cell = + this->begin(0); + cell != this->end(0); + ++cell) + { + // skip coarse cells that are not ours + if (tree_exists_locally( + parallel_forest, + coarse_cell_to_p4est_tree_permutation[cell->index()]) == false) + continue; + + // initialize auxiliary top level p4est quadrant + typename dealii::internal::p4est::types::quadrant + p4est_coarse_cell; + dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + + // determine tree to start recursion on + typename dealii::internal::p4est::types::tree *tree = + init_tree(cell->index()); + + update_quadrant_cell_relations_recursively( + local_quadrant_cell_relations, *tree, cell, p4est_coarse_cell); + } + } + + + template void Triangulation::attach_mesh_data() @@ -3952,6 +3810,12 @@ namespace parallel return; } + // check if local_quadrant_cell_relations have been previously gathered + // correctly + Assert(local_quadrant_cell_relations.size() == + static_cast(parallel_forest->local_num_quadrants), + ExcInternalError()); + // reallocate user_data in p4est void *userptr = parallel_forest->user_pointer; dealii::internal::p4est::functions::reset_data( @@ -3961,30 +3825,63 @@ namespace parallel nullptr); parallel_forest->user_pointer = userptr; - - // Recurse over p4est and Triangulation - // to find refined/coarsened/kept - // cells. Then query and attach the data. - for (typename Triangulation::cell_iterator cell = - this->begin(0); - cell != this->end(0); - ++cell) + for (const auto &it_rel : local_quadrant_cell_relations) { - // skip coarse cells, that are not ours - if (tree_exists_locally( - parallel_forest, - coarse_cell_to_p4est_tree_permutation[cell->index()]) == false) - continue; + const auto &q = std::get<0>(it_rel); + const auto &cell_status = std::get<1>(it_rel); + const auto &cell_it = std::get<2>(it_rel); - typename dealii::internal::p4est::types::quadrant - p4est_coarse_cell; - typename dealii::internal::p4est::types::tree *tree = - init_tree(cell->index()); + *static_cast::CellStatus *>( + q->p.user_data) = cell_status; - dealii::internal::p4est::init_coarse_quadrant(p4est_coarse_cell); + switch (cell_status) + { + case parallel::distributed::Triangulation::CELL_PERSIST: + case parallel::distributed::Triangulation::CELL_COARSEN: + case parallel::distributed::Triangulation::CELL_REFINE: +# ifdef DEBUG + if (cell_status == parallel::distributed:: + Triangulation::CELL_COARSEN) + { + // double check the condition that we will only ever attach + // data to cells with children when we get here. however, we + // can only tolerate one level of coarsening at a time, so + // check that the children are all active + Assert(cell_it->active() == false, ExcInternalError()); + for (unsigned int c = 0; + c < GeometryInfo::max_children_per_cell; + ++c) + Assert(cell_it->child(c)->active(), ExcInternalError()); + } + else + { + // double check the condition that we will only ever attach + // data to active cells when we get here + Assert(cell_it->active(), ExcInternalError()); + } +# endif + // call all callbacks + for (const auto &it_pack : cell_attached_data.pack_callbacks) + { + void *ptr = static_cast(q->p.user_data) + + it_pack.first; // add offset + (it_pack.second)(cell_it, cell_status, ptr); + } + break; - attach_mesh_data_recursively( - *tree, cell, p4est_coarse_cell, cell_attached_data.pack_callbacks); + case parallel::distributed::Triangulation::CELL_INVALID: + // do nothing on these cells + break; + + default: + Assert(false, ExcInternalError()); + break; + } } } @@ -3992,8 +3889,14 @@ namespace parallel template std::vector - Triangulation::get_cell_weights() + Triangulation::get_cell_weights() const { + // check if local_quadrant_cell_relations have been previously gathered + // correctly + Assert(local_quadrant_cell_relations.size() == + static_cast(parallel_forest->local_num_quadrants), + ExcInternalError()); + // Allocate the space for the weights. In fact we do not know yet, how // many cells we own after the refinement (only p4est knows that // at this point). We simply reserve n_active_cells space and if many @@ -4002,40 +3905,58 @@ namespace parallel std::vector weights; weights.reserve(this->n_active_cells()); - // Recurse over p4est and Triangulation + // Iterate over p4est and Triangulation relations // to find refined/coarsened/kept // cells. Then append cell_weight. // Note that we need to follow the p4est ordering // instead of the deal.II ordering to get the cell_weights // in the same order p4est will encounter them during repartitioning. - for (unsigned int c = 0; c < this->n_cells(0); ++c) + for (const auto &it_rel : local_quadrant_cell_relations) { - // skip coarse cells, that are not ours - if (tree_exists_locally(parallel_forest, c) == false) - continue; + const auto &cell_status = std::get<1>(it_rel); + const auto &cell_it = std::get<2>(it_rel); - const unsigned int coarse_cell_index = - p4est_tree_to_coarse_cell_permutation[c]; + switch (cell_status) + { + case parallel::distributed::Triangulation::CELL_PERSIST: + weights.push_back(1000); + weights.back() += this->signals.cell_weight( + cell_it, + parallel::distributed::Triangulation::CELL_PERSIST); + break; - const typename Triangulation::cell_iterator - dealii_coarse_cell(this, 0, coarse_cell_index); + case parallel::distributed::Triangulation::CELL_REFINE: + case parallel::distributed::Triangulation::CELL_INVALID: + { + // calculate weight of parent cell + unsigned int parent_weight = 1000; + parent_weight += this->signals.cell_weight( + cell_it, + parallel::distributed::Triangulation:: + CELL_REFINE); + // assign the weight of the parent cell equally to all + // children + weights.push_back(parent_weight); + break; + } - typename dealii::internal::p4est::types::quadrant - p4est_coarse_cell; - dealii::internal::p4est::functions::quadrant_set_morton( - &p4est_coarse_cell, - /*level=*/0, - /*index=*/0); - p4est_coarse_cell.p.which_tree = c; - - const typename dealii::internal::p4est::types::tree *tree = - init_tree(coarse_cell_index); - - get_cell_weights_recursively(*tree, - dealii_coarse_cell, - p4est_coarse_cell, - this->signals, - weights); + case parallel::distributed::Triangulation::CELL_COARSEN: + weights.push_back(1000); + weights.back() += this->signals.cell_weight( + cell_it, + parallel::distributed::Triangulation::CELL_COARSEN); + break; + + default: + Assert(false, ExcInternalError()); + break; + } } return weights; -- 2.39.5