1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2015 - 2020 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 #include <deal.II/base/mpi.h> 17 #include <deal.II/base/utilities.h> 18 19 #include <deal.II/distributed/shared_tria.h> 20 #include <deal.II/distributed/tria.h> 21 22 #include <deal.II/grid/filtered_iterator.h> 23 #include <deal.II/grid/grid_tools.h> 24 #include <deal.II/grid/tria.h> 25 #include <deal.II/grid/tria_accessor.h> 26 #include <deal.II/grid/tria_iterator.h> 27 28 #include <deal.II/lac/sparsity_tools.h> 29 30 31 DEAL_II_NAMESPACE_OPEN 32 33 #ifdef DEAL_II_WITH_MPI 34 namespace parallel 35 { 36 namespace shared 37 { 38 template <int dim, int spacedim> Triangulation(MPI_Comm mpi_communicator,const typename dealii::Triangulation<dim,spacedim>::MeshSmoothing smooth_grid,const bool allow_artificial_cells,const Settings settings)39 Triangulation<dim, spacedim>::Triangulation( 40 MPI_Comm mpi_communicator, 41 const typename dealii::Triangulation<dim, spacedim>::MeshSmoothing 42 smooth_grid, 43 const bool allow_artificial_cells, 44 const Settings settings) 45 : dealii::parallel::TriangulationBase<dim, spacedim>(mpi_communicator, 46 smooth_grid, 47 false) 48 , settings(settings) 49 , allow_artificial_cells(allow_artificial_cells) 50 { 51 const auto partition_settings = 52 (partition_zoltan | partition_metis | partition_zorder | 53 partition_custom_signal) & 54 settings; 55 (void)partition_settings; 56 Assert(partition_settings == partition_auto || 57 partition_settings == partition_metis || 58 partition_settings == partition_zoltan || 59 partition_settings == partition_zorder || 60 partition_settings == partition_custom_signal, 61 ExcMessage("Settings must contain exactly one type of the active " 62 "cell partitioning scheme.")); 63 64 if (settings & construct_multigrid_hierarchy) 65 Assert(allow_artificial_cells, 66 ExcMessage("construct_multigrid_hierarchy requires " 67 "allow_artificial_cells to be set to true.")); 68 } 69 70 71 72 template <int dim, int spacedim> 73 bool is_multilevel_hierarchy_constructed() const74 Triangulation<dim, spacedim>::is_multilevel_hierarchy_constructed() const 75 { 76 return (settings & construct_multigrid_hierarchy); 77 } 78 79 80 81 template <int dim, int spacedim> 82 void partition()83 Triangulation<dim, spacedim>::partition() 84 { 85 # ifdef DEBUG 86 // Check that all meshes are the same (or at least have the same 87 // total number of active cells): 88 const unsigned int max_active_cells = 89 Utilities::MPI::max(this->n_active_cells(), this->get_communicator()); 90 Assert( 91 max_active_cells == this->n_active_cells(), 92 ExcMessage( 93 "A parallel::shared::Triangulation needs to be refined in the same " 94 "way on all processors, but the participating processors don't " 95 "agree on the number of active cells.")); 96 # endif 97 98 auto partition_settings = (partition_zoltan | partition_metis | 99 partition_zorder | partition_custom_signal) & 100 settings; 101 if (partition_settings == partition_auto) 102 # ifdef DEAL_II_TRILINOS_WITH_ZOLTAN 103 partition_settings = partition_zoltan; 104 # elif defined DEAL_II_WITH_METIS 105 partition_settings = partition_metis; 106 # else 107 partition_settings = partition_zorder; 108 # endif 109 110 if (partition_settings == partition_zoltan) 111 { 112 # ifndef DEAL_II_TRILINOS_WITH_ZOLTAN 113 AssertThrow(false, 114 ExcMessage( 115 "Choosing 'partition_zoltan' requires the library " 116 "to be compiled with support for Zoltan! " 117 "Instead, you might use 'partition_auto' to select " 118 "a partitioning algorithm that is supported " 119 "by your current configuration.")); 120 # else 121 GridTools::partition_triangulation( 122 this->n_subdomains, *this, SparsityTools::Partitioner::zoltan); 123 # endif 124 } 125 else if (partition_settings == partition_metis) 126 { 127 # ifndef DEAL_II_WITH_METIS 128 AssertThrow(false, 129 ExcMessage( 130 "Choosing 'partition_metis' requires the library " 131 "to be compiled with support for METIS! " 132 "Instead, you might use 'partition_auto' to select " 133 "a partitioning algorithm that is supported " 134 "by your current configuration.")); 135 # else 136 GridTools::partition_triangulation(this->n_subdomains, 137 *this, 138 SparsityTools::Partitioner::metis); 139 # endif 140 } 141 else if (partition_settings == partition_zorder) 142 { 143 GridTools::partition_triangulation_zorder(this->n_subdomains, *this); 144 } 145 else if (partition_settings == partition_custom_signal) 146 { 147 // User partitions mesh manually 148 } 149 else 150 { 151 AssertThrow(false, ExcInternalError()) 152 } 153 154 // do not partition multigrid levels if user is 155 // defining a custom partition 156 if ((settings & construct_multigrid_hierarchy) && 157 !(settings & partition_custom_signal)) 158 dealii::GridTools::partition_multigrid_levels(*this); 159 160 true_subdomain_ids_of_cells.resize(this->n_active_cells()); 161 162 // loop over all cells and mark artificial: 163 typename parallel::shared::Triangulation<dim, 164 spacedim>::active_cell_iterator 165 cell = this->begin_active(), 166 endc = this->end(); 167 168 if (allow_artificial_cells) 169 { 170 // get active halo layer of (ghost) cells 171 // parallel::shared::Triangulation<dim>:: 172 std::function<bool( 173 const typename parallel::shared::Triangulation<dim, spacedim>:: 174 active_cell_iterator &)> 175 predicate = IteratorFilters::SubdomainEqualTo(this->my_subdomain); 176 177 const std::vector<typename parallel::shared::Triangulation< 178 dim, 179 spacedim>::active_cell_iterator> 180 active_halo_layer_vector = 181 dealii::GridTools::compute_active_cell_halo_layer(*this, 182 predicate); 183 std::set<typename parallel::shared::Triangulation<dim, spacedim>:: 184 active_cell_iterator> 185 active_halo_layer(active_halo_layer_vector.begin(), 186 active_halo_layer_vector.end()); 187 188 for (unsigned int index = 0; cell != endc; cell++, index++) 189 { 190 // store original/true subdomain ids: 191 true_subdomain_ids_of_cells[index] = cell->subdomain_id(); 192 193 if (cell->is_locally_owned() == false && 194 active_halo_layer.find(cell) == active_halo_layer.end()) 195 cell->set_subdomain_id(numbers::artificial_subdomain_id); 196 } 197 198 // loop over all cells in multigrid hierarchy and mark artificial: 199 if (settings & construct_multigrid_hierarchy) 200 { 201 true_level_subdomain_ids_of_cells.resize(this->n_levels()); 202 203 std::function<bool( 204 const typename parallel::shared::Triangulation<dim, spacedim>:: 205 cell_iterator &)> 206 predicate = IteratorFilters::LocallyOwnedLevelCell(); 207 for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl) 208 { 209 true_level_subdomain_ids_of_cells[lvl].resize( 210 this->n_cells(lvl)); 211 212 const std::vector<typename parallel::shared::Triangulation< 213 dim, 214 spacedim>::cell_iterator> 215 level_halo_layer_vector = 216 dealii::GridTools::compute_cell_halo_layer_on_level( 217 *this, predicate, lvl); 218 std::set<typename parallel::shared:: 219 Triangulation<dim, spacedim>::cell_iterator> 220 level_halo_layer(level_halo_layer_vector.begin(), 221 level_halo_layer_vector.end()); 222 223 typename parallel::shared::Triangulation<dim, spacedim>:: 224 cell_iterator cell = this->begin(lvl), 225 endc = this->end(lvl); 226 for (unsigned int index = 0; cell != endc; cell++, index++) 227 { 228 // Store true level subdomain IDs before setting 229 // artificial 230 true_level_subdomain_ids_of_cells[lvl][index] = 231 cell->level_subdomain_id(); 232 233 // for active cells, we must have knowledge of level 234 // subdomain ids of all neighbors to our subdomain, not 235 // just neighbors on the same level. if the cells 236 // subdomain id was not set to artitficial above, we will 237 // also keep its level subdomain id since it is either 238 // owned by this processor or in the ghost layer of the 239 // active mesh. 240 if (cell->is_active() && 241 cell->subdomain_id() != 242 numbers::artificial_subdomain_id) 243 continue; 244 245 // we must have knowledge of our parent in the hierarchy 246 if (cell->has_children()) 247 { 248 bool keep_cell = false; 249 for (unsigned int c = 0; 250 c < GeometryInfo<dim>::max_children_per_cell; 251 ++c) 252 if (cell->child(c)->level_subdomain_id() == 253 this->my_subdomain) 254 { 255 keep_cell = true; 256 break; 257 } 258 if (keep_cell) 259 continue; 260 } 261 262 // we must have knowledge of our neighbors on the same 263 // level 264 if (!cell->is_locally_owned_on_level() && 265 level_halo_layer.find(cell) != level_halo_layer.end()) 266 continue; 267 268 // mark all other cells to artificial 269 cell->set_level_subdomain_id( 270 numbers::artificial_subdomain_id); 271 } 272 } 273 } 274 } 275 else 276 { 277 // just store true subdomain ids 278 for (unsigned int index = 0; cell != endc; cell++, index++) 279 true_subdomain_ids_of_cells[index] = cell->subdomain_id(); 280 } 281 282 # ifdef DEBUG 283 { 284 // Assert that each cell is owned by a processor 285 unsigned int n_my_cells = 0; 286 typename parallel::shared::Triangulation<dim, 287 spacedim>::active_cell_iterator 288 cell = this->begin_active(), 289 endc = this->end(); 290 for (; cell != endc; ++cell) 291 if (cell->is_locally_owned()) 292 n_my_cells += 1; 293 294 const unsigned int total_cells = 295 Utilities::MPI::sum(n_my_cells, this->get_communicator()); 296 Assert(total_cells == this->n_active_cells(), 297 ExcMessage("Not all cells are assigned to a processor.")); 298 } 299 300 // If running with multigrid, assert that each level 301 // cell is owned by a processor 302 if (settings & construct_multigrid_hierarchy) 303 { 304 unsigned int n_my_cells = 0; 305 typename parallel::shared::Triangulation<dim, spacedim>::cell_iterator 306 cell = this->begin(), 307 endc = this->end(); 308 for (; cell != endc; ++cell) 309 if (cell->is_locally_owned_on_level()) 310 n_my_cells += 1; 311 312 const unsigned int total_cells = 313 Utilities::MPI::sum(n_my_cells, this->get_communicator()); 314 Assert(total_cells == this->n_cells(), 315 ExcMessage("Not all cells are assigned to a processor.")); 316 } 317 # endif 318 } 319 320 321 322 template <int dim, int spacedim> 323 bool with_artificial_cells() const324 Triangulation<dim, spacedim>::with_artificial_cells() const 325 { 326 return allow_artificial_cells; 327 } 328 329 330 331 template <int dim, int spacedim> 332 const std::vector<types::subdomain_id> & get_true_subdomain_ids_of_cells() const333 Triangulation<dim, spacedim>::get_true_subdomain_ids_of_cells() const 334 { 335 return true_subdomain_ids_of_cells; 336 } 337 338 339 340 template <int dim, int spacedim> 341 const std::vector<types::subdomain_id> & get_true_level_subdomain_ids_of_cells(const unsigned int level) const342 Triangulation<dim, spacedim>::get_true_level_subdomain_ids_of_cells( 343 const unsigned int level) const 344 { 345 Assert(level < true_level_subdomain_ids_of_cells.size(), 346 ExcInternalError()); 347 Assert(true_level_subdomain_ids_of_cells[level].size() == 348 this->n_cells(level), 349 ExcInternalError()); 350 return true_level_subdomain_ids_of_cells[level]; 351 } 352 353 354 355 template <int dim, int spacedim> 356 void execute_coarsening_and_refinement()357 Triangulation<dim, spacedim>::execute_coarsening_and_refinement() 358 { 359 dealii::Triangulation<dim, spacedim>::execute_coarsening_and_refinement(); 360 partition(); 361 this->update_number_cache(); 362 } 363 364 365 366 template <int dim, int spacedim> 367 void create_triangulation(const std::vector<Point<spacedim>> & vertices,const std::vector<CellData<dim>> & cells,const SubCellData & subcelldata)368 Triangulation<dim, spacedim>::create_triangulation( 369 const std::vector<Point<spacedim>> &vertices, 370 const std::vector<CellData<dim>> & cells, 371 const SubCellData & subcelldata) 372 { 373 try 374 { 375 dealii::Triangulation<dim, spacedim>::create_triangulation( 376 vertices, cells, subcelldata); 377 } 378 catch ( 379 const typename dealii::Triangulation<dim, spacedim>::DistortedCellList 380 &) 381 { 382 // the underlying triangulation should not be checking for distorted 383 // cells 384 Assert(false, ExcInternalError()); 385 } 386 partition(); 387 this->update_number_cache(); 388 } 389 390 391 392 template <int dim, int spacedim> 393 void create_triangulation(const TriangulationDescription::Description<dim,spacedim> & construction_data)394 Triangulation<dim, spacedim>::create_triangulation( 395 const TriangulationDescription::Description<dim, spacedim> 396 &construction_data) 397 { 398 (void)construction_data; 399 400 Assert(false, ExcInternalError()); 401 } 402 403 404 405 template <int dim, int spacedim> 406 void copy_triangulation(const dealii::Triangulation<dim,spacedim> & other_tria)407 Triangulation<dim, spacedim>::copy_triangulation( 408 const dealii::Triangulation<dim, spacedim> &other_tria) 409 { 410 Assert( 411 (dynamic_cast< 412 const dealii::parallel::distributed::Triangulation<dim, spacedim> *>( 413 &other_tria) == nullptr), 414 ExcMessage( 415 "Cannot use this function on parallel::distributed::Triangulation.")); 416 417 dealii::parallel::TriangulationBase<dim, spacedim>::copy_triangulation( 418 other_tria); 419 partition(); 420 this->update_number_cache(); 421 } 422 } // namespace shared 423 } // namespace parallel 424 425 #else 426 427 namespace parallel 428 { 429 namespace shared 430 { 431 template <int dim, int spacedim> 432 bool 433 Triangulation<dim, spacedim>::with_artificial_cells() const 434 { 435 Assert(false, ExcNotImplemented()); 436 return true; 437 } 438 439 template <int dim, int spacedim> 440 bool 441 Triangulation<dim, spacedim>::is_multilevel_hierarchy_constructed() const 442 { 443 return false; 444 } 445 446 template <int dim, int spacedim> 447 const std::vector<unsigned int> & 448 Triangulation<dim, spacedim>::get_true_subdomain_ids_of_cells() const 449 { 450 Assert(false, ExcNotImplemented()); 451 return true_subdomain_ids_of_cells; 452 } 453 454 455 456 template <int dim, int spacedim> 457 const std::vector<unsigned int> & 458 Triangulation<dim, spacedim>::get_true_level_subdomain_ids_of_cells( 459 const unsigned int) const 460 { 461 Assert(false, ExcNotImplemented()); 462 return true_level_subdomain_ids_of_cells; 463 } 464 } // namespace shared 465 } // namespace parallel 466 467 468 #endif 469 470 471 /*-------------- Explicit Instantiations -------------------------------*/ 472 #include "shared_tria.inst" 473 474 DEAL_II_NAMESPACE_CLOSE 475