1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2000 - 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 #ifndef dealii_fe_tools_extrapolate_templates_H 17 #define dealii_fe_tools_extrapolate_templates_H 18 19 20 #include <deal.II/base/config.h> 21 22 #include <deal.II/distributed/p4est_wrappers.h> 23 #include <deal.II/distributed/tria.h> 24 25 #include <deal.II/dofs/dof_handler.h> 26 #include <deal.II/dofs/dof_tools.h> 27 28 #include <deal.II/fe/fe_tools.h> 29 #include <deal.II/fe/fe_tools_interpolate.templates.h> 30 31 #include <deal.II/grid/tria.h> 32 33 #include <deal.II/lac/affine_constraints.h> 34 #include <deal.II/lac/block_vector.h> 35 #include <deal.II/lac/la_parallel_block_vector.h> 36 #include <deal.II/lac/la_parallel_vector.h> 37 #include <deal.II/lac/la_vector.h> 38 #include <deal.II/lac/petsc_block_vector.h> 39 #include <deal.II/lac/petsc_vector.h> 40 #include <deal.II/lac/trilinos_parallel_block_vector.h> 41 #include <deal.II/lac/trilinos_vector.h> 42 43 #include <queue> 44 45 DEAL_II_NAMESPACE_OPEN 46 47 namespace FETools 48 { 49 namespace internal 50 { 51 #ifndef DEAL_II_WITH_P4EST 52 // Dummy implementation in case p4est is not available. 53 template <int dim, int spacedim, class OutVector> 54 class ExtrapolateImplementation 55 { 56 public: ExtrapolateImplementation()57 ExtrapolateImplementation() 58 { 59 Assert(false, ExcNotImplemented()); 60 }; 61 62 template <class InVector> 63 void extrapolate_parallel(const InVector &,const DoFHandler<dim,spacedim> &,OutVector &)64 extrapolate_parallel(const InVector & /*u2_relevant*/, 65 const DoFHandler<dim, spacedim> & /*dof2*/, 66 OutVector & /*u2*/) 67 { 68 Assert(false, ExcNotImplemented()) 69 } 70 }; 71 #else 72 // Implementation of the @p extrapolate function 73 // on parallel distributed grids. 74 template <int dim, int spacedim, class OutVector> 75 class ExtrapolateImplementation 76 { 77 public: 78 ExtrapolateImplementation(); 79 80 template <class InVector> 81 void 82 extrapolate_parallel(const InVector & u2_relevant, 83 const DoFHandler<dim, spacedim> &dof2, 84 OutVector & u2); 85 86 private: 87 // A shortcut for the type of the OutVector 88 using value_type = typename OutVector::value_type; 89 90 // A structure holding all data to 91 // set dofs recursively on cells of arbitrary level 92 struct WorkPackage 93 { 94 const typename dealii::internal::p4est::types<dim>::forest forest; 95 const typename dealii::internal::p4est::types<dim>::tree tree; 96 const typename dealii::internal::p4est::types<dim>::locidx tree_index; 97 const typename DoFHandler<dim, spacedim>::cell_iterator dealii_cell; 98 const typename dealii::internal::p4est::types<dim>::quadrant p4est_cell; 99 100 WorkPackage( 101 const typename dealii::internal::p4est::types<dim>::forest &forest_, 102 const typename dealii::internal::p4est::types<dim>::tree & tree_, 103 const typename dealii::internal::p4est::types<dim>::locidx 104 & tree_index_, 105 const typename DoFHandler<dim, spacedim>::cell_iterator &dealii_cell_, 106 const typename dealii::internal::p4est::types<dim>::quadrant 107 &p4est_cell_) 108 : forest(forest_) 109 , tree(tree_) 110 , tree_index(tree_index_) 111 , dealii_cell(dealii_cell_) 112 , p4est_cell(p4est_cell_) 113 {} 114 }; 115 116 117 // A structure holding all data 118 // of cells needed from other processes 119 // for the extrapolate algorithm. 120 struct CellData 121 { 122 CellData(); 123 124 CellData(const unsigned int dofs_per_cell); 125 126 Vector<value_type> dof_values; 127 128 unsigned int tree_index; 129 130 typename dealii::internal::p4est::types<dim>::quadrant quadrant; 131 132 int receiver; 133 134 bool 135 operator<(const CellData &rhs) const 136 { 137 if (dealii::internal::p4est::functions<dim>::quadrant_compare( 138 &quadrant, &rhs.quadrant) < 0) 139 return true; 140 141 return false; 142 } 143 144 unsigned int 145 bytes_for_buffer() const 146 { 147 return (sizeof(unsigned int) + // dofs_per_cell 148 dof_values.size() * sizeof(value_type) + // dof_values 149 sizeof(unsigned int) + // tree_index 150 sizeof(typename dealii::internal::p4est::types< 151 dim>::quadrant)); // quadrant 152 } 153 154 void 155 pack_data(std::vector<char> &buffer) const 156 { 157 buffer.resize(bytes_for_buffer()); 158 159 char *ptr = buffer.data(); 160 161 unsigned int n_dofs = dof_values.size(); 162 std::memcpy(ptr, &n_dofs, sizeof(unsigned int)); 163 ptr += sizeof(unsigned int); 164 165 std::memcpy(ptr, dof_values.begin(), n_dofs * sizeof(value_type)); 166 ptr += n_dofs * sizeof(value_type); 167 168 std::memcpy(ptr, &tree_index, sizeof(unsigned int)); 169 ptr += sizeof(unsigned int); 170 171 std::memcpy( 172 ptr, 173 &quadrant, 174 sizeof(typename dealii::internal::p4est::types<dim>::quadrant)); 175 ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant); 176 177 Assert(ptr == buffer.data() + buffer.size(), ExcInternalError()); 178 } 179 180 void 181 unpack_data(const std::vector<char> &buffer) 182 { 183 const char * ptr = buffer.data(); 184 unsigned int n_dofs; 185 memcpy(&n_dofs, ptr, sizeof(unsigned int)); 186 ptr += sizeof(unsigned int); 187 188 dof_values.reinit(n_dofs); 189 std::memcpy(dof_values.begin(), ptr, n_dofs * sizeof(value_type)); 190 ptr += n_dofs * sizeof(value_type); 191 192 std::memcpy(&tree_index, ptr, sizeof(unsigned int)); 193 ptr += sizeof(unsigned int); 194 195 std::memcpy( 196 &quadrant, 197 ptr, 198 sizeof(typename dealii::internal::p4est::types<dim>::quadrant)); 199 ptr += sizeof(typename dealii::internal::p4est::types<dim>::quadrant); 200 201 Assert(ptr == buffer.data() + buffer.size(), ExcInternalError()); 202 } 203 }; 204 205 // Problem: The function extrapolates a polynomial 206 // function from a finer mesh of size $h$ to a polynmial 207 // function of higher degree but on a coarser mesh of 208 // size $2h$. Therefore the mesh has to consist of patches 209 // of four (2d) or eight (3d) cells and the function requires 210 // that the mesh is refined globally at least once. 211 // The algorithm starts on the coarsest level of the grid, 212 // loops over all cells and if a cell has at least one active child, 213 // dof values are set via 214 // cell->get_interpolated_dof_values(u_input, dof_values) 215 // cell->set_dof_values_by_interpolation(dof_values, u_output) 216 // both *_interpolation_* functions traverse recursively 217 // over all children down to the active cells 218 // and get/set dof values by interpolation. 219 // 220 // On distributed parallel grids the problem arises, that 221 // if a cell has at least one active child, there is no guarantee 222 // that all children of this cell belong to this process. 223 // There might be children which are owned by other processes 224 // and the algorithm needs to find and has to get the 225 // dof values from these processes and so on. 226 // 227 // Algorithm: 228 // 1) Loop over all coarse cells 229 // From each coarse cell traverse down the tree 230 // of refined cells and search for active children 231 // If there is an active child, check, if all 232 // other children down to the finest level are part 233 // of this process, if not, add the cell to the list 234 // of needs. 235 // 2) Send the list of needs to other processes 236 // Each process has a list of needs and after 237 // the loop over all coarse cells (all trees) 238 // is finished, this list can be send to other processes. 239 // 3) Compute the needs required by other processes 240 // While computing the values required from other 241 // processes there can arise new needs and they are 242 // stored in a list again. 243 // 4) Send the computed values and the list of new needs around 244 // 245 // This procedure has to be repeated until no process needs any 246 // new need and all needs are computed, but there are at most the 247 // "number of grid levels" rounds of sending/receiving cell data. 248 // 249 // Then each process has all data needed for extrapolation. 250 251 // driver function sending/receiving all values from 252 // different processes 253 template <class InVector> 254 void 255 compute_all_non_local_data(const DoFHandler<dim, spacedim> &dof2, 256 const InVector & u); 257 258 // traverse recursively over 259 // the cells of this tree and 260 // interpolate over patches which 261 // are part of this process 262 template <class InVector> 263 void 264 interpolate_recursively( 265 const typename dealii::internal::p4est::types<dim>::forest &forest, 266 const typename dealii::internal::p4est::types<dim>::tree & tree, 267 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 268 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 269 const typename dealii::internal::p4est::types<dim>::quadrant 270 & p4est_cell, 271 const InVector &u1, 272 OutVector & u2); 273 274 // get dof values for this 275 // cell by interpolation 276 // if a child is reached which 277 // is not part of this process 278 // a new need is created 279 template <class InVector> 280 void 281 get_interpolated_dof_values( 282 const typename dealii::internal::p4est::types<dim>::forest &forest, 283 const typename dealii::internal::p4est::types<dim>::tree & tree, 284 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 285 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 286 const typename dealii::internal::p4est::types<dim>::quadrant 287 & p4est_cell, 288 const InVector & u, 289 Vector<value_type> & interpolated_values, 290 std::vector<CellData> &new_needs); 291 292 // set dof values for this 293 // cell by interpolation 294 void 295 set_dof_values_by_interpolation( 296 const typename DoFHandler<dim, spacedim>::cell_iterator &dealii_cell, 297 const typename dealii::internal::p4est::types<dim>::quadrant 298 & p4est_cell, 299 const Vector<value_type> &interpolated_values, 300 OutVector & u); 301 302 // compute all cell_data 303 // needed from other processes 304 // to interpolate the part of 305 // this process 306 void 307 compute_needs(const DoFHandler<dim, spacedim> &dof2, 308 std::vector<CellData> & new_needs); 309 310 // traverse over the tree 311 // and look for patches this 312 // process has to work on 313 void 314 traverse_tree_recursively( 315 const typename dealii::internal::p4est::types<dim>::forest &forest, 316 const typename dealii::internal::p4est::types<dim>::tree & tree, 317 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 318 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 319 const typename dealii::internal::p4est::types<dim>::quadrant 320 & p4est_cell, 321 std::vector<CellData> &new_needs); 322 323 // traverse recursively 324 // over a patch and look 325 // for cells needed from 326 // other processes for interpolation 327 void 328 traverse_patch_recursively( 329 const typename dealii::internal::p4est::types<dim>::forest &forest, 330 const typename dealii::internal::p4est::types<dim>::tree & tree, 331 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 332 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 333 const typename dealii::internal::p4est::types<dim>::quadrant 334 & p4est_cell, 335 std::vector<CellData> &new_needs); 336 337 // compute dof values of all 338 // cells collected in cells_to_compute 339 // computed cells are deleted 340 // from cells_to_compute and 341 // stored in computed_cells 342 // during computation there can 343 // be new cells needed from 344 // other processes, these cells 345 // are stored in new_needs 346 template <class InVector> 347 void 348 compute_cells(const DoFHandler<dim, spacedim> &dof2, 349 const InVector & u, 350 std::vector<CellData> & cells_to_compute, 351 std::vector<CellData> & computed_cells, 352 std::vector<CellData> & new_needs); 353 354 // traverse over the tree 355 // and compute cells 356 template <class InVector> 357 void 358 compute_cells_in_tree_recursively( 359 const typename dealii::internal::p4est::types<dim>::forest &forest, 360 const typename dealii::internal::p4est::types<dim>::tree & tree, 361 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 362 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 363 const typename dealii::internal::p4est::types<dim>::quadrant 364 & p4est_cell, 365 const InVector & u, 366 std::vector<CellData> &cells_to_compute, 367 std::vector<CellData> &computed_cells, 368 std::vector<CellData> &new_needs); 369 370 // send all cell_data in the vector 371 // cells_to_send to their receivers 372 // and receives a vector of cell_data 373 void 374 send_cells(const std::vector<CellData> &cells_to_send, 375 std::vector<CellData> & received_cells) const; 376 377 // add new cell_data to 378 // the ordered list new_needs 379 // uses cell_data_insert 380 static void 381 add_new_need( 382 const typename dealii::internal::p4est::types<dim>::forest &forest, 383 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 384 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 385 const typename dealii::internal::p4est::types<dim>::quadrant 386 & p4est_cell, 387 std::vector<CellData> &new_needs); 388 389 // binary search in cells_list 390 // assume that cells_list 391 // is ordered 392 static int 393 cell_data_search(const CellData & cell_data, 394 const std::vector<CellData> &cells_list); 395 396 // insert cell_data into a sorted 397 // vector cells_list at the correct 398 // position if cell_data 399 // not exists already in cells_list 400 static void 401 cell_data_insert(const CellData & cell_data, 402 std::vector<CellData> &cells_list); 403 404 MPI_Comm communicator; 405 406 // a vector of all cells this process has 407 // computed or received data 408 std::vector<CellData> available_cells; 409 410 // stores the indices of dofs on more refined ghosted cells along 411 // with the maximum level 412 std::map<types::global_dof_index, int> dofs_on_refined_neighbors; 413 414 // counts the send/receive round we are in 415 unsigned int round; 416 }; 417 418 template <class OutVector> 419 class ExtrapolateImplementation<1, 1, OutVector> 420 { 421 public: 422 ExtrapolateImplementation() 423 { 424 AssertThrow(false, ExcNotImplemented()) 425 } 426 427 template <class InVector> 428 void 429 extrapolate_parallel(const InVector & /*u2_relevant*/, 430 const DoFHandler<1, 1> & /*dof2*/, 431 OutVector & /*u2*/) 432 {} 433 }; 434 435 template <class OutVector> 436 class ExtrapolateImplementation<1, 2, OutVector> 437 { 438 public: 439 ExtrapolateImplementation() 440 { 441 AssertThrow(false, ExcNotImplemented()) 442 } 443 444 template <class InVector> 445 void 446 extrapolate_parallel(const InVector & /*u2_relevant*/, 447 const DoFHandler<1, 2> & /*dof2*/, 448 OutVector & /*u2*/) 449 {} 450 }; 451 452 template <class OutVector> 453 class ExtrapolateImplementation<1, 3, OutVector> 454 { 455 public: 456 ExtrapolateImplementation() 457 { 458 AssertThrow(false, ExcNotImplemented()) 459 } 460 461 template <class InVector> 462 void 463 extrapolate_parallel(const InVector & /*u2_relevant*/, 464 const DoFHandler<1, 3> & /*dof2*/, 465 OutVector & /*u2*/) 466 {} 467 }; 468 469 470 471 template <int dim, int spacedim, class OutVector> 472 ExtrapolateImplementation<dim, spacedim, OutVector>:: 473 ExtrapolateImplementation() 474 : round(0) 475 {} 476 477 478 479 template <int dim, int spacedim, class OutVector> 480 ExtrapolateImplementation<dim, spacedim, OutVector>::CellData::CellData() 481 : tree_index(0) 482 , receiver(0) 483 {} 484 485 486 487 template <int dim, int spacedim, class OutVector> 488 ExtrapolateImplementation<dim, spacedim, OutVector>::CellData::CellData( 489 const unsigned int dofs_per_cell) 490 : tree_index(0) 491 , receiver(0) 492 { 493 dof_values.reinit(dofs_per_cell); 494 } 495 496 497 498 template <int dim, int spacedim, class OutVector> 499 template <class InVector> 500 void 501 ExtrapolateImplementation<dim, spacedim, OutVector>:: 502 interpolate_recursively( 503 const typename dealii::internal::p4est::types<dim>::forest &forest, 504 const typename dealii::internal::p4est::types<dim>::tree & tree, 505 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 506 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 507 const typename dealii::internal::p4est::types<dim>::quadrant 508 & p4est_cell, 509 const InVector &u1, 510 OutVector & u2) 511 { 512 // check if this cell exists in the local p4est 513 int idx = sc_array_bsearch( 514 const_cast<sc_array_t *>(&tree.quadrants), 515 &p4est_cell, 516 dealii::internal::p4est::functions<dim>::quadrant_compare); 517 518 // if neither this cell nor one of it's children belongs to us, don't do 519 // anything 520 if (idx == -1 && 521 (dealii::internal::p4est::functions<dim>::quadrant_overlaps_tree( 522 const_cast<typename dealii::internal::p4est::types<dim>::tree *>( 523 &tree), 524 &p4est_cell) == false)) 525 return; 526 527 bool p4est_has_children = (idx == -1); 528 529 bool locally_owned_children = false; 530 if (p4est_has_children) 531 { 532 Assert(dealii_cell->has_children(), ExcInternalError()); 533 534 // check if at least one child is locally owned on our process 535 for (unsigned int child_n = 0; child_n < dealii_cell->n_children(); 536 ++child_n) 537 if (dealii_cell->child(child_n)->is_active()) 538 if (dealii_cell->child(child_n)->is_locally_owned()) 539 { 540 locally_owned_children = true; 541 break; 542 } 543 } 544 545 if (locally_owned_children) 546 { 547 const FiniteElement<dim, spacedim> &fe = 548 dealii_cell->get_dof_handler().get_fe(); 549 const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); 550 551 Vector<typename OutVector::value_type> interpolated_values( 552 dofs_per_cell); 553 554 std::vector<CellData> new_needs; 555 get_interpolated_dof_values(forest, 556 tree, 557 tree_index, 558 dealii_cell, 559 p4est_cell, 560 u1, 561 interpolated_values, 562 new_needs); 563 564 // at this point of 565 // the procedure no new 566 // needs should come up 567 Assert(new_needs.size() == 0, ExcInternalError()); 568 569 set_dof_values_by_interpolation(dealii_cell, 570 p4est_cell, 571 interpolated_values, 572 u2); 573 } 574 } 575 576 577 578 template <int dim, int spacedim, class OutVector> 579 template <class InVector> 580 void 581 ExtrapolateImplementation<dim, spacedim, OutVector>:: 582 get_interpolated_dof_values( 583 const typename dealii::internal::p4est::types<dim>::forest &forest, 584 const typename dealii::internal::p4est::types<dim>::tree & tree, 585 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 586 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 587 const typename dealii::internal::p4est::types<dim>::quadrant 588 & p4est_cell, 589 const InVector & u, 590 Vector<value_type> & interpolated_values, 591 std::vector<CellData> &new_needs) 592 { 593 if (dealii_cell->is_active()) 594 { 595 if (dealii_cell->is_locally_owned()) 596 { 597 // if this is one of our cells, 598 // get dof values from input vector 599 dealii_cell->get_dof_values(u, interpolated_values); 600 } 601 else 602 { 603 add_new_need( 604 forest, tree_index, dealii_cell, p4est_cell, new_needs); 605 } 606 } 607 else 608 { 609 const FiniteElement<dim, spacedim> &fe = 610 dealii_cell->get_dof_handler().get_fe(); 611 const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); 612 613 Assert(interpolated_values.size() == dofs_per_cell, 614 ExcDimensionMismatch(interpolated_values.size(), 615 dofs_per_cell)); 616 Assert(u.size() == dealii_cell->get_dof_handler().n_dofs(), 617 ExcDimensionMismatch(u.size(), 618 dealii_cell->get_dof_handler().n_dofs())); 619 620 Vector<value_type> tmp1(dofs_per_cell); 621 Vector<value_type> tmp2(dofs_per_cell); 622 623 interpolated_values = 0; 624 std::vector<bool> restriction_is_additive(dofs_per_cell); 625 for (unsigned int i = 0; i < dofs_per_cell; ++i) 626 restriction_is_additive[i] = fe.restriction_is_additive(i); 627 628 typename dealii::internal::p4est::types<dim>::quadrant 629 p4est_child[GeometryInfo<dim>::max_children_per_cell]; 630 631 dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, 632 p4est_child); 633 634 bool found_child = true; 635 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; 636 ++c) 637 { 638 if (dealii::internal::p4est::functions<dim>:: 639 quadrant_overlaps_tree( 640 const_cast< 641 typename dealii::internal::p4est::types<dim>::tree *>( 642 &tree), 643 &p4est_child[c]) == false) 644 { 645 // this is a cell this process needs 646 // data from another process 647 648 // check if this cell is 649 // available from other 650 // computed patches 651 CellData cell_data; 652 cell_data.quadrant = p4est_child[c]; 653 int pos = cell_data_search(cell_data, available_cells); 654 655 if (pos == -1) 656 { 657 // data is not available 658 // create a new need 659 found_child = false; 660 661 add_new_need(forest, 662 tree_index, 663 dealii_cell->child(c), 664 p4est_child[c], 665 new_needs); 666 } 667 else 668 { 669 Assert(available_cells[pos].dof_values.size() == 670 dofs_per_cell, 671 ExcDimensionMismatch( 672 available_cells[pos].dof_values.size(), 673 dofs_per_cell)); 674 675 tmp1 = available_cells[pos].dof_values; 676 } 677 } 678 else 679 { 680 // get the values from the present child, if necessary by 681 // interpolation either from its own children or 682 // by interpolating from the finite element on an active 683 // child to the finite element space requested here 684 get_interpolated_dof_values(forest, 685 tree, 686 tree_index, 687 dealii_cell->child(c), 688 p4est_child[c], 689 u, 690 tmp1, 691 new_needs); 692 } 693 694 if (found_child) 695 { 696 // interpolate these to the mother cell 697 fe.get_restriction_matrix(c, dealii_cell->refinement_case()) 698 .vmult(tmp2, tmp1); 699 700 // and add up or set them in the output vector 701 for (unsigned int i = 0; i < dofs_per_cell; ++i) 702 if (tmp2(i) != value_type()) 703 { 704 if (restriction_is_additive[i]) 705 interpolated_values(i) += tmp2(i); 706 else 707 interpolated_values(i) = tmp2(i); 708 } 709 } 710 } 711 712 if (found_child == false) 713 interpolated_values = 0; 714 } 715 } 716 717 718 719 template <int dim, int spacedim, class OutVector> 720 void 721 ExtrapolateImplementation<dim, spacedim, OutVector>:: 722 set_dof_values_by_interpolation( 723 const typename DoFHandler<dim, spacedim>::cell_iterator &dealii_cell, 724 const typename dealii::internal::p4est::types<dim>::quadrant 725 & p4est_cell, 726 const Vector<value_type> &local_values, 727 OutVector & u) 728 { 729 const FiniteElement<dim, spacedim> &fe = 730 dealii_cell->get_dof_handler().get_fe(); 731 const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); 732 733 if (dealii_cell->is_active()) 734 { 735 if (dealii_cell->is_locally_owned()) 736 { 737 std::vector<types::global_dof_index> indices(dofs_per_cell); 738 dealii_cell->get_dof_indices(indices); 739 for (unsigned int j = 0; j < dofs_per_cell; ++j) 740 { 741 // don't set this dof if there is a more refined ghosted 742 // neighbor setting this dof. 743 const bool on_refined_neighbor = 744 (dofs_on_refined_neighbors.find(indices[j]) != 745 dofs_on_refined_neighbors.end()); 746 if (!(on_refined_neighbor && 747 dofs_on_refined_neighbors[indices[j]] > 748 dealii_cell->level())) 749 ::dealii::internal::ElementAccess<OutVector>::set( 750 local_values(j), indices[j], u); 751 } 752 } 753 } 754 else 755 { 756 Assert(local_values.size() == dofs_per_cell, 757 ExcDimensionMismatch(local_values.size(), dofs_per_cell)); 758 Assert(u.size() == dealii_cell->get_dof_handler().n_dofs(), 759 ExcDimensionMismatch(u.size(), 760 dealii_cell->get_dof_handler().n_dofs())); 761 762 Vector<value_type> tmp(dofs_per_cell); 763 764 typename dealii::internal::p4est::types<dim>::quadrant 765 p4est_child[GeometryInfo<dim>::max_children_per_cell]; 766 767 dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, 768 p4est_child); 769 770 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; 771 ++c) 772 { 773 if (tmp.size() > 0) 774 fe.get_prolongation_matrix(c, dealii_cell->refinement_case()) 775 .vmult(tmp, local_values); 776 777 set_dof_values_by_interpolation(dealii_cell->child(c), 778 p4est_child[c], 779 tmp, 780 u); 781 } 782 } 783 } 784 785 786 787 template <int dim, int spacedim, class OutVector> 788 void 789 ExtrapolateImplementation<dim, spacedim, OutVector>::compute_needs( 790 const DoFHandler<dim, spacedim> &dof2, 791 std::vector<CellData> & new_needs) 792 { 793 const parallel::distributed::Triangulation<dim, spacedim> *tr = 794 (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> 795 *>(&dof2.get_triangulation())); 796 797 Assert(tr != nullptr, ExcInternalError()); 798 799 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof2.begin(0), 800 endc = dof2.end(0); 801 for (; cell != endc; ++cell) 802 { 803 if (dealii::internal::p4est::tree_exists_locally<dim>( 804 tr->parallel_forest, 805 tr->coarse_cell_to_p4est_tree_permutation[cell->index()]) == 806 false) 807 continue; 808 809 typename dealii::internal::p4est::types<dim>::quadrant 810 p4est_coarse_cell; 811 const unsigned int tree_index = 812 tr->coarse_cell_to_p4est_tree_permutation[cell->index()]; 813 typename dealii::internal::p4est::types<dim>::tree *tree = 814 tr->init_tree(cell->index()); 815 816 dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell); 817 818 // make sure that each cell on the 819 // coarsest level is at least once 820 // refined, otherwise, these cells 821 // can't be treated and would 822 // generate a bogus result 823 { 824 int idx = sc_array_bsearch( 825 const_cast<sc_array_t *>(&tree->quadrants), 826 &p4est_coarse_cell, 827 dealii::internal::p4est::functions<dim>::quadrant_compare); 828 829 AssertThrow(idx == -1, ExcGridNotRefinedAtLeastOnce()); 830 } 831 832 traverse_tree_recursively(*tr->parallel_forest, 833 *tree, 834 tree_index, 835 cell, 836 p4est_coarse_cell, 837 new_needs); 838 } 839 } 840 841 842 843 template <int dim, int spacedim, class OutVector> 844 void 845 ExtrapolateImplementation<dim, spacedim, OutVector>:: 846 traverse_tree_recursively( 847 const typename dealii::internal::p4est::types<dim>::forest &forest, 848 const typename dealii::internal::p4est::types<dim>::tree & tree, 849 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 850 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 851 const typename dealii::internal::p4est::types<dim>::quadrant 852 & p4est_cell, 853 std::vector<CellData> &new_needs) 854 { 855 // check if this cell exists in the local p4est 856 int idx = sc_array_bsearch( 857 const_cast<sc_array_t *>(&tree.quadrants), 858 &p4est_cell, 859 dealii::internal::p4est::functions<dim>::quadrant_compare); 860 861 // if neither this cell nor one of it's children belongs to us, don't do 862 // anything 863 if (idx == -1 && 864 (dealii::internal::p4est::functions<dim>::quadrant_overlaps_tree( 865 const_cast<typename dealii::internal::p4est::types<dim>::tree *>( 866 &tree), 867 &p4est_cell) == false)) 868 return; 869 870 bool p4est_has_children = (idx == -1); 871 872 // this cell is part of a patch 873 // this process has to interpolate on 874 // if there is at least one locally 875 // owned child 876 bool locally_owned_children = false; 877 if (p4est_has_children) 878 { 879 Assert(dealii_cell->has_children(), ExcInternalError()); 880 881 for (unsigned int child_n = 0; child_n < dealii_cell->n_children(); 882 ++child_n) 883 if (dealii_cell->child(child_n)->is_active()) 884 if (dealii_cell->child(child_n)->is_locally_owned()) 885 { 886 locally_owned_children = true; 887 break; 888 } 889 } 890 891 // traverse the patch recursively 892 // to find cells needed from 893 // other processes 894 if (locally_owned_children) 895 { 896 traverse_patch_recursively( 897 forest, tree, tree_index, dealii_cell, p4est_cell, new_needs); 898 } 899 900 // traverse recursively over this tree 901 if (p4est_has_children) 902 { 903 typename dealii::internal::p4est::types<dim>::quadrant 904 p4est_child[GeometryInfo<dim>::max_children_per_cell]; 905 906 dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, 907 p4est_child); 908 909 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; 910 ++c) 911 { 912 traverse_tree_recursively(forest, 913 tree, 914 tree_index, 915 dealii_cell->child(c), 916 p4est_child[c], 917 new_needs); 918 } 919 } 920 } 921 922 923 924 template <int dim, int spacedim, class OutVector> 925 void 926 ExtrapolateImplementation<dim, spacedim, OutVector>:: 927 traverse_patch_recursively( 928 const typename dealii::internal::p4est::types<dim>::forest &forest, 929 const typename dealii::internal::p4est::types<dim>::tree & tree, 930 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 931 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 932 const typename dealii::internal::p4est::types<dim>::quadrant 933 & p4est_cell, 934 std::vector<CellData> &new_needs) 935 { 936 if (dealii_cell->has_children()) 937 { 938 Assert(dealii_cell->has_children(), ExcInternalError()); 939 940 typename dealii::internal::p4est::types<dim>::quadrant 941 p4est_child[GeometryInfo<dim>::max_children_per_cell]; 942 943 dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, 944 p4est_child); 945 946 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; 947 ++c) 948 { 949 // check if this child 950 // is locally available 951 // in the p4est 952 if (dealii::internal::p4est::functions<dim>:: 953 quadrant_overlaps_tree( 954 const_cast< 955 typename dealii::internal::p4est::types<dim>::tree *>( 956 &tree), 957 &p4est_child[c]) == false) 958 { 959 // this is a cell for which this process 960 // needs data from another process 961 // so add the cell to the list 962 add_new_need(forest, 963 tree_index, 964 dealii_cell->child(c), 965 p4est_child[c], 966 new_needs); 967 } 968 else 969 { 970 // at least some part of 971 // the tree rooted in this 972 // child is locally available 973 traverse_patch_recursively(forest, 974 tree, 975 tree_index, 976 dealii_cell->child(c), 977 p4est_child[c], 978 new_needs); 979 } 980 } 981 } 982 } 983 984 985 986 template <int dim, int spacedim, class OutVector> 987 template <class InVector> 988 void 989 ExtrapolateImplementation<dim, spacedim, OutVector>::compute_cells( 990 const DoFHandler<dim, spacedim> &dof2, 991 const InVector & u, 992 std::vector<CellData> & cells_to_compute, 993 std::vector<CellData> & computed_cells, 994 std::vector<CellData> & new_needs) 995 { 996 const parallel::distributed::Triangulation<dim, spacedim> *tr = 997 (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> 998 *>(&dof2.get_triangulation())); 999 1000 Assert(tr != nullptr, ExcInternalError()); 1001 1002 // collect in a set all trees this 1003 // process has to compute cells on 1004 std::set<unsigned int> trees; 1005 for (typename std::vector<CellData>::const_iterator it = 1006 cells_to_compute.begin(); 1007 it != cells_to_compute.end(); 1008 ++it) 1009 trees.insert(it->tree_index); 1010 1011 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof2.begin(0), 1012 endc = dof2.end(0); 1013 for (; cell != endc; ++cell) 1014 { 1015 // check if this is a tree this process has to 1016 // work on and that this tree is in the p4est 1017 const unsigned int tree_index = 1018 tr->coarse_cell_to_p4est_tree_permutation[cell->index()]; 1019 1020 if ((trees.find(tree_index) == trees.end()) || 1021 (dealii::internal::p4est::tree_exists_locally<dim>( 1022 tr->parallel_forest, tree_index) == false)) 1023 continue; 1024 1025 typename dealii::internal::p4est::types<dim>::quadrant 1026 p4est_coarse_cell; 1027 typename dealii::internal::p4est::types<dim>::tree *tree = 1028 tr->init_tree(cell->index()); 1029 1030 dealii::internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell); 1031 1032 compute_cells_in_tree_recursively(*tr->parallel_forest, 1033 *tree, 1034 tree_index, 1035 cell, 1036 p4est_coarse_cell, 1037 u, 1038 cells_to_compute, 1039 computed_cells, 1040 new_needs); 1041 } 1042 } 1043 1044 1045 1046 template <int dim, int spacedim, class OutVector> 1047 template <class InVector> 1048 void 1049 ExtrapolateImplementation<dim, spacedim, OutVector>:: 1050 compute_cells_in_tree_recursively( 1051 const typename dealii::internal::p4est::types<dim>::forest &forest, 1052 const typename dealii::internal::p4est::types<dim>::tree & tree, 1053 const typename dealii::internal::p4est::types<dim>::locidx &tree_index, 1054 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 1055 const typename dealii::internal::p4est::types<dim>::quadrant 1056 & p4est_cell, 1057 const InVector & u, 1058 std::vector<CellData> &cells_to_compute, 1059 std::vector<CellData> &computed_cells, 1060 std::vector<CellData> &new_needs) 1061 { 1062 if (cells_to_compute.size() == 0) 1063 return; 1064 1065 // check if this cell exists in the local p4est 1066 int idx = sc_array_bsearch( 1067 const_cast<sc_array_t *>(&tree.quadrants), 1068 &p4est_cell, 1069 dealii::internal::p4est::functions<dim>::quadrant_compare); 1070 1071 // if neither this cell nor one of it's children belongs to us, don't do 1072 // anything 1073 if (idx == -1 && 1074 (dealii::internal::p4est::functions<dim>::quadrant_overlaps_tree( 1075 const_cast<typename dealii::internal::p4est::types<dim>::tree *>( 1076 &tree), 1077 &p4est_cell) == false)) 1078 return; 1079 1080 bool p4est_has_children = (idx == -1); 1081 1082 // check if this quadrant is in the list 1083 CellData cell_data; 1084 cell_data.quadrant = p4est_cell; 1085 int pos = cell_data_search(cell_data, cells_to_compute); 1086 if (pos != -1) 1087 { 1088 std::vector<CellData> tmp; 1089 // compute dof values 1090 get_interpolated_dof_values(forest, 1091 tree, 1092 tree_index, 1093 dealii_cell, 1094 p4est_cell, 1095 u, 1096 cells_to_compute[pos].dof_values, 1097 tmp); 1098 // there is no new cell_data 1099 // this process needs for computing this cell, 1100 // store cell_data in the list of 1101 // computed cells and erase this cell 1102 // from the list of cells to compute 1103 if (tmp.size() == 0) 1104 { 1105 cell_data_insert(cells_to_compute[pos], computed_cells); 1106 cells_to_compute.erase(cells_to_compute.begin() + pos); 1107 } 1108 else 1109 { 1110 for (unsigned int i = 0; i < tmp.size(); ++i) 1111 cell_data_insert(tmp[i], new_needs); 1112 } 1113 } 1114 1115 // search recursively over this tree 1116 if (p4est_has_children) 1117 { 1118 typename dealii::internal::p4est::types<dim>::quadrant 1119 p4est_child[GeometryInfo<dim>::max_children_per_cell]; 1120 1121 dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, 1122 p4est_child); 1123 1124 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; 1125 ++c) 1126 { 1127 compute_cells_in_tree_recursively(forest, 1128 tree, 1129 tree_index, 1130 dealii_cell->child(c), 1131 p4est_child[c], 1132 u, 1133 cells_to_compute, 1134 computed_cells, 1135 new_needs); 1136 } 1137 } 1138 } 1139 1140 1141 1142 template <int dim, int spacedim, class OutVector> 1143 void 1144 ExtrapolateImplementation<dim, spacedim, OutVector>::send_cells( 1145 const std::vector<CellData> &cells_to_send, 1146 std::vector<CellData> & received_cells) const 1147 { 1148 std::vector<std::vector<char>> sendbuffers(cells_to_send.size()); 1149 std::vector<std::vector<char>>::iterator buffer = sendbuffers.begin(); 1150 std::vector<MPI_Request> requests(cells_to_send.size()); 1151 std::vector<unsigned int> destinations; 1152 1153 // Protect the communication below: 1154 static Utilities::MPI::CollectiveMutex mutex; 1155 Utilities::MPI::CollectiveMutex::ScopedLock lock(mutex, communicator); 1156 1157 // We pick a new tag in each round. Wrap around after 10 rounds: 1158 const int mpi_tag = 1159 Utilities::MPI::internal::Tags::fe_tools_extrapolate + round % 10; 1160 1161 // send data 1162 unsigned int idx = 0; 1163 for (typename std::vector<CellData>::const_iterator it = 1164 cells_to_send.begin(); 1165 it != cells_to_send.end(); 1166 ++it, ++idx) 1167 { 1168 destinations.push_back(it->receiver); 1169 1170 it->pack_data(*buffer); 1171 const int ierr = MPI_Isend(buffer->data(), 1172 buffer->size(), 1173 MPI_BYTE, 1174 it->receiver, 1175 mpi_tag, 1176 communicator, 1177 &requests[idx]); 1178 AssertThrowMPI(ierr); 1179 } 1180 1181 Assert(destinations.size() == cells_to_send.size(), ExcInternalError()); 1182 1183 const unsigned int n_senders = 1184 Utilities::MPI::compute_n_point_to_point_communications(communicator, 1185 destinations); 1186 1187 // receive data 1188 std::vector<char> receive; 1189 CellData cell_data; 1190 for (unsigned int index = 0; index < n_senders; ++index) 1191 { 1192 MPI_Status status; 1193 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag, communicator, &status); 1194 AssertThrowMPI(ierr); 1195 1196 int len; 1197 ierr = MPI_Get_count(&status, MPI_BYTE, &len); 1198 AssertThrowMPI(ierr); 1199 receive.resize(len); 1200 1201 char *buf = receive.data(); 1202 ierr = MPI_Recv(buf, 1203 len, 1204 MPI_BYTE, 1205 status.MPI_SOURCE, 1206 status.MPI_TAG, 1207 communicator, 1208 &status); 1209 AssertThrowMPI(ierr); 1210 1211 cell_data.unpack_data(receive); 1212 1213 // this process has to send this 1214 // cell back to the sender 1215 // the receiver is the old sender 1216 cell_data.receiver = status.MPI_SOURCE; 1217 1218 received_cells.push_back(cell_data); 1219 } 1220 1221 if (requests.size() > 0) 1222 { 1223 const int ierr = 1224 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); 1225 AssertThrowMPI(ierr); 1226 } 1227 1228 // finally sort the list of cells 1229 std::sort(received_cells.begin(), received_cells.end()); 1230 } 1231 1232 1233 1234 template <int dim, int spacedim, class OutVector> 1235 void 1236 ExtrapolateImplementation<dim, spacedim, OutVector>::add_new_need( 1237 const typename dealii::internal::p4est::types<dim>::forest & forest, 1238 const typename dealii::internal::p4est::types<dim>::locidx & tree_index, 1239 const typename DoFHandler<dim, spacedim>::cell_iterator & dealii_cell, 1240 const typename dealii::internal::p4est::types<dim>::quadrant &p4est_cell, 1241 std::vector<CellData> & new_needs) 1242 { 1243 const FiniteElement<dim, spacedim> &fe = 1244 dealii_cell->get_dof_handler().get_fe(); 1245 const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); 1246 1247 CellData cell_data(dofs_per_cell); 1248 cell_data.quadrant = p4est_cell; 1249 cell_data.tree_index = tree_index; 1250 cell_data.receiver = 1251 dealii::internal::p4est::functions<dim>::comm_find_owner( 1252 const_cast<typename dealii::internal::p4est::types<dim>::forest *>( 1253 &forest), 1254 tree_index, 1255 &p4est_cell, 1256 dealii_cell->level_subdomain_id()); 1257 1258 cell_data_insert(cell_data, new_needs); 1259 } 1260 1261 1262 1263 template <int dim, int spacedim, class OutVector> 1264 int 1265 ExtrapolateImplementation<dim, spacedim, OutVector>::cell_data_search( 1266 const CellData & cell_data, 1267 const std::vector<CellData> &cells_list) 1268 { 1269 typename std::vector<CellData>::const_iterator bound = 1270 std::lower_bound(cells_list.begin(), cells_list.end(), cell_data); 1271 1272 if ((bound != cells_list.end()) && !(cell_data < *bound)) 1273 return static_cast<int>(bound - cells_list.begin()); 1274 1275 return -1; 1276 } 1277 1278 1279 1280 template <int dim, int spacedim, class OutVector> 1281 void 1282 ExtrapolateImplementation<dim, spacedim, OutVector>::cell_data_insert( 1283 const CellData & cell_data, 1284 std::vector<CellData> &cells_list) 1285 { 1286 typename std::vector<CellData>::iterator bound = 1287 std::lower_bound(cells_list.begin(), cells_list.end(), cell_data); 1288 1289 if ((bound == cells_list.end()) || (cell_data < *bound)) 1290 cells_list.insert(bound, 1, cell_data); 1291 } 1292 1293 1294 1295 template <int dim, int spacedim, class OutVector> 1296 template <class InVector> 1297 void 1298 ExtrapolateImplementation<dim, spacedim, OutVector>:: 1299 compute_all_non_local_data(const DoFHandler<dim, spacedim> &dof2, 1300 const InVector & u) 1301 { 1302 std::vector<CellData> cells_we_need, cells_to_compute, received_cells, 1303 received_needs, new_needs, computed_cells, cells_to_send; 1304 1305 // reset the round count we will use in send_cells 1306 round = 0; 1307 1308 // Compute all the cells needed from other processes. 1309 compute_needs(dof2, cells_we_need); 1310 1311 // Send the cells needed to their owners and receive 1312 // a list of cells other processes need from us. 1313 send_cells(cells_we_need, received_needs); 1314 1315 // The list of received needs can contain some cells more than once 1316 // because different processes may need data from the same cell. 1317 // To compute data only once, generate a vector with unique entries 1318 // and distribute the computed data afterwards back to a vector with 1319 // correct receivers. 1320 // Computing cell_data can cause a need for data from some new cells. 1321 // If a cell is computed, send it back to their senders, maybe receive 1322 // new needs and compute again, do not wait that all cells are computed 1323 // or all needs are collected. 1324 // Otherwise we can run into a deadlock if a cell needed from another 1325 // process itself needs some data from us. 1326 unsigned int ready = 0; 1327 do 1328 { 1329 for (unsigned int i = 0; i < received_needs.size(); ++i) 1330 cell_data_insert(received_needs[i], cells_to_compute); 1331 1332 compute_cells(dof2, u, cells_to_compute, computed_cells, new_needs); 1333 1334 // if there are no cells to compute and no new needs, stop 1335 ready = 1336 Utilities::MPI::sum(new_needs.size() + cells_to_compute.size(), 1337 communicator); 1338 1339 for (typename std::vector<CellData>::const_iterator comp = 1340 computed_cells.begin(); 1341 comp != computed_cells.end(); 1342 ++comp) 1343 { 1344 // store computed cells... 1345 cell_data_insert(*comp, available_cells); 1346 1347 // ...and generate a vector of computed cells with correct 1348 // receivers, then delete this received need from the list 1349 typename std::vector<CellData>::iterator recv = 1350 received_needs.begin(); 1351 while (recv != received_needs.end()) 1352 { 1353 if (dealii::internal::p4est::quadrant_is_equal<dim>( 1354 recv->quadrant, comp->quadrant)) 1355 { 1356 recv->dof_values = comp->dof_values; 1357 cells_to_send.push_back(*recv); 1358 received_needs.erase(recv); 1359 recv = received_needs.begin(); 1360 } 1361 else 1362 ++recv; 1363 } 1364 } 1365 1366 // increase the round counter, such that we are sure to only send 1367 // and receive data from the correct call 1368 ++round; 1369 1370 send_cells(cells_to_send, received_cells); 1371 1372 // store received cell_data 1373 for (typename std::vector<CellData>::const_iterator recv = 1374 received_cells.begin(); 1375 recv != received_cells.end(); 1376 ++recv) 1377 { 1378 cell_data_insert(*recv, available_cells); 1379 } 1380 1381 // increase the round counter, such that we are sure to only send 1382 // and receive data from the correct call 1383 ++round; 1384 1385 // finally send and receive new needs and start a new round 1386 send_cells(new_needs, received_needs); 1387 } 1388 while (ready != 0); 1389 } 1390 1391 1392 1393 template <int dim, int spacedim, class OutVector> 1394 template <class InVector> 1395 void 1396 ExtrapolateImplementation<dim, spacedim, OutVector>::extrapolate_parallel( 1397 const InVector & u2_relevant, 1398 const DoFHandler<dim, spacedim> &dof2, 1399 OutVector & u2) 1400 { 1401 const parallel::distributed::Triangulation<dim, spacedim> *tr = 1402 (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> 1403 *>(&dof2.get_triangulation())); 1404 1405 Assert( 1406 tr != nullptr, 1407 ExcMessage( 1408 "Extrapolate in parallel only works for parallel distributed triangulations!")); 1409 1410 communicator = tr->get_communicator(); 1411 1412 compute_all_non_local_data(dof2, u2_relevant); 1413 1414 // exclude dofs on more refined ghosted cells 1415 const FiniteElement<dim, spacedim> &fe = dof2.get_fe(); 1416 if (fe.max_dofs_per_face() > 0) 1417 { 1418 const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); 1419 std::vector<types::global_dof_index> indices(dofs_per_cell); 1420 typename DoFHandler<dim, spacedim>::active_cell_iterator 1421 cell = dof2.begin_active(), 1422 endc = dof2.end(); 1423 for (; cell != endc; ++cell) 1424 if (cell->is_ghost()) 1425 { 1426 cell->get_dof_indices(indices); 1427 for (const unsigned int face : 1428 GeometryInfo<dim>::face_indices()) 1429 if (!cell->at_boundary(face)) 1430 { 1431 const typename DoFHandler<dim, spacedim>::cell_iterator 1432 neighbor = cell->neighbor(face); 1433 if (neighbor->level() != cell->level()) 1434 for (unsigned int i = 0; i < fe.n_dofs_per_face(face); 1435 ++i) 1436 { 1437 const types::global_dof_index index = 1438 indices[fe.face_to_cell_index(i, face)]; 1439 ; 1440 const bool index_stored = 1441 (dofs_on_refined_neighbors.find(index) != 1442 dofs_on_refined_neighbors.end()); 1443 const unsigned int level = 1444 index_stored ? 1445 std::max(cell->level(), 1446 dofs_on_refined_neighbors[index]) : 1447 cell->level(); 1448 dofs_on_refined_neighbors[index] = level; 1449 } 1450 } 1451 } 1452 } 1453 1454 // after all dof values on 1455 // not owned patch cells 1456 // are computed, start 1457 // the interpolation 1458 u2 = typename OutVector::value_type(0.); 1459 1460 std::queue<WorkPackage> queue; 1461 { 1462 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof2.begin(0), 1463 endc = dof2.end(0); 1464 for (; cell != endc; ++cell) 1465 { 1466 if (dealii::internal::p4est::tree_exists_locally<dim>( 1467 tr->parallel_forest, 1468 tr->coarse_cell_to_p4est_tree_permutation[cell->index()]) == 1469 false) 1470 continue; 1471 1472 typename dealii::internal::p4est::types<dim>::quadrant 1473 p4est_coarse_cell; 1474 const unsigned int tree_index = 1475 tr->coarse_cell_to_p4est_tree_permutation[cell->index()]; 1476 typename dealii::internal::p4est::types<dim>::tree *tree = 1477 tr->init_tree(cell->index()); 1478 1479 dealii::internal::p4est::init_coarse_quadrant<dim>( 1480 p4est_coarse_cell); 1481 1482 const WorkPackage data( 1483 *tr->parallel_forest, *tree, tree_index, cell, p4est_coarse_cell); 1484 1485 queue.push(data); 1486 } 1487 } 1488 1489 while (!queue.empty()) 1490 { 1491 const WorkPackage &data = queue.front(); 1492 1493 const typename dealii::internal::p4est::types<dim>::forest &forest = 1494 data.forest; 1495 const typename dealii::internal::p4est::types<dim>::tree &tree = 1496 data.tree; 1497 const typename dealii::internal::p4est::types<dim>::locidx 1498 &tree_index = data.tree_index; 1499 const typename DoFHandler<dim, spacedim>::cell_iterator &dealii_cell = 1500 data.dealii_cell; 1501 const typename dealii::internal::p4est::types<dim>::quadrant 1502 &p4est_cell = data.p4est_cell; 1503 1504 interpolate_recursively( 1505 forest, tree, tree_index, dealii_cell, p4est_cell, u2_relevant, u2); 1506 1507 // traverse recursively over this tree 1508 if (dealii_cell->has_children()) 1509 { 1510 Assert(dealii_cell->has_children(), ExcInternalError()); 1511 typename dealii::internal::p4est::types<dim>::quadrant 1512 p4est_child[GeometryInfo<dim>::max_children_per_cell]; 1513 1514 dealii::internal::p4est::init_quadrant_children<dim>(p4est_cell, 1515 p4est_child); 1516 1517 for (unsigned int c = 0; 1518 c < GeometryInfo<dim>::max_children_per_cell; 1519 ++c) 1520 queue.push(WorkPackage(forest, 1521 tree, 1522 tree_index, 1523 dealii_cell->child(c), 1524 p4est_child[c])); 1525 } 1526 queue.pop(); 1527 } 1528 1529 1530 u2.compress(VectorOperation::insert); 1531 } 1532 #endif // DEAL_II_WITH_P4EST 1533 1534 template <class VectorType, typename dummy = void> 1535 struct BlockTypeHelper 1536 { 1537 using type = VectorType; 1538 }; 1539 1540 template <class VectorType> 1541 struct BlockTypeHelper< 1542 VectorType, 1543 typename std::enable_if<IsBlockVector<VectorType>::value>::type> 1544 { 1545 using type = typename VectorType::BlockType; 1546 }; 1547 1548 template <class VectorType> 1549 using BlockType = typename BlockTypeHelper<VectorType>::type; 1550 1551 template <class VectorType, class DH> 1552 void 1553 reinit_distributed(const DH &dh, VectorType &vector) 1554 { 1555 vector.reinit(dh.n_dofs()); 1556 } 1557 1558 #ifdef DEAL_II_WITH_PETSC 1559 template <int dim, int spacedim> 1560 void 1561 reinit_distributed(const DoFHandler<dim, spacedim> &dh, 1562 PETScWrappers::MPI::Vector & vector) 1563 { 1564 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1565 dynamic_cast< 1566 const parallel::distributed::Triangulation<dim, spacedim> *>( 1567 &dh.get_triangulation()); 1568 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1569 1570 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1571 vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); 1572 } 1573 #endif // DEAL_II_WITH_PETSC 1574 1575 #ifdef DEAL_II_WITH_TRILINOS 1576 template <int dim, int spacedim> 1577 void 1578 reinit_distributed(const DoFHandler<dim, spacedim> &dh, 1579 TrilinosWrappers::MPI::Vector & vector) 1580 { 1581 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1582 dynamic_cast< 1583 const parallel::distributed::Triangulation<dim, spacedim> *>( 1584 &dh.get_triangulation()); 1585 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1586 1587 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1588 vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); 1589 } 1590 1591 1592 1593 # ifdef DEAL_II_WITH_MPI 1594 # ifdef DEAL_II_TRILINOS_WITH_TPETRA 1595 template <int dim, int spacedim, typename Number> 1596 void 1597 reinit_distributed(const DoFHandler<dim, spacedim> & dh, 1598 LinearAlgebra::TpetraWrappers::Vector<Number> &vector) 1599 { 1600 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1601 dynamic_cast< 1602 const parallel::distributed::Triangulation<dim, spacedim> *>( 1603 &dh.get_triangulation()); 1604 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1605 1606 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1607 vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); 1608 } 1609 # endif 1610 1611 template <int dim, int spacedim> 1612 void 1613 reinit_distributed(const DoFHandler<dim, spacedim> & dh, 1614 LinearAlgebra::EpetraWrappers::Vector &vector) 1615 { 1616 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1617 dynamic_cast< 1618 const parallel::distributed::Triangulation<dim, spacedim> *>( 1619 &dh.get_triangulation()); 1620 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1621 1622 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1623 vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); 1624 } 1625 # endif 1626 #endif // DEAL_II_WITH_TRILINOS 1627 1628 template <int dim, int spacedim, typename Number> 1629 void 1630 reinit_distributed(const DoFHandler<dim, spacedim> & dh, 1631 LinearAlgebra::distributed::Vector<Number> &vector) 1632 { 1633 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1634 dynamic_cast< 1635 const parallel::distributed::Triangulation<dim, spacedim> *>( 1636 &dh.get_triangulation()); 1637 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1638 1639 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1640 vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); 1641 } 1642 1643 1644 1645 template <class VectorType, class DH> 1646 void 1647 reinit_ghosted(const DH & /*dh*/, VectorType & /*vector*/) 1648 { 1649 Assert(false, ExcNotImplemented()); 1650 } 1651 1652 #ifdef DEAL_II_WITH_PETSC 1653 template <int dim, int spacedim> 1654 void 1655 reinit_ghosted(const DoFHandler<dim, spacedim> &dh, 1656 PETScWrappers::MPI::Vector & vector) 1657 { 1658 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1659 dynamic_cast< 1660 const parallel::distributed::Triangulation<dim, spacedim> *>( 1661 &dh.get_triangulation()); 1662 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1663 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1664 IndexSet locally_relevant_dofs; 1665 DoFTools::extract_locally_relevant_dofs(dh, locally_relevant_dofs); 1666 vector.reinit(locally_owned_dofs, 1667 locally_relevant_dofs, 1668 parallel_tria->get_communicator()); 1669 } 1670 #endif // DEAL_II_WITH_PETSC 1671 1672 #ifdef DEAL_II_WITH_TRILINOS 1673 template <int dim, int spacedim> 1674 void 1675 reinit_ghosted(const DoFHandler<dim, spacedim> &dh, 1676 TrilinosWrappers::MPI::Vector & vector) 1677 { 1678 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1679 dynamic_cast< 1680 const parallel::distributed::Triangulation<dim, spacedim> *>( 1681 &dh.get_triangulation()); 1682 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1683 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1684 IndexSet locally_relevant_dofs; 1685 DoFTools::extract_locally_relevant_dofs(dh, locally_relevant_dofs); 1686 vector.reinit(locally_owned_dofs, 1687 locally_relevant_dofs, 1688 parallel_tria->get_communicator()); 1689 } 1690 #endif // DEAL_II_WITH_TRILINOS 1691 1692 template <int dim, int spacedim, typename Number> 1693 void 1694 reinit_ghosted(const DoFHandler<dim, spacedim> & dh, 1695 LinearAlgebra::distributed::Vector<Number> &vector) 1696 { 1697 const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria = 1698 dynamic_cast< 1699 const parallel::distributed::Triangulation<dim, spacedim> *>( 1700 &dh.get_triangulation()); 1701 Assert(parallel_tria != nullptr, ExcNotImplemented()); 1702 const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); 1703 IndexSet locally_relevant_dofs; 1704 DoFTools::extract_locally_relevant_dofs(dh, locally_relevant_dofs); 1705 vector.reinit(locally_owned_dofs, 1706 locally_relevant_dofs, 1707 parallel_tria->get_communicator()); 1708 } 1709 1710 1711 1712 template <int dim, class InVector, class OutVector, int spacedim> 1713 void 1714 extrapolate_serial(const InVector & u3, 1715 const DoFHandler<dim, spacedim> &dof2, 1716 OutVector & u2) 1717 { 1718 const unsigned int dofs_per_cell = dof2.get_fe().n_dofs_per_cell(); 1719 Vector<typename OutVector::value_type> dof_values(dofs_per_cell); 1720 1721 // then traverse grid bottom up 1722 for (unsigned int level = 0; 1723 level < dof2.get_triangulation().n_levels() - 1; 1724 ++level) 1725 { 1726 typename DoFHandler<dim, spacedim>::cell_iterator cell = 1727 dof2.begin(level), 1728 endc = 1729 dof2.end(level); 1730 1731 for (; cell != endc; ++cell) 1732 if (!cell->is_active()) 1733 { 1734 // check whether this 1735 // cell has active 1736 // children 1737 bool active_children = false; 1738 for (unsigned int child_n = 0; child_n < cell->n_children(); 1739 ++child_n) 1740 if (cell->child(child_n)->is_active()) 1741 { 1742 active_children = true; 1743 break; 1744 } 1745 1746 // if there are active 1747 // children, this process 1748 // has to work on this 1749 // cell. get the data 1750 // from the one vector 1751 // and set it on the 1752 // other 1753 if (active_children) 1754 { 1755 cell->get_interpolated_dof_values(u3, dof_values); 1756 cell->set_dof_values_by_interpolation(dof_values, u2); 1757 } 1758 } 1759 } 1760 } 1761 } // namespace internal 1762 1763 template <int dim, class InVector, class OutVector, int spacedim> 1764 void 1765 extrapolate(const DoFHandler<dim, spacedim> &dof1, 1766 const InVector & u1, 1767 const DoFHandler<dim, spacedim> &dof2, 1768 OutVector & u2) 1769 { 1770 AffineConstraints<typename OutVector::value_type> dummy; 1771 dummy.close(); 1772 extrapolate(dof1, u1, dof2, dummy, u2); 1773 } 1774 1775 1776 1777 template <int dim, class InVector, class OutVector, int spacedim> 1778 void 1779 extrapolate( 1780 const DoFHandler<dim, spacedim> & dof1, 1781 const InVector & u1, 1782 const DoFHandler<dim, spacedim> & dof2, 1783 const AffineConstraints<typename OutVector::value_type> &constraints, 1784 OutVector & u2) 1785 { 1786 Assert(dof1.get_fe(0).n_components() == dof2.get_fe(0).n_components(), 1787 ExcDimensionMismatch(dof1.get_fe(0).n_components(), 1788 dof2.get_fe(0).n_components())); 1789 Assert(&dof1.get_triangulation() == &dof2.get_triangulation(), 1790 ExcTriangulationMismatch()); 1791 Assert(u1.size() == dof1.n_dofs(), 1792 ExcDimensionMismatch(u1.size(), dof1.n_dofs())); 1793 Assert(u2.size() == dof2.n_dofs(), 1794 ExcDimensionMismatch(u2.size(), dof2.n_dofs())); 1795 1796 // make sure that each cell on the coarsest level is at least once refined, 1797 // otherwise, these cells can't be treated and would generate a bogus result 1798 { 1799 typename DoFHandler<dim, spacedim>::cell_iterator cell = dof2.begin(0), 1800 endc = dof2.end(0); 1801 for (; cell != endc; ++cell) 1802 Assert(cell->has_children() || cell->is_artificial(), 1803 ExcGridNotRefinedAtLeastOnce()); 1804 } 1805 1806 1807 internal::BlockType<OutVector> u3; 1808 internal::reinit_distributed(dof2, u3); 1809 if (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> 1810 *>(&dof2.get_triangulation()) != nullptr) 1811 { 1812 interpolate(dof1, u1, dof2, constraints, u3); 1813 1814 internal::BlockType<OutVector> u3_relevant; 1815 internal::reinit_ghosted(dof2, u3_relevant); 1816 u3_relevant = u3; 1817 1818 internal::ExtrapolateImplementation<dim, spacedim, OutVector> 1819 implementation; 1820 implementation.extrapolate_parallel(u3_relevant, dof2, u2); 1821 } 1822 else 1823 { 1824 interpolate(dof1, u1, dof2, constraints, u3); 1825 1826 internal::extrapolate_serial(u3, dof2, u2); 1827 } 1828 1829 constraints.distribute(u2); 1830 } 1831 } // end of namespace FETools 1832 1833 DEAL_II_NAMESPACE_CLOSE 1834 1835 /*-------------------- fe_tools_extrapolate_templates.h -------------------*/ 1836 #endif // dealii_fe_tools_extrapolate_templates_H 1837