From 8a514aabf43de16b6f578ce8ea26d5eba65827da Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 27 Aug 2021 08:12:40 +0200 Subject: [PATCH] Restructure TriangulationDescription::Utilities::create_description_from_triangulation() 2 --- source/grid/tria_description.cc | 1223 ++++++++++++++++--------------- 1 file changed, 626 insertions(+), 597 deletions(-) diff --git a/source/grid/tria_description.cc b/source/grid/tria_description.cc index fa6df29355..37210bbb8c 100644 --- a/source/grid/tria_description.cc +++ b/source/grid/tria_description.cc @@ -82,6 +82,244 @@ namespace TriangulationDescription { namespace { + /** + * A temporal class similar to TriangulationDescription::Description. + */ + template + struct DescriptionTemp + { + /** + * Serialization function for packing and unpacking the content of this + * class. + */ + template + void + serialize(Archive &ar, const unsigned int /*version*/) + { + ar &coarse_cells; + ar &coarse_cell_vertices; + ar &coarse_cell_index_to_coarse_cell_id; + ar &cell_infos; + } + + /** + * Send the DescriptionTemp objects to the right process. + */ + void + collect( + const std::vector & relevant_processes, + const std::vector> &description_temp, + const MPI_Comm & comm) + { + dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess + process( + [&]() { return relevant_processes; }, + [&](const unsigned int other_rank, + std::vector &send_buffer) { + const auto ptr = std::find(relevant_processes.begin(), + relevant_processes.end(), + other_rank); + + Assert(ptr != relevant_processes.end(), ExcInternalError()); + + const auto other_rank_index = + std::distance(relevant_processes.begin(), ptr); + + send_buffer = + dealii::Utilities::pack(description_temp[other_rank_index], + false); + }, + [&](const unsigned int &, + const std::vector &recv_buffer, + std::vector &) { + this->merge( + dealii::Utilities::unpack>( + recv_buffer, false)); + }); + + dealii::Utilities::MPI::ConsensusAlgorithms::Selector( + process, comm) + .run(); + } + + /** + * Merge the given DescriptionTemp object into the current one. Note + * that no duplicate cells, vertices, ... are removed at this stage. + * This done by the reduce() function. + */ + void + merge(const DescriptionTemp &other) + { + this->cell_infos.resize( + std::max(other.cell_infos.size(), this->cell_infos.size())); + + this->coarse_cells.insert(this->coarse_cells.end(), + other.coarse_cells.begin(), + other.coarse_cells.end()); + this->coarse_cell_vertices.insert(this->coarse_cell_vertices.end(), + other.coarse_cell_vertices.begin(), + other.coarse_cell_vertices.end()); + this->coarse_cell_index_to_coarse_cell_id.insert( + this->coarse_cell_index_to_coarse_cell_id.end(), + other.coarse_cell_index_to_coarse_cell_id.begin(), + other.coarse_cell_index_to_coarse_cell_id.end()); + + for (unsigned int i = 0; i < this->cell_infos.size(); ++i) + this->cell_infos[i].insert(this->cell_infos[i].end(), + other.cell_infos[i].begin(), + other.cell_infos[i].end()); + } + + /** + * Remove all duplicate information obtained during merge(). + */ + void + reduce() + { + // make coarse cells unique + { + std::vector, + unsigned int>> + temp; + + for (unsigned int i = 0; i < this->coarse_cells.size(); ++i) + temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i], + this->coarse_cells[i], + i); + + std::sort(temp.begin(), + temp.end(), + [](const auto &a, const auto &b) { + return std::get<0>(a) < std::get<0>(b); + }); + temp.erase(std::unique(temp.begin(), + temp.end(), + [](const auto &a, const auto &b) { + return std::get<0>(a) == std::get<0>(b); + }), + temp.end()); + std::sort(temp.begin(), + temp.end(), + [](const auto &a, const auto &b) { + return std::get<2>(a) < std::get<2>(b); + }); + + this->coarse_cell_index_to_coarse_cell_id.resize(temp.size()); + this->coarse_cells.resize(temp.size()); + + for (unsigned int i = 0; i < temp.size(); ++i) + { + this->coarse_cell_index_to_coarse_cell_id[i] = + std::get<0>(temp[i]); + this->coarse_cells[i] = std::get<1>(temp[i]); + } + } + + // make coarse cell vertices unique + { + std::sort(this->coarse_cell_vertices.begin(), + this->coarse_cell_vertices.end(), + [](const auto &a, const auto &b) { + return a.first < b.first; + }); + this->coarse_cell_vertices.erase( + std::unique(this->coarse_cell_vertices.begin(), + this->coarse_cell_vertices.end(), + [](const auto &a, const auto &b) { + return a.first == b.first; + }), + this->coarse_cell_vertices.end()); + } + + // make cells unique + for (unsigned int i = 0; i < this->cell_infos.size(); ++i) + { + if (this->cell_infos[i].size() == 0) + continue; + + std::sort(this->cell_infos[i].begin(), + this->cell_infos[i].end(), + [](const auto &a, const auto &b) { + return a.id < b.id; + }); + + std::vector> temp; + temp.push_back(this->cell_infos[i][0]); + + for (unsigned int j = 1; j < this->cell_infos[i].size(); ++j) + if (temp.back().id == cell_infos[i][j].id) + { + temp.back().subdomain_id = + std::min(temp.back().subdomain_id, + this->cell_infos[i][j].subdomain_id); + temp.back().level_subdomain_id = + std::min(temp.back().level_subdomain_id, + this->cell_infos[i][j].level_subdomain_id); + } + else + { + temp.push_back(this->cell_infos[i][j]); + } + + this->cell_infos[i] = temp; + } + } + + /** + * Convert the content of this class to a + * TriangulationDescription::Description, which can be used during + * Triangulation::create_triangulation(). + */ + Description + convert(const MPI_Comm comm, + const typename Triangulation::MeshSmoothing + mesh_smoothing, + const TriangulationDescription::Settings settings) + { + Description description; + + // copy communicator + description.comm = comm; + + description.settings = settings; + + // use mesh smoothing from base triangulation + description.smoothing = mesh_smoothing; + + std::map map; + + for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i) + { + description.coarse_cell_vertices.push_back( + this->coarse_cell_vertices[i].second); + map[this->coarse_cell_vertices[i].first] = i; + } + + description.coarse_cells = this->coarse_cells; + + for (auto &cell : description.coarse_cells) + for (unsigned int v = 0; v < cell.vertices.size(); ++v) + cell.vertices[v] = map[cell.vertices[v]]; + + description.coarse_cell_index_to_coarse_cell_id = + this->coarse_cell_index_to_coarse_cell_id; + description.cell_infos = this->cell_infos; + + return description; + } + + std::vector> coarse_cells; + + std::vector>> + coarse_cell_vertices; + + std::vector coarse_cell_index_to_coarse_cell_id; + + std::vector>> cell_infos; + }; + /** * Set the user_flag of a cell and of all its parent cells. */ @@ -94,6 +332,345 @@ namespace TriangulationDescription if (cell->level() != 0) set_user_flag_and_of_its_parents(cell->parent()); } + + /** + * Helper calls for the + * TriangulationDescription::Utilities::create_description_from_triangulation() + * functions. + */ + template + class CreateDescriptionFromTriangulationHelper + { + public: + CreateDescriptionFromTriangulationHelper( + const dealii::Triangulation &tria, + const std::function::cell_iterator + &)> & subdomain_id_function, + const std::function::cell_iterator + &)> & level_subdomain_id_function, + const MPI_Comm & comm, + const TriangulationDescription::Settings settings) + : tria(tria) + , subdomain_id_function(subdomain_id_function) + , level_subdomain_id_function(level_subdomain_id_function) + , comm(comm) + , settings(settings) + , construct_multigrid( + settings & + TriangulationDescription::Settings::construct_multigrid_hierarchy) + { + Assert( + !(settings & TriangulationDescription::Settings:: + construct_multigrid_hierarchy) || + (tria.get_mesh_smoothing() & + Triangulation::limit_level_difference_at_vertices), + ExcMessage( + "Source triangulation has to be set up with " + "limit_level_difference_at_vertices if the construction of the " + "multigrid hierarchy is requested!")); + + GridTools::collect_coinciding_vertices( + tria, coinciding_vertex_groups, vertex_to_coinciding_vertex_group); + } + + template + T + create_description_for_rank(const unsigned int my_rank) const + { + T construction_data; + + set_additional_data(construction_data); + + // helper function, which collects all vertices belonging to a cell + // (also taking into account periodicity) + const auto + add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells = + [this](const auto & cell, + std::vector &vertices_owned_by_locally_owned_cells) { + // add local vertices + for (const auto v : cell->vertex_indices()) + { + vertices_owned_by_locally_owned_cells[cell->vertex_index( + v)] = true; + const auto coinciding_vertex_group = + vertex_to_coinciding_vertex_group.find( + cell->vertex_index(v)); + if (coinciding_vertex_group != + vertex_to_coinciding_vertex_group.end()) + for (const auto &co_vertex : coinciding_vertex_groups.at( + coinciding_vertex_group->second)) + vertices_owned_by_locally_owned_cells[co_vertex] = true; + } + }; + + // 1) collect locally relevant cells (set user_flag) + std::vector old_user_flags; + tria.save_user_flags(old_user_flags); + + // 1a) clear user_flags + const_cast &>(tria) + .clear_user_flags(); + + // 1b) loop over levels (from fine to coarse) and mark on each level + // the locally relevant cells + for (int level = tria.get_triangulation().n_global_levels() - 1; + level >= 0; + --level) + { + // collect vertices connected to a (on any level) locally owned + // cell + std::vector vertices_owned_by_locally_owned_cells_on_level( + tria.n_vertices()); + for (const auto &cell : tria.cell_iterators_on_level(level)) + if (construct_multigrid && + (level_subdomain_id_function(cell) == my_rank)) + add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( + cell, vertices_owned_by_locally_owned_cells_on_level); + + for (const auto &cell : tria.active_cell_iterators()) + if (subdomain_id_function(cell) == my_rank) + add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( + cell, vertices_owned_by_locally_owned_cells_on_level); + + // helper function to determine if cell is locally relevant + // (i.e. a cell which is connected to a vertex via a locally + // owned cell) + const auto is_locally_relevant_on_level = [&](const auto &cell) { + for (const auto v : cell->vertex_indices()) + if (vertices_owned_by_locally_owned_cells_on_level + [cell->vertex_index(v)]) + return true; + return false; + }; + + // mark all locally relevant cells + for (const auto &cell : tria.cell_iterators_on_level(level)) + if (is_locally_relevant_on_level(cell)) + set_user_flag_and_of_its_parents(cell); + } + + // 2) set_up coarse-grid triangulation + { + std::vector vertices_locally_relevant(tria.n_vertices(), + false); + + // a) loop over all cells + for (const auto &cell : tria.cell_iterators_on_level(0)) + { + if (!cell->user_flag_set()) + continue; + + // extract cell definition (with old numbering of vertices) + dealii::CellData cell_data(cell->n_vertices()); + cell_data.material_id = cell->material_id(); + cell_data.manifold_id = cell->manifold_id(); + for (const auto v : cell->vertex_indices()) + cell_data.vertices[v] = cell->vertex_index(v); + construction_data.coarse_cells.push_back(cell_data); + + // save indices of each vertex of this cell + for (const auto v : cell->vertex_indices()) + vertices_locally_relevant[cell->vertex_index(v)] = true; + + // save translation for corase grid: lid -> gid + construction_data.coarse_cell_index_to_coarse_cell_id.push_back( + cell->id().get_coarse_cell_id()); + } + + add_vertices(construction_data, vertices_locally_relevant); + } + + + // 3) collect info of each cell + construction_data.cell_infos.resize( + tria.get_triangulation().n_global_levels()); + + // collect local vertices on active level + std::vector vertices_owned_by_locally_owned_active_cells( + tria.n_vertices()); + for (const auto &cell : tria.active_cell_iterators()) + if (subdomain_id_function(cell) == my_rank) + add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( + cell, vertices_owned_by_locally_owned_active_cells); + + // helper function to determine if cell is locally relevant + // on active level + const auto is_locally_relevant_on_active_level = + [&](const auto &cell) { + if (cell->is_active()) + for (const auto v : cell->vertex_indices()) + if (vertices_owned_by_locally_owned_active_cells + [cell->vertex_index(v)]) + return true; + return false; + }; + + for (unsigned int level = 0; + level < tria.get_triangulation().n_global_levels(); + ++level) + { + // collect local vertices on level + std::vector vertices_owned_by_locally_owned_cells_on_level( + tria.n_vertices()); + for (const auto &cell : tria.cell_iterators_on_level(level)) + if ((construct_multigrid && + (level_subdomain_id_function(cell) == my_rank)) || + (cell->is_active() && + subdomain_id_function(cell) == my_rank)) + add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( + cell, vertices_owned_by_locally_owned_cells_on_level); + + // helper function to determine if cell is locally relevant + // on level + const auto is_locally_relevant_on_level = [&](const auto &cell) { + for (const auto v : cell->vertex_indices()) + if (vertices_owned_by_locally_owned_cells_on_level + [cell->vertex_index(v)]) + return true; + return false; + }; + + auto &level_cell_infos = construction_data.cell_infos[level]; + for (const auto &cell : tria.cell_iterators_on_level(level)) + { + // check if cell is locally relevant + if (!(cell->user_flag_set())) + continue; + + CellData cell_info; + + // save coarse-cell id + cell_info.id = cell->id().template to_binary(); + + // save boundary_ids of each face of this cell + for (const auto f : cell->face_indices()) + { + types::boundary_id boundary_ind = + cell->face(f)->boundary_id(); + if (boundary_ind != numbers::internal_face_boundary_id) + cell_info.boundary_ids.emplace_back(f, boundary_ind); + } + + // save manifold id + { + // ... of cell + cell_info.manifold_id = cell->manifold_id(); + + // ... of lines + if (dim >= 2) + for (const auto line : cell->line_indices()) + cell_info.manifold_line_ids[line] = + cell->line(line)->manifold_id(); + + // ... of quads + if (dim == 3) + for (const auto f : cell->face_indices()) + cell_info.manifold_quad_ids[f] = + cell->quad(f)->manifold_id(); + } + + // subdomain and level subdomain id + cell_info.subdomain_id = numbers::artificial_subdomain_id; + cell_info.level_subdomain_id = + numbers::artificial_subdomain_id; + + if (is_locally_relevant_on_active_level(cell)) + { + cell_info.subdomain_id = subdomain_id_function(cell); + + cell_info.level_subdomain_id = + level_subdomain_id_function(cell); + } + else if (is_locally_relevant_on_level(cell)) + { + cell_info.level_subdomain_id = + level_subdomain_id_function(cell); + } + else + { + // cell is locally relevant but an artificial cell + } + + level_cell_infos.emplace_back(cell_info); + } + } + + const_cast &>(tria) + .load_user_flags(old_user_flags); + + return construction_data; + } + + private: + void + set_additional_data(Description &construction_data) const + { + construction_data.comm = comm; + construction_data.smoothing = tria.get_mesh_smoothing(); + construction_data.settings = settings; + } + + void + set_additional_data(DescriptionTemp &) const + { + // nothing to do + } + + void + add_vertices(Description &construction_data, + const std::vector &vertices_locally_relevant) const + { + std::vector vertices_locally_relevant_indices( + vertices_locally_relevant.size()); + + // enumerate locally relevant vertices + unsigned int vertex_counter = 0; + for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i) + if (vertices_locally_relevant[i]) + { + construction_data.coarse_cell_vertices.push_back( + tria.get_vertices()[i]); + vertices_locally_relevant_indices[i] = vertex_counter++; + } + + // correct vertices of cells (make them local) + for (auto &cell : construction_data.coarse_cells) + for (unsigned int v = 0; v < cell.vertices.size(); ++v) + cell.vertices[v] = + vertices_locally_relevant_indices[cell.vertices[v]]; + } + + void + add_vertices(DescriptionTemp &construction_data, + const std::vector &vertices_locally_relevant) const + { + for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i) + if (vertices_locally_relevant[i]) + construction_data.coarse_cell_vertices.emplace_back( + i, tria.get_vertices()[i]); + } + + + const dealii::Triangulation &tria; + const std::function::cell_iterator &)> + subdomain_id_function; + const std::function::cell_iterator &)> + level_subdomain_id_function; + + const MPI_Comm & comm; + const TriangulationDescription::Settings settings; + const bool construct_multigrid; + + std::map> + coinciding_vertex_groups; + std::map vertex_to_coinciding_vertex_group; + }; + } // namespace @@ -105,7 +682,7 @@ namespace TriangulationDescription const TriangulationDescription::Settings settings, const unsigned int my_rank_in) { - if (auto tria_pdt = dynamic_cast< + if (const auto tria_pdt = dynamic_cast< const parallel::distributed::Triangulation *>(&tria)) Assert(comm == tria_pdt->get_communicator(), ExcMessage("MPI communicators do not match.")); @@ -141,255 +718,22 @@ namespace TriangulationDescription ExcMessage("This type of triangulation is not supported!")); } - Description construction_data; - - // store the communicator - construction_data.comm = comm; - - std::map> - coinciding_vertex_groups; - std::map vertex_to_coinciding_vertex_group; - - GridTools::collect_coinciding_vertices(tria, - coinciding_vertex_groups, - vertex_to_coinciding_vertex_group); - - // helper function, which collects all vertices belonging to a cell - // (also taking into account periodicity) - const auto add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells = - [&coinciding_vertex_groups, &vertex_to_coinciding_vertex_group]( - const auto & cell, - std::vector &vertices_owned_by_locally_owned_cells) { - // add local vertices - for (const auto v : cell->vertex_indices()) - { - vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] = - true; - const auto coinciding_vertex_group = - vertex_to_coinciding_vertex_group.find(cell->vertex_index(v)); - if (coinciding_vertex_group != - vertex_to_coinciding_vertex_group.end()) - for (const auto &co_vertex : coinciding_vertex_groups.at( - coinciding_vertex_group->second)) - vertices_owned_by_locally_owned_cells[co_vertex] = true; - } - }; - - construction_data.smoothing = tria.get_mesh_smoothing(); - construction_data.settings = settings; - - const bool construct_multigrid = - settings & - TriangulationDescription::Settings::construct_multigrid_hierarchy; - - Assert( - !(settings & - TriangulationDescription::Settings::construct_multigrid_hierarchy) || - (tria.get_mesh_smoothing() & - Triangulation::limit_level_difference_at_vertices), - ExcMessage( - "Source triangulation has to be set up with " - "limit_level_difference_at_vertices if the construction of the " - "multigrid hierarchy is requested!")); - - // 1) collect locally relevant cells (set user_flag) - std::vector old_user_flags; - tria.save_user_flags(old_user_flags); - - // 1a) clear user_flags - const_cast &>(tria) - .clear_user_flags(); - - // 1b) loop over levels (from fine to coarse) and mark on each level - // the locally relevant cells - for (int level = tria.get_triangulation().n_global_levels() - 1; - level >= 0; - --level) - { - // collect vertices connected to a (on any level) locally owned - // cell - std::vector vertices_owned_by_locally_owned_cells_on_level( - tria.n_vertices()); - for (const auto &cell : tria.cell_iterators_on_level(level)) - if (construct_multigrid && (cell->level_subdomain_id() == my_rank)) - add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( - cell, vertices_owned_by_locally_owned_cells_on_level); - - for (const auto &cell : tria.active_cell_iterators()) - if (cell->subdomain_id() == my_rank) - add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( - cell, vertices_owned_by_locally_owned_cells_on_level); - - // helper function to determine if cell is locally relevant - // (i.e. a cell which is connected to a vertex via a locally - // owned cell) - const auto is_locally_relevant_on_level = [&](const auto &cell) { - for (const auto v : cell->vertex_indices()) - if (vertices_owned_by_locally_owned_cells_on_level - [cell->vertex_index(v)]) - return true; - return false; - }; - - // mark all locally relevant cells - for (const auto &cell : tria.cell_iterators_on_level(level)) - if (is_locally_relevant_on_level(cell)) - set_user_flag_and_of_its_parents(cell); - } - - // 2) set_up coarse-grid triangulation - { - std::map vertices_locally_relevant; - - // a) loop over all cells - for (const auto &cell : tria.cell_iterators_on_level(0)) - { - if (!cell->user_flag_set()) - continue; - - // extract cell definition (with old numbering of vertices) - dealii::CellData cell_data(cell->n_vertices()); - cell_data.material_id = cell->material_id(); - cell_data.manifold_id = cell->manifold_id(); - for (const auto v : cell->vertex_indices()) - cell_data.vertices[v] = cell->vertex_index(v); - construction_data.coarse_cells.push_back(cell_data); - - // save indices of each vertex of this cell - for (const auto v : cell->vertex_indices()) - vertices_locally_relevant[cell->vertex_index(v)] = - numbers::invalid_unsigned_int; - - // save translation for corase grid: lid -> gid - construction_data.coarse_cell_index_to_coarse_cell_id.push_back( - cell->id().get_coarse_cell_id()); - } - - // b) enumerate locally relevant vertices - unsigned int vertex_counter = 0; - for (auto &vertex : vertices_locally_relevant) - { - construction_data.coarse_cell_vertices.push_back( - tria.get_vertices()[vertex.first]); - vertex.second = vertex_counter++; - } - - // c) correct vertices of cells (make them local) - for (auto &cell : construction_data.coarse_cells) - for (unsigned int v = 0; v < cell.vertices.size(); ++v) - cell.vertices[v] = vertices_locally_relevant[cell.vertices[v]]; - } - - - // 3) collect info of each cell - construction_data.cell_infos.resize( - tria.get_triangulation().n_global_levels()); - - // collect local vertices on active level - std::vector vertices_owned_by_locally_owned_active_cells( - tria.n_vertices()); - for (const auto &cell : tria.active_cell_iterators()) - if (cell->subdomain_id() == my_rank) - add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( - cell, vertices_owned_by_locally_owned_active_cells); - - // helper function to determine if cell is locally relevant - // on active level - const auto is_locally_relevant_on_active_level = [&](const auto &cell) { - if (cell->is_active()) - for (const auto v : cell->vertex_indices()) - if (vertices_owned_by_locally_owned_active_cells[cell->vertex_index( - v)]) - return true; - return false; + const auto subdomain_id_function = [](const auto &cell) { + return cell->subdomain_id(); }; - for (unsigned int level = 0; - level < tria.get_triangulation().n_global_levels(); - ++level) - { - // collect local vertices on level - std::vector vertices_owned_by_locally_owned_cells_on_level( - tria.n_vertices()); - for (const auto &cell : tria.cell_iterators_on_level(level)) - if ((construct_multigrid && - (cell->level_subdomain_id() == my_rank)) || - (cell->is_active() && cell->subdomain_id() == my_rank)) - add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( - cell, vertices_owned_by_locally_owned_cells_on_level); - - // helper function to determine if cell is locally relevant - // on level - const auto is_locally_relevant_on_level = [&](const auto &cell) { - for (const auto v : cell->vertex_indices()) - if (vertices_owned_by_locally_owned_cells_on_level - [cell->vertex_index(v)]) - return true; - return false; - }; - - auto &level_cell_infos = construction_data.cell_infos[level]; - for (const auto &cell : tria.cell_iterators_on_level(level)) - { - // check if cell is locally relevant - if (!(cell->user_flag_set())) - continue; - - CellData cell_info; - - // save coarse-cell id - cell_info.id = cell->id().template to_binary(); - - // save boundary_ids of each face of this cell - for (const auto f : cell->face_indices()) - { - types::boundary_id boundary_ind = - cell->face(f)->boundary_id(); - if (boundary_ind != numbers::internal_face_boundary_id) - cell_info.boundary_ids.emplace_back(f, boundary_ind); - } - - // save manifold id - { - // ... of cell - cell_info.manifold_id = cell->manifold_id(); - - // ... of lines - if (dim >= 2) - for (const auto line : cell->line_indices()) - cell_info.manifold_line_ids[line] = - cell->line(line)->manifold_id(); - - // ... of quads - if (dim == 3) - for (const auto f : cell->face_indices()) - cell_info.manifold_quad_ids[f] = - cell->quad(f)->manifold_id(); - } - - // subdomain and level subdomain id - cell_info.subdomain_id = numbers::artificial_subdomain_id; - cell_info.level_subdomain_id = numbers::artificial_subdomain_id; - - if (is_locally_relevant_on_active_level(cell)) - { - cell_info.level_subdomain_id = cell->level_subdomain_id(); - cell_info.subdomain_id = cell->subdomain_id(); - } - else if (is_locally_relevant_on_level(cell)) - { - cell_info.level_subdomain_id = cell->level_subdomain_id(); - } - // else: cell is locally relevant but an artificial cell - - level_cell_infos.emplace_back(cell_info); - } - } - - const_cast &>(tria).load_user_flags( - old_user_flags); + const auto level_subdomain_id_function = [](const auto &cell) { + return cell->level_subdomain_id(); + }; - return construction_data; + return CreateDescriptionFromTriangulationHelper( + tria, + subdomain_id_function, + level_subdomain_id_function, + comm, + settings) + .template create_description_for_rank>( + my_rank); } @@ -522,211 +866,6 @@ namespace TriangulationDescription #endif } - template - struct DescriptionTemp - { - template - void - serialize(Archive &ar, const unsigned int /*version*/) - { - ar &coarse_cells; - ar &coarse_cell_vertices; - ar &coarse_cell_index_to_coarse_cell_id; - ar &cell_infos; - } - - void - merge(const DescriptionTemp &other) - { - this->cell_infos.resize( - std::max(other.cell_infos.size(), this->cell_infos.size())); - - this->coarse_cells.insert(this->coarse_cells.end(), - other.coarse_cells.begin(), - other.coarse_cells.end()); - this->coarse_cell_vertices.insert(this->coarse_cell_vertices.end(), - other.coarse_cell_vertices.begin(), - other.coarse_cell_vertices.end()); - this->coarse_cell_index_to_coarse_cell_id.insert( - this->coarse_cell_index_to_coarse_cell_id.end(), - other.coarse_cell_index_to_coarse_cell_id.begin(), - other.coarse_cell_index_to_coarse_cell_id.end()); - - for (unsigned int i = 0; i < this->cell_infos.size(); ++i) - this->cell_infos[i].insert(this->cell_infos[i].end(), - other.cell_infos[i].begin(), - other.cell_infos[i].end()); - } - - void - reduce() - { - // make coarse cells unique - { - std::vector, - unsigned int>> - temp; - - for (unsigned int i = 0; i < this->coarse_cells.size(); ++i) - temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i], - this->coarse_cells[i], - i); - - std::sort(temp.begin(), temp.end(), [](const auto &a, const auto &b) { - return std::get<0>(a) < std::get<0>(b); - }); - temp.erase(std::unique(temp.begin(), - temp.end(), - [](const auto &a, const auto &b) { - return std::get<0>(a) == std::get<0>(b); - }), - temp.end()); - std::sort(temp.begin(), temp.end(), [](const auto &a, const auto &b) { - return std::get<2>(a) < std::get<2>(b); - }); - - this->coarse_cell_index_to_coarse_cell_id.resize(temp.size()); - this->coarse_cells.resize(temp.size()); - - for (unsigned int i = 0; i < temp.size(); ++i) - { - this->coarse_cell_index_to_coarse_cell_id[i] = - std::get<0>(temp[i]); - this->coarse_cells[i] = std::get<1>(temp[i]); - } - } - - // make coarse cell vertices unique - { - std::sort(this->coarse_cell_vertices.begin(), - this->coarse_cell_vertices.end(), - [](const auto &a, const auto &b) { - return a.first < b.first; - }); - this->coarse_cell_vertices.erase( - std::unique(this->coarse_cell_vertices.begin(), - this->coarse_cell_vertices.end(), - [](const auto &a, const auto &b) { - return a.first == b.first; - }), - this->coarse_cell_vertices.end()); - } - - // make cells unique - for (unsigned int i = 0; i < this->cell_infos.size(); ++i) - { - std::sort(this->cell_infos[i].begin(), - this->cell_infos[i].end(), - [](const auto &a, const auto &b) { return a.id < b.id; }); - this->cell_infos[i].erase(std::unique(this->cell_infos[i].begin(), - this->cell_infos[i].end(), - [](const auto &a, - const auto &b) { - return a.id == b.id; - }), - this->cell_infos[i].end()); - } - } - - Description - convert(const MPI_Comm comm, - const typename Triangulation::MeshSmoothing - mesh_smoothing) - { - Description description; - - // copy communicator - description.comm = comm; - - // set settings (no multigrid levels are set for now) - description.settings = - TriangulationDescription::Settings::default_setting; - - // use mesh smoothing from base triangulation - description.smoothing = mesh_smoothing; - - std::map map; - - for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i) - { - description.coarse_cell_vertices.push_back( - this->coarse_cell_vertices[i].second); - map[this->coarse_cell_vertices[i].first] = i; - } - - description.coarse_cells = this->coarse_cells; - - for (auto &cell : description.coarse_cells) - for (unsigned int v = 0; v < cell.vertices.size(); ++v) - cell.vertices[v] = map[cell.vertices[v]]; - - description.coarse_cell_index_to_coarse_cell_id = - this->coarse_cell_index_to_coarse_cell_id; - description.cell_infos = this->cell_infos; - - return description; - } - - std::vector> coarse_cells; - - std::vector>> - coarse_cell_vertices; - - std::vector coarse_cell_index_to_coarse_cell_id; - - std::vector>> cell_infos; - }; - - - - template - void - fill_cell_infos(const TriaIterator> &cell, - std::vector>> & cell_infos, - const LinearAlgebra::distributed::Vector &partition) - { - CellData cell_info; - - // save coarse-cell id - cell_info.id = cell->id().template to_binary(); - - // save boundary_ids of each face of this cell - for (const auto f : cell->face_indices()) - { - types::boundary_id boundary_ind = cell->face(f)->boundary_id(); - if (boundary_ind != numbers::internal_face_boundary_id) - cell_info.boundary_ids.emplace_back(f, boundary_ind); - } - - // save manifold id - { - // ... of cell - cell_info.manifold_id = cell->manifold_id(); - - // ... of lines - if (dim >= 2) - for (const auto line : cell->line_indices()) - cell_info.manifold_line_ids[line] = cell->line(line)->manifold_id(); - - // ... of quads - if (dim == 3) - for (const auto f : cell->face_indices()) - cell_info.manifold_quad_ids[f] = cell->quad(f)->manifold_id(); - } - - // subdomain and level subdomain id - if (cell->is_active() && (cell->is_artificial() == false)) - cell_info.subdomain_id = static_cast( - partition[cell->global_active_cell_index()]); - else - cell_info.subdomain_id = numbers::artificial_subdomain_id; - - cell_info.level_subdomain_id = numbers::artificial_subdomain_id; - - cell_infos[cell->level()].emplace_back(cell_info); - } - template @@ -741,8 +880,10 @@ namespace TriangulationDescription #endif if (partition.size() == 0) - return create_description_from_triangulation(tria, - tria.get_communicator()); + { + return create_description_from_triangulation(tria, + tria.get_communicator()); + } partition.update_ghost_values(); @@ -758,164 +899,52 @@ namespace TriangulationDescription relevant_processes.end()); }(); - // 2) determine coinciding vertices (important for periodicity) - std::map> - coinciding_vertex_groups; - std::map vertex_to_coinciding_vertex_group; - - GridTools::collect_coinciding_vertices(tria, - coinciding_vertex_groups, - vertex_to_coinciding_vertex_group); - - const auto add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells = - [&coinciding_vertex_groups, &vertex_to_coinciding_vertex_group]( - const auto & cell, - std::vector &vertices_owned_by_locally_owned_cells) { - for (const auto v : cell->vertex_indices()) - { - vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] = - true; - const auto coinciding_vertex_group = - vertex_to_coinciding_vertex_group.find(cell->vertex_index(v)); - if (coinciding_vertex_group != - vertex_to_coinciding_vertex_group.end()) - for (const auto &co_vertex : coinciding_vertex_groups.at( - coinciding_vertex_group->second)) - vertices_owned_by_locally_owned_cells[co_vertex] = true; - } - }; - - // 3) create a description (locally owned cell and a layer of ghost cells - // and all their parents) - std::vector> description_temp( - relevant_processes.size()); + TriangulationDescription::Settings settings = + TriangulationDescription::Settings::default_setting; - std::vector old_user_flags; - if (description_temp.size() > 0) - tria.save_user_flags(old_user_flags); - - for (unsigned int i = 0; i < description_temp.size(); ++i) - { - const unsigned int proc = relevant_processes[i]; - auto & description_temp_i = description_temp[i]; - description_temp_i.cell_infos.resize( - tria.get_triangulation().n_global_levels()); - - // clear user_flags - const_cast &>(tria) - .clear_user_flags(); - - for (int level = tria.get_triangulation().n_global_levels() - 1; - level >= 0; - --level) - { - // mark all vertices attached to locally owned cells - std::vector vertices_owned_by_locally_owned_cells_on_level( - tria.n_vertices()); - for (const auto &cell : tria.active_cell_iterators()) - if (cell->is_active() && cell->is_locally_owned() && - static_cast( - partition[cell->global_active_cell_index()]) == proc) - add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells( - cell, vertices_owned_by_locally_owned_cells_on_level); - - // helper function if a cell is locally relevant (active and - // connected to cell via a vertex) - const auto is_locally_relevant_on_level = [&](const auto &cell) { - for (const auto v : cell->vertex_indices()) - if (vertices_owned_by_locally_owned_cells_on_level - [cell->vertex_index(v)]) - return true; - return false; - }; - - // collect locally relevant cells (including their parents) - for (const auto &cell : tria.cell_iterators_on_level(level)) - if (is_locally_relevant_on_level(cell)) - set_user_flag_and_of_its_parents(cell); - } - - for (unsigned int level = 0; level < tria.n_global_levels(); ++level) - for (const auto &cell : tria.cell_iterators_on_level(level)) - if (cell->user_flag_set()) - fill_cell_infos(cell, description_temp_i.cell_infos, partition); - - // collect coarse-grid cells - std::vector vertices_locally_relevant(tria.n_vertices(), false); + const auto subdomain_id_function = [&partition](const auto &cell) { + if ((cell->is_active() && (cell->is_artificial() == false))) + return static_cast( + partition[cell->global_active_cell_index()]); + else + return numbers::artificial_subdomain_id; + }; - for (const auto &cell : tria.cell_iterators_on_level(0)) - { - if (!cell->user_flag_set()) - continue; + const auto level_subdomain_id_function = [](const auto & /*cell*/) { + return numbers::artificial_subdomain_id; + }; - // extract cell definition - dealii::CellData cell_data(cell->n_vertices()); - cell_data.material_id = cell->material_id(); - cell_data.manifold_id = cell->manifold_id(); - for (const auto v : cell->vertex_indices()) - cell_data.vertices[v] = cell->vertex_index(v); - description_temp_i.coarse_cells.push_back(cell_data); - - // mark cell vertices as relevant - for (const auto v : cell->vertex_indices()) - vertices_locally_relevant[cell->vertex_index(v)] = true; - - // save translation for corase grid: lid -> gid - description_temp_i.coarse_cell_index_to_coarse_cell_id.push_back( - cell->id().get_coarse_cell_id()); - } + CreateDescriptionFromTriangulationHelper helper( + tria, + subdomain_id_function, + level_subdomain_id_function, + tria.get_communicator(), + settings); - // collect coarse-grid vertices - for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i) - if (vertices_locally_relevant[i]) - description_temp_i.coarse_cell_vertices.emplace_back( - i, tria.get_vertices()[i]); - } + // create a description (locally owned cell and a layer of ghost cells + // and all their parents) + std::vector> descriptions_per_rank; + descriptions_per_rank.reserve(relevant_processes.size()); - // restore flags - if (description_temp.size() > 0) - const_cast &>(tria) - .load_user_flags(old_user_flags); + for (const auto rank : relevant_processes) + descriptions_per_rank.emplace_back( + helper.template create_description_for_rank< + DescriptionTemp>(rank)); // collect description from all processes that used to own locally-owned // active cells of this process in a single description DescriptionTemp description_merged; - - dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess - process([&]() { return relevant_processes; }, - [&](const unsigned int other_rank, - std::vector &send_buffer) { - const auto ptr = std::find(relevant_processes.begin(), - relevant_processes.end(), - other_rank); - - Assert(ptr != relevant_processes.end(), ExcInternalError()); - - const auto other_rank_index = - std::distance(relevant_processes.begin(), ptr); - - send_buffer = - dealii::Utilities::pack(description_temp[other_rank_index], - false); - }, - [&](const unsigned int &, - const std::vector &recv_buffer, - std::vector &) { - description_merged.merge( - dealii::Utilities::unpack>( - recv_buffer, false)); - }); - - dealii::Utilities::MPI::ConsensusAlgorithms::Selector( - process, partition.get_mpi_communicator()) - .run(); + description_merged.collect(relevant_processes, + descriptions_per_rank, + partition.get_mpi_communicator()); // remove redundant entries description_merged.reduce(); // convert to actual description return description_merged.convert(partition.get_mpi_communicator(), - tria.get_mesh_smoothing()); + tria.get_mesh_smoothing(), + settings); } } // namespace Utilities -- 2.39.5