#include <deal.II/base/geometry_info.h>
#include <deal.II/base/mpi.h>
+#include <deal.II/base/mpi_consensus_algorithms.h>
#include <deal.II/distributed/tria.h>
#endif
}
+ template <int dim, int spacedim>
+ struct DescriptionTemp
+ {
+ template <class Archive>
+ 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<dim, spacedim> &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<std::tuple<types::coarse_cell_id,
+ dealii::CellData<dim>,
+ 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<dim, spacedim>
+ convert(const MPI_Comm comm,
+ const typename Triangulation<dim, spacedim>::MeshSmoothing
+ mesh_smoothing)
+ {
+ Description<dim, spacedim> 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<unsigned int, unsigned int> 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<dealii::CellData<dim>> coarse_cells;
+
+ std::vector<std::pair<unsigned int, Point<spacedim>>>
+ coarse_cell_vertices;
+
+ std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
+
+ std::vector<std::vector<CellData<dim>>> cell_infos;
+ };
+
+
+
+ template <int dim, int spacedim>
+ void
+ fill_cell_infos(const TriaIterator<CellAccessor<dim, spacedim>> &cell,
+ std::vector<std::vector<CellData<dim>>> & cell_infos,
+ const LinearAlgebra::distributed::Vector<double> &partition)
+ {
+ if (cell->user_flag_set())
+ return; // cell has been already added -> nothing to do
+
+ cell->set_user_flag();
+
+ CellData<dim> cell_info;
+
+ // save coarse-cell id
+ cell_info.id = cell->id().template to_binary<dim>();
+
+ // 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_info.subdomain_id = static_cast<unsigned int>(
+ 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);
+
+ if (cell->level() != 0) // proceed with parent
+ fill_cell_infos(cell->parent(), cell_infos, partition);
+ }
+
+
+
+ template <int dim, int spacedim>
+ Description<dim, spacedim>
+ create_description_from_triangulation(
+ const Triangulation<dim, spacedim> & tria,
+ const LinearAlgebra::distributed::Vector<double> &partition)
+ {
+ // 1) determine processes owning locally owned cells
+ const std::vector<unsigned int> relevant_processes = [&]() {
+ std::set<unsigned int> relevant_processes;
+
+ for (unsigned int i = 0; i < partition.local_size(); ++i)
+ relevant_processes.insert(
+ static_cast<unsigned int>(partition.local_element(i)));
+
+ return std::vector<unsigned int>(relevant_processes.begin(),
+ relevant_processes.end());
+ }();
+
+ // 2) determine coinciding vertices (important for periodicity)
+ std::map<unsigned int, std::vector<unsigned int>>
+ coinciding_vertex_groups;
+ std::map<unsigned int, unsigned int> 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<bool> &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<DescriptionTemp<dim, spacedim>> description_temp(
+ relevant_processes.size());
+
+ 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
+ std::vector<bool> old_user_flags;
+ tria.save_user_flags(old_user_flags);
+
+ const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
+ .clear_user_flags();
+
+ // mark all vertices attached to locally owned cells
+ std::vector<bool> vertices_owned_by_locally_owned_cells(
+ tria.n_vertices());
+ for (const auto &cell : tria.active_cell_iterators())
+ if (cell->is_locally_owned() &&
+ static_cast<unsigned int>(
+ 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);
+
+ // 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[cell->vertex_index(v)])
+ return true;
+ return false;
+ };
+
+ // collect locally relevant cells (including their parents)
+ for (const auto &cell : tria.active_cell_iterators())
+ if (is_locally_relevant_on_level(cell))
+ fill_cell_infos(cell, description_temp_i.cell_infos, partition);
+
+ // collect coarse-grid cells
+ std::vector<bool> vertices_locally_relevant(tria.n_vertices(), false);
+
+ for (const auto &cell : tria.cell_iterators_on_level(0))
+ {
+ if (!cell->user_flag_set())
+ continue;
+
+ // extract cell definition
+ dealii::CellData<dim> 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());
+ }
+
+ // 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]);
+
+ // restore flags
+ const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
+ .load_user_flags(old_user_flags);
+ }
+
+ // collect description from all processes that used to own locally-owned
+ // active cells of this process in a single description
+ DescriptionTemp<dim, spacedim> description_merged;
+
+ dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess<char, char>
+ process([&]() { return relevant_processes; },
+ [&](const unsigned int other_rank,
+ std::vector<char> &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<char> &recv_buffer,
+ std::vector<char> &) {
+ description_merged.merge(
+ dealii::Utilities::unpack<DescriptionTemp<dim, spacedim>>(
+ recv_buffer, false));
+ });
+
+ dealii::Utilities::MPI::ConsensusAlgorithms::Selector<char, char>(
+ process, tria.get_communicator())
+ .run();
+
+ // remove redundant entries
+ description_merged.reduce();
+
+ // convert to actual description
+ return description_merged.convert(tria.get_communicator(),
+ tria.get_mesh_smoothing());
+ }
+
} // namespace Utilities
} // namespace TriangulationDescription
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test
+// TriangulationDescription::Utilities::create_description_from_triangulation()
+// with repartitioning capabilities.
+
+#include <deal.II/base/mpi_consensus_algorithms.h>
+
+#include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include "./tests.h"
+
+using namespace dealii;
+
+
+template <int dim, int spacedim>
+LinearAlgebra::distributed::Vector<double>
+partition_distributed_triangulation(const Triangulation<dim, spacedim> &tria_in,
+ const unsigned int n_partitions)
+{
+ const auto tria =
+ dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria_in);
+
+ Assert(tria, ExcNotImplemented());
+
+ LinearAlgebra::distributed::Vector<double> partition(
+ tria->global_active_cell_index_partitioner().lock());
+
+ for (const auto &cell : tria_in.active_cell_iterators())
+ if (cell->is_locally_owned())
+ partition[cell->global_active_cell_index()] =
+ std::floor(cell->center()[0] * n_partitions);
+
+ partition.update_ghost_values();
+
+ return partition;
+}
+
+
+template <int dim>
+void
+test(const MPI_Comm comm, const unsigned int n_partitions)
+{
+ parallel::distributed::Triangulation<dim> tria(comm);
+ GridGenerator::subdivided_hyper_cube(tria, 4);
+ tria.refine_global(3);
+
+ const auto partition_new =
+ partition_distributed_triangulation(tria, n_partitions);
+
+ // repartition triangulation so that it has strided partitioning
+ const auto construction_data =
+ TriangulationDescription::Utilities::create_description_from_triangulation(
+ tria, partition_new);
+
+ parallel::fullydistributed::Triangulation<dim> tria_pft(comm);
+ tria_pft.create_triangulation(construction_data);
+
+ FE_Q<dim> fe(2);
+ DoFHandler<dim> dof_handler(tria_pft);
+ dof_handler.distribute_dofs(fe);
+
+ // print statistics
+ print_statistics(tria_pft);
+ print_statistics(dof_handler);
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+ MPILogInitAll all;
+
+ MPI_Comm comm = MPI_COMM_WORLD;
+
+ deallog.push("all");
+ test<2>(comm, Utilities::MPI::n_mpi_processes(comm));
+ deallog.pop();
+
+ // test that we can eliminate processes
+ deallog.push("reduced");
+ test<2>(comm,
+ std::max<unsigned int>(1, Utilities::MPI::n_mpi_processes(comm) / 2));
+ deallog.pop();
+}
projects / dealii.git / commitdiff