1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 1998 - 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 17 #include <deal.II/base/geometry_info.h> 18 #include <deal.II/base/memory_consumption.h> 19 #include <deal.II/base/partitioner.h> 20 #include <deal.II/base/thread_management.h> 21 #include <deal.II/base/utilities.h> 22 #include <deal.II/base/work_stream.h> 23 24 #include <deal.II/distributed/shared_tria.h> 25 #include <deal.II/distributed/tria.h> 26 27 #include <deal.II/dofs/dof_accessor.h> 28 #include <deal.II/dofs/dof_handler.h> 29 #include <deal.II/dofs/dof_handler_policy.h> 30 31 #include <deal.II/fe/fe.h> 32 33 #include <deal.II/grid/grid_tools.h> 34 #include <deal.II/grid/tria.h> 35 #include <deal.II/grid/tria_iterator.h> 36 37 #include <algorithm> 38 #include <memory> 39 #include <numeric> 40 #include <set> 41 42 DEAL_II_NAMESPACE_OPEN 43 44 45 namespace internal 46 { 47 namespace DoFHandlerImplementation 48 { 49 namespace Policy 50 { 51 // use class dealii::DoFHandler instead 52 // of namespace internal::DoFHandler in 53 // the following 54 using dealii::DoFHandler; 55 56 namespace 57 { 58 /** 59 * We'll use this constant to mark unnumbered degrees of freedom 60 * during their enumeration in multiple parts of the code. 61 * This constant is necessary to distinguish between valid and 62 * invalid degrees of freedom. 63 */ 64 const types::global_dof_index enumeration_dof_index = 65 numbers::invalid_dof_index - 1; 66 67 /** 68 * Update the cache used for cell dof indices on all (non-artificial) 69 * active cells of the given DoFHandler. 70 */ 71 template <int dim, int spacedim> 72 void update_all_active_cell_dof_indices_caches(const DoFHandler<dim,spacedim> & dof_handler)73 update_all_active_cell_dof_indices_caches( 74 const DoFHandler<dim, spacedim> &dof_handler) 75 { 76 const auto worker = [](const auto &cell, void *, void *) { 77 if (!cell->is_artificial()) 78 cell->update_cell_dof_indices_cache(); 79 }; 80 81 // parallelize filling all of the cell caches. by using 82 // WorkStream, we make sure that we only run through the 83 // range of iterators once, whereas a parallel_for loop 84 // for example has to split the range multiple times, 85 // which is expensive because cell iterators are not 86 // random access iterators with a cheap operator- 87 WorkStream::run(dof_handler.begin_active(), 88 dof_handler.end(), 89 worker, 90 /* copier */ std::function<void(void *)>(), 91 /* scratch_data */ nullptr, 92 /* copy_data */ nullptr, 93 2 * MultithreadInfo::n_threads(), 94 /* chunk_size = */ 32); 95 } 96 97 98 /** 99 * Update the cache used for cell dof indices on all (non-artificial) 100 * level (multigrid) cells of the given DoFHandler. 101 */ 102 template <int dim, int spacedim> 103 void update_all_level_cell_dof_indices_caches(const DoFHandler<dim,spacedim> & dof_handler)104 update_all_level_cell_dof_indices_caches( 105 const DoFHandler<dim, spacedim> &dof_handler) 106 { 107 const auto worker = [](const auto &cell, void *, void *) { 108 if (cell->has_children() || !cell->is_artificial()) 109 cell->update_cell_dof_indices_cache(); 110 }; 111 112 // parallelize filling all of the cell caches. by using 113 // WorkStream, we make sure that we only run through the 114 // range of iterators once, whereas a parallel_for loop 115 // for example has to split the range multiple times, 116 // which is expensive because cell iterators are not 117 // random access iterators with a cheap operator- 118 WorkStream::run(dof_handler.begin(), 119 dof_handler.end(), 120 worker, 121 /* copier */ std::function<void(void *)>(), 122 /* scratch_data */ nullptr, 123 /* copy_data */ nullptr, 124 2 * MultithreadInfo::n_threads(), 125 /* chunk_size = */ 32); 126 } 127 128 129 using DoFIdentities = 130 std::vector<std::pair<unsigned int, unsigned int>>; 131 132 133 /** 134 * Make sure that the given @p identities pointer points to a 135 * valid array. If the pointer is zero beforehand, create an 136 * entry with the correct data. If it is nonzero, don't touch 137 * it. 138 * 139 * @p structdim denotes the dimension of the objects on which 140 * identities are to be represented, i.e. zero for vertices, 141 * one for lines, etc. 142 */ 143 template <int structdim, int dim, int spacedim> 144 const std::unique_ptr<DoFIdentities> & ensure_existence_and_return_dof_identities(const FiniteElement<dim,spacedim> & fe1,const FiniteElement<dim,spacedim> & fe2,std::unique_ptr<DoFIdentities> & identities)145 ensure_existence_and_return_dof_identities( 146 const FiniteElement<dim, spacedim> &fe1, 147 const FiniteElement<dim, spacedim> &fe2, 148 std::unique_ptr<DoFIdentities> & identities) 149 { 150 // see if we need to fill this entry, or whether it already 151 // exists 152 if (identities.get() == nullptr) 153 { 154 switch (structdim) 155 { 156 case 0: 157 { 158 identities = std::make_unique<DoFIdentities>( 159 fe1.hp_vertex_dof_identities(fe2)); 160 break; 161 } 162 163 case 1: 164 { 165 identities = std::make_unique<DoFIdentities>( 166 fe1.hp_line_dof_identities(fe2)); 167 break; 168 } 169 170 case 2: 171 { 172 identities = std::make_unique<DoFIdentities>( 173 fe1.hp_quad_dof_identities(fe2)); 174 break; 175 } 176 177 default: 178 Assert(false, ExcNotImplemented()); 179 } 180 181 // double check whether the newly created entries make 182 // any sense at all 183 for (unsigned int i = 0; i < identities->size(); ++i) 184 { 185 Assert((*identities)[i].first < 186 fe1.template n_dofs_per_object<structdim>(), 187 ExcInternalError()); 188 Assert((*identities)[i].second < 189 fe2.template n_dofs_per_object<structdim>(), 190 ExcInternalError()); 191 } 192 } 193 194 return identities; 195 } 196 } // namespace 197 198 199 200 struct Implementation 201 { 202 /* -------------- distribute_dofs functionality ------------- */ 203 204 /** 205 * Compute identities between DoFs located on vertices. Called from 206 * distribute_dofs(). 207 */ 208 template <int dim, int spacedim> 209 static std::map<types::global_dof_index, types::global_dof_index> compute_vertex_dof_identitiesinternal::DoFHandlerImplementation::Policy::Implementation210 compute_vertex_dof_identities( 211 const DoFHandler<dim, spacedim> &dof_handler) 212 { 213 Assert( 214 dof_handler.hp_capability_enabled == true, 215 (typename DoFHandler<dim, spacedim>::ExcNotAvailableWithoutHP())); 216 217 std::map<types::global_dof_index, types::global_dof_index> 218 dof_identities; 219 220 // Note: we may wish to have something here similar to what 221 // we do for lines and quads, namely that we only identify 222 // dofs for any fe towards the most dominating one. however, 223 // it is not clear whether this is actually necessary for 224 // vertices at all, I can't think of a finite element that 225 // would make that necessary... 226 dealii::Table<2, std::unique_ptr<DoFIdentities>> 227 vertex_dof_identities(dof_handler.get_fe_collection().size(), 228 dof_handler.get_fe_collection().size()); 229 230 // loop over all vertices and see which one we need to work on 231 for (unsigned int vertex_index = 0; 232 vertex_index < dof_handler.get_triangulation().n_vertices(); 233 ++vertex_index) 234 if (dof_handler.get_triangulation() 235 .get_used_vertices()[vertex_index] == true) 236 { 237 const unsigned int n_active_fe_indices = 238 dealii::internal::DoFAccessorImplementation::Implementation:: 239 n_active_fe_indices(dof_handler, 240 0, 241 vertex_index, 242 std::integral_constant<int, 0>()); 243 244 if (n_active_fe_indices > 1) 245 { 246 const std::set<unsigned int> fe_indices = 247 dealii::internal::DoFAccessorImplementation:: 248 Implementation::get_active_fe_indices( 249 dof_handler, 250 0, 251 vertex_index, 252 std::integral_constant<int, 0>()); 253 254 // find out which is the most dominating finite 255 // element of the ones that are used on this vertex 256 unsigned int most_dominating_fe_index = 257 dof_handler.get_fe_collection().find_dominating_fe( 258 fe_indices, 259 /*codim*/ dim); 260 261 // if we haven't found a dominating finite element, 262 // choose the very first one to be dominant 263 if (most_dominating_fe_index == 264 numbers::invalid_unsigned_int) 265 most_dominating_fe_index = 266 dealii::internal::DoFAccessorImplementation:: 267 Implementation::nth_active_fe_index( 268 dof_handler, 269 0, 270 vertex_index, 271 0, 272 std::integral_constant<int, 0>()); 273 274 // loop over the indices of all the finite 275 // elements that are not dominating, and 276 // identify their dofs to the most dominating 277 // one 278 for (const auto &other_fe_index : fe_indices) 279 if (other_fe_index != most_dominating_fe_index) 280 { 281 // make sure the entry in the equivalence 282 // table exists 283 const auto &identities = 284 *ensure_existence_and_return_dof_identities<0>( 285 dof_handler.get_fe(most_dominating_fe_index), 286 dof_handler.get_fe(other_fe_index), 287 vertex_dof_identities[most_dominating_fe_index] 288 [other_fe_index]); 289 290 // then loop through the identities we 291 // have. first get the global numbers of the 292 // dofs we want to identify and make sure they 293 // are not yet constrained to anything else, 294 // except for to each other. use the rule that 295 // we will always constrain the dof with the 296 // higher fe index to the one with the lower, 297 // to avoid circular reasoning. 298 for (const auto &identity : identities) 299 { 300 const types::global_dof_index primary_dof_index = 301 dealii::internal::DoFAccessorImplementation:: 302 Implementation::get_dof_index( 303 dof_handler, 304 0, 305 vertex_index, 306 most_dominating_fe_index, 307 identity.first, 308 std::integral_constant<int, 0>()); 309 const types::global_dof_index 310 dependent_dof_index = 311 dealii::internal::DoFAccessorImplementation:: 312 Implementation::get_dof_index( 313 dof_handler, 314 0, 315 vertex_index, 316 other_fe_index, 317 identity.second, 318 std::integral_constant<int, 0>()); 319 320 // on subdomain boundaries, we will 321 // encounter invalid DoFs on ghost cells, 322 // for which we have not yet distributed 323 // valid indices. depending on which finte 324 // element is dominating the other on this 325 // interface, we either have to constrain 326 // the valid to the invalid indices, or vice 327 // versa. 328 // 329 // we only store an identity if we are about 330 // to overwrite a valid DoF. we will skip 331 // constraining invalid DoFs for now, and 332 // consider them later in Phase 5. 333 if (dependent_dof_index != 334 numbers::invalid_dof_index) 335 { 336 // if the DoF indices of both elements 337 // are already distributed, i.e., both 338 // of these 'fe_indices' are associated 339 // with a locally owned cell, then we 340 // should either not have a dof_identity 341 // yet, or it must come out here to be 342 // exactly as we had computed before 343 if (primary_dof_index != 344 numbers::invalid_dof_index) 345 Assert( 346 (dof_identities.find(primary_dof_index) == 347 dof_identities.end()) || 348 (dof_identities[dependent_dof_index] == 349 primary_dof_index), 350 ExcInternalError()); 351 352 dof_identities[dependent_dof_index] = 353 primary_dof_index; 354 } 355 } 356 } 357 } 358 } 359 360 return dof_identities; 361 } 362 363 364 /** 365 * Compute identities between DoFs located on lines. Called from 366 * distribute_dofs(). 367 */ 368 template <int spacedim> 369 static std::map<types::global_dof_index, types::global_dof_index> compute_line_dof_identitiesinternal::DoFHandlerImplementation::Policy::Implementation370 compute_line_dof_identities(const DoFHandler<1, spacedim> &dof_handler) 371 { 372 (void)dof_handler; 373 Assert( 374 dof_handler.hp_capability_enabled == true, 375 (typename DoFHandler<1, spacedim>::ExcNotAvailableWithoutHP())); 376 377 return std::map<types::global_dof_index, types::global_dof_index>(); 378 } 379 380 381 template <int dim, int spacedim> 382 static std::map<types::global_dof_index, types::global_dof_index> compute_line_dof_identitiesinternal::DoFHandlerImplementation::Policy::Implementation383 compute_line_dof_identities( 384 const DoFHandler<dim, spacedim> &dof_handler) 385 { 386 Assert( 387 dof_handler.hp_capability_enabled == true, 388 (typename DoFHandler<dim, spacedim>::ExcNotAvailableWithoutHP())); 389 390 std::map<types::global_dof_index, types::global_dof_index> 391 dof_identities; 392 393 // we will mark lines that we have already treated, so first save and 394 // clear the user flags on lines and later restore them 395 std::vector<bool> user_flags; 396 dof_handler.get_triangulation().save_user_flags_line(user_flags); 397 const_cast<dealii::Triangulation<dim, spacedim> &>( 398 dof_handler.get_triangulation()) 399 .clear_user_flags_line(); 400 401 // An implementation of the algorithm described in the hp paper, 402 // including the modification mentioned later in the "complications in 403 // 3-d" subsections 404 // 405 // as explained there, we do something only if there are exactly 2 406 // finite elements associated with an object. if there is only one, 407 // then there is nothing to do anyway, and if there are 3 or more, 408 // then we can get into trouble. note that this only happens for lines 409 // in 3d and higher, and for quads only in 4d and higher, so this 410 // isn't a particularly frequent case 411 // 412 // there is one case, however, that we would like to handle (see, for 413 // example, the hp/crash_15 testcase): if we have 414 // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i', 415 // then we should be able to handle this because we can simply unify 416 // *all* dofs, not only a some. so what we do is to first treat all 417 // pairs of finite elements that have *identical* dofs, and then only 418 // deal with those that are not identical of which we can handle at 419 // most 2 420 dealii::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities( 421 dof_handler.fe_collection.size(), dof_handler.fe_collection.size()); 422 423 for (const auto &cell : dof_handler.active_cell_iterators()) 424 for (const auto l : cell->line_indices()) 425 if (cell->line(l)->user_flag_set() == false) 426 { 427 const auto line = cell->line(l); 428 line->set_user_flag(); 429 430 unsigned int unique_sets_of_dofs = 431 line->n_active_fe_indices(); 432 433 // do a first loop over all sets of dofs and do identity 434 // uniquification 435 const unsigned int n_active_fe_indices = 436 line->n_active_fe_indices(); 437 for (unsigned int f = 0; f < n_active_fe_indices; ++f) 438 for (unsigned int g = f + 1; g < n_active_fe_indices; ++g) 439 { 440 const unsigned int fe_index_1 = 441 line->nth_active_fe_index(f), 442 fe_index_2 = 443 line->nth_active_fe_index(g); 444 445 // as described in the hp paper, we only unify on lines 446 // when there are at most two different FE objects 447 // assigned on it. 448 // however, more than two 'active_fe_indices' can be 449 // attached that still fulfill the above criterion, 450 // i.e. when two different FiniteElement objects are 451 // assigned to neighboring cells that map their degrees 452 // of freedom one-to-one. 453 // we cannot verify with certainty if two dofs each of 454 // separate FiniteElement objects actually map 455 // one-to-one. however, checking for the number of 456 // 'dofs_per_line' turned out to be a reasonable 457 // approach, that also works for e.g. two different 458 // FE_Q objects of the same order, from which one is 459 // enhanced by a bubble function that is zero on the 460 // boundary. 461 if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() == 462 dof_handler.get_fe(fe_index_2) 463 .n_dofs_per_line()) && 464 (dof_handler.get_fe(fe_index_1).n_dofs_per_line() > 465 0)) 466 { 467 // the number of dofs per line is identical 468 const unsigned int dofs_per_line = 469 dof_handler.get_fe(fe_index_1).n_dofs_per_line(); 470 471 const auto &identities = 472 *ensure_existence_and_return_dof_identities<1>( 473 dof_handler.get_fe(fe_index_1), 474 dof_handler.get_fe(fe_index_2), 475 line_dof_identities[fe_index_1][fe_index_2]); 476 // see if these sets of dofs are identical. the 477 // first condition for this is that indeed there are 478 // n identities 479 if (identities.size() == dofs_per_line) 480 { 481 unsigned int i = 0; 482 for (; i < dofs_per_line; ++i) 483 if ((identities[i].first != i) && 484 (identities[i].second != i)) 485 // not an identity 486 break; 487 488 if (i == dofs_per_line) 489 { 490 // The line dofs (i.e., the ones interior to 491 // a line) of these two finite elements are 492 // identical. Note that there could be 493 // situations when one element still 494 // dominates another, e.g.: FE_Q(2) x 495 // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1) 496 497 --unique_sets_of_dofs; 498 499 // determine which one of both finite 500 // elements is the dominating one. 501 const std::set<unsigned int> fe_indices{ 502 fe_index_1, fe_index_2}; 503 504 unsigned int dominating_fe_index = 505 dof_handler.get_fe_collection() 506 .find_dominating_fe(fe_indices, 507 /*codim=*/dim - 1); 508 unsigned int other_fe_index = 509 numbers::invalid_unsigned_int; 510 511 if (dominating_fe_index != 512 numbers::invalid_unsigned_int) 513 other_fe_index = 514 (dominating_fe_index == fe_index_1) ? 515 fe_index_2 : 516 fe_index_1; 517 else 518 { 519 // if we haven't found a dominating 520 // finite element, choose the one with 521 // the lower index to be dominating 522 dominating_fe_index = fe_index_1; 523 other_fe_index = fe_index_2; 524 } 525 526 for (unsigned int j = 0; j < dofs_per_line; 527 ++j) 528 { 529 const types::global_dof_index 530 primary_dof_index = line->dof_index( 531 j, dominating_fe_index); 532 const types::global_dof_index 533 dependent_dof_index = 534 line->dof_index(j, other_fe_index); 535 536 // on subdomain boundaries, we will 537 // encounter invalid DoFs on ghost 538 // cells, for which we have not yet 539 // distributed valid indices. depending 540 // on which finte element is dominating 541 // the other on this interface, we 542 // either have to constrain the valid to 543 // the invalid indices, or vice versa. 544 // 545 // we only store an identity if we are 546 // about to overwrite a valid DoF. we 547 // will skip constraining invalid DoFs 548 // for now, and consider them later in 549 // Phase 5. 550 if (dependent_dof_index != 551 numbers::invalid_dof_index) 552 { 553 if (primary_dof_index != 554 numbers::invalid_dof_index) 555 { 556 // if primary dof was already 557 // constrained, constrain to 558 // that one, otherwise constrain 559 // dependent to primary 560 if (dof_identities.find( 561 primary_dof_index) != 562 dof_identities.end()) 563 { 564 // if the DoF indices of 565 // both elements are already 566 // distributed, i.e., both 567 // of these 'fe_indices' are 568 // associated with a locally 569 // owned cell, then we 570 // should either not have a 571 // dof_identity yet, or it 572 // must come out here to be 573 // exactly as we had 574 // computed before 575 Assert( 576 dof_identities.find( 577 dof_identities 578 [primary_dof_index]) == 579 dof_identities.end(), 580 ExcInternalError()); 581 582 dof_identities 583 [dependent_dof_index] = 584 dof_identities 585 [primary_dof_index]; 586 } 587 else 588 { 589 // see comment above for an 590 // explanation of this 591 // assertion 592 Assert( 593 (dof_identities.find( 594 primary_dof_index) == 595 dof_identities.end()) || 596 (dof_identities 597 [dependent_dof_index] == 598 primary_dof_index), 599 ExcInternalError()); 600 601 dof_identities 602 [dependent_dof_index] = 603 primary_dof_index; 604 } 605 } 606 else 607 { 608 // set dependent_dof to 609 // primary_dof_index, which is 610 // invalid 611 dof_identities 612 [dependent_dof_index] = 613 numbers::invalid_dof_index; 614 } 615 } 616 } 617 } 618 } 619 } 620 } 621 622 // if at this point, there is only one unique set of dofs 623 // left, then we have taken care of everything above. if there 624 // are two, then we need to deal with them here. if there are 625 // more, then we punt, as described in the paper (and 626 // mentioned above) 627 // TODO: The check for 'dim==2' was inserted by intuition. It 628 // fixes 629 // the previous problems with step-27 in 3D. But an 630 // explanation for this is still required, and what we do here 631 // is not what we describe in the paper!. 632 if ((unique_sets_of_dofs == 2) && (dim == 2)) 633 { 634 const std::set<unsigned int> fe_indices = 635 line->get_active_fe_indices(); 636 637 // find out which is the most dominating finite element of 638 // the ones that are used on this line 639 const unsigned int most_dominating_fe_index = 640 dof_handler.get_fe_collection().find_dominating_fe( 641 fe_indices, 642 /*codim=*/dim - 1); 643 644 // if we found the most dominating element, then use this 645 // to eliminate some of the degrees of freedom by 646 // identification. otherwise, the code that computes 647 // hanging node constraints will have to deal with it by 648 // computing appropriate constraints along this face/edge 649 if (most_dominating_fe_index != 650 numbers::invalid_unsigned_int) 651 { 652 // loop over the indices of all the finite elements 653 // that are not dominating, and identify their dofs to 654 // the most dominating one 655 for (const auto &other_fe_index : fe_indices) 656 if (other_fe_index != most_dominating_fe_index) 657 { 658 const auto &identities = 659 *ensure_existence_and_return_dof_identities< 660 1>(dof_handler.get_fe( 661 most_dominating_fe_index), 662 dof_handler.get_fe(other_fe_index), 663 line_dof_identities 664 [most_dominating_fe_index] 665 [other_fe_index]); 666 667 for (const auto &identity : identities) 668 { 669 const types::global_dof_index 670 primary_dof_index = line->dof_index( 671 identity.first, 672 most_dominating_fe_index); 673 const types::global_dof_index 674 dependent_dof_index = 675 line->dof_index(identity.second, 676 other_fe_index); 677 678 // on subdomain boundaries, we will 679 // encounter invalid DoFs on ghost cells, 680 // for which we have not yet distributed 681 // valid indices. depending on which finte 682 // element is dominating the other on this 683 // interface, we either have to constrain 684 // the valid to the invalid indices, or vice 685 // versa. 686 // 687 // we only store an identity if we are about 688 // to overwrite a valid DoF. we will skip 689 // constraining invalid DoFs for now, and 690 // consider them later in Phase 5. 691 if (dependent_dof_index != 692 numbers::invalid_dof_index) 693 { 694 // if the DoF indices of both elements 695 // are already distributed, i.e., both 696 // of these 'fe_indices' are associated 697 // with a locally owned cell, then we 698 // should either not have a dof_identity 699 // yet, or it must come out here to be 700 // exactly as we had computed before 701 if (primary_dof_index != 702 numbers::invalid_dof_index) 703 Assert((dof_identities.find( 704 primary_dof_index) == 705 dof_identities.end()) || 706 (dof_identities 707 [dependent_dof_index] == 708 primary_dof_index), 709 ExcInternalError()); 710 711 dof_identities[dependent_dof_index] = 712 primary_dof_index; 713 } 714 } 715 } 716 } 717 } 718 } 719 720 // finally restore the user flags 721 const_cast<dealii::Triangulation<dim, spacedim> &>( 722 dof_handler.get_triangulation()) 723 .load_user_flags_line(user_flags); 724 725 return dof_identities; 726 } 727 728 729 730 /** 731 * Compute identities between DoFs located on quads. Called from 732 * distribute_dofs(). 733 */ 734 template <int dim, int spacedim> 735 static std::map<types::global_dof_index, types::global_dof_index> compute_quad_dof_identitiesinternal::DoFHandlerImplementation::Policy::Implementation736 compute_quad_dof_identities( 737 const DoFHandler<dim, spacedim> &dof_handler) 738 { 739 (void)dof_handler; 740 Assert( 741 dof_handler.hp_capability_enabled == true, 742 (typename DoFHandler<dim, spacedim>::ExcNotAvailableWithoutHP())); 743 744 // this function should only be called for dim<3 where there are 745 // no quad dof identies. for dim==3, the specialization below should 746 // take care of it 747 Assert(dim < 3, ExcInternalError()); 748 749 return std::map<types::global_dof_index, types::global_dof_index>(); 750 } 751 752 753 template <int spacedim> 754 static std::map<types::global_dof_index, types::global_dof_index> compute_quad_dof_identitiesinternal::DoFHandlerImplementation::Policy::Implementation755 compute_quad_dof_identities(const DoFHandler<3, spacedim> &dof_handler) 756 { 757 Assert( 758 dof_handler.hp_capability_enabled == true, 759 (typename DoFHandler<3, spacedim>::ExcNotAvailableWithoutHP())); 760 761 const int dim = 3; 762 763 std::map<types::global_dof_index, types::global_dof_index> 764 dof_identities; 765 766 767 // we will mark quads that we have already treated, so first 768 // save and clear the user flags on quads and later restore 769 // them 770 std::vector<bool> user_flags; 771 dof_handler.get_triangulation().save_user_flags_quad(user_flags); 772 const_cast<dealii::Triangulation<dim, spacedim> &>( 773 dof_handler.get_triangulation()) 774 .clear_user_flags_quad(); 775 776 // An implementation of the algorithm described in the hp 777 // paper, including the modification mentioned later in the 778 // "complications in 3-d" subsections 779 // 780 // as explained there, we do something only if there are 781 // exactly 2 finite elements associated with an object. if 782 // there is only one, then there is nothing to do anyway, 783 // and if there are 3 or more, then we can get into 784 // trouble. note that this only happens for lines in 3d and 785 // higher, and for quads only in 4d and higher, so this 786 // isn't a particularly frequent case 787 dealii::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities( 788 dof_handler.fe_collection.size(), 789 dof_handler.fe_collection.size(), 790 2 /*triangle (0) or quadrilateral (1)*/); 791 792 for (const auto &cell : dof_handler.active_cell_iterators()) 793 for (const auto q : cell->face_indices()) 794 if ((cell->quad(q)->user_flag_set() == false) && 795 (cell->quad(q)->n_active_fe_indices() == 2)) 796 { 797 const auto quad = cell->quad(q); 798 quad->set_user_flag(); 799 800 const std::set<unsigned int> fe_indices = 801 quad->get_active_fe_indices(); 802 803 // find out which is the most dominating finite 804 // element of the ones that are used on this quad 805 const unsigned int most_dominating_fe_index = 806 dof_handler.get_fe_collection().find_dominating_fe( 807 fe_indices, 808 /*codim=*/dim - 2); 809 810 // if we found the most dominating element, then use 811 // this to eliminate some of the degrees of freedom 812 // by identification. otherwise, the code that 813 // computes hanging node constraints will have to 814 // deal with it by computing appropriate constraints 815 // along this face/edge 816 if (most_dominating_fe_index != numbers::invalid_unsigned_int) 817 { 818 // loop over the indices of all the finite 819 // elements that are not dominating, and 820 // identify their dofs to the most dominating 821 // one 822 for (const auto &other_fe_index : fe_indices) 823 if (other_fe_index != most_dominating_fe_index) 824 { 825 const auto &identities = 826 *ensure_existence_and_return_dof_identities<2>( 827 dof_handler.get_fe(most_dominating_fe_index), 828 dof_handler.get_fe(other_fe_index), 829 quad_dof_identities 830 [most_dominating_fe_index][other_fe_index] 831 [cell->quad(q)->reference_cell_type() == 832 ReferenceCell::Type::Quad]); 833 834 for (const auto &identity : identities) 835 { 836 const types::global_dof_index 837 primary_dof_index = 838 quad->dof_index(identity.first, 839 most_dominating_fe_index); 840 const types::global_dof_index 841 dependent_dof_index = 842 quad->dof_index(identity.second, 843 other_fe_index); 844 845 // we only store an identity if we are about to 846 // overwrite a valid degree of freedom. we will 847 // skip invalid degrees of freedom (that are 848 // associated with ghost cells) for now, and 849 // consider them later in phase 5. 850 if (dependent_dof_index != 851 numbers::invalid_dof_index) 852 { 853 // if the DoF indices of both elements are 854 // already distributed, i.e., both of these 855 // 'fe_indices' are associated with a 856 // locally owned cell, then we should either 857 // not have a dof_identity yet, or it must 858 // come out here to be exactly as we had 859 // computed before 860 if (primary_dof_index != 861 numbers::invalid_dof_index) 862 Assert((dof_identities.find( 863 primary_dof_index) == 864 dof_identities.end()) || 865 (dof_identities 866 [dependent_dof_index] == 867 primary_dof_index), 868 ExcInternalError()); 869 870 dof_identities[dependent_dof_index] = 871 primary_dof_index; 872 } 873 } 874 } 875 } 876 } 877 878 // finally restore the user flags 879 const_cast<dealii::Triangulation<dim, spacedim> &>( 880 dof_handler.get_triangulation()) 881 .load_user_flags_quad(user_flags); 882 883 return dof_identities; 884 } 885 886 887 888 /** 889 * Compute the constraints that correspond to unifying DoF indices 890 * on vertices, lines, and quads. Do so in parallel. 891 */ 892 template <int dim, int spacedim> 893 static void compute_dof_identitiesinternal::DoFHandlerImplementation::Policy::Implementation894 compute_dof_identities(std::vector<std::map<types::global_dof_index, 895 types::global_dof_index>> 896 &all_constrained_indices, 897 const DoFHandler<dim, spacedim> &dof_handler) 898 { 899 if (dof_handler.hp_capability_enabled == false) 900 return; 901 902 Assert(all_constrained_indices.size() == dim, ExcInternalError()); 903 904 Threads::TaskGroup<> tasks; 905 906 unsigned int i = 0; 907 tasks += Threads::new_task([&, i]() { 908 all_constrained_indices[i] = 909 compute_vertex_dof_identities(dof_handler); 910 }); 911 912 if (dim > 1) 913 { 914 ++i; 915 tasks += Threads::new_task([&, i]() { 916 all_constrained_indices[i] = 917 compute_line_dof_identities(dof_handler); 918 }); 919 } 920 921 if (dim > 2) 922 { 923 ++i; 924 tasks += Threads::new_task([&, i]() { 925 all_constrained_indices[i] = 926 compute_quad_dof_identities(dof_handler); 927 }); 928 } 929 930 tasks.join_all(); 931 } 932 933 934 935 /** 936 * Once degrees of freedom have been distributed on all cells, we may 937 * want to eliminate duplicates, and enumerate the remaining ones 938 * consecutively. This particular function is responsible for the 939 * latter part. 940 * 941 * This function stores the new indices of all DoFs in ascending order 942 * in @p new_dof_indices, which can be used to renumber all DoFs with 943 * the renumber_dofs() function later. 944 * 945 * This vector will contain enumerated indices, skipping invalid indices 946 * previously stored in it. Additionally, if a 947 * @p all_constrained_indices parameter is provided, DoF identity 948 * relations will be considered as well during the enumeration process 949 * by identifying similar DoFs on vertices, lines and quads. 950 * 951 * Returns the final number of degrees of freedom, which is the number 952 * of all valid DoF indices in @p new_dof_indices. 953 */ 954 template <int dim, int spacedim> 955 static types::global_dof_index enumerate_dof_indices_for_renumberinginternal::DoFHandlerImplementation::Policy::Implementation956 enumerate_dof_indices_for_renumbering( 957 std::vector<types::global_dof_index> &new_dof_indices, 958 const std::vector< 959 std::map<types::global_dof_index, types::global_dof_index>> 960 &all_constrained_indices, 961 const DoFHandler<dim, spacedim> &) 962 { 963 Assert(all_constrained_indices.size() == dim, ExcInternalError()); 964 965 // first preset the new DoF indices that are identities 966 for (const auto &constrained_dof_indices : all_constrained_indices) 967 for (const auto &p : constrained_dof_indices) 968 if (new_dof_indices[p.first] != numbers::invalid_dof_index) 969 { 970 Assert(new_dof_indices[p.first] == enumeration_dof_index, 971 ExcInternalError()); 972 973 new_dof_indices[p.first] = p.second; 974 } 975 976 // then enumerate the rest 977 types::global_dof_index next_free_dof = 0; 978 for (auto &new_dof_index : new_dof_indices) 979 if (new_dof_index == enumeration_dof_index) 980 new_dof_index = next_free_dof++; 981 982 // then loop over all those that are constrained and record the 983 // new dof number for those 984 for (const auto &constrained_dof_indices : all_constrained_indices) 985 for (const auto &p : constrained_dof_indices) 986 if (new_dof_indices[p.first] != numbers::invalid_dof_index) 987 { 988 Assert(new_dof_indices[p.first] != enumeration_dof_index, 989 ExcInternalError()); 990 991 if (p.second != numbers::invalid_dof_index) 992 new_dof_indices[p.first] = new_dof_indices[p.second]; 993 } 994 995 for (const types::global_dof_index new_dof_index : new_dof_indices) 996 { 997 (void)new_dof_index; 998 Assert(new_dof_index != enumeration_dof_index, 999 ExcInternalError()); 1000 Assert(new_dof_index < next_free_dof || 1001 new_dof_index == numbers::invalid_dof_index, 1002 ExcInternalError()); 1003 } 1004 1005 return next_free_dof; 1006 } 1007 1008 1009 1010 /** 1011 * Once degrees of freedom have been distributed on all cells, see if 1012 * we can identify DoFs on neighboring cells. This function does 1013 * nothing unless the DoFHandler has hp capabilities. 1014 * 1015 * Return the final number of degrees of freedom, which is the old one 1016 * minus however many were identified. 1017 */ 1018 template <int dim, int spacedim> 1019 static types::global_dof_index unify_dof_indicesinternal::DoFHandlerImplementation::Policy::Implementation1020 unify_dof_indices(const DoFHandler<dim, spacedim> &dof_handler, 1021 const unsigned int n_dofs_before_identification, 1022 const bool check_validity) 1023 { 1024 if (dof_handler.hp_capability_enabled == false) 1025 return n_dofs_before_identification; 1026 1027 std::vector< 1028 std::map<types::global_dof_index, types::global_dof_index>> 1029 all_constrained_indices(dim); 1030 compute_dof_identities(all_constrained_indices, dof_handler); 1031 1032 std::vector<dealii::types::global_dof_index> renumbering( 1033 n_dofs_before_identification, enumeration_dof_index); 1034 const types::global_dof_index n_dofs = 1035 enumerate_dof_indices_for_renumbering(renumbering, 1036 all_constrained_indices, 1037 dof_handler); 1038 1039 renumber_dofs(renumbering, IndexSet(0), dof_handler, check_validity); 1040 1041 update_all_active_cell_dof_indices_caches(dof_handler); 1042 1043 return n_dofs; 1044 } 1045 1046 1047 1048 /** 1049 * Merge invalid DoF indices on vertices located on ghost interfaces 1050 * by a dominating valid one. 1051 */ 1052 template <int dim, int spacedim> 1053 static void merge_invalid_vertex_dofs_on_ghost_interfacesinternal::DoFHandlerImplementation::Policy::Implementation1054 merge_invalid_vertex_dofs_on_ghost_interfaces( 1055 DoFHandler<dim, spacedim> &dof_handler) 1056 { 1057 Assert( 1058 dof_handler.hp_capability_enabled == true, 1059 (typename DoFHandler<dim, spacedim>::ExcNotAvailableWithoutHP())); 1060 1061 // Note: we may wish to have something here similar to what 1062 // we do for lines and quads, namely that we only identify 1063 // dofs for any fe towards the most dominating one. however, 1064 // it is not clear whether this is actually necessary for 1065 // vertices at all, I can't think of a finite element that 1066 // would make that necessary... 1067 dealii::Table<2, std::unique_ptr<DoFIdentities>> 1068 vertex_dof_identities(dof_handler.get_fe_collection().size(), 1069 dof_handler.get_fe_collection().size()); 1070 1071 // mark all vertices on ghost cells 1072 std::vector<bool> include_vertex( 1073 dof_handler.get_triangulation().n_vertices(), false); 1074 if (dynamic_cast<const dealii::parallel:: 1075 DistributedTriangulationBase<dim, spacedim> *>( 1076 &dof_handler.get_triangulation()) != nullptr) 1077 for (const auto &cell : dof_handler.active_cell_iterators()) 1078 if (cell->is_ghost()) 1079 for (const unsigned int v : cell->vertex_indices()) 1080 include_vertex[cell->vertex_index(v)] = true; 1081 1082 // loop over all vertices and see which one we need to work on 1083 for (unsigned int vertex_index = 0; 1084 vertex_index < dof_handler.get_triangulation().n_vertices(); 1085 ++vertex_index) 1086 if ((dof_handler.get_triangulation() 1087 .get_used_vertices()[vertex_index] == true) && 1088 (include_vertex[vertex_index] == true)) 1089 { 1090 const unsigned int n_active_fe_indices = 1091 dealii::internal::DoFAccessorImplementation::Implementation:: 1092 n_active_fe_indices(dof_handler, 1093 0, 1094 vertex_index, 1095 std::integral_constant<int, 0>()); 1096 1097 if (n_active_fe_indices > 1) 1098 { 1099 const std::set<unsigned int> fe_indices = 1100 dealii::internal::DoFAccessorImplementation:: 1101 Implementation::get_active_fe_indices( 1102 dof_handler, 1103 0, 1104 vertex_index, 1105 std::integral_constant<int, 0>()); 1106 1107 // find out which is the most dominating finite 1108 // element of the ones that are used on this vertex 1109 unsigned int most_dominating_fe_index = 1110 dof_handler.get_fe_collection().find_dominating_fe( 1111 fe_indices, 1112 /*codim=*/dim); 1113 1114 // if we haven't found a dominating finite element, 1115 // choose the very first one to be dominant similar 1116 // to compute_vertex_dof_identities() 1117 if (most_dominating_fe_index == 1118 numbers::invalid_unsigned_int) 1119 most_dominating_fe_index = 1120 dealii::internal::DoFAccessorImplementation:: 1121 Implementation::nth_active_fe_index( 1122 dof_handler, 1123 0, 1124 vertex_index, 1125 0, 1126 std::integral_constant<int, 0>()); 1127 1128 // loop over the indices of all the finite 1129 // elements that are not dominating, and 1130 // identify their dofs to the most dominating 1131 // one 1132 for (const auto &other_fe_index : fe_indices) 1133 if (other_fe_index != most_dominating_fe_index) 1134 { 1135 // make sure the entry in the equivalence 1136 // table exists 1137 const auto &identities = 1138 *ensure_existence_and_return_dof_identities<0>( 1139 dof_handler.get_fe(most_dominating_fe_index), 1140 dof_handler.get_fe(other_fe_index), 1141 vertex_dof_identities[most_dominating_fe_index] 1142 [other_fe_index]); 1143 1144 // then loop through the identities we 1145 // have. first get the global numbers of the 1146 // dofs we want to identify and make sure they 1147 // are not yet constrained to anything else, 1148 // except for to each other. use the rule that 1149 // we will always constrain the dof with the 1150 // higher fe index to the one with the lower, 1151 // to avoid circular reasoning. 1152 for (const auto &identity : identities) 1153 { 1154 const types::global_dof_index primary_dof_index = 1155 dealii::internal::DoFAccessorImplementation:: 1156 Implementation::get_dof_index( 1157 dof_handler, 1158 0, 1159 vertex_index, 1160 most_dominating_fe_index, 1161 identity.first, 1162 std::integral_constant<int, 0>()); 1163 const types::global_dof_index 1164 dependent_dof_index = 1165 dealii::internal::DoFAccessorImplementation:: 1166 Implementation::get_dof_index( 1167 dof_handler, 1168 0, 1169 vertex_index, 1170 other_fe_index, 1171 identity.second, 1172 std::integral_constant<int, 0>()); 1173 1174 // check if we are on an interface between 1175 // a locally owned and a ghost cell on which 1176 // we need to work on. 1177 // 1178 // all degrees of freedom belonging to 1179 // dominating fe indices or to a processor 1180 // with a higher rank have been set at this 1181 // point (either in Phase 2, or after the 1182 // first ghost exchange in Phase 5). thus, 1183 // we only have to set the indices of 1184 // degrees of freedom that have been 1185 // previously flagged invalid. 1186 if ((dependent_dof_index == 1187 numbers::invalid_dof_index) && 1188 (primary_dof_index != 1189 numbers::invalid_dof_index)) 1190 dealii::internal::DoFAccessorImplementation:: 1191 Implementation::set_dof_index( 1192 dof_handler, 1193 0, 1194 vertex_index, 1195 other_fe_index, 1196 identity.second, 1197 std::integral_constant<int, 0>(), 1198 primary_dof_index); 1199 } 1200 } 1201 } 1202 } 1203 } 1204 1205 1206 1207 /** 1208 * Merge invalid DoF indices on lines located on ghost interfaces 1209 * by a dominating valid one. 1210 */ 1211 template <int spacedim> merge_invalid_line_dofs_on_ghost_interfacesinternal::DoFHandlerImplementation::Policy::Implementation1212 static void merge_invalid_line_dofs_on_ghost_interfaces( 1213 DoFHandler<1, spacedim> &dof_handler) 1214 { 1215 (void)dof_handler; 1216 Assert( 1217 dof_handler.hp_capability_enabled == true, 1218 (typename DoFHandler<1, spacedim>::ExcNotAvailableWithoutHP())); 1219 } 1220 1221 1222 template <int dim, int spacedim> 1223 static void merge_invalid_line_dofs_on_ghost_interfacesinternal::DoFHandlerImplementation::Policy::Implementation1224 merge_invalid_line_dofs_on_ghost_interfaces( 1225 DoFHandler<dim, spacedim> &dof_handler) 1226 { 1227 Assert( 1228 dof_handler.hp_capability_enabled == true, 1229 (typename DoFHandler<dim, spacedim>::ExcNotAvailableWithoutHP())); 1230 1231 // we will mark lines that we have already treated, so first save and 1232 // clear the user flags on lines and later restore them 1233 std::vector<bool> user_flags; 1234 dof_handler.get_triangulation().save_user_flags_line(user_flags); 1235 const_cast<dealii::Triangulation<dim, spacedim> &>( 1236 dof_handler.get_triangulation()) 1237 .clear_user_flags_line(); 1238 1239 // mark all lines on ghost cells 1240 for (const auto &cell : dof_handler.active_cell_iterators()) 1241 if (cell->is_ghost()) 1242 for (const auto l : cell->line_indices()) 1243 cell->line(l)->set_user_flag(); 1244 1245 // An implementation of the algorithm described in the hp paper, 1246 // including the modification mentioned later in the "complications in 1247 // 3-d" subsections 1248 // 1249 // as explained there, we do something only if there are exactly 2 1250 // finite elements associated with an object. if there is only one, 1251 // then there is nothing to do anyway, and if there are 3 or more, 1252 // then we can get into trouble. note that this only happens for lines 1253 // in 3d and higher, and for quads only in 4d and higher, so this 1254 // isn't a particularly frequent case 1255 // 1256 // there is one case, however, that we would like to handle (see, for 1257 // example, the hp/crash_15 testcase): if we have 1258 // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i', 1259 // then we should be able to handle this because we can simply unify 1260 // *all* dofs, not only a some. so what we do is to first treat all 1261 // pairs of finite elements that have *identical* dofs, and then only 1262 // deal with those that are not identical of which we can handle at 1263 // most 2 1264 dealii::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities( 1265 dof_handler.fe_collection.size(), dof_handler.fe_collection.size()); 1266 1267 for (const auto &cell : dof_handler.active_cell_iterators()) 1268 for (const auto l : cell->line_indices()) 1269 if ((cell->is_locally_owned()) && 1270 (cell->line(l)->user_flag_set() == true)) 1271 { 1272 const auto line = cell->line(l); 1273 line->clear_user_flag(); 1274 1275 unsigned int unique_sets_of_dofs = 1276 line->n_active_fe_indices(); 1277 1278 // do a first loop over all sets of dofs and do identity 1279 // uniquification 1280 const unsigned int n_active_fe_indices = 1281 line->n_active_fe_indices(); 1282 for (unsigned int f = 0; f < n_active_fe_indices; ++f) 1283 for (unsigned int g = f + 1; g < n_active_fe_indices; ++g) 1284 { 1285 const unsigned int fe_index_1 = 1286 line->nth_active_fe_index(f), 1287 fe_index_2 = 1288 line->nth_active_fe_index(g); 1289 1290 if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() == 1291 dof_handler.get_fe(fe_index_2) 1292 .n_dofs_per_line()) && 1293 (dof_handler.get_fe(fe_index_1).n_dofs_per_line() > 1294 0)) 1295 { 1296 // the number of dofs per line is identical 1297 const unsigned int dofs_per_line = 1298 dof_handler.get_fe(fe_index_1).n_dofs_per_line(); 1299 1300 const auto &identities = 1301 *ensure_existence_and_return_dof_identities<1>( 1302 dof_handler.get_fe(fe_index_1), 1303 dof_handler.get_fe(fe_index_2), 1304 line_dof_identities[fe_index_1][fe_index_2]); 1305 // see if these sets of dofs are identical. the 1306 // first condition for this is that indeed there are 1307 // n identities 1308 if (identities.size() == dofs_per_line) 1309 { 1310 unsigned int i = 0; 1311 for (; i < dofs_per_line; ++i) 1312 if ((identities[i].first != i) && 1313 (identities[i].second != i)) 1314 // not an identity 1315 break; 1316 1317 if (i == dofs_per_line) 1318 { 1319 // The line dofs (i.e., the ones interior to 1320 // a line) of these two finite elements are 1321 // identical. Note that there could be 1322 // situations when one element still 1323 // dominates another, e.g.: FE_Q(2) x 1324 // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1) 1325 1326 --unique_sets_of_dofs; 1327 1328 // determine which one of both finite 1329 // elements is the dominating one. 1330 const std::set<unsigned int> fe_indices{ 1331 fe_index_1, fe_index_2}; 1332 1333 unsigned int dominating_fe_index = 1334 dof_handler.get_fe_collection() 1335 .find_dominating_fe(fe_indices, 1336 /*codim*/ dim - 1); 1337 unsigned int other_fe_index = 1338 numbers::invalid_unsigned_int; 1339 1340 if (dominating_fe_index != 1341 numbers::invalid_unsigned_int) 1342 other_fe_index = 1343 (dominating_fe_index == fe_index_1) ? 1344 fe_index_2 : 1345 fe_index_1; 1346 else 1347 { 1348 // if we haven't found a dominating 1349 // finite element, choose the one with 1350 // the lower index to be dominating 1351 dominating_fe_index = fe_index_1; 1352 other_fe_index = fe_index_2; 1353 } 1354 1355 for (unsigned int j = 0; j < dofs_per_line; 1356 ++j) 1357 { 1358 const types::global_dof_index 1359 primary_dof_index = line->dof_index( 1360 j, dominating_fe_index); 1361 const types::global_dof_index 1362 dependent_dof_index = 1363 line->dof_index(j, other_fe_index); 1364 1365 // check if we are on an interface 1366 // between a locally owned and a ghost 1367 // cell on which we need to work on. 1368 // 1369 // all degrees of freedom belonging to 1370 // dominating fe_indices or to a 1371 // processor with a higher rank have 1372 // been set at this point (either in 1373 // Phase 2, or after the first ghost 1374 // exchange in Phase 5). thus, we only 1375 // have to set the indices of degrees 1376 // of freedom that have been previously 1377 // flagged invalid. 1378 if ((dependent_dof_index == 1379 numbers::invalid_dof_index) && 1380 (primary_dof_index != 1381 numbers::invalid_dof_index)) 1382 line->set_dof_index(j, 1383 primary_dof_index, 1384 fe_index_2); 1385 } 1386 } 1387 } 1388 } 1389 } 1390 1391 // if at this point, there is only one unique set of dofs 1392 // left, then we have taken care of everything above. if there 1393 // are two, then we need to deal with them here. if there are 1394 // more, then we punt, as described in the paper (and 1395 // mentioned above) 1396 // TODO: The check for 'dim==2' was inserted by intuition. It 1397 // fixes 1398 // the previous problems with step-27 in 3D. But an 1399 // explanation for this is still required, and what we do here 1400 // is not what we describe in the paper!. 1401 if ((unique_sets_of_dofs == 2) && (dim == 2)) 1402 { 1403 const std::set<unsigned int> fe_indices = 1404 line->get_active_fe_indices(); 1405 1406 // find out which is the most dominating finite element of 1407 // the ones that are used on this line 1408 const unsigned int most_dominating_fe_index = 1409 dof_handler.get_fe_collection().find_dominating_fe( 1410 fe_indices, 1411 /*codim=*/dim - 1); 1412 1413 // if we found the most dominating element, then use this 1414 // to eliminate some of the degrees of freedom by 1415 // identification. otherwise, the code that computes 1416 // hanging node constraints will have to deal with it by 1417 // computing appropriate constraints along this face/edge 1418 if (most_dominating_fe_index != 1419 numbers::invalid_unsigned_int) 1420 { 1421 // loop over the indices of all the finite elements 1422 // that are not dominating, and identify their dofs to 1423 // the most dominating one 1424 for (const auto &other_fe_index : fe_indices) 1425 if (other_fe_index != most_dominating_fe_index) 1426 { 1427 const auto &identities = 1428 *ensure_existence_and_return_dof_identities< 1429 1>(dof_handler.get_fe( 1430 most_dominating_fe_index), 1431 dof_handler.get_fe(other_fe_index), 1432 line_dof_identities 1433 [most_dominating_fe_index] 1434 [other_fe_index]); 1435 1436 for (const auto &identity : identities) 1437 { 1438 const types::global_dof_index 1439 primary_dof_index = line->dof_index( 1440 identity.first, 1441 most_dominating_fe_index); 1442 const types::global_dof_index 1443 dependent_dof_index = 1444 line->dof_index(identity.second, 1445 other_fe_index); 1446 1447 // check if we are on an interface between 1448 // a locally owned and a ghost cell on which 1449 // we need to work on. 1450 // 1451 // all degrees of freedom belonging to 1452 // dominating fe indices or to a processor 1453 // with a higher rank have been set at this 1454 // point (either in Phase 2, or after the 1455 // first ghost exchange in Phase 5). thus, 1456 // we only have to set the indices of 1457 // degrees of freedom that have been 1458 // previously flagged invalid. 1459 if ((dependent_dof_index == 1460 numbers::invalid_dof_index) && 1461 (primary_dof_index != 1462 numbers::invalid_dof_index)) 1463 line->set_dof_index(identity.second, 1464 primary_dof_index, 1465 other_fe_index); 1466 } 1467 } 1468 } 1469 } 1470 } 1471 1472 // finally restore the user flags 1473 const_cast<dealii::Triangulation<dim, spacedim> &>( 1474 dof_handler.get_triangulation()) 1475 .load_user_flags_line(user_flags); 1476 } 1477 1478 1479 1480 /** 1481 * Merge invalid DoF indices on quads located on ghost interfaces 1482 * by a dominating valid one. 1483 */ 1484 template <int dim, int spacedim> 1485 static void merge_invalid_quad_dofs_on_ghost_interfacesinternal::DoFHandlerImplementation::Policy::Implementation1486 merge_invalid_quad_dofs_on_ghost_interfaces( 1487 DoFHandler<dim, spacedim> &dof_handler) 1488 { 1489 (void)dof_handler; 1490 Assert( 1491 dof_handler.hp_capability_enabled == true, 1492 (typename DoFHandler<dim, spacedim>::ExcNotAvailableWithoutHP())); 1493 1494 // this function should only be called for dim<3 where there are 1495 // no quad dof identies. for dim>=3, the specialization below should 1496 // take care of it 1497 Assert(dim < 3, ExcInternalError()); 1498 } 1499 1500 1501 template <int spacedim> merge_invalid_quad_dofs_on_ghost_interfacesinternal::DoFHandlerImplementation::Policy::Implementation1502 static void merge_invalid_quad_dofs_on_ghost_interfaces( 1503 DoFHandler<3, spacedim> &dof_handler) 1504 { 1505 Assert( 1506 dof_handler.hp_capability_enabled == true, 1507 (typename DoFHandler<3, spacedim>::ExcNotAvailableWithoutHP())); 1508 1509 const int dim = 3; 1510 1511 // we will mark quads that we have already treated, so first 1512 // save and clear the user flags on quads and later restore 1513 // them 1514 std::vector<bool> user_flags; 1515 dof_handler.get_triangulation().save_user_flags_quad(user_flags); 1516 const_cast<dealii::Triangulation<dim, spacedim> &>( 1517 dof_handler.get_triangulation()) 1518 .clear_user_flags_quad(); 1519 1520 // mark all quads on ghost cells 1521 for (const auto &cell : dof_handler.active_cell_iterators()) 1522 if (cell->is_ghost()) 1523 for (const auto q : cell->face_indices()) 1524 cell->quad(q)->set_user_flag(); 1525 1526 // An implementation of the algorithm described in the hp 1527 // paper, including the modification mentioned later in the 1528 // "complications in 3-d" subsections 1529 // 1530 // as explained there, we do something only if there are 1531 // exactly 2 finite elements associated with an object. if 1532 // there is only one, then there is nothing to do anyway, 1533 // and if there are 3 or more, then we can get into 1534 // trouble. note that this only happens for lines in 3d and 1535 // higher, and for quads only in 4d and higher, so this 1536 // isn't a particularly frequent case 1537 dealii::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities( 1538 dof_handler.fe_collection.size(), 1539 dof_handler.fe_collection.size(), 1540 2 /*triangle (0) or quadrilateral (1)*/); 1541 1542 for (const auto &cell : dof_handler.active_cell_iterators()) 1543 for (const auto q : cell->face_indices()) 1544 if ((cell->is_locally_owned()) && 1545 (cell->quad(q)->user_flag_set() == true) && 1546 (cell->quad(q)->n_active_fe_indices() == 2)) 1547 { 1548 const auto quad = cell->quad(q); 1549 quad->clear_user_flag(); 1550 1551 const std::set<unsigned int> fe_indices = 1552 quad->get_active_fe_indices(); 1553 1554 // find out which is the most dominating finite 1555 // element of the ones that are used on this quad 1556 const unsigned int most_dominating_fe_index = 1557 dof_handler.get_fe_collection().find_dominating_fe( 1558 fe_indices, 1559 /*codim=*/dim - 2); 1560 1561 // if we found the most dominating element, then use 1562 // this to eliminate some of the degrees of freedom 1563 // by identification. otherwise, the code that 1564 // computes hanging node constraints will have to 1565 // deal with it by computing appropriate constraints 1566 // along this face/edge 1567 if (most_dominating_fe_index != numbers::invalid_unsigned_int) 1568 { 1569 // loop over the indices of all the finite 1570 // elements that are not dominating, and 1571 // identify their dofs to the most dominating 1572 // one 1573 for (const auto &other_fe_index : fe_indices) 1574 if (other_fe_index != most_dominating_fe_index) 1575 { 1576 const auto &identities = 1577 *ensure_existence_and_return_dof_identities<2>( 1578 dof_handler.get_fe(most_dominating_fe_index), 1579 dof_handler.get_fe(other_fe_index), 1580 quad_dof_identities 1581 [most_dominating_fe_index][other_fe_index] 1582 [cell->quad(q)->reference_cell_type() == 1583 ReferenceCell::Type::Quad]); 1584 1585 for (const auto &identity : identities) 1586 { 1587 const types::global_dof_index 1588 primary_dof_index = 1589 quad->dof_index(identity.first, 1590 most_dominating_fe_index); 1591 const types::global_dof_index 1592 dependent_dof_index = 1593 quad->dof_index(identity.second, 1594 other_fe_index); 1595 1596 // check if we are on an interface between 1597 // a locally owned and a ghost cell on which 1598 // we need to work on. 1599 // 1600 // all degrees of freedom belonging to 1601 // dominating fe indices or to a processor with 1602 // a higher rank have been set at this point 1603 // (either in Phase 2, or after the first ghost 1604 // exchange in Phase 5). thus, we only have to 1605 // set the indices of degrees of freedom that 1606 // have been previously flagged invalid. 1607 if ((dependent_dof_index == 1608 numbers::invalid_dof_index) && 1609 (primary_dof_index != 1610 numbers::invalid_dof_index)) 1611 quad->set_dof_index(identity.second, 1612 primary_dof_index, 1613 other_fe_index); 1614 } 1615 } 1616 } 1617 } 1618 1619 // finally restore the user flags 1620 const_cast<dealii::Triangulation<dim, spacedim> &>( 1621 dof_handler.get_triangulation()) 1622 .load_user_flags_quad(user_flags); 1623 } 1624 1625 1626 1627 /** 1628 * After DoF unification in Phase 2, we may still have invalid DoF 1629 * indices left on ghost interfaces that are dominated by a 1630 * FiniteElement object of an adjacent cell, which could be either a 1631 * ghost cell or locally owned, whose DoF indices we did not know at the 1632 * moment of unification. After the first ghost exchange in Phase 5, we 1633 * know the indices of all dominating DoFs, and we have to assign those 1634 * invalid entries to their corresponding global value. 1635 * 1636 * This function does nothing unless the DoFHandler has hp 1637 * capabilities. 1638 */ 1639 template <int dim, int spacedim> 1640 static void merge_invalid_dof_indices_on_ghost_interfacesinternal::DoFHandlerImplementation::Policy::Implementation1641 merge_invalid_dof_indices_on_ghost_interfaces( 1642 DoFHandler<dim, spacedim> &dof_handler) 1643 { 1644 if (dof_handler.hp_capability_enabled == false) 1645 return; 1646 1647 { 1648 Threads::TaskGroup<> tasks; 1649 1650 tasks += Threads::new_task([&]() { 1651 merge_invalid_vertex_dofs_on_ghost_interfaces(dof_handler); 1652 }); 1653 1654 if (dim > 1) 1655 { 1656 tasks += Threads::new_task([&]() { 1657 merge_invalid_line_dofs_on_ghost_interfaces(dof_handler); 1658 }); 1659 } 1660 1661 if (dim > 2) 1662 { 1663 tasks += Threads::new_task([&]() { 1664 merge_invalid_quad_dofs_on_ghost_interfaces(dof_handler); 1665 }); 1666 } 1667 1668 tasks.join_all(); 1669 } 1670 1671 update_all_active_cell_dof_indices_caches(dof_handler); 1672 } 1673 1674 1675 1676 /** 1677 * Distribute degrees of freedom on all cells, or on cells with the 1678 * correct subdomain_id if the corresponding argument is not equal to 1679 * numbers::invalid_subdomain_id. Return the total number of dofs 1680 * distributed. 1681 */ 1682 template <int dim, int spacedim> 1683 static types::global_dof_index distribute_dofsinternal::DoFHandlerImplementation::Policy::Implementation1684 distribute_dofs(const types::subdomain_id subdomain_id, 1685 DoFHandler<dim, spacedim> &dof_handler) 1686 { 1687 Assert(dof_handler.get_triangulation().n_levels() > 0, 1688 ExcMessage("Empty triangulation")); 1689 1690 // Step 1: distribute dofs on all cells, but definitely 1691 // exclude artificial cells 1692 types::global_dof_index next_free_dof = 0; 1693 1694 std::vector<types::global_dof_index> dof_indices; 1695 1696 for (auto cell : dof_handler.active_cell_iterators()) 1697 if (!cell->is_artificial()) 1698 if ((subdomain_id == numbers::invalid_subdomain_id) || 1699 (cell->subdomain_id() == subdomain_id)) 1700 { 1701 dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1702 1703 // circumvent cache 1704 internal::DoFAccessorImplementation::Implementation:: 1705 get_dof_indices(*cell, 1706 dof_indices, 1707 cell->active_fe_index()); 1708 1709 for (auto &dof_index : dof_indices) 1710 if (dof_index == numbers::invalid_dof_index) 1711 dof_index = next_free_dof++; 1712 1713 cell->set_dof_indices(dof_indices); 1714 } 1715 1716 update_all_active_cell_dof_indices_caches(dof_handler); 1717 1718 return next_free_dof; 1719 } 1720 1721 1722 1723 /** 1724 * During DoF distribution, DoFs on ghost interfaces get different 1725 * indices assigned by each adjacent subdomain. We need to clarify 1726 * ownership of those DoFs by imposing a criterion, which is that 1727 * the adjacent subdomain belonging to the processor of lowest rank 1728 * prescribes its index. 1729 * 1730 * Thus, we have to invalidate all those DoFs on ghost cells that 1731 * belong to processors of lower rank than the current one, which is 1732 * the @p subdomain_id parameter. Later during ghost exchange in 1733 * Phase 5, these values will be overwritten by the correct one. The 1734 * invalidated indices will be stored in the @p renumbering parameter. 1735 */ 1736 template <int dim, int spacedim> 1737 static void invalidate_dof_indices_on_weaker_ghost_cells_for_renumberinginternal::DoFHandlerImplementation::Policy::Implementation1738 invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering( 1739 std::vector<types::global_dof_index> &renumbering, 1740 const types::subdomain_id subdomain_id, 1741 const DoFHandler<dim, spacedim> & dof_handler) 1742 { 1743 std::vector<types::global_dof_index> local_dof_indices; 1744 1745 for (const auto &cell : dof_handler.active_cell_iterators()) 1746 if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id)) 1747 { 1748 // we found a neighboring ghost cell whose subdomain 1749 // is "stronger" than our own subdomain 1750 1751 // delete all dofs that live there and that we have 1752 // previously assigned a number to (i.e. the ones on 1753 // the interface) 1754 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1755 cell->get_dof_indices(local_dof_indices); 1756 for (const auto &local_dof_index : local_dof_indices) 1757 if (local_dof_index != numbers::invalid_dof_index) 1758 renumbering[local_dof_index] = numbers::invalid_dof_index; 1759 } 1760 } 1761 1762 1763 1764 /* -------------- distribute_mg_dofs functionality ------------- */ 1765 1766 1767 1768 template <int dim, int spacedim> 1769 static types::global_dof_index distribute_dofs_on_levelinternal::DoFHandlerImplementation::Policy::Implementation1770 distribute_dofs_on_level(const types::subdomain_id level_subdomain_id, 1771 DoFHandler<dim, spacedim> &dof_handler, 1772 const unsigned int level) 1773 { 1774 Assert(dof_handler.hp_capability_enabled == false, 1775 ExcInternalError()); 1776 1777 const dealii::Triangulation<dim, spacedim> &tria = 1778 dof_handler.get_triangulation(); 1779 Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation")); 1780 if (level >= tria.n_levels()) 1781 return 0; // this is allowed for multigrid 1782 1783 types::global_dof_index next_free_dof = 0; 1784 1785 std::vector<types::global_dof_index> dof_indices; 1786 1787 for (auto cell : dof_handler.cell_iterators_on_level(level)) 1788 if ((level_subdomain_id == numbers::invalid_subdomain_id) || 1789 (cell->level_subdomain_id() == level_subdomain_id)) 1790 { 1791 dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 1792 1793 cell->get_mg_dof_indices(dof_indices); 1794 1795 for (auto &dof_index : dof_indices) 1796 if (dof_index == numbers::invalid_dof_index) 1797 dof_index = next_free_dof++; 1798 1799 cell->set_mg_dof_indices(dof_indices); 1800 } 1801 1802 return next_free_dof; 1803 } 1804 1805 1806 1807 /* --------------------- renumber_dofs functionality ---------------- */ 1808 1809 1810 /** 1811 * The part of the renumber_dofs() functionality that operates on faces. 1812 * This part is dimension dependent and so needs to be implemented in 1813 * three separate specializations of the function. 1814 * 1815 * See renumber_dofs() for the meaning of the arguments. 1816 */ 1817 template <int dim, int spacedim> 1818 static void renumber_face_dofsinternal::DoFHandlerImplementation::Policy::Implementation1819 renumber_face_dofs( 1820 const std::vector<types::global_dof_index> &new_numbers, 1821 const IndexSet & indices_we_care_about, 1822 DoFHandler<dim, spacedim> & dof_handler) 1823 { 1824 for (unsigned int d = 1; d < dim; d++) 1825 for (auto &i : dof_handler.object_dof_indices[0][d]) 1826 if (i != numbers::invalid_dof_index) 1827 i = ((indices_we_care_about.size() == 0) ? 1828 new_numbers[i] : 1829 new_numbers[indices_we_care_about.index_within_set(i)]); 1830 } 1831 1832 1833 1834 template <int dim, int spacedim> 1835 static void renumber_vertex_dofsinternal::DoFHandlerImplementation::Policy::Implementation1836 renumber_vertex_dofs( 1837 const std::vector<types::global_dof_index> &new_numbers, 1838 const IndexSet & indices_we_care_about, 1839 DoFHandler<dim, spacedim> & dof_handler, 1840 const bool check_validity) 1841 { 1842 if (dof_handler.hp_capability_enabled == false) 1843 { 1844 // we can not use cell iterators in this function since then 1845 // we would renumber the dofs on the interface of two cells 1846 // more than once. Anyway, this way it's not only more 1847 // correct but also faster; note, however, that dof numbers 1848 // may be invalid_dof_index, namely when the appropriate 1849 // vertex/line/etc is unused 1850 for (std::vector<types::global_dof_index>::iterator i = 1851 dof_handler.object_dof_indices[0][0].begin(); 1852 i != dof_handler.object_dof_indices[0][0].end(); 1853 ++i) 1854 if (*i != numbers::invalid_dof_index) 1855 *i = 1856 (indices_we_care_about.size() == 0) ? 1857 (new_numbers[*i]) : 1858 (new_numbers[indices_we_care_about.index_within_set(*i)]); 1859 else if (check_validity) 1860 // if index is invalid_dof_index: check if this one 1861 // really is unused 1862 Assert(dof_handler.get_triangulation().vertex_used( 1863 (i - dof_handler.object_dof_indices[0][0].begin()) / 1864 dof_handler.get_fe().n_dofs_per_vertex()) == false, 1865 ExcInternalError()); 1866 return; 1867 } 1868 1869 1870 for (unsigned int vertex_index = 0; 1871 vertex_index < dof_handler.get_triangulation().n_vertices(); 1872 ++vertex_index) 1873 { 1874 const unsigned int n_active_fe_indices = 1875 dealii::internal::DoFAccessorImplementation::Implementation:: 1876 n_active_fe_indices(dof_handler, 1877 0, 1878 vertex_index, 1879 std::integral_constant<int, 0>()); 1880 1881 // if this vertex is unused, then we really ought not to have 1882 // allocated any space for it, i.e., n_active_fe_indices should be 1883 // zero, and there is no space to actually store dof indices for 1884 // this vertex 1885 if (dof_handler.get_triangulation().vertex_used(vertex_index) == 1886 false) 1887 Assert(n_active_fe_indices == 0, ExcInternalError()); 1888 1889 // otherwise the vertex is used; it may still not hold any dof 1890 // indices if it is located on an artificial cell and not adjacent 1891 // to a ghost cell, but in that case there is simply nothing for 1892 // us to do 1893 for (unsigned int f = 0; f < n_active_fe_indices; ++f) 1894 { 1895 const unsigned int fe_index = 1896 dealii::internal::DoFAccessorImplementation:: 1897 Implementation::nth_active_fe_index( 1898 dof_handler, 1899 0, 1900 vertex_index, 1901 f, 1902 std::integral_constant<int, 0>()); 1903 1904 for (unsigned int d = 0; 1905 d < dof_handler.get_fe(fe_index).n_dofs_per_vertex(); 1906 ++d) 1907 { 1908 const types::global_dof_index old_dof_index = 1909 dealii::internal::DoFAccessorImplementation:: 1910 Implementation::get_dof_index( 1911 dof_handler, 1912 0, 1913 vertex_index, 1914 fe_index, 1915 d, 1916 std::integral_constant<int, 0>()); 1917 1918 // if check_validity was set, then we are to verify that 1919 // the previous indices were all valid. this really should 1920 // be the case: we allocated space for these vertex dofs, 1921 // i.e., at least one adjacent cell has a valid 1922 // active_fe_index, so there are DoFs that really live 1923 // on this vertex. if check_validity is set, then we 1924 // must make sure that they have been set to something 1925 // useful 1926 if (check_validity) 1927 Assert(old_dof_index != numbers::invalid_dof_index, 1928 ExcInternalError()); 1929 1930 if (old_dof_index != numbers::invalid_dof_index) 1931 { 1932 // In the following blocks, we first check whether 1933 // we were given an IndexSet of DoFs to touch. If not 1934 // (the first 'if' case here), then we are in the 1935 // sequential case and are allowed to touch all DoFs. 1936 // 1937 // If yes (the 'else' case), then we need to 1938 // distinguish whether the DoF whose number we want to 1939 // touch is in fact locally owned (i.e., is in the 1940 // index set) and then we can actually assign it a new 1941 // number; otherwise, we have encountered a 1942 // non-locally owned DoF for which we don't know the 1943 // new number yet and so set it to an invalid index. 1944 // This will later be fixed up after the first ghost 1945 // exchange phase when we unify hp DoFs on neighboring 1946 // cells. 1947 if (indices_we_care_about.size() == 0) 1948 dealii::internal::DoFAccessorImplementation:: 1949 Implementation::set_dof_index( 1950 dof_handler, 1951 0, 1952 vertex_index, 1953 fe_index, 1954 d, 1955 std::integral_constant<int, 0>(), 1956 new_numbers[old_dof_index]); 1957 else 1958 { 1959 if (indices_we_care_about.is_element( 1960 old_dof_index)) 1961 dealii::internal::DoFAccessorImplementation:: 1962 Implementation::set_dof_index( 1963 dof_handler, 1964 0, 1965 vertex_index, 1966 fe_index, 1967 d, 1968 std::integral_constant<int, 0>(), 1969 new_numbers[indices_we_care_about 1970 .index_within_set( 1971 old_dof_index)]); 1972 else 1973 dealii::internal::DoFAccessorImplementation:: 1974 Implementation::set_dof_index( 1975 dof_handler, 1976 0, 1977 vertex_index, 1978 fe_index, 1979 d, 1980 std::integral_constant<int, 0>(), 1981 numbers::invalid_dof_index); 1982 } 1983 } 1984 } 1985 } 1986 } 1987 } 1988 1989 1990 1991 template <int dim, int spacedim> 1992 static void renumber_cell_dofsinternal::DoFHandlerImplementation::Policy::Implementation1993 renumber_cell_dofs( 1994 const std::vector<types::global_dof_index> &new_numbers, 1995 const IndexSet & indices_we_care_about, 1996 DoFHandler<dim, spacedim> & dof_handler) 1997 { 1998 if (dof_handler.hp_capability_enabled == false) 1999 { 2000 for (unsigned int level = 0; 2001 level < dof_handler.object_dof_indices.size(); 2002 ++level) 2003 for (auto &i : dof_handler.object_dof_indices[level][dim]) 2004 if (i != numbers::invalid_dof_index) 2005 i = ((indices_we_care_about.size() == 0) ? 2006 new_numbers[i] : 2007 new_numbers[indices_we_care_about.index_within_set( 2008 i)]); 2009 return; 2010 } 2011 2012 for (const auto &cell : dof_handler.active_cell_iterators()) 2013 if (!cell->is_artificial()) 2014 { 2015 const unsigned int fe_index = cell->active_fe_index(); 2016 2017 for (unsigned int d = 0; 2018 d < dof_handler.get_fe(fe_index) 2019 .template n_dofs_per_object<dim>(); 2020 ++d) 2021 { 2022 const types::global_dof_index old_dof_index = 2023 cell->dof_index(d, fe_index); 2024 if (old_dof_index != numbers::invalid_dof_index) 2025 { 2026 // In the following blocks, we first check whether 2027 // we were given an IndexSet of DoFs to touch. If not 2028 // (the first 'if' case here), then we are in the 2029 // sequential case and are allowed to touch all DoFs. 2030 // 2031 // If yes (the 'else' case), then we need to distinguish 2032 // whether the DoF whose number we want to touch is in 2033 // fact locally owned (i.e., is in the index set) and 2034 // then we can actually assign it a new number; 2035 // otherwise, we have encountered a non-locally owned 2036 // DoF for which we don't know the new number yet and so 2037 // set it to an invalid index. This will later be fixed 2038 // up after the first ghost exchange phase when we unify 2039 // hp DoFs on neighboring cells. 2040 if (indices_we_care_about.size() == 0) 2041 cell->set_dof_index(d, 2042 new_numbers[old_dof_index], 2043 fe_index); 2044 else 2045 { 2046 if (indices_we_care_about.is_element(old_dof_index)) 2047 cell->set_dof_index( 2048 d, 2049 new_numbers[indices_we_care_about 2050 .index_within_set(old_dof_index)], 2051 fe_index); 2052 else 2053 cell->set_dof_index(d, 2054 numbers::invalid_dof_index, 2055 fe_index); 2056 } 2057 } 2058 } 2059 } 2060 } 2061 2062 2063 2064 template <int spacedim> 2065 static void renumber_face_dofsinternal::DoFHandlerImplementation::Policy::Implementation2066 renumber_face_dofs( 2067 const std::vector<types::global_dof_index> & /*new_numbers*/, 2068 const IndexSet & /*indices_we_care_about*/, 2069 DoFHandler<1, spacedim> & /*dof_handler*/) 2070 { 2071 // nothing to do in 1d since there are no separate faces -- we've 2072 // already taken care of this when dealing with the vertices 2073 } 2074 2075 2076 2077 template <int spacedim> 2078 static void renumber_face_dofsinternal::DoFHandlerImplementation::Policy::Implementation2079 renumber_face_dofs( 2080 const std::vector<types::global_dof_index> &new_numbers, 2081 const IndexSet & indices_we_care_about, 2082 DoFHandler<2, spacedim> & dof_handler) 2083 { 2084 const unsigned int dim = 2; 2085 2086 if (dof_handler.hp_capability_enabled == false) 2087 { 2088 for (unsigned int d = 1; d < dim; d++) 2089 for (auto &i : dof_handler.object_dof_indices[0][d]) 2090 if (i != numbers::invalid_dof_index) 2091 i = ((indices_we_care_about.size() == 0) ? 2092 new_numbers[i] : 2093 new_numbers[indices_we_care_about.index_within_set( 2094 i)]); 2095 return; 2096 } 2097 2098 // deal with DoFs on lines 2099 { 2100 // save user flags on lines so we can use them to mark lines 2101 // we've already treated 2102 std::vector<bool> saved_line_user_flags; 2103 const_cast<dealii::Triangulation<dim, spacedim> &>( 2104 dof_handler.get_triangulation()) 2105 .save_user_flags_line(saved_line_user_flags); 2106 const_cast<dealii::Triangulation<dim, spacedim> &>( 2107 dof_handler.get_triangulation()) 2108 .clear_user_flags_line(); 2109 2110 for (const auto &cell : dof_handler.active_cell_iterators()) 2111 if (!cell->is_artificial()) 2112 for (const auto l : cell->line_indices()) 2113 if (cell->line(l)->user_flag_set() == false) 2114 { 2115 const auto line = cell->line(l); 2116 line->set_user_flag(); 2117 2118 const unsigned int n_active_fe_indices = 2119 line->n_active_fe_indices(); 2120 2121 for (unsigned int f = 0; f < n_active_fe_indices; ++f) 2122 { 2123 const unsigned int fe_index = 2124 line->nth_active_fe_index(f); 2125 2126 for (unsigned int d = 0; 2127 d < 2128 dof_handler.get_fe(fe_index).n_dofs_per_line(); 2129 ++d) 2130 { 2131 const types::global_dof_index old_dof_index = 2132 line->dof_index(d, fe_index); 2133 if (old_dof_index != numbers::invalid_dof_index) 2134 { 2135 // In the following blocks, we first check 2136 // whether we were given an IndexSet of DoFs 2137 // to touch. If not (the first 'if' case 2138 // here), then we are in the sequential case 2139 // and are allowed to touch all DoFs. 2140 // 2141 // If yes (the 'else' case), then we need to 2142 // distinguish whether the DoF whose number we 2143 // want to touch is in fact locally owned 2144 // (i.e., is in the index set) and then we can 2145 // actually assign it a new number; otherwise, 2146 // we have encountered a non-locally owned DoF 2147 // for which we don't know the new number yet 2148 // and so set it to an invalid index. This 2149 // will later be fixed up after the first 2150 // ghost exchange phase when we unify hp DoFs 2151 // on neighboring cells. 2152 if (indices_we_care_about.size() == 0) 2153 line->set_dof_index( 2154 d, new_numbers[old_dof_index], fe_index); 2155 else 2156 { 2157 if (indices_we_care_about.is_element( 2158 old_dof_index)) 2159 line->set_dof_index( 2160 d, 2161 new_numbers[indices_we_care_about 2162 .index_within_set( 2163 old_dof_index)], 2164 fe_index); 2165 else 2166 line->set_dof_index( 2167 d, 2168 numbers::invalid_dof_index, 2169 fe_index); 2170 } 2171 } 2172 } 2173 } 2174 } 2175 2176 // at the end, restore the user 2177 // flags for the lines 2178 const_cast<dealii::Triangulation<dim, spacedim> &>( 2179 dof_handler.get_triangulation()) 2180 .load_user_flags_line(saved_line_user_flags); 2181 } 2182 } 2183 2184 2185 2186 template <int spacedim> 2187 static void renumber_face_dofsinternal::DoFHandlerImplementation::Policy::Implementation2188 renumber_face_dofs( 2189 const std::vector<types::global_dof_index> &new_numbers, 2190 const IndexSet & indices_we_care_about, 2191 DoFHandler<3, spacedim> & dof_handler) 2192 { 2193 const unsigned int dim = 3; 2194 2195 if (dof_handler.hp_capability_enabled == false) 2196 { 2197 for (unsigned int d = 1; d < dim; d++) 2198 for (auto &i : dof_handler.object_dof_indices[0][d]) 2199 if (i != numbers::invalid_dof_index) 2200 i = ((indices_we_care_about.size() == 0) ? 2201 new_numbers[i] : 2202 new_numbers[indices_we_care_about.index_within_set( 2203 i)]); 2204 return; 2205 } 2206 2207 // deal with DoFs on lines 2208 { 2209 // save user flags on lines so we can use them to mark lines 2210 // we've already treated 2211 std::vector<bool> saved_line_user_flags; 2212 const_cast<dealii::Triangulation<dim, spacedim> &>( 2213 dof_handler.get_triangulation()) 2214 .save_user_flags_line(saved_line_user_flags); 2215 const_cast<dealii::Triangulation<dim, spacedim> &>( 2216 dof_handler.get_triangulation()) 2217 .clear_user_flags_line(); 2218 2219 for (const auto &cell : dof_handler.active_cell_iterators()) 2220 if (!cell->is_artificial()) 2221 for (const auto l : cell->line_indices()) 2222 if (cell->line(l)->user_flag_set() == false) 2223 { 2224 const auto line = cell->line(l); 2225 line->set_user_flag(); 2226 2227 const unsigned int n_active_fe_indices = 2228 line->n_active_fe_indices(); 2229 2230 for (unsigned int f = 0; f < n_active_fe_indices; ++f) 2231 { 2232 const unsigned int fe_index = 2233 line->nth_active_fe_index(f); 2234 2235 for (unsigned int d = 0; 2236 d < 2237 dof_handler.get_fe(fe_index).n_dofs_per_line(); 2238 ++d) 2239 { 2240 const types::global_dof_index old_dof_index = 2241 line->dof_index(d, fe_index); 2242 if (old_dof_index != numbers::invalid_dof_index) 2243 { 2244 // In the following blocks, we first check 2245 // whether we were given an IndexSet of DoFs 2246 // to touch. If not (the first 'if' case 2247 // here), then we are in the sequential case 2248 // and are allowed to touch all DoFs. 2249 // 2250 // If yes (the 'else' case), then we need to 2251 // distinguish whether the DoF whose number we 2252 // want to touch is in fact locally owned 2253 // (i.e., is in the index set) and then we can 2254 // actually assign it a new number; otherwise, 2255 // we have encountered a non-locally owned DoF 2256 // for which we don't know the new number yet 2257 // and so set it to an invalid index. This 2258 // will later be fixed up after the first 2259 // ghost exchange phase when we unify hp DoFs 2260 // on neighboring cells. 2261 if (indices_we_care_about.size() == 0) 2262 line->set_dof_index( 2263 d, new_numbers[old_dof_index], fe_index); 2264 else if (indices_we_care_about.is_element( 2265 old_dof_index)) 2266 line->set_dof_index( 2267 d, 2268 new_numbers[indices_we_care_about 2269 .index_within_set( 2270 old_dof_index)], 2271 fe_index); 2272 else 2273 line->set_dof_index( 2274 d, numbers::invalid_dof_index, fe_index); 2275 } 2276 } 2277 } 2278 } 2279 2280 // at the end, restore the user 2281 // flags for the lines 2282 const_cast<dealii::Triangulation<dim, spacedim> &>( 2283 dof_handler.get_triangulation()) 2284 .load_user_flags_line(saved_line_user_flags); 2285 } 2286 2287 // then deal with dofs on quads 2288 { 2289 std::vector<bool> saved_quad_user_flags; 2290 const_cast<dealii::Triangulation<dim, spacedim> &>( 2291 dof_handler.get_triangulation()) 2292 .save_user_flags_quad(saved_quad_user_flags); 2293 const_cast<dealii::Triangulation<dim, spacedim> &>( 2294 dof_handler.get_triangulation()) 2295 .clear_user_flags_quad(); 2296 2297 for (const auto &cell : dof_handler.active_cell_iterators()) 2298 if (!cell->is_artificial()) 2299 for (const auto q : cell->face_indices()) 2300 if (cell->quad(q)->user_flag_set() == false) 2301 { 2302 const auto quad = cell->quad(q); 2303 quad->set_user_flag(); 2304 2305 const unsigned int n_active_fe_indices = 2306 quad->n_active_fe_indices(); 2307 2308 for (unsigned int f = 0; f < n_active_fe_indices; ++f) 2309 { 2310 const unsigned int fe_index = 2311 quad->nth_active_fe_index(f); 2312 2313 for (unsigned int d = 0; 2314 d < 2315 dof_handler.get_fe(fe_index).n_dofs_per_quad(); 2316 ++d) 2317 { 2318 const types::global_dof_index old_dof_index = 2319 quad->dof_index(d, fe_index); 2320 if (old_dof_index != numbers::invalid_dof_index) 2321 { 2322 // In the following blocks, we first check 2323 // whether we were given an IndexSet of DoFs 2324 // to touch. If not (the first 'if' case 2325 // here), then we are in the sequential case 2326 // and are allowed to touch all DoFs. 2327 // 2328 // If yes (the 'else' case), then we need to 2329 // distinguish whether the DoF whose number we 2330 // want to touch is in fact locally owned 2331 // (i.e., is in the index set) and then we can 2332 // actually assign it a new number; otherwise, 2333 // we have encountered a non-locally owned DoF 2334 // for which we don't know the new number yet 2335 // and so set it to an invalid index. This 2336 // will later be fixed up after the first 2337 // ghost exchange phase when we unify hp DoFs 2338 // on neighboring cells. 2339 if (indices_we_care_about.size() == 0) 2340 quad->set_dof_index( 2341 d, new_numbers[old_dof_index], fe_index); 2342 else 2343 { 2344 if (indices_we_care_about.is_element( 2345 old_dof_index)) 2346 quad->set_dof_index( 2347 d, 2348 new_numbers[indices_we_care_about 2349 .index_within_set( 2350 old_dof_index)], 2351 fe_index); 2352 else 2353 quad->set_dof_index( 2354 d, 2355 numbers::invalid_dof_index, 2356 fe_index); 2357 } 2358 } 2359 } 2360 } 2361 } 2362 2363 // at the end, restore the user flags for the quads 2364 const_cast<dealii::Triangulation<dim, spacedim> &>( 2365 dof_handler.get_triangulation()) 2366 .load_user_flags_quad(saved_quad_user_flags); 2367 } 2368 } 2369 2370 2371 2372 /** 2373 * Implementation of DoFHandler::renumber_dofs() 2374 * 2375 * If the second argument has any elements set, elements of 2376 * the then the vector of new numbers do not relate to the old 2377 * DoF number but instead to the index of the old DoF number 2378 * within the set of locally owned DoFs. 2379 * 2380 * (The IndexSet argument is not used in 1d because we only need 2381 * it for parallel meshes and 1d doesn't support that right now.) 2382 */ 2383 template <int dim, int space_dim> 2384 static void renumber_dofsinternal::DoFHandlerImplementation::Policy::Implementation2385 renumber_dofs(const std::vector<types::global_dof_index> &new_numbers, 2386 const IndexSet & indices_we_care_about, 2387 const DoFHandler<dim, space_dim> &dof_handler, 2388 const bool check_validity) 2389 { 2390 if (dim == 1) 2391 Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented()); 2392 2393 // renumber DoF indices on vertices, cells, and faces. this 2394 // can be done in parallel because the respective functions 2395 // work on separate data structures 2396 Threads::TaskGroup<> tasks; 2397 tasks += Threads::new_task([&]() { 2398 renumber_vertex_dofs(new_numbers, 2399 indices_we_care_about, 2400 const_cast<DoFHandler<dim, space_dim> &>( 2401 dof_handler), 2402 check_validity); 2403 }); 2404 tasks += Threads::new_task([&]() { 2405 renumber_face_dofs(new_numbers, 2406 indices_we_care_about, 2407 const_cast<DoFHandler<dim, space_dim> &>( 2408 dof_handler)); 2409 }); 2410 tasks += Threads::new_task([&]() { 2411 renumber_cell_dofs(new_numbers, 2412 indices_we_care_about, 2413 const_cast<DoFHandler<dim, space_dim> &>( 2414 dof_handler)); 2415 }); 2416 tasks.join_all(); 2417 2418 // update the cache used for cell dof indices 2419 update_all_active_cell_dof_indices_caches( 2420 const_cast<DoFHandler<dim, space_dim> &>(dof_handler)); 2421 } 2422 2423 2424 2425 /* --------------------- renumber_mg_dofs functionality ---------------- 2426 */ 2427 2428 /** 2429 * The part of the renumber_mg_dofs() functionality that is dimension 2430 * independent because it renumbers the DoF indices on vertex interiors 2431 * (which exist for all dimensions). 2432 * 2433 * See renumber_mg_dofs() for the meaning of the arguments. 2434 */ 2435 template <int dim, int spacedim> 2436 static void renumber_vertex_mg_dofsinternal::DoFHandlerImplementation::Policy::Implementation2437 renumber_vertex_mg_dofs( 2438 const std::vector<dealii::types::global_dof_index> &new_numbers, 2439 const IndexSet & indices_we_care_about, 2440 DoFHandler<dim, spacedim> &dof_handler, 2441 const unsigned int level, 2442 const bool check_validity) 2443 { 2444 (void)check_validity; 2445 Assert(level < dof_handler.get_triangulation().n_levels(), 2446 ExcInternalError()); 2447 2448 for (typename std::vector< 2449 typename DoFHandler<dim, spacedim>::MGVertexDoFs>::iterator i = 2450 dof_handler.mg_vertex_dofs.begin(); 2451 i != dof_handler.mg_vertex_dofs.end(); 2452 ++i) 2453 // if the present vertex lives on the current level 2454 if ((i->get_coarsest_level() <= level) && 2455 (i->get_finest_level() >= level)) 2456 for (unsigned int d = 0; 2457 d < dof_handler.get_fe().n_dofs_per_vertex(); 2458 ++d) 2459 { 2460 const dealii::types::global_dof_index idx = 2461 i->get_index(level, 2462 d, 2463 dof_handler.get_fe().n_dofs_per_vertex()); 2464 2465 if (idx != numbers::invalid_dof_index) 2466 { 2467 Assert(check_validity == false || 2468 (indices_we_care_about.size() > 0 ? 2469 indices_we_care_about.is_element(idx) : 2470 (idx < new_numbers.size())), 2471 ExcInternalError()); 2472 i->set_index(level, 2473 d, 2474 dof_handler.get_fe().n_dofs_per_vertex(), 2475 (indices_we_care_about.size() == 0) ? 2476 (new_numbers[idx]) : 2477 (new_numbers[indices_we_care_about 2478 .index_within_set(idx)])); 2479 } 2480 } 2481 } 2482 2483 2484 2485 /** 2486 * The part of the renumber_dofs() functionality that is dimension 2487 * independent because it renumbers the DoF indices on cell interiors 2488 * (which exist for all dimensions). 2489 * 2490 * See renumber_mg_dofs() for the meaning of the arguments. 2491 */ 2492 template <int dim, int spacedim> 2493 static void renumber_cell_mg_dofsinternal::DoFHandlerImplementation::Policy::Implementation2494 renumber_cell_mg_dofs( 2495 const std::vector<dealii::types::global_dof_index> &new_numbers, 2496 const IndexSet & indices_we_care_about, 2497 DoFHandler<dim, spacedim> &dof_handler, 2498 const unsigned int level) 2499 { 2500 for (std::vector<types::global_dof_index>::iterator i = 2501 dof_handler.mg_levels[level]->dof_object.dofs.begin(); 2502 i != dof_handler.mg_levels[level]->dof_object.dofs.end(); 2503 ++i) 2504 { 2505 if (*i != numbers::invalid_dof_index) 2506 { 2507 Assert((indices_we_care_about.size() > 0 ? 2508 indices_we_care_about.is_element(*i) : 2509 (*i < new_numbers.size())), 2510 ExcInternalError()); 2511 *i = 2512 (indices_we_care_about.size() == 0) ? 2513 (new_numbers[*i]) : 2514 (new_numbers[indices_we_care_about.index_within_set(*i)]); 2515 } 2516 } 2517 } 2518 2519 2520 2521 /** 2522 * The part of the renumber_mg_dofs() functionality that operates on 2523 * faces. This part is dimension dependent and so needs to be 2524 * implemented in three separate specializations of the function. 2525 * 2526 * See renumber_mg_dofs() for the meaning of the arguments. 2527 */ 2528 template <int spacedim> 2529 static void renumber_face_mg_dofsinternal::DoFHandlerImplementation::Policy::Implementation2530 renumber_face_mg_dofs( 2531 const std::vector<types::global_dof_index> & /*new_numbers*/, 2532 const IndexSet & /*indices_we_care_about*/, 2533 DoFHandler<1, spacedim> & /*dof_handler*/, 2534 const unsigned int /*level*/, 2535 const bool /*check_validity*/) 2536 { 2537 // nothing to do in 1d because there are no separate faces 2538 } 2539 2540 2541 2542 template <int spacedim> 2543 static void renumber_face_mg_dofsinternal::DoFHandlerImplementation::Policy::Implementation2544 renumber_face_mg_dofs( 2545 const std::vector<dealii::types::global_dof_index> &new_numbers, 2546 const IndexSet & indices_we_care_about, 2547 DoFHandler<2, spacedim> &dof_handler, 2548 const unsigned int level, 2549 const bool check_validity) 2550 { 2551 if (dof_handler.get_fe().n_dofs_per_line() > 0) 2552 { 2553 // save user flags as they will be modified 2554 std::vector<bool> user_flags; 2555 dof_handler.get_triangulation().save_user_flags(user_flags); 2556 const_cast<dealii::Triangulation<2, spacedim> &>( 2557 dof_handler.get_triangulation()) 2558 .clear_user_flags(); 2559 2560 // flag all lines adjacent to cells of the current 2561 // level, as those lines logically belong to the same 2562 // level as the cell, at least for for isotropic 2563 // refinement 2564 for (const auto &cell : 2565 dof_handler.cell_iterators_on_level(level)) 2566 if (cell->level_subdomain_id() != 2567 numbers::artificial_subdomain_id) 2568 for (const unsigned int line : cell->face_indices()) 2569 cell->face(line)->set_user_flag(); 2570 2571 for (const auto &cell : 2572 dof_handler.cell_iterators_on_level(level)) 2573 for (const auto l : cell->line_indices()) 2574 if (cell->line(l)->user_flag_set()) 2575 { 2576 for (unsigned int d = 0; 2577 d < dof_handler.get_fe().n_dofs_per_line(); 2578 ++d) 2579 { 2580 const dealii::types::global_dof_index idx = 2581 cell->line(l)->mg_dof_index(level, d); 2582 if (check_validity) 2583 Assert(idx != numbers::invalid_dof_index, 2584 ExcInternalError()); 2585 2586 if (idx != numbers::invalid_dof_index) 2587 cell->line(l)->set_mg_dof_index( 2588 level, 2589 d, 2590 ((indices_we_care_about.size() == 0) ? 2591 new_numbers[idx] : 2592 new_numbers[indices_we_care_about 2593 .index_within_set(idx)])); 2594 } 2595 cell->line(l)->clear_user_flag(); 2596 } 2597 // finally, restore user flags 2598 const_cast<dealii::Triangulation<2, spacedim> &>( 2599 dof_handler.get_triangulation()) 2600 .load_user_flags(user_flags); 2601 } 2602 } 2603 2604 2605 2606 template <int spacedim> 2607 static void renumber_face_mg_dofsinternal::DoFHandlerImplementation::Policy::Implementation2608 renumber_face_mg_dofs( 2609 const std::vector<dealii::types::global_dof_index> &new_numbers, 2610 const IndexSet & indices_we_care_about, 2611 DoFHandler<3, spacedim> &dof_handler, 2612 const unsigned int level, 2613 const bool check_validity) 2614 { 2615 if (dof_handler.get_fe().n_dofs_per_line() > 0 || 2616 dof_handler.get_fe().max_dofs_per_quad() > 0) 2617 { 2618 // save user flags as they will be modified 2619 std::vector<bool> user_flags; 2620 dof_handler.get_triangulation().save_user_flags(user_flags); 2621 const_cast<dealii::Triangulation<3, spacedim> &>( 2622 dof_handler.get_triangulation()) 2623 .clear_user_flags(); 2624 2625 // flag all lines adjacent to cells of the current 2626 // level, as those lines logically belong to the same 2627 // level as the cell, at least for isotropic refinement 2628 for (const auto &cell : 2629 dof_handler.cell_iterators_on_level(level)) 2630 if (cell->level_subdomain_id() != 2631 numbers::artificial_subdomain_id) 2632 for (const auto line : cell->line_indices()) 2633 cell->line(line)->set_user_flag(); 2634 2635 for (const auto &cell : 2636 dof_handler.cell_iterators_on_level(level)) 2637 for (const auto l : cell->line_indices()) 2638 if (cell->line(l)->user_flag_set()) 2639 { 2640 for (unsigned int d = 0; 2641 d < dof_handler.get_fe().n_dofs_per_line(); 2642 ++d) 2643 { 2644 const dealii::types::global_dof_index idx = 2645 cell->line(l)->mg_dof_index(level, d); 2646 if (check_validity) 2647 Assert(idx != numbers::invalid_dof_index, 2648 ExcInternalError()); 2649 2650 if (idx != numbers::invalid_dof_index) 2651 cell->line(l)->set_mg_dof_index( 2652 level, 2653 d, 2654 ((indices_we_care_about.size() == 0) ? 2655 new_numbers[idx] : 2656 new_numbers[indices_we_care_about 2657 .index_within_set(idx)])); 2658 } 2659 cell->line(l)->clear_user_flag(); 2660 } 2661 2662 // flag all quads adjacent to cells of the current level, as 2663 // those quads logically belong to the same level as the cell, 2664 // at least for isotropic refinement 2665 for (const auto &cell : 2666 dof_handler.cell_iterators_on_level(level)) 2667 if (cell->level_subdomain_id() != 2668 numbers::artificial_subdomain_id) 2669 for (const auto quad : cell->face_indices()) 2670 cell->quad(quad)->set_user_flag(); 2671 2672 for (const auto &cell : dof_handler.cell_iterators()) 2673 for (const auto l : cell->face_indices()) 2674 if (cell->quad(l)->user_flag_set()) 2675 { 2676 for (unsigned int d = 0; 2677 d < dof_handler.get_fe().n_dofs_per_quad(); 2678 ++d) 2679 { 2680 const dealii::types::global_dof_index idx = 2681 cell->quad(l)->mg_dof_index(level, d); 2682 if (check_validity) 2683 Assert(idx != numbers::invalid_dof_index, 2684 ExcInternalError()); 2685 2686 if (idx != numbers::invalid_dof_index) 2687 cell->quad(l)->set_mg_dof_index( 2688 level, 2689 d, 2690 ((indices_we_care_about.size() == 0) ? 2691 new_numbers[idx] : 2692 new_numbers[indices_we_care_about 2693 .index_within_set(idx)])); 2694 } 2695 cell->quad(l)->clear_user_flag(); 2696 } 2697 2698 // finally, restore user flags 2699 const_cast<dealii::Triangulation<3, spacedim> &>( 2700 dof_handler.get_triangulation()) 2701 .load_user_flags(user_flags); 2702 } 2703 } 2704 2705 2706 2707 template <int dim, int spacedim> 2708 static void renumber_mg_dofsinternal::DoFHandlerImplementation::Policy::Implementation2709 renumber_mg_dofs( 2710 const std::vector<dealii::types::global_dof_index> &new_numbers, 2711 const IndexSet & indices_we_care_about, 2712 DoFHandler<dim, spacedim> &dof_handler, 2713 const unsigned int level, 2714 const bool check_validity) 2715 { 2716 Assert( 2717 dof_handler.hp_capability_enabled == false, 2718 (typename DoFHandler<dim, spacedim>::ExcNotImplementedWithHP())); 2719 2720 Assert(level < dof_handler.get_triangulation().n_global_levels(), 2721 ExcInternalError()); 2722 2723 // renumber DoF indices on vertices, cells, and faces. this 2724 // can be done in parallel because the respective functions 2725 // work on separate data structures 2726 Threads::TaskGroup<> tasks; 2727 tasks += Threads::new_task([&]() { 2728 renumber_vertex_mg_dofs(new_numbers, 2729 indices_we_care_about, 2730 dof_handler, 2731 level, 2732 check_validity); 2733 }); 2734 tasks += Threads::new_task([&]() { 2735 renumber_face_mg_dofs(new_numbers, 2736 indices_we_care_about, 2737 dof_handler, 2738 level, 2739 check_validity); 2740 }); 2741 tasks += Threads::new_task([&]() { 2742 renumber_cell_mg_dofs(new_numbers, 2743 indices_we_care_about, 2744 dof_handler, 2745 level); 2746 }); 2747 tasks.join_all(); 2748 } 2749 }; 2750 2751 2752 2753 /* --------------------- class Sequential ---------------- */ 2754 2755 2756 2757 template <int dim, int spacedim> Sequential(DoFHandler<dim,spacedim> & dof_handler)2758 Sequential<dim, spacedim>::Sequential( 2759 DoFHandler<dim, spacedim> &dof_handler) 2760 : dof_handler(&dof_handler) 2761 {} 2762 2763 2764 2765 template <int dim, int spacedim> 2766 NumberCache distribute_dofs() const2767 Sequential<dim, spacedim>::distribute_dofs() const 2768 { 2769 const types::global_dof_index n_initial_dofs = 2770 Implementation::distribute_dofs(numbers::invalid_subdomain_id, 2771 *dof_handler); 2772 2773 const types::global_dof_index n_dofs = 2774 Implementation::unify_dof_indices(*dof_handler, 2775 n_initial_dofs, 2776 /*check_validity=*/true); 2777 2778 // return a sequential, complete index set 2779 return NumberCache(n_dofs); 2780 } 2781 2782 2783 2784 template <int dim, int spacedim> 2785 std::vector<NumberCache> distribute_mg_dofs() const2786 Sequential<dim, spacedim>::distribute_mg_dofs() const 2787 { 2788 std::vector<bool> user_flags; 2789 dof_handler->get_triangulation().save_user_flags(user_flags); 2790 2791 const_cast<dealii::Triangulation<dim, spacedim> &>( 2792 dof_handler->get_triangulation()) 2793 .clear_user_flags(); 2794 2795 std::vector<NumberCache> number_caches; 2796 number_caches.reserve(dof_handler->get_triangulation().n_levels()); 2797 for (unsigned int level = 0; 2798 level < dof_handler->get_triangulation().n_levels(); 2799 ++level) 2800 { 2801 // first distribute dofs on this level 2802 const types::global_dof_index n_level_dofs = 2803 Implementation::distribute_dofs_on_level( 2804 numbers::invalid_subdomain_id, *dof_handler, level); 2805 2806 // then add a complete, sequential index set 2807 number_caches.emplace_back(n_level_dofs); 2808 } 2809 2810 const_cast<dealii::Triangulation<dim, spacedim> &>( 2811 dof_handler->get_triangulation()) 2812 .load_user_flags(user_flags); 2813 2814 return number_caches; 2815 } 2816 2817 2818 2819 template <int dim, int spacedim> 2820 NumberCache renumber_dofs(const std::vector<types::global_dof_index> & new_numbers) const2821 Sequential<dim, spacedim>::renumber_dofs( 2822 const std::vector<types::global_dof_index> &new_numbers) const 2823 { 2824 Implementation::renumber_dofs(new_numbers, 2825 IndexSet(0), 2826 *dof_handler, 2827 /*check_validity=*/true); 2828 2829 // return a sequential, complete index set. take into account that the 2830 // number of DoF indices may in fact be smaller than there were before 2831 // if some previously separately numbered dofs have been identified. 2832 // this is, for example, what we do when the DoFHandler has hp 2833 // capabilities enabled: it first enumerates all DoFs on cells 2834 // independently, and then unifies some located at vertices or faces; 2835 // this leaves us with fewer DoFs than there were before, so use the 2836 // largest index as the one to determine the size of the index space 2837 return NumberCache( 2838 *std::max_element(new_numbers.begin(), new_numbers.end()) + 1); 2839 } 2840 2841 2842 2843 template <int dim, int spacedim> 2844 NumberCache renumber_mg_dofs(const unsigned int level,const std::vector<types::global_dof_index> & new_numbers) const2845 Sequential<dim, spacedim>::renumber_mg_dofs( 2846 const unsigned int level, 2847 const std::vector<types::global_dof_index> &new_numbers) const 2848 { 2849 Implementation::renumber_mg_dofs( 2850 new_numbers, IndexSet(0), *dof_handler, level, true); 2851 2852 // return a sequential, complete index set 2853 return NumberCache(new_numbers.size()); 2854 } 2855 2856 2857 /* --------------------- class ParallelShared ---------------- */ 2858 2859 2860 template <int dim, int spacedim> ParallelShared(DoFHandler<dim,spacedim> & dof_handler)2861 ParallelShared<dim, spacedim>::ParallelShared( 2862 DoFHandler<dim, spacedim> &dof_handler) 2863 : dof_handler(&dof_handler) 2864 {} 2865 2866 2867 2868 namespace 2869 { 2870 /** 2871 * This function is a variation of 2872 * DoFTools::get_dof_subdomain_association() with the exception that it 2873 * is (i) specialized for parallel::shared::Triangulation objects, and 2874 * (ii) does not assume that the internal number cache of the DoFHandler 2875 * has already been set up. In can, consequently, be called at a point 2876 * when we are still distributing degrees of freedom. 2877 */ 2878 template <int dim, int spacedim> 2879 std::vector<types::subdomain_id> get_dof_subdomain_association(const DoFHandler<dim,spacedim> & dof_handler,const types::global_dof_index n_dofs,const unsigned int n_procs)2880 get_dof_subdomain_association( 2881 const DoFHandler<dim, spacedim> &dof_handler, 2882 const types::global_dof_index n_dofs, 2883 const unsigned int n_procs) 2884 { 2885 (void)n_procs; 2886 std::vector<types::subdomain_id> subdomain_association( 2887 n_dofs, numbers::invalid_subdomain_id); 2888 std::vector<types::global_dof_index> local_dof_indices; 2889 local_dof_indices.reserve( 2890 dof_handler.get_fe_collection().max_dofs_per_cell()); 2891 2892 // loop over all cells and record which subdomain a DoF belongs to. 2893 // give to the smaller subdomain_id in case it is on an interface 2894 for (const auto &cell : dof_handler.active_cell_iterators()) 2895 { 2896 // get the owner of the cell; note that we have made sure above 2897 // that all cells are either locally owned or ghosts (not 2898 // artificial), so this call will always yield the true owner 2899 const types::subdomain_id subdomain_id = cell->subdomain_id(); 2900 const unsigned int dofs_per_cell = 2901 cell->get_fe().n_dofs_per_cell(); 2902 local_dof_indices.resize(dofs_per_cell); 2903 cell->get_dof_indices(local_dof_indices); 2904 2905 // set subdomain ids. if dofs already have their values set then 2906 // they must be on partition interfaces. In that case assign them 2907 // to the processor with the smaller subdomain id. 2908 for (unsigned int i = 0; i < dofs_per_cell; ++i) 2909 if (subdomain_association[local_dof_indices[i]] == 2910 numbers::invalid_subdomain_id) 2911 subdomain_association[local_dof_indices[i]] = subdomain_id; 2912 else if (subdomain_association[local_dof_indices[i]] > 2913 subdomain_id) 2914 { 2915 subdomain_association[local_dof_indices[i]] = subdomain_id; 2916 } 2917 } 2918 2919 Assert(std::find(subdomain_association.begin(), 2920 subdomain_association.end(), 2921 numbers::invalid_subdomain_id) == 2922 subdomain_association.end(), 2923 ExcInternalError()); 2924 2925 Assert(*std::max_element(subdomain_association.begin(), 2926 subdomain_association.end()) < n_procs, 2927 ExcInternalError()); 2928 2929 return subdomain_association; 2930 } 2931 2932 2933 /** 2934 * level subdomain association. Similar to the above function only 2935 * for level meshes. This function assigns boundary dofs in 2936 * the same way as parallel::DistributedTriangulationBase (proc with 2937 * smallest index) instead of the coin flip method above. 2938 */ 2939 template <int dim, int spacedim> 2940 std::vector<types::subdomain_id> get_dof_level_subdomain_association(const DoFHandler<dim,spacedim> & dof_handler,const types::global_dof_index n_dofs_on_level,const unsigned int n_procs,const unsigned int level)2941 get_dof_level_subdomain_association( 2942 const DoFHandler<dim, spacedim> &dof_handler, 2943 const types::global_dof_index n_dofs_on_level, 2944 const unsigned int n_procs, 2945 const unsigned int level) 2946 { 2947 (void)n_procs; 2948 std::vector<types::subdomain_id> level_subdomain_association( 2949 n_dofs_on_level, numbers::invalid_subdomain_id); 2950 std::vector<types::global_dof_index> local_dof_indices; 2951 local_dof_indices.reserve( 2952 dof_handler.get_fe_collection().max_dofs_per_cell()); 2953 2954 // loop over all cells and record which subdomain a DoF belongs to. 2955 // interface goes to proccessor with smaller subdomain id 2956 for (const auto &cell : dof_handler.cell_iterators_on_level(level)) 2957 { 2958 // get the owner of the cell; note that we have made sure above 2959 // that all cells are either locally owned or ghosts (not 2960 // artificial), so this call will always yield the true owner 2961 const types::subdomain_id level_subdomain_id = 2962 cell->level_subdomain_id(); 2963 const unsigned int dofs_per_cell = 2964 cell->get_fe().n_dofs_per_cell(); 2965 local_dof_indices.resize(dofs_per_cell); 2966 cell->get_mg_dof_indices(local_dof_indices); 2967 2968 // set level subdomain ids. if dofs already have their values set 2969 // then they must be on partition interfaces. In that case assign 2970 // them to the processor with the smaller subdomain id. 2971 for (unsigned int i = 0; i < dofs_per_cell; ++i) 2972 if (level_subdomain_association[local_dof_indices[i]] == 2973 numbers::invalid_subdomain_id) 2974 level_subdomain_association[local_dof_indices[i]] = 2975 level_subdomain_id; 2976 else if (level_subdomain_association[local_dof_indices[i]] > 2977 level_subdomain_id) 2978 { 2979 level_subdomain_association[local_dof_indices[i]] = 2980 level_subdomain_id; 2981 } 2982 } 2983 2984 Assert(std::find(level_subdomain_association.begin(), 2985 level_subdomain_association.end(), 2986 numbers::invalid_subdomain_id) == 2987 level_subdomain_association.end(), 2988 ExcInternalError()); 2989 2990 Assert(*std::max_element(level_subdomain_association.begin(), 2991 level_subdomain_association.end()) < n_procs, 2992 ExcInternalError()); 2993 2994 return level_subdomain_association; 2995 } 2996 } // namespace 2997 2998 2999 3000 template <int dim, int spacedim> 3001 NumberCache distribute_dofs() const3002 ParallelShared<dim, spacedim>::distribute_dofs() const 3003 { 3004 const dealii::parallel::shared::Triangulation<dim, spacedim> *tr = 3005 (dynamic_cast< 3006 const dealii::parallel::shared::Triangulation<dim, spacedim> *>( 3007 &this->dof_handler->get_triangulation())); 3008 Assert(tr != nullptr, ExcInternalError()); 3009 3010 const unsigned int n_procs = 3011 Utilities::MPI::n_mpi_processes(tr->get_communicator()); 3012 3013 // If the underlying shared::Tria allows artificial cells, 3014 // then save the current set of subdomain ids, and set 3015 // subdomain ids to the "true" owner of each cell. we later 3016 // restore these flags 3017 std::vector<types::subdomain_id> saved_subdomain_ids; 3018 if (tr->with_artificial_cells()) 3019 { 3020 saved_subdomain_ids.resize(tr->n_active_cells()); 3021 3022 const std::vector<types::subdomain_id> &true_subdomain_ids = 3023 tr->get_true_subdomain_ids_of_cells(); 3024 3025 for (const auto &cell : tr->active_cell_iterators()) 3026 { 3027 const unsigned int index = cell->active_cell_index(); 3028 saved_subdomain_ids[index] = cell->subdomain_id(); 3029 cell->set_subdomain_id(true_subdomain_ids[index]); 3030 } 3031 } 3032 3033 // first let the sequential algorithm do its magic. it is going to 3034 // enumerate DoFs on all cells, regardless of owner 3035 const types::global_dof_index n_initial_dofs = 3036 Implementation::distribute_dofs(numbers::invalid_subdomain_id, 3037 *this->dof_handler); 3038 3039 const types::global_dof_index n_dofs = 3040 Implementation::unify_dof_indices(*this->dof_handler, 3041 n_initial_dofs, 3042 /*check_validity=*/true); 3043 3044 // then re-enumerate them based on their subdomain association. 3045 // for this, we first have to identify for each current DoF 3046 // index which subdomain they belong to. ideally, we would 3047 // like to call DoFRenumbering::subdomain_wise(), but 3048 // because the NumberCache of the current DoFHandler is not 3049 // fully set up yet, we can't quite do that. also, that 3050 // function has to deal with other kinds of triangulations as 3051 // well, whereas we here know what kind of triangulation 3052 // we have and can simplify the code accordingly 3053 std::vector<types::global_dof_index> new_dof_indices( 3054 n_dofs, enumeration_dof_index); 3055 { 3056 // first get the association of each dof with a subdomain and 3057 // determine the total number of subdomain ids used 3058 const std::vector<types::subdomain_id> subdomain_association = 3059 get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs); 3060 3061 // then renumber the subdomains by first looking at those belonging 3062 // to subdomain 0, then those of subdomain 1, etc. note that the 3063 // algorithm is stable, i.e. if two dofs i,j have i<j and belong to 3064 // the same subdomain, then they will be in this order also after 3065 // reordering 3066 types::global_dof_index next_free_index = 0; 3067 for (types::subdomain_id subdomain = 0; subdomain < n_procs; 3068 ++subdomain) 3069 for (types::global_dof_index i = 0; i < n_dofs; ++i) 3070 if (subdomain_association[i] == subdomain) 3071 { 3072 Assert(new_dof_indices[i] == enumeration_dof_index, 3073 ExcInternalError()); 3074 new_dof_indices[i] = next_free_index; 3075 ++next_free_index; 3076 } 3077 3078 // we should have numbered all dofs 3079 Assert(next_free_index == n_dofs, ExcInternalError()); 3080 Assert(std::find(new_dof_indices.begin(), 3081 new_dof_indices.end(), 3082 enumeration_dof_index) == new_dof_indices.end(), 3083 ExcInternalError()); 3084 } 3085 // finally do the renumbering. we can use the sequential 3086 // version of the function because we do things on all 3087 // cells and all cells have their subdomain ids and DoFs 3088 // correctly set 3089 Implementation::renumber_dofs(new_dof_indices, 3090 IndexSet(0), 3091 *this->dof_handler, 3092 /*check_validity=*/true); 3093 3094 // update the number cache. for this, we first have to find the 3095 // subdomain association for each DoF again following renumbering, from 3096 // which we can then compute the IndexSets of locally owned DoFs for all 3097 // processors. all other fields then follow from this 3098 // 3099 // given the way we enumerate degrees of freedom, the locally owned 3100 // ranges must all be contiguous and consecutive. this makes filling 3101 // the IndexSets cheap. an assertion at the top verifies that this 3102 // assumption is true 3103 const std::vector<types::subdomain_id> subdomain_association = 3104 get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs); 3105 3106 for (unsigned int i = 1; i < n_dofs; ++i) 3107 Assert(subdomain_association[i] >= subdomain_association[i - 1], 3108 ExcInternalError()); 3109 3110 std::vector<IndexSet> locally_owned_dofs_per_processor( 3111 n_procs, IndexSet(n_dofs)); 3112 { 3113 // we know that the set of subdomain indices is contiguous from 3114 // the assertion above; find the start and end index for each 3115 // processor, taking into account that sometimes a processor 3116 // may not in fact have any DoFs at all. we do the latter 3117 // by just identifying contiguous ranges of subdomain_ids 3118 // and filling IndexSets for those subdomains; subdomains 3119 // that don't appear will lead to IndexSets that are simply 3120 // never touched and remain empty as initialized above. 3121 unsigned int start_index = 0; 3122 unsigned int end_index = 0; 3123 while (start_index < n_dofs) 3124 { 3125 while ((end_index) < n_dofs && 3126 (subdomain_association[end_index] == 3127 subdomain_association[start_index])) 3128 ++end_index; 3129 3130 // we've now identified a range of same indices. set that 3131 // range in the corresponding IndexSet 3132 if (end_index > start_index) 3133 { 3134 const unsigned int subdomain_owner = 3135 subdomain_association[start_index]; 3136 locally_owned_dofs_per_processor[subdomain_owner].add_range( 3137 start_index, end_index); 3138 } 3139 3140 // then move on to thinking about the next range 3141 start_index = end_index; 3142 } 3143 } 3144 3145 // finally, restore current subdomain ids 3146 if (tr->with_artificial_cells()) 3147 for (const auto &cell : tr->active_cell_iterators()) 3148 cell->set_subdomain_id( 3149 saved_subdomain_ids[cell->active_cell_index()]); 3150 3151 // return a NumberCache object made up from the sets of locally 3152 // owned DoFs 3153 return NumberCache( 3154 locally_owned_dofs_per_processor, 3155 this->dof_handler->get_triangulation().locally_owned_subdomain()); 3156 } 3157 3158 3159 3160 template <int dim, int spacedim> 3161 std::vector<NumberCache> distribute_mg_dofs() const3162 ParallelShared<dim, spacedim>::distribute_mg_dofs() const 3163 { 3164 const dealii::parallel::shared::Triangulation<dim, spacedim> *tr = 3165 (dynamic_cast< 3166 const dealii::parallel::shared::Triangulation<dim, spacedim> *>( 3167 &this->dof_handler->get_triangulation())); 3168 Assert(tr != nullptr, ExcInternalError()); 3169 3170 const unsigned int n_procs = 3171 Utilities::MPI::n_mpi_processes(tr->get_communicator()); 3172 const unsigned int n_levels = tr->n_global_levels(); 3173 3174 std::vector<NumberCache> number_caches; 3175 number_caches.reserve(n_levels); 3176 3177 // We create an index set for each level 3178 for (unsigned int lvl = 0; lvl < n_levels; ++lvl) 3179 { 3180 // If the underlying shared::Tria allows artificial cells, 3181 // then save the current set of level subdomain ids, and set 3182 // subdomain ids to the "true" owner of each cell. we later 3183 // restore these flags 3184 // Note: "allows_artificial_cells" is currently enforced for 3185 // MG computations. 3186 std::vector<types::subdomain_id> saved_level_subdomain_ids; 3187 saved_level_subdomain_ids.resize(tr->n_cells(lvl)); 3188 { 3189 typename dealii::parallel::shared::Triangulation<dim, spacedim>:: 3190 cell_iterator cell = 3191 this->dof_handler->get_triangulation().begin( 3192 lvl), 3193 endc = 3194 this->dof_handler->get_triangulation().end(lvl); 3195 3196 const std::vector<types::subdomain_id> &true_level_subdomain_ids = 3197 tr->get_true_level_subdomain_ids_of_cells(lvl); 3198 3199 for (unsigned int index = 0; cell != endc; ++cell, ++index) 3200 { 3201 saved_level_subdomain_ids[index] = cell->level_subdomain_id(); 3202 cell->set_level_subdomain_id(true_level_subdomain_ids[index]); 3203 } 3204 } 3205 3206 // Next let the sequential algorithm do its magic. it is going to 3207 // enumerate DoFs on all cells on the given level, regardless of 3208 // owner 3209 const types::global_dof_index n_dofs_on_level = 3210 Implementation::distribute_dofs_on_level( 3211 numbers::invalid_subdomain_id, *this->dof_handler, lvl); 3212 3213 // then re-enumerate them based on their level subdomain 3214 // association. for this, we first have to identify for each current 3215 // DoF index which subdomain they belong to. ideally, we would like 3216 // to call DoFRenumbering::subdomain_wise(), but because the 3217 // NumberCache of the current DoFHandler is not fully set up yet, we 3218 // can't quite do that. also, that function has to deal with other 3219 // kinds of triangulations as well, whereas we here know what kind 3220 // of triangulation we have and can simplify the code accordingly 3221 std::vector<types::global_dof_index> new_dof_indices( 3222 n_dofs_on_level, numbers::invalid_dof_index); 3223 { 3224 // first get the association of each dof with a subdomain and 3225 // determine the total number of subdomain ids used 3226 const std::vector<types::subdomain_id> 3227 level_subdomain_association = 3228 get_dof_level_subdomain_association(*this->dof_handler, 3229 n_dofs_on_level, 3230 n_procs, 3231 lvl); 3232 3233 // then renumber the subdomains by first looking at those 3234 // belonging to subdomain 0, then those of subdomain 1, etc. note 3235 // that the algorithm is stable, i.e. if two dofs i,j have i<j and 3236 // belong to the same subdomain, then they will be in this order 3237 // also after reordering 3238 types::global_dof_index next_free_index = 0; 3239 for (types::subdomain_id level_subdomain = 0; 3240 level_subdomain < n_procs; 3241 ++level_subdomain) 3242 for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i) 3243 if (level_subdomain_association[i] == level_subdomain) 3244 { 3245 Assert(new_dof_indices[i] == numbers::invalid_dof_index, 3246 ExcInternalError()); 3247 new_dof_indices[i] = next_free_index; 3248 ++next_free_index; 3249 } 3250 3251 // we should have numbered all dofs 3252 Assert(next_free_index == n_dofs_on_level, ExcInternalError()); 3253 Assert(std::find(new_dof_indices.begin(), 3254 new_dof_indices.end(), 3255 numbers::invalid_dof_index) == 3256 new_dof_indices.end(), 3257 ExcInternalError()); 3258 } 3259 3260 // finally do the renumbering. we can use the sequential 3261 // version of the function because we do things on all 3262 // cells and all cells have their subdomain ids and DoFs 3263 // correctly set 3264 Implementation::renumber_mg_dofs( 3265 new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true); 3266 3267 // update the number cache. for this, we first have to find the 3268 // level subdomain association for each DoF again following 3269 // renumbering, from which we can then compute the IndexSets of 3270 // locally owned DoFs for all processors. all other fields then 3271 // follow from this 3272 // 3273 // given the way we enumerate degrees of freedom, the locally owned 3274 // ranges must all be contiguous and consecutive. this makes filling 3275 // the IndexSets cheap. an assertion at the top verifies that this 3276 // assumption is true 3277 const std::vector<types::subdomain_id> level_subdomain_association = 3278 get_dof_level_subdomain_association(*this->dof_handler, 3279 n_dofs_on_level, 3280 n_procs, 3281 lvl); 3282 3283 for (unsigned int i = 1; i < n_dofs_on_level; ++i) 3284 Assert(level_subdomain_association[i] >= 3285 level_subdomain_association[i - 1], 3286 ExcInternalError()); 3287 3288 std::vector<IndexSet> locally_owned_dofs_per_processor( 3289 n_procs, IndexSet(n_dofs_on_level)); 3290 { 3291 // we know that the set of subdomain indices is contiguous from 3292 // the assertion above; find the start and end index for each 3293 // processor, taking into account that sometimes a processor 3294 // may not in fact have any DoFs at all. we do the latter 3295 // by just identifying contiguous ranges of level_subdomain_ids 3296 // and filling IndexSets for those subdomains; subdomains 3297 // that don't appear will lead to IndexSets that are simply 3298 // never touched and remain empty as initialized above. 3299 unsigned int start_index = 0; 3300 unsigned int end_index = 0; 3301 while (start_index < n_dofs_on_level) 3302 { 3303 while ((end_index) < n_dofs_on_level && 3304 (level_subdomain_association[end_index] == 3305 level_subdomain_association[start_index])) 3306 ++end_index; 3307 3308 // we've now identified a range of same indices. set that 3309 // range in the corresponding IndexSet 3310 if (end_index > start_index) 3311 { 3312 const unsigned int level_subdomain_owner = 3313 level_subdomain_association[start_index]; 3314 locally_owned_dofs_per_processor[level_subdomain_owner] 3315 .add_range(start_index, end_index); 3316 } 3317 3318 // then move on to thinking about the next range 3319 start_index = end_index; 3320 } 3321 } 3322 3323 // finally, restore current level subdomain ids 3324 { 3325 typename dealii::parallel::shared::Triangulation<dim, spacedim>:: 3326 cell_iterator cell = 3327 this->dof_handler->get_triangulation().begin( 3328 lvl), 3329 endc = 3330 this->dof_handler->get_triangulation().end(lvl); 3331 3332 for (unsigned int index = 0; cell != endc; ++cell, ++index) 3333 cell->set_level_subdomain_id(saved_level_subdomain_ids[index]); 3334 3335 // add NumberCache for current level 3336 number_caches.emplace_back( 3337 NumberCache(locally_owned_dofs_per_processor, 3338 this->dof_handler->get_triangulation() 3339 .locally_owned_subdomain())); 3340 } 3341 } 3342 3343 return number_caches; 3344 } 3345 3346 3347 3348 template <int dim, int spacedim> 3349 NumberCache renumber_dofs(const std::vector<types::global_dof_index> & new_numbers) const3350 ParallelShared<dim, spacedim>::renumber_dofs( 3351 const std::vector<types::global_dof_index> &new_numbers) const 3352 { 3353 #ifndef DEAL_II_WITH_MPI 3354 (void)new_numbers; 3355 Assert(false, ExcNotImplemented()); 3356 return NumberCache(); 3357 #else 3358 // Similar to distribute_dofs() we need to have a special treatment in 3359 // case artificial cells are present. 3360 const dealii::parallel::shared::Triangulation<dim, spacedim> *tr = 3361 (dynamic_cast< 3362 const dealii::parallel::shared::Triangulation<dim, spacedim> *>( 3363 &this->dof_handler->get_triangulation())); 3364 Assert(tr != nullptr, ExcInternalError()); 3365 3366 typename dealii::parallel::shared::Triangulation<dim, spacedim>:: 3367 active_cell_iterator 3368 cell = this->dof_handler->get_triangulation().begin_active(), 3369 endc = this->dof_handler->get_triangulation().end(); 3370 std::vector<types::subdomain_id> current_subdomain_ids( 3371 tr->n_active_cells()); 3372 const std::vector<types::subdomain_id> &true_subdomain_ids = 3373 tr->get_true_subdomain_ids_of_cells(); 3374 if (tr->with_artificial_cells()) 3375 for (unsigned int index = 0; cell != endc; cell++, index++) 3376 { 3377 current_subdomain_ids[index] = cell->subdomain_id(); 3378 cell->set_subdomain_id(true_subdomain_ids[index]); 3379 } 3380 3381 std::vector<types::global_dof_index> global_gathered_numbers( 3382 this->dof_handler->n_dofs(), 0); 3383 // as we call DoFRenumbering::subdomain_wise (*dof_handler) from 3384 // distribute_dofs(), we need to support sequential-like input. 3385 // Distributed-like input from, for example, component_wise renumbering 3386 // is also supported. 3387 if (new_numbers.size() == this->dof_handler->n_dofs()) 3388 { 3389 global_gathered_numbers = new_numbers; 3390 } 3391 else 3392 { 3393 Assert(new_numbers.size() == 3394 this->dof_handler->locally_owned_dofs().n_elements(), 3395 ExcInternalError()); 3396 const unsigned int n_cpu = 3397 Utilities::MPI::n_mpi_processes(tr->get_communicator()); 3398 std::vector<types::global_dof_index> gathered_new_numbers( 3399 this->dof_handler->n_dofs(), 0); 3400 Assert(Utilities::MPI::this_mpi_process(tr->get_communicator()) == 3401 this->dof_handler->get_triangulation() 3402 .locally_owned_subdomain(), 3403 ExcInternalError()) 3404 3405 // gather new numbers among processors into one vector 3406 { 3407 std::vector<types::global_dof_index> new_numbers_copy( 3408 new_numbers); 3409 3410 // store the number of elements that are to be received from each 3411 // process 3412 std::vector<int> rcounts(n_cpu); 3413 3414 types::global_dof_index shift = 0; 3415 // set rcounts based on new_numbers: 3416 int cur_count = new_numbers_copy.size(); 3417 int ierr = MPI_Allgather(&cur_count, 3418 1, 3419 MPI_INT, 3420 rcounts.data(), 3421 1, 3422 MPI_INT, 3423 tr->get_communicator()); 3424 AssertThrowMPI(ierr); 3425 3426 // compute the displacements (relative to recvbuf) 3427 // at which to place the incoming data from process i 3428 std::vector<int> displacements(n_cpu); 3429 for (unsigned int i = 0; i < n_cpu; i++) 3430 { 3431 displacements[i] = shift; 3432 shift += rcounts[i]; 3433 } 3434 Assert(new_numbers_copy.size() == 3435 static_cast<unsigned int>( 3436 rcounts[Utilities::MPI::this_mpi_process( 3437 tr->get_communicator())]), 3438 ExcInternalError()); 3439 ierr = MPI_Allgatherv(new_numbers_copy.data(), 3440 new_numbers_copy.size(), 3441 DEAL_II_DOF_INDEX_MPI_TYPE, 3442 gathered_new_numbers.data(), 3443 rcounts.data(), 3444 displacements.data(), 3445 DEAL_II_DOF_INDEX_MPI_TYPE, 3446 tr->get_communicator()); 3447 AssertThrowMPI(ierr); 3448 } 3449 3450 // put new numbers according to the current 3451 // locally_owned_dofs_per_processor IndexSets 3452 types::global_dof_index shift = 0; 3453 // flag_1 and flag_2 are 3454 // used to control that there is a 3455 // one-to-one relation between old and new DoFs. 3456 std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0); 3457 std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0); 3458 for (unsigned int i = 0; i < n_cpu; i++) 3459 { 3460 const IndexSet iset = 3461 this->dof_handler->locally_owned_dofs_per_processor()[i]; 3462 for (types::global_dof_index ind = 0; ind < iset.n_elements(); 3463 ind++) 3464 { 3465 const types::global_dof_index target = 3466 iset.nth_index_in_set(ind); 3467 const types::global_dof_index value = 3468 gathered_new_numbers[shift + ind]; 3469 Assert(target < this->dof_handler->n_dofs(), 3470 ExcInternalError()); 3471 Assert(value < this->dof_handler->n_dofs(), 3472 ExcInternalError()); 3473 global_gathered_numbers[target] = value; 3474 flag_1[target]++; 3475 flag_2[value]++; 3476 } 3477 shift += iset.n_elements(); 3478 } 3479 3480 Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1, 3481 ExcInternalError()); 3482 Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1, 3483 ExcInternalError()); 3484 Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1, 3485 ExcInternalError()); 3486 Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1, 3487 ExcInternalError()); 3488 } 3489 3490 // let the sequential algorithm do its magic; ignore the 3491 // return type, but reconstruct the number cache based on 3492 // which DoFs each process owns 3493 Implementation::renumber_dofs(global_gathered_numbers, 3494 IndexSet(0), 3495 *this->dof_handler, 3496 /*check_validity=*/true); 3497 3498 const NumberCache number_cache( 3499 DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler), 3500 this->dof_handler->get_triangulation().locally_owned_subdomain()); 3501 3502 // restore artificial cells 3503 cell = tr->begin_active(); 3504 if (tr->with_artificial_cells()) 3505 for (unsigned int index = 0; cell != endc; cell++, index++) 3506 cell->set_subdomain_id(current_subdomain_ids[index]); 3507 3508 return number_cache; 3509 #endif 3510 } 3511 3512 3513 3514 template <int dim, int spacedim> 3515 NumberCache renumber_mg_dofs(const unsigned int,const std::vector<types::global_dof_index> &) const3516 ParallelShared<dim, spacedim>::renumber_mg_dofs( 3517 const unsigned int /*level*/, 3518 const std::vector<types::global_dof_index> & /*new_numbers*/) const 3519 { 3520 // multigrid is not currently implemented for shared triangulations 3521 Assert(false, ExcNotImplemented()); 3522 3523 return {}; 3524 } 3525 3526 3527 3528 /* --------------------- class ParallelDistributed ---------------- */ 3529 3530 #ifdef DEAL_II_WITH_MPI 3531 3532 namespace 3533 { 3534 template <int dim, int spacedim> 3535 void communicate_mg_ghost_cells(const typename dealii::parallel::DistributedTriangulationBase<dim,spacedim> & tria,DoFHandler<dim,spacedim> & dof_handler)3536 communicate_mg_ghost_cells( 3537 const typename dealii::parallel:: 3538 DistributedTriangulationBase<dim, spacedim> &tria, 3539 DoFHandler<dim, spacedim> & dof_handler) 3540 { 3541 (void)tria; 3542 const auto pack = [&](const auto &cell) { 3543 // why would somebody request a cell that is not ours? 3544 Assert(cell->level_subdomain_id() == tria.locally_owned_subdomain(), 3545 ExcInternalError()); 3546 3547 std::vector<dealii::types::global_dof_index> data( 3548 cell->get_fe().n_dofs_per_cell()); 3549 cell->get_mg_dof_indices(data); 3550 3551 return data; 3552 }; 3553 3554 const auto unpack = [](const auto &cell, const auto &dofs) { 3555 Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(), 3556 ExcInternalError()); 3557 3558 Assert(cell->level_subdomain_id() != 3559 dealii::numbers::artificial_subdomain_id, 3560 ExcInternalError()); 3561 3562 std::vector<dealii::types::global_dof_index> dof_indices( 3563 cell->get_fe().n_dofs_per_cell()); 3564 cell->get_mg_dof_indices(dof_indices); 3565 3566 bool complete = true; 3567 for (unsigned int i = 0; i < dof_indices.size(); ++i) 3568 if (dofs[i] != numbers::invalid_dof_index) 3569 { 3570 Assert((dof_indices[i] == (numbers::invalid_dof_index)) || 3571 (dof_indices[i] == dofs[i]), 3572 ExcInternalError()); 3573 dof_indices[i] = dofs[i]; 3574 } 3575 else 3576 complete = false; 3577 3578 if (!complete) 3579 const_cast< 3580 typename DoFHandler<dim, spacedim>::level_cell_iterator &>(cell) 3581 ->set_user_flag(); 3582 else 3583 const_cast< 3584 typename DoFHandler<dim, spacedim>::level_cell_iterator &>(cell) 3585 ->clear_user_flag(); 3586 3587 const_cast< 3588 typename DoFHandler<dim, spacedim>::level_cell_iterator &>(cell) 3589 ->set_mg_dof_indices(dof_indices); 3590 }; 3591 3592 const auto filter = [](const auto &cell) { 3593 return cell->user_flag_set(); 3594 }; 3595 3596 GridTools::exchange_cell_data_to_level_ghosts< 3597 std::vector<types::global_dof_index>, 3598 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter); 3599 } 3600 3601 3602 3603 template <int spacedim> 3604 void communicate_mg_ghost_cells(const typename dealii::parallel::distributed::Triangulation<1,spacedim> &,DoFHandler<1,spacedim> &)3605 communicate_mg_ghost_cells(const typename dealii::parallel:: 3606 distributed::Triangulation<1, spacedim> &, 3607 DoFHandler<1, spacedim> &) 3608 { 3609 Assert(false, ExcNotImplemented()); 3610 } 3611 3612 3613 3614 /** 3615 * A function that communicates the DoF indices from that subset of 3616 * locally owned cells that have their user indices set to the 3617 * corresponding ghost cells on other processors. 3618 * 3619 * This function makes use of the user flags in the following 3620 * way: 3621 * - On locally owned cells, the flag indicates whether we still 3622 * need to send the DoF indices to other processors on which 3623 * the current cell is a ghost. In phase 1, this is true for 3624 * all locally owned cells that are adjacent to ghost cells 3625 * in some way. In phase 2, this is only true if before phase 3626 * 1 we did not know all dof indices yet 3627 * - On ghost cells, the flag indicates whether we still expect 3628 * information to be sent to us. In phase 1, this is true for 3629 * all ghost cells. In phase 2, this is only true if we 3630 * did not receive a complete set of DoF indices in phase 1. 3631 */ 3632 template <int dim, int spacedim> 3633 void communicate_dof_indices_on_marked_cells(const DoFHandler<dim,spacedim> & dof_handler)3634 communicate_dof_indices_on_marked_cells( 3635 const DoFHandler<dim, spacedim> &dof_handler) 3636 { 3637 # ifndef DEAL_II_WITH_MPI 3638 (void)dof_handler; 3639 Assert(false, ExcNotImplemented()); 3640 # else 3641 3642 // define functions that pack data on cells that are ghost cells 3643 // somewhere else, and unpack data on cells where we get information 3644 // from elsewhere 3645 const auto pack = [](const auto &cell) { 3646 Assert(cell->is_locally_owned(), ExcInternalError()); 3647 3648 std::vector<dealii::types::global_dof_index> data( 3649 cell->get_fe().n_dofs_per_cell()); 3650 cell->get_dof_indices(data); 3651 3652 return data; 3653 }; 3654 3655 const auto unpack = [](const auto &cell, const auto &dofs) { 3656 Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(), 3657 ExcInternalError()); 3658 3659 Assert(cell->is_ghost(), ExcInternalError()); 3660 3661 std::vector<dealii::types::global_dof_index> dof_indices( 3662 cell->get_fe().n_dofs_per_cell()); 3663 cell->update_cell_dof_indices_cache(); 3664 cell->get_dof_indices(dof_indices); 3665 3666 bool complete = true; 3667 for (unsigned int i = 0; i < dof_indices.size(); ++i) 3668 if (dofs[i] != numbers::invalid_dof_index) 3669 { 3670 Assert((dof_indices[i] == (numbers::invalid_dof_index)) || 3671 (dof_indices[i] == dofs[i]), 3672 ExcInternalError()); 3673 dof_indices[i] = dofs[i]; 3674 } 3675 else 3676 complete = false; 3677 3678 if (!complete) 3679 const_cast< 3680 typename DoFHandler<dim, spacedim>::active_cell_iterator &>( 3681 cell) 3682 ->set_user_flag(); 3683 else 3684 const_cast< 3685 typename DoFHandler<dim, spacedim>::active_cell_iterator &>( 3686 cell) 3687 ->clear_user_flag(); 3688 3689 const_cast< 3690 typename DoFHandler<dim, spacedim>::active_cell_iterator &>(cell) 3691 ->set_dof_indices(dof_indices); 3692 }; 3693 3694 const auto filter = [](const auto &cell) { 3695 return cell->user_flag_set(); 3696 }; 3697 3698 GridTools::exchange_cell_data_to_ghosts< 3699 std::vector<types::global_dof_index>, 3700 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter); 3701 3702 // finally update the cell DoF indices caches to make sure 3703 // our internal data structures are consistent 3704 update_all_active_cell_dof_indices_caches(dof_handler); 3705 3706 3707 // have a barrier so that sends between two calls to this 3708 // function are not mixed up. 3709 // 3710 // this is necessary because above we just see if there are 3711 // messages and then receive them, without discriminating 3712 // where they come from and whether they were sent in phase 3713 // 1 or 2 (the function is called twice in a row). the need 3714 // for a global communication step like this barrier could 3715 // be avoided by receiving messages specifically from those 3716 // processors from which we expect messages, and by using 3717 // different tags for phase 1 and 2, but the cost of a 3718 // barrier is negligible compared to everything else we do 3719 // here 3720 if (const auto *triangulation = 3721 dynamic_cast<const dealii::parallel:: 3722 DistributedTriangulationBase<dim, spacedim> *>( 3723 &dof_handler.get_triangulation())) 3724 { 3725 const int ierr = MPI_Barrier(triangulation->get_communicator()); 3726 AssertThrowMPI(ierr); 3727 } 3728 else 3729 { 3730 Assert(false, 3731 ExcMessage( 3732 "The function communicate_dof_indices_on_marked_cells() " 3733 "only works with parallel distributed triangulations.")); 3734 } 3735 # endif 3736 } 3737 3738 3739 3740 } // namespace 3741 3742 #endif // DEAL_II_WITH_MPI 3743 3744 3745 3746 template <int dim, int spacedim> ParallelDistributed(DoFHandler<dim,spacedim> & dof_handler)3747 ParallelDistributed<dim, spacedim>::ParallelDistributed( 3748 DoFHandler<dim, spacedim> &dof_handler) 3749 : dof_handler(&dof_handler) 3750 {} 3751 3752 3753 3754 template <int dim, int spacedim> 3755 NumberCache distribute_dofs() const3756 ParallelDistributed<dim, spacedim>::distribute_dofs() const 3757 { 3758 #ifndef DEAL_II_WITH_MPI 3759 Assert(false, ExcNotImplemented()); 3760 return NumberCache(); 3761 #else 3762 3763 dealii::parallel::DistributedTriangulationBase<dim, spacedim> 3764 *triangulation = 3765 (dynamic_cast< 3766 dealii::parallel::DistributedTriangulationBase<dim, spacedim> *>( 3767 const_cast<dealii::Triangulation<dim, spacedim> *>( 3768 &dof_handler->get_triangulation()))); 3769 Assert(triangulation != nullptr, ExcInternalError()); 3770 3771 const types::subdomain_id subdomain_id = 3772 triangulation->locally_owned_subdomain(); 3773 3774 3775 /* 3776 The following algorithm has a number of stages that are all 3777 documented in the paper that describes the parallel::distributed 3778 functionality: 3779 3780 1/ locally enumerate dofs on locally owned cells 3781 2/ eliminate dof duplicates on all cells. 3782 un-numerate those that are on interfaces with ghost 3783 cells and that we don't own based on the tie-breaking 3784 criterion. unify dofs afterwards. 3785 3/ unify dofs and re-enumerate the remaining valid ones. 3786 the end result is that we only enumerate locally owned 3787 DoFs 3788 4/ shift indices so that each processor has a unique 3789 range of indices 3790 5/ for all locally owned cells that are ghost 3791 cells somewhere else, send our own DoF indices 3792 to the appropriate set of other processors. 3793 overwrite invalid DoF indices on ghost interfaces 3794 with the corresponding valid ones that we now know. 3795 6/ send DoF indices again to get the correct indices 3796 on ghost cells that we may not have known earlier 3797 */ 3798 3799 // --------- Phase 1: enumerate dofs on locally owned cells 3800 const types::global_dof_index n_initial_local_dofs = 3801 Implementation::distribute_dofs(subdomain_id, *dof_handler); 3802 3803 // --------- Phase 2: eliminate dof duplicates on all cells: 3804 // - un-numerate dofs on interfaces to ghost cells 3805 // that we don't own 3806 // - in case of hp support, unify dofs 3807 std::vector<dealii::types::global_dof_index> renumbering( 3808 n_initial_local_dofs, enumeration_dof_index); 3809 3810 // first, we invalidate degrees of freedom that belong to processors 3811 // of a lower rank, from which we will receive the final (and lower) 3812 // degrees of freedom later. 3813 Implementation:: 3814 invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering( 3815 renumbering, subdomain_id, *dof_handler); 3816 3817 // then, we identify DoF duplicates if the DoFHandler has hp 3818 // capabilities 3819 std::vector<std::map<types::global_dof_index, types::global_dof_index>> 3820 all_constrained_indices(dim); 3821 Implementation::compute_dof_identities(all_constrained_indices, 3822 *dof_handler); 3823 3824 // --------- Phase 3: re-enumerate the valid degrees of freedom 3825 // consecutively. thus, we finally receive the 3826 // correct number of locally owned DoFs after 3827 // this step. 3828 // 3829 // the order in which we handle Phases 2 and 3 is important, 3830 // since we want to clarify ownership of degrees of freedom before 3831 // we actually unify and enumerate their indices. otherwise, we could 3832 // end up having a degee of freedom to which only invalid indices will 3833 // be assigned. 3834 const types::global_dof_index n_locally_owned_dofs = 3835 Implementation::enumerate_dof_indices_for_renumbering( 3836 renumbering, all_constrained_indices, *dof_handler); 3837 3838 // --------- Phase 4: shift indices so that each processor has a unique 3839 // range of indices 3840 dealii::types::global_dof_index my_shift = 0; 3841 const int ierr = 3842 MPI_Exscan(DEAL_II_MPI_CONST_CAST(&n_locally_owned_dofs), 3843 &my_shift, 3844 1, 3845 DEAL_II_DOF_INDEX_MPI_TYPE, 3846 MPI_SUM, 3847 triangulation->get_communicator()); 3848 AssertThrowMPI(ierr); 3849 3850 // make dof indices globally consecutive 3851 for (auto &new_index : renumbering) 3852 if (new_index != numbers::invalid_dof_index) 3853 new_index += my_shift; 3854 3855 // now re-enumerate all dofs to this shifted and condensed 3856 // numbering form. we renumber some dofs as invalid, so 3857 // choose the nocheck-version. 3858 Implementation::renumber_dofs(renumbering, 3859 IndexSet(0), 3860 *dof_handler, 3861 /*check_validity=*/false); 3862 3863 // now a little bit of housekeeping 3864 const dealii::types::global_dof_index n_global_dofs = 3865 Utilities::MPI::sum(n_locally_owned_dofs, 3866 triangulation->get_communicator()); 3867 3868 NumberCache number_cache; 3869 number_cache.n_global_dofs = n_global_dofs; 3870 number_cache.n_locally_owned_dofs = n_locally_owned_dofs; 3871 number_cache.locally_owned_dofs = IndexSet(n_global_dofs); 3872 number_cache.locally_owned_dofs.add_range(my_shift, 3873 my_shift + 3874 n_locally_owned_dofs); 3875 number_cache.locally_owned_dofs.compress(); 3876 3877 // this ends the phase where we enumerate degrees of freedom on 3878 // each processor. what is missing is communicating DoF indices 3879 // on ghost cells 3880 3881 // --------- Phase 5: for all locally owned cells that are ghost 3882 // cells somewhere else, send our own DoF indices 3883 // to the appropriate set of other processors 3884 { 3885 std::vector<bool> user_flags; 3886 triangulation->save_user_flags(user_flags); 3887 triangulation->clear_user_flags(); 3888 3889 // mark all cells that either have to send data (locally 3890 // owned cells that are adjacent to ghost neighbors in some 3891 // way) or receive data (all ghost cells) via the user flags 3892 for (const auto &cell : dof_handler->active_cell_iterators()) 3893 if (cell->is_ghost()) 3894 cell->set_user_flag(); 3895 3896 // Send and receive cells. After this, only the local cells 3897 // are marked, that received new data. This has to be 3898 // communicated in a second communication step. 3899 // 3900 // as explained in the 'distributed' paper, this has to be 3901 // done twice 3902 communicate_dof_indices_on_marked_cells(*dof_handler); 3903 3904 // If the DoFHandler has hp capabilities enabled, then we may have 3905 // received valid indices of degrees of freedom that are dominated 3906 // by a fe object adjacent to a ghost interface. thus, we overwrite 3907 // the remaining invalid indices with the valid ones in this step. 3908 Implementation::merge_invalid_dof_indices_on_ghost_interfaces( 3909 *dof_handler); 3910 3911 // --------- Phase 6: all locally owned cells have their correct 3912 // DoF indices set. however, some ghost cells 3913 // may still have invalid ones. thus, exchange 3914 // one more time. 3915 communicate_dof_indices_on_marked_cells(*dof_handler); 3916 3917 // at this point, we must have taken care of the data transfer 3918 // on all cells we had previously marked. verify this 3919 # ifdef DEBUG 3920 for (const auto &cell : dof_handler->active_cell_iterators()) 3921 Assert(cell->user_flag_set() == false, ExcInternalError()); 3922 # endif 3923 3924 triangulation->load_user_flags(user_flags); 3925 } 3926 3927 # ifdef DEBUG 3928 // check that we are really done 3929 { 3930 std::vector<dealii::types::global_dof_index> local_dof_indices; 3931 3932 for (const auto &cell : dof_handler->active_cell_iterators()) 3933 if (!cell->is_artificial()) 3934 { 3935 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 3936 cell->get_dof_indices(local_dof_indices); 3937 if (local_dof_indices.end() != 3938 std::find(local_dof_indices.begin(), 3939 local_dof_indices.end(), 3940 numbers::invalid_dof_index)) 3941 { 3942 if (cell->is_ghost()) 3943 { 3944 Assert(false, 3945 ExcMessage( 3946 "A ghost cell ended up with incomplete " 3947 "DoF index information. This should not " 3948 "have happened!")); 3949 } 3950 else 3951 { 3952 Assert( 3953 false, 3954 ExcMessage( 3955 "A locally owned cell ended up with incomplete " 3956 "DoF index information. This should not " 3957 "have happened!")); 3958 } 3959 } 3960 } 3961 } 3962 # endif // DEBUG 3963 return number_cache; 3964 #endif // DEAL_II_WITH_MPI 3965 } 3966 3967 3968 3969 template <int dim, int spacedim> 3970 std::vector<NumberCache> distribute_mg_dofs() const3971 ParallelDistributed<dim, spacedim>::distribute_mg_dofs() const 3972 { 3973 #ifndef DEAL_II_WITH_MPI 3974 Assert(false, ExcNotImplemented()); 3975 return std::vector<NumberCache>(); 3976 #else 3977 3978 dealii::parallel::DistributedTriangulationBase<dim, spacedim> 3979 *triangulation = 3980 (dynamic_cast< 3981 dealii::parallel::DistributedTriangulationBase<dim, spacedim> *>( 3982 const_cast<dealii::Triangulation<dim, spacedim> *>( 3983 &dof_handler->get_triangulation()))); 3984 Assert(triangulation != nullptr, ExcInternalError()); 3985 3986 AssertThrow((triangulation->is_multilevel_hierarchy_constructed()), 3987 ExcMessage( 3988 "Multigrid DoFs can only be distributed on a parallel " 3989 "Triangulation if the flag construct_multigrid_hierarchy " 3990 "is set in the constructor.")); 3991 3992 // loop over all levels that exist globally (across all 3993 // processors), even if the current processor does not in fact 3994 // have any cells on that level or if the local part of the 3995 // Triangulation has fewer levels. we need to do this because 3996 // we need to communicate across all processors on all levels 3997 const unsigned int n_levels = triangulation->n_global_levels(); 3998 std::vector<NumberCache> number_caches; 3999 number_caches.reserve(n_levels); 4000 for (unsigned int level = 0; level < n_levels; ++level) 4001 { 4002 NumberCache level_number_cache; 4003 4004 //* 1. distribute on own subdomain 4005 const unsigned int n_initial_local_dofs = 4006 Implementation::distribute_dofs_on_level( 4007 triangulation->locally_owned_subdomain(), *dof_handler, level); 4008 4009 //* 2. iterate over ghostcells and kill dofs that are not 4010 // owned by us 4011 std::vector<dealii::types::global_dof_index> renumbering( 4012 n_initial_local_dofs); 4013 for (dealii::types::global_dof_index i = 0; i < renumbering.size(); 4014 ++i) 4015 renumbering[i] = i; 4016 4017 if (level < triangulation->n_levels()) 4018 { 4019 std::vector<dealii::types::global_dof_index> local_dof_indices; 4020 4021 for (const auto &cell : 4022 dof_handler->cell_iterators_on_level(level)) 4023 if (cell->level_subdomain_id() != 4024 numbers::artificial_subdomain_id && 4025 (cell->level_subdomain_id() < 4026 triangulation->locally_owned_subdomain())) 4027 { 4028 // we found a neighboring ghost cell whose 4029 // subdomain is "stronger" than our own 4030 // subdomain 4031 4032 // delete all dofs that live there and that we 4033 // have previously assigned a number to 4034 // (i.e. the ones on the interface) 4035 local_dof_indices.resize( 4036 cell->get_fe().n_dofs_per_cell()); 4037 cell->get_mg_dof_indices(local_dof_indices); 4038 for (unsigned int i = 0; 4039 i < cell->get_fe().n_dofs_per_cell(); 4040 ++i) 4041 if (local_dof_indices[i] != numbers::invalid_dof_index) 4042 renumbering[local_dof_indices[i]] = 4043 numbers::invalid_dof_index; 4044 } 4045 } 4046 4047 level_number_cache.n_locally_owned_dofs = 0; 4048 for (types::global_dof_index &index : renumbering) 4049 if (index != numbers::invalid_dof_index) 4050 index = level_number_cache.n_locally_owned_dofs++; 4051 4052 //* 3. communicate local dofcount and shift ids to make 4053 // them unique 4054 dealii::types::global_dof_index my_shift = 0; 4055 int ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST( 4056 &level_number_cache.n_locally_owned_dofs), 4057 &my_shift, 4058 1, 4059 DEAL_II_DOF_INDEX_MPI_TYPE, 4060 MPI_SUM, 4061 triangulation->get_communicator()); 4062 AssertThrowMPI(ierr); 4063 4064 // The last processor knows about the total number of dofs, so we 4065 // can use a cheaper broadcast rather than an MPI_Allreduce via 4066 // MPI::sum(). 4067 level_number_cache.n_global_dofs = 4068 my_shift + level_number_cache.n_locally_owned_dofs; 4069 ierr = MPI_Bcast(&level_number_cache.n_global_dofs, 4070 1, 4071 DEAL_II_DOF_INDEX_MPI_TYPE, 4072 Utilities::MPI::n_mpi_processes( 4073 triangulation->get_communicator()) - 4074 1, 4075 triangulation->get_communicator()); 4076 4077 // shift indices 4078 for (types::global_dof_index &index : renumbering) 4079 if (index != numbers::invalid_dof_index) 4080 index += my_shift; 4081 4082 // now re-enumerate all dofs to this shifted and condensed 4083 // numbering form. we renumber some dofs as invalid, so 4084 // choose the nocheck-version of the function 4085 // 4086 // of course there is nothing for us to renumber if the 4087 // level we are currently dealing with doesn't even exist 4088 // within the current triangulation, so skip renumbering 4089 // in that case 4090 if (level < triangulation->n_levels()) 4091 Implementation::renumber_mg_dofs( 4092 renumbering, IndexSet(0), *dof_handler, level, false); 4093 4094 // now a little bit of housekeeping 4095 level_number_cache.locally_owned_dofs = 4096 IndexSet(level_number_cache.n_global_dofs); 4097 level_number_cache.locally_owned_dofs.add_range( 4098 my_shift, my_shift + level_number_cache.n_locally_owned_dofs); 4099 level_number_cache.locally_owned_dofs.compress(); 4100 4101 number_caches.emplace_back(level_number_cache); 4102 } 4103 4104 4105 //* communicate ghost DoFs 4106 // We mark all ghost cells by setting the user_flag and then request 4107 // these cells from the corresponding owners. As this information 4108 // can be incomplete, 4109 { 4110 std::vector<bool> user_flags; 4111 triangulation->save_user_flags(user_flags); 4112 triangulation->clear_user_flags(); 4113 4114 // mark all ghost cells for transfer 4115 { 4116 for (const auto &cell : dof_handler->cell_iterators()) 4117 if (cell->level_subdomain_id() != 4118 dealii::numbers::artificial_subdomain_id && 4119 !cell->is_locally_owned_on_level()) 4120 cell->set_user_flag(); 4121 } 4122 4123 // Phase 1. Request all marked cells from corresponding owners. If we 4124 // managed to get every DoF, remove the user_flag, otherwise we 4125 // will request them again in the step below. 4126 communicate_mg_ghost_cells(*triangulation, *dof_handler); 4127 4128 // have a barrier so that sends from above and below this 4129 // place are not mixed up. 4130 // 4131 // this is necessary because above we just see if there are 4132 // messages and then receive them, without discriminating 4133 // where they come from and whether they were sent in phase 4134 // 1 or 2 in communicate_mg_ghost_cells() on another 4135 // processor. the need for a global communication step like 4136 // this barrier could be avoided by receiving messages 4137 // specifically from those processors from which we expect 4138 // messages, and by using different tags for phase 1 and 2, 4139 // but the cost of a barrier is negligible compared to 4140 // everything else we do here 4141 const int ierr = MPI_Barrier(triangulation->get_communicator()); 4142 AssertThrowMPI(ierr); 4143 4144 // Phase 2, only request the cells that were not completed 4145 // in Phase 1. 4146 communicate_mg_ghost_cells(*triangulation, *dof_handler); 4147 4148 # ifdef DEBUG 4149 // make sure we have removed all flags: 4150 { 4151 for (const auto &cell : dof_handler->cell_iterators()) 4152 if (cell->level_subdomain_id() != 4153 dealii::numbers::artificial_subdomain_id && 4154 !cell->is_locally_owned_on_level()) 4155 Assert(cell->user_flag_set() == false, ExcInternalError()); 4156 } 4157 # endif 4158 4159 triangulation->load_user_flags(user_flags); 4160 } 4161 4162 4163 4164 # ifdef DEBUG 4165 // check that we are really done 4166 { 4167 std::vector<dealii::types::global_dof_index> local_dof_indices; 4168 for (const auto &cell : dof_handler->cell_iterators()) 4169 if (cell->level_subdomain_id() != 4170 dealii::numbers::artificial_subdomain_id) 4171 { 4172 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 4173 cell->get_mg_dof_indices(local_dof_indices); 4174 if (local_dof_indices.end() != 4175 std::find(local_dof_indices.begin(), 4176 local_dof_indices.end(), 4177 numbers::invalid_dof_index)) 4178 { 4179 Assert(false, ExcMessage("not all DoFs got distributed!")); 4180 } 4181 } 4182 } 4183 # endif // DEBUG 4184 4185 return number_caches; 4186 4187 #endif // DEAL_II_WITH_MPI 4188 } 4189 4190 4191 template <int dim, int spacedim> 4192 NumberCache renumber_dofs(const std::vector<dealii::types::global_dof_index> & new_numbers) const4193 ParallelDistributed<dim, spacedim>::renumber_dofs( 4194 const std::vector<dealii::types::global_dof_index> &new_numbers) const 4195 { 4196 (void)new_numbers; 4197 4198 Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(), 4199 ExcInternalError()); 4200 4201 #ifndef DEAL_II_WITH_MPI 4202 Assert(false, ExcNotImplemented()); 4203 return NumberCache(); 4204 #else 4205 4206 dealii::parallel::DistributedTriangulationBase<dim, spacedim> 4207 *triangulation = 4208 (dynamic_cast< 4209 dealii::parallel::DistributedTriangulationBase<dim, spacedim> *>( 4210 const_cast<dealii::Triangulation<dim, spacedim> *>( 4211 &dof_handler->get_triangulation()))); 4212 Assert(triangulation != nullptr, ExcInternalError()); 4213 4214 4215 // We start by checking whether only the numbering within the MPI 4216 // ranks changed. In that case, we can apply the renumbering with some 4217 // local renumbering only (this is similar to the renumber_mg_dofs() 4218 // function below) 4219 const bool locally_owned_set_changes = 4220 std::any_of(new_numbers.cbegin(), 4221 new_numbers.cend(), 4222 [this](const types::global_dof_index i) { 4223 return dof_handler->locally_owned_dofs().is_element( 4224 i) == false; 4225 }); 4226 4227 if (Utilities::MPI::sum(static_cast<unsigned int>( 4228 locally_owned_set_changes), 4229 triangulation->get_communicator()) == 0) 4230 { 4231 // Since only the order within the local subdomains has changed, 4232 // all we need to do is to propagate the knowledge about the 4233 // numbers from the locally owned dofs (given by the new_numbers 4234 // array) to all ghosted dofs on neighboring processors. We can do 4235 // this by ghost layer exchange routines as in parallel vectors: 4236 // We create an IndexSet for the relevant dofs and then export 4237 // into an array of those values via Utilities::MPI::Partitioner. 4238 IndexSet relevant_dofs; 4239 DoFTools::extract_locally_relevant_dofs(*dof_handler, 4240 relevant_dofs); 4241 std::vector<types::global_dof_index> ghosted_new_numbers( 4242 relevant_dofs.n_elements()); 4243 { 4244 Utilities::MPI::Partitioner partitioner( 4245 dof_handler->locally_owned_dofs(), 4246 relevant_dofs, 4247 triangulation->get_communicator()); 4248 4249 // choose some number that makes it unlikely to get conflicts 4250 // with other ongoing non-blocking communication (there 4251 // shouldn't be any at this place in most programs). 4252 const unsigned int communication_channel = 19; 4253 std::vector<types::global_dof_index> temp_array( 4254 partitioner.n_import_indices()); 4255 std::vector<MPI_Request> requests; 4256 partitioner.export_to_ghosted_array_start( 4257 communication_channel, 4258 make_array_view(new_numbers), 4259 make_array_view(temp_array), 4260 ArrayView<types::global_dof_index>( 4261 ghosted_new_numbers.data() + new_numbers.size(), 4262 partitioner.n_ghost_indices()), 4263 requests); 4264 partitioner.export_to_ghosted_array_finish( 4265 ArrayView<types::global_dof_index>( 4266 ghosted_new_numbers.data() + new_numbers.size(), 4267 partitioner.n_ghost_indices()), 4268 requests); 4269 4270 // we need to fill the indices of the locally owned part into 4271 // the new numbers array, which is not provided by the parallel 4272 // partitioner. their right position is somewhere in the middle 4273 // of the array, so we first copy the ghosted part from smaller 4274 // ranks to the front, then insert the data in the middle. 4275 unsigned int n_ghosts_on_smaller_ranks = 0; 4276 for (std::pair<unsigned int, unsigned int> t : 4277 partitioner.ghost_targets()) 4278 { 4279 if (t.first > partitioner.this_mpi_process()) 4280 break; 4281 n_ghosts_on_smaller_ranks += t.second; 4282 } 4283 if (n_ghosts_on_smaller_ranks > 0) 4284 { 4285 Assert(ghosted_new_numbers.data() != nullptr, 4286 ExcInternalError()); 4287 std::memmove(ghosted_new_numbers.data(), 4288 ghosted_new_numbers.data() + new_numbers.size(), 4289 sizeof(types::global_dof_index) * 4290 n_ghosts_on_smaller_ranks); 4291 } 4292 if (new_numbers.size() > 0) 4293 { 4294 Assert(new_numbers.data() != nullptr, ExcInternalError()); 4295 std::memcpy(ghosted_new_numbers.data() + 4296 n_ghosts_on_smaller_ranks, 4297 new_numbers.data(), 4298 sizeof(types::global_dof_index) * 4299 new_numbers.size()); 4300 } 4301 } 4302 4303 // In case we do not carry any relevant dof (but only some remote 4304 // processor), we do not need to call the renumbering. We call the 4305 // version without validity check because vertex dofs will be 4306 // set already in the artificial region. 4307 if (relevant_dofs.n_elements() > 0) 4308 Implementation::renumber_dofs(ghosted_new_numbers, 4309 relevant_dofs, 4310 *dof_handler, 4311 /*check_validity=*/false); 4312 4313 NumberCache number_cache; 4314 number_cache.locally_owned_dofs = dof_handler->locally_owned_dofs(); 4315 number_cache.n_global_dofs = dof_handler->n_dofs(); 4316 number_cache.n_locally_owned_dofs = 4317 number_cache.locally_owned_dofs.n_elements(); 4318 return number_cache; 4319 } 4320 else 4321 { 4322 // Now back to the more complicated case 4323 // 4324 // First figure out the new set of locally owned DoF indices. 4325 // If we own no DoFs, we still need to go through this function, 4326 // but we can skip this calculation. 4327 // 4328 // The IndexSet::add_indices() function is substantially more 4329 // efficient if the set of indices is already sorted because 4330 // it can then insert ranges instead of individual elements. 4331 // consequently, pre-sort the array of new indices 4332 IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs()); 4333 if (dof_handler->n_locally_owned_dofs() > 0) 4334 { 4335 std::vector<dealii::types::global_dof_index> 4336 new_numbers_sorted = new_numbers; 4337 std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end()); 4338 4339 my_locally_owned_new_dof_indices.add_indices( 4340 new_numbers_sorted.begin(), new_numbers_sorted.end()); 4341 my_locally_owned_new_dof_indices.compress(); 4342 4343 Assert(my_locally_owned_new_dof_indices.n_elements() == 4344 new_numbers.size(), 4345 ExcInternalError()); 4346 } 4347 4348 // delete all knowledge of DoF indices that are not locally 4349 // owned. we do so by getting DoF indices on cells, checking 4350 // whether they are locally owned, if not, setting them to 4351 // an invalid value, and then setting them again on the current 4352 // cell 4353 // 4354 // DoFs we (i) know about, and (ii) don't own locally must be 4355 // located either on ghost cells, or on the interface between a 4356 // locally owned cell and a ghost cell. In any case, it is 4357 // sufficient to kill them only from the ghost side cell, so loop 4358 // only over ghost cells 4359 { 4360 std::vector<dealii::types::global_dof_index> local_dof_indices; 4361 4362 for (auto cell : dof_handler->active_cell_iterators()) 4363 if (cell->is_ghost()) 4364 { 4365 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); 4366 cell->get_dof_indices(local_dof_indices); 4367 4368 for (unsigned int i = 0; 4369 i < cell->get_fe().n_dofs_per_cell(); 4370 ++i) 4371 // delete a DoF index if it has not already been deleted 4372 // (e.g., by visiting a neighboring cell, if it is on the 4373 // boundary), and if we don't own it 4374 if ((local_dof_indices[i] != 4375 numbers::invalid_dof_index) && 4376 (!dof_handler->locally_owned_dofs().is_element( 4377 local_dof_indices[i]))) 4378 local_dof_indices[i] = numbers::invalid_dof_index; 4379 4380 cell->set_dof_indices(local_dof_indices); 4381 } 4382 } 4383 4384 4385 // renumber. Skip when there is nothing to do because we own no DoF. 4386 if (dof_handler->locally_owned_dofs().n_elements() > 0) 4387 Implementation::renumber_dofs(new_numbers, 4388 dof_handler->locally_owned_dofs(), 4389 *dof_handler, 4390 /*check_validity=*/false); 4391 4392 // Communicate newly assigned DoF indices to other processors 4393 // and get the same information for our own ghost cells. 4394 // 4395 // This is the same as phase 5+6 in the distribute_dofs() algorithm, 4396 // taking into account that we have to unify a few DoFs in between 4397 // then communication phases if we do hp numbering 4398 { 4399 std::vector<bool> user_flags; 4400 triangulation->save_user_flags(user_flags); 4401 triangulation->clear_user_flags(); 4402 4403 // mark all own cells for transfer 4404 for (const auto &cell : dof_handler->active_cell_iterators()) 4405 if (cell->is_ghost()) 4406 cell->set_user_flag(); 4407 4408 4409 // Send and receive cells. After this, only the local cells 4410 // are marked, that received new data. This has to be 4411 // communicated in a second communication step. 4412 // 4413 // as explained in the 'distributed' paper, this has to be 4414 // done twice 4415 communicate_dof_indices_on_marked_cells(*dof_handler); 4416 4417 // if the DoFHandler has hp capabilities then we may have 4418 // received valid indices of degrees of freedom that are 4419 // dominated by a fe object adjacent to a ghost interface. 4420 // thus, we overwrite the remaining invalid indices with the 4421 // valid ones in this step. 4422 Implementation::merge_invalid_dof_indices_on_ghost_interfaces( 4423 *dof_handler); 4424 4425 communicate_dof_indices_on_marked_cells(*dof_handler); 4426 4427 triangulation->load_user_flags(user_flags); 4428 } 4429 4430 NumberCache number_cache; 4431 number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices; 4432 number_cache.n_global_dofs = dof_handler->n_dofs(); 4433 number_cache.n_locally_owned_dofs = 4434 number_cache.locally_owned_dofs.n_elements(); 4435 return number_cache; 4436 } 4437 #endif 4438 } 4439 4440 4441 4442 template <int dim, int spacedim> 4443 NumberCache renumber_mg_dofs(const unsigned int level,const std::vector<types::global_dof_index> & new_numbers) const4444 ParallelDistributed<dim, spacedim>::renumber_mg_dofs( 4445 const unsigned int level, 4446 const std::vector<types::global_dof_index> &new_numbers) const 4447 { 4448 // we only implement the case where the multigrid numbers are 4449 // renumbered within the processor's partition, rather than the most 4450 // general case 4451 const IndexSet index_set = dof_handler->locally_owned_mg_dofs(level); 4452 4453 #ifdef DEAL_II_WITH_MPI 4454 4455 const dealii::parallel::TriangulationBase<dim, spacedim> *tr = 4456 (dynamic_cast<const dealii::parallel::TriangulationBase<dim, spacedim> 4457 *>(&this->dof_handler->get_triangulation())); 4458 Assert(tr != nullptr, ExcInternalError()); 4459 4460 const unsigned int my_rank = 4461 Utilities::MPI::this_mpi_process(tr->get_communicator()); 4462 4463 # ifdef DEBUG 4464 for (types::global_dof_index i : new_numbers) 4465 { 4466 Assert(index_set.is_element(i), 4467 ExcNotImplemented( 4468 "Renumberings that change the locally owned mg dofs " 4469 "partitioning are currently not implemented for " 4470 "the multigrid levels")); 4471 } 4472 # endif 4473 4474 // we need to access all locally relevant degrees of freedom. we 4475 // use Utilities::MPI::Partitioner for handling the data exchange 4476 // of the new numbers, which is simply the extraction of ghost data 4477 IndexSet relevant_dofs; 4478 DoFTools::extract_locally_relevant_level_dofs(*dof_handler, 4479 level, 4480 relevant_dofs); 4481 std::vector<types::global_dof_index> ghosted_new_numbers( 4482 relevant_dofs.n_elements()); 4483 { 4484 Utilities::MPI::Partitioner partitioner(index_set, 4485 relevant_dofs, 4486 tr->get_communicator()); 4487 std::vector<types::global_dof_index> temp_array( 4488 partitioner.n_import_indices()); 4489 const unsigned int communication_channel = 17; 4490 std::vector<MPI_Request> requests; 4491 partitioner.export_to_ghosted_array_start( 4492 communication_channel, 4493 make_array_view(new_numbers), 4494 make_array_view(temp_array), 4495 ArrayView<types::global_dof_index>(ghosted_new_numbers.data() + 4496 new_numbers.size(), 4497 partitioner.n_ghost_indices()), 4498 requests); 4499 partitioner.export_to_ghosted_array_finish( 4500 ArrayView<types::global_dof_index>(ghosted_new_numbers.data() + 4501 new_numbers.size(), 4502 partitioner.n_ghost_indices()), 4503 requests); 4504 4505 // we need to fill the indices of the locally owned part into the 4506 // new numbers array. their right position is somewhere in the 4507 // middle of the array, so we first copy the ghosted part from 4508 // smaller ranks to the front, then insert the data in the middle. 4509 unsigned int n_ghosts_on_smaller_ranks = 0; 4510 for (std::pair<unsigned int, unsigned int> t : 4511 partitioner.ghost_targets()) 4512 { 4513 if (t.first > my_rank) 4514 break; 4515 n_ghosts_on_smaller_ranks += t.second; 4516 } 4517 if (n_ghosts_on_smaller_ranks > 0) 4518 { 4519 Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError()); 4520 std::memmove(ghosted_new_numbers.data(), 4521 ghosted_new_numbers.data() + new_numbers.size(), 4522 sizeof(types::global_dof_index) * 4523 n_ghosts_on_smaller_ranks); 4524 } 4525 if (new_numbers.size() > 0) 4526 { 4527 Assert(new_numbers.data() != nullptr, ExcInternalError()); 4528 std::memcpy(ghosted_new_numbers.data() + 4529 n_ghosts_on_smaller_ranks, 4530 new_numbers.data(), 4531 sizeof(types::global_dof_index) * new_numbers.size()); 4532 } 4533 } 4534 4535 // in case we do not own any of the given level (but only some remote 4536 // processor), we do not need to call the renumbering 4537 if (level < this->dof_handler->get_triangulation().n_levels() && 4538 relevant_dofs.n_elements() > 0) 4539 Implementation::renumber_mg_dofs( 4540 ghosted_new_numbers, relevant_dofs, *dof_handler, level, true); 4541 #else 4542 (void)new_numbers; 4543 Assert(false, ExcNotImplemented()); 4544 #endif 4545 4546 NumberCache number_cache; 4547 number_cache.locally_owned_dofs = index_set; 4548 number_cache.n_global_dofs = dof_handler->n_dofs(); 4549 number_cache.n_locally_owned_dofs = 4550 number_cache.locally_owned_dofs.n_elements(); 4551 return number_cache; 4552 } 4553 } // namespace Policy 4554 } // namespace DoFHandlerImplementation 4555 } // namespace internal 4556 4557 4558 4559 /*-------------- Explicit Instantiations -------------------------------*/ 4560 #include "dof_handler_policy.inst" 4561 4562 4563 DEAL_II_NAMESPACE_CLOSE 4564