1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2008 - 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/distributed/p4est_wrappers.h> 18 #include <deal.II/distributed/tria.h> 19 20 DEAL_II_NAMESPACE_OPEN 21 22 #ifdef DEAL_II_WITH_P4EST 23 24 namespace internal 25 { 26 namespace p4est 27 { 28 namespace 29 { 30 template <int dim, int spacedim> 31 typename dealii::Triangulation<dim, spacedim>::cell_iterator cell_from_quad(const dealii::parallel::distributed::Triangulation<dim,spacedim> * triangulation,const typename dealii::internal::p4est::types<dim>::topidx treeidx,const typename dealii::internal::p4est::types<dim>::quadrant & quad)32 cell_from_quad( 33 const dealii::parallel::distributed::Triangulation<dim, spacedim> 34 *triangulation, 35 const typename dealii::internal::p4est::types<dim>::topidx treeidx, 36 const typename dealii::internal::p4est::types<dim>::quadrant &quad) 37 { 38 int i, l = quad.level; 39 dealii::types::global_dof_index dealii_index = 40 triangulation->get_p4est_tree_to_coarse_cell_permutation()[treeidx]; 41 42 for (i = 0; i < l; i++) 43 { 44 typename dealii::Triangulation<dim, spacedim>::cell_iterator cell( 45 triangulation, i, dealii_index); 46 const int child_id = 47 dealii::internal::p4est::functions<dim>::quadrant_ancestor_id( 48 &quad, i + 1); 49 Assert(cell->has_children(), 50 ExcMessage("p4est quadrant does not correspond to a cell!")); 51 dealii_index = cell->child_index(child_id); 52 } 53 54 typename dealii::Triangulation<dim, spacedim>::cell_iterator out_cell( 55 triangulation, l, dealii_index); 56 57 return out_cell; 58 } 59 60 /** 61 * This is the callback data structure used to fill 62 * vertices_with_ghost_neighbors via the p4est_iterate tool 63 */ 64 template <int dim, int spacedim> 65 struct FindGhosts 66 { 67 const typename dealii::parallel::distributed::Triangulation<dim, 68 spacedim> 69 * triangulation; 70 sc_array_t *subids; 71 std::map<unsigned int, std::set<dealii::types::subdomain_id>> 72 *vertices_with_ghost_neighbors; 73 }; 74 75 76 /** At a corner (vertex), determine if any of the neighboring cells are 77 * ghosts. If there are, find out their subdomain ids, and if this is a 78 * local vertex, then add these subdomain ids to the map 79 * vertices_with_ghost_neighbors of that index 80 */ 81 template <int dim, int spacedim> 82 void find_ghosts_corner(typename dealii::internal::p4est::iter<dim>::corner_info * info,void * user_data)83 find_ghosts_corner( 84 typename dealii::internal::p4est::iter<dim>::corner_info *info, 85 void * user_data) 86 { 87 int i, j; 88 int nsides = info->sides.elem_count; 89 auto *sides = reinterpret_cast< 90 typename dealii::internal::p4est::iter<dim>::corner_side *>( 91 info->sides.array); 92 FindGhosts<dim, spacedim> *fg = 93 static_cast<FindGhosts<dim, spacedim> *>(user_data); 94 sc_array_t *subids = fg->subids; 95 const dealii::parallel::distributed::Triangulation<dim, spacedim> 96 * triangulation = fg->triangulation; 97 int nsubs; 98 dealii::types::subdomain_id *subdomain_ids; 99 std::map<unsigned int, std::set<dealii::types::subdomain_id>> 100 *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; 101 102 subids->elem_count = 0; 103 for (i = 0; i < nsides; i++) 104 { 105 if (sides[i].is_ghost) 106 { 107 typename dealii::parallel::distributed:: 108 Triangulation<dim, spacedim>::cell_iterator cell = 109 cell_from_quad(triangulation, 110 sides[i].treeid, 111 *(sides[i].quad)); 112 Assert(cell->is_ghost(), 113 ExcMessage("ghost quad did not find ghost cell")); 114 dealii::types::subdomain_id *subid = 115 static_cast<dealii::types::subdomain_id *>( 116 sc_array_push(subids)); 117 *subid = cell->subdomain_id(); 118 } 119 } 120 121 if (!subids->elem_count) 122 { 123 return; 124 } 125 126 nsubs = static_cast<int>(subids->elem_count); 127 subdomain_ids = 128 reinterpret_cast<dealii::types::subdomain_id *>(subids->array); 129 130 for (i = 0; i < nsides; i++) 131 { 132 if (!sides[i].is_ghost) 133 { 134 typename dealii::parallel::distributed:: 135 Triangulation<dim, spacedim>::cell_iterator cell = 136 cell_from_quad(triangulation, 137 sides[i].treeid, 138 *(sides[i].quad)); 139 140 Assert(!cell->is_ghost(), 141 ExcMessage("local quad found ghost cell")); 142 143 for (j = 0; j < nsubs; j++) 144 { 145 (*vertices_with_ghost_neighbors)[cell->vertex_index( 146 sides[i].corner)] 147 .insert(subdomain_ids[j]); 148 } 149 } 150 } 151 152 subids->elem_count = 0; 153 } 154 155 /** Similar to find_ghosts_corner, but for the hanging vertex in the 156 * middle of an edge 157 */ 158 template <int dim, int spacedim> 159 void find_ghosts_edge(typename dealii::internal::p4est::iter<dim>::edge_info * info,void * user_data)160 find_ghosts_edge( 161 typename dealii::internal::p4est::iter<dim>::edge_info *info, 162 void * user_data) 163 { 164 int i, j, k; 165 int nsides = info->sides.elem_count; 166 auto *sides = reinterpret_cast< 167 typename dealii::internal::p4est::iter<dim>::edge_side *>( 168 info->sides.array); 169 auto * fg = static_cast<FindGhosts<dim, spacedim> *>(user_data); 170 sc_array_t *subids = fg->subids; 171 const dealii::parallel::distributed::Triangulation<dim, spacedim> 172 * triangulation = fg->triangulation; 173 int nsubs; 174 dealii::types::subdomain_id *subdomain_ids; 175 std::map<unsigned int, std::set<dealii::types::subdomain_id>> 176 *vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; 177 178 subids->elem_count = 0; 179 for (i = 0; i < nsides; i++) 180 { 181 if (sides[i].is_hanging) 182 { 183 for (j = 0; j < 2; j++) 184 { 185 if (sides[i].is.hanging.is_ghost[j]) 186 { 187 typename dealii::parallel::distributed:: 188 Triangulation<dim, spacedim>::cell_iterator cell = 189 cell_from_quad(triangulation, 190 sides[i].treeid, 191 *(sides[i].is.hanging.quad[j])); 192 dealii::types::subdomain_id *subid = 193 static_cast<dealii::types::subdomain_id *>( 194 sc_array_push(subids)); 195 *subid = cell->subdomain_id(); 196 } 197 } 198 } 199 } 200 201 if (!subids->elem_count) 202 { 203 return; 204 } 205 206 nsubs = static_cast<int>(subids->elem_count); 207 subdomain_ids = 208 reinterpret_cast<dealii::types::subdomain_id *>(subids->array); 209 210 for (i = 0; i < nsides; i++) 211 { 212 if (sides[i].is_hanging) 213 { 214 for (j = 0; j < 2; j++) 215 { 216 if (!sides[i].is.hanging.is_ghost[j]) 217 { 218 typename dealii::parallel::distributed:: 219 Triangulation<dim, spacedim>::cell_iterator cell = 220 cell_from_quad(triangulation, 221 sides[i].treeid, 222 *(sides[i].is.hanging.quad[j])); 223 224 for (k = 0; k < nsubs; k++) 225 { 226 (*vertices_with_ghost_neighbors) 227 [cell->vertex_index( 228 p8est_edge_corners[sides[i].edge][1 ^ j])] 229 .insert(subdomain_ids[k]); 230 } 231 } 232 } 233 } 234 } 235 236 subids->elem_count = 0; 237 } 238 239 /** Similar to find_ghosts_corner, but for the hanging vertex in the 240 * middle of a face 241 */ 242 template <int dim, int spacedim> 243 void find_ghosts_face(typename dealii::internal::p4est::iter<dim>::face_info * info,void * user_data)244 find_ghosts_face( 245 typename dealii::internal::p4est::iter<dim>::face_info *info, 246 void * user_data) 247 { 248 int i, j, k; 249 int nsides = info->sides.elem_count; 250 auto *sides = reinterpret_cast< 251 typename dealii::internal::p4est::iter<dim>::face_side *>( 252 info->sides.array); 253 FindGhosts<dim, spacedim> *fg = 254 static_cast<FindGhosts<dim, spacedim> *>(user_data); 255 sc_array_t *subids = fg->subids; 256 const dealii::parallel::distributed::Triangulation<dim, spacedim> 257 * triangulation = fg->triangulation; 258 int nsubs; 259 dealii::types::subdomain_id *subdomain_ids; 260 std::map<unsigned int, std::set<dealii::types::subdomain_id>> 261 * vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors; 262 int limit = (dim == 2) ? 2 : 4; 263 264 subids->elem_count = 0; 265 for (i = 0; i < nsides; i++) 266 { 267 if (sides[i].is_hanging) 268 { 269 for (j = 0; j < limit; j++) 270 { 271 if (sides[i].is.hanging.is_ghost[j]) 272 { 273 typename dealii::parallel::distributed:: 274 Triangulation<dim, spacedim>::cell_iterator cell = 275 cell_from_quad(triangulation, 276 sides[i].treeid, 277 *(sides[i].is.hanging.quad[j])); 278 dealii::types::subdomain_id *subid = 279 static_cast<dealii::types::subdomain_id *>( 280 sc_array_push(subids)); 281 *subid = cell->subdomain_id(); 282 } 283 } 284 } 285 } 286 287 if (!subids->elem_count) 288 { 289 return; 290 } 291 292 nsubs = static_cast<int>(subids->elem_count); 293 subdomain_ids = 294 reinterpret_cast<dealii::types::subdomain_id *>(subids->array); 295 296 for (i = 0; i < nsides; i++) 297 { 298 if (sides[i].is_hanging) 299 { 300 for (j = 0; j < limit; j++) 301 { 302 if (!sides[i].is.hanging.is_ghost[j]) 303 { 304 typename dealii::parallel::distributed:: 305 Triangulation<dim, spacedim>::cell_iterator cell = 306 cell_from_quad(triangulation, 307 sides[i].treeid, 308 *(sides[i].is.hanging.quad[j])); 309 310 for (k = 0; k < nsubs; k++) 311 { 312 if (dim == 2) 313 { 314 (*vertices_with_ghost_neighbors) 315 [cell->vertex_index( 316 p4est_face_corners[sides[i].face] 317 [(limit - 1) ^ j])] 318 .insert(subdomain_ids[k]); 319 } 320 else 321 { 322 (*vertices_with_ghost_neighbors) 323 [cell->vertex_index( 324 p8est_face_corners[sides[i].face] 325 [(limit - 1) ^ j])] 326 .insert(subdomain_ids[k]); 327 } 328 } 329 } 330 } 331 } 332 } 333 334 subids->elem_count = 0; 335 } 336 } // namespace 337 338 339 int (&functions<2>::quadrant_compare)(const void *v1, const void *v2) = 340 p4est_quadrant_compare; 341 342 void (&functions<2>::quadrant_childrenv)(const types<2>::quadrant *q, 343 types<2>::quadrant c[]) = 344 p4est_quadrant_childrenv; 345 346 int (&functions<2>::quadrant_overlaps_tree)(types<2>::tree * tree, 347 const types<2>::quadrant *q) = 348 p4est_quadrant_overlaps_tree; 349 350 void (&functions<2>::quadrant_set_morton)(types<2>::quadrant *quadrant, 351 int level, 352 uint64_t id) = 353 p4est_quadrant_set_morton; 354 355 int (&functions<2>::quadrant_is_equal)(const types<2>::quadrant *q1, 356 const types<2>::quadrant *q2) = 357 p4est_quadrant_is_equal; 358 359 int (&functions<2>::quadrant_is_sibling)(const types<2>::quadrant *q1, 360 const types<2>::quadrant *q2) = 361 p4est_quadrant_is_sibling; 362 363 int (&functions<2>::quadrant_is_ancestor)(const types<2>::quadrant *q1, 364 const types<2>::quadrant *q2) = 365 p4est_quadrant_is_ancestor; 366 367 int (&functions<2>::quadrant_ancestor_id)(const types<2>::quadrant *q, 368 int level) = 369 p4est_quadrant_ancestor_id; 370 371 int (&functions<2>::comm_find_owner)(types<2>::forest * p4est, 372 const types<2>::locidx which_tree, 373 const types<2>::quadrant *q, 374 const int guess) = 375 p4est_comm_find_owner; 376 377 types<2>::connectivity *(&functions<2>::connectivity_new)( 378 types<2>::topidx num_vertices, 379 types<2>::topidx num_trees, 380 types<2>::topidx num_corners, 381 types<2>::topidx num_vtt) = p4est_connectivity_new; 382 383 void (&functions<2>::connectivity_join_faces)(types<2>::connectivity *conn, 384 types<2>::topidx tree_left, 385 types<2>::topidx tree_right, 386 int face_left, 387 int face_right, 388 int orientation) = 389 p4est_connectivity_join_faces; 390 391 void (&functions<2>::connectivity_destroy)( 392 p4est_connectivity_t *connectivity) = p4est_connectivity_destroy; 393 394 types<2>::forest *(&functions<2>::new_forest)( 395 MPI_Comm mpicomm, 396 types<2>::connectivity *connectivity, 397 types<2>::locidx min_quadrants, 398 int min_level, 399 int fill_uniform, 400 std::size_t data_size, 401 p4est_init_t init_fn, 402 void * user_pointer) = p4est_new_ext; 403 404 void (&functions<2>::destroy)(types<2>::forest *p4est) = p4est_destroy; 405 406 void (&functions<2>::refine)(types<2>::forest *p4est, 407 int refine_recursive, 408 p4est_refine_t refine_fn, 409 p4est_init_t init_fn) = p4est_refine; 410 411 void (&functions<2>::coarsen)(types<2>::forest *p4est, 412 int coarsen_recursive, 413 p4est_coarsen_t coarsen_fn, 414 p4est_init_t init_fn) = p4est_coarsen; 415 416 void (&functions<2>::balance)(types<2>::forest * p4est, 417 types<2>::balance_type btype, 418 p4est_init_t init_fn) = p4est_balance; 419 420 types<2>::gloidx (&functions<2>::partition)(types<2>::forest *p4est, 421 int partition_for_coarsening, 422 p4est_weight_t weight_fn) = 423 p4est_partition_ext; 424 425 void (&functions<2>::save)(const char * filename, 426 types<2>::forest *p4est, 427 int save_data) = p4est_save; 428 429 types<2>::forest *(&functions<2>::load_ext)( 430 const char * filename, 431 MPI_Comm mpicomm, 432 std::size_t data_size, 433 int load_data, 434 int autopartition, 435 int broadcasthead, 436 void * user_pointer, 437 types<2>::connectivity **p4est) = p4est_load_ext; 438 439 int (&functions<2>::connectivity_save)( 440 const char * filename, 441 types<2>::connectivity *connectivity) = p4est_connectivity_save; 442 443 int (&functions<2>::connectivity_is_valid)( 444 types<2>::connectivity *connectivity) = p4est_connectivity_is_valid; 445 446 types<2>::connectivity *(&functions<2>::connectivity_load)( 447 const char * filename, 448 std::size_t *length) = p4est_connectivity_load; 449 450 unsigned int (&functions<2>::checksum)(types<2>::forest *p4est) = 451 p4est_checksum; 452 453 void (&functions<2>::vtk_write_file)(types<2>::forest *p4est, 454 p4est_geometry_t *, 455 const char *baseName) = 456 p4est_vtk_write_file; 457 458 types<2>::ghost *(&functions<2>::ghost_new)(types<2>::forest * p4est, 459 types<2>::balance_type btype) = 460 p4est_ghost_new; 461 462 void (&functions<2>::ghost_destroy)(types<2>::ghost *ghost) = 463 p4est_ghost_destroy; 464 465 void (&functions<2>::reset_data)(types<2>::forest *p4est, 466 std::size_t data_size, 467 p4est_init_t init_fn, 468 void *user_pointer) = p4est_reset_data; 469 470 std::size_t (&functions<2>::forest_memory_used)(types<2>::forest *p4est) = 471 p4est_memory_used; 472 473 std::size_t (&functions<2>::connectivity_memory_used)( 474 types<2>::connectivity *p4est) = p4est_connectivity_memory_used; 475 476 constexpr unsigned int functions<2>::max_level; 477 478 void (&functions<2>::transfer_fixed)(const types<2>::gloidx *dest_gfq, 479 const types<2>::gloidx *src_gfq, 480 MPI_Comm mpicomm, 481 int tag, 482 void * dest_data, 483 const void * src_data, 484 std::size_t data_size) = 485 p4est_transfer_fixed; 486 487 types<2>::transfer_context *(&functions<2>::transfer_fixed_begin)( 488 const types<2>::gloidx *dest_gfq, 489 const types<2>::gloidx *src_gfq, 490 MPI_Comm mpicomm, 491 int tag, 492 void * dest_data, 493 const void * src_data, 494 std::size_t data_size) = p4est_transfer_fixed_begin; 495 496 void (&functions<2>::transfer_fixed_end)(types<2>::transfer_context *tc) = 497 p4est_transfer_fixed_end; 498 499 void (&functions<2>::transfer_custom)(const types<2>::gloidx *dest_gfq, 500 const types<2>::gloidx *src_gfq, 501 MPI_Comm mpicomm, 502 int tag, 503 void * dest_data, 504 const int * dest_sizes, 505 const void * src_data, 506 const int * src_sizes) = 507 p4est_transfer_custom; 508 509 types<2>::transfer_context *(&functions<2>::transfer_custom_begin)( 510 const types<2>::gloidx *dest_gfq, 511 const types<2>::gloidx *src_gfq, 512 MPI_Comm mpicomm, 513 int tag, 514 void * dest_data, 515 const int * dest_sizes, 516 const void * src_data, 517 const int * src_sizes) = p4est_transfer_custom_begin; 518 519 void (&functions<2>::transfer_custom_end)(types<2>::transfer_context *tc) = 520 p4est_transfer_custom_end; 521 522 523 524 int (&functions<3>::quadrant_compare)(const void *v1, const void *v2) = 525 p8est_quadrant_compare; 526 527 void (&functions<3>::quadrant_childrenv)(const types<3>::quadrant *q, 528 types<3>::quadrant c[]) = 529 p8est_quadrant_childrenv; 530 531 int (&functions<3>::quadrant_overlaps_tree)(types<3>::tree * tree, 532 const types<3>::quadrant *q) = 533 p8est_quadrant_overlaps_tree; 534 535 void (&functions<3>::quadrant_set_morton)(types<3>::quadrant *quadrant, 536 int level, 537 uint64_t id) = 538 p8est_quadrant_set_morton; 539 540 int (&functions<3>::quadrant_is_equal)(const types<3>::quadrant *q1, 541 const types<3>::quadrant *q2) = 542 p8est_quadrant_is_equal; 543 544 int (&functions<3>::quadrant_is_sibling)(const types<3>::quadrant *q1, 545 const types<3>::quadrant *q2) = 546 p8est_quadrant_is_sibling; 547 548 int (&functions<3>::quadrant_is_ancestor)(const types<3>::quadrant *q1, 549 const types<3>::quadrant *q2) = 550 p8est_quadrant_is_ancestor; 551 552 int (&functions<3>::quadrant_ancestor_id)(const types<3>::quadrant *q, 553 int level) = 554 p8est_quadrant_ancestor_id; 555 556 int (&functions<3>::comm_find_owner)(types<3>::forest * p4est, 557 const types<3>::locidx which_tree, 558 const types<3>::quadrant *q, 559 const int guess) = 560 p8est_comm_find_owner; 561 562 types<3>::connectivity *(&functions<3>::connectivity_new)( 563 types<3>::topidx num_vertices, 564 types<3>::topidx num_trees, 565 types<3>::topidx num_edges, 566 types<3>::topidx num_ett, 567 types<3>::topidx num_corners, 568 types<3>::topidx num_ctt) = p8est_connectivity_new; 569 570 void (&functions<3>::connectivity_destroy)( 571 p8est_connectivity_t *connectivity) = p8est_connectivity_destroy; 572 573 void (&functions<3>::connectivity_join_faces)(types<3>::connectivity *conn, 574 types<3>::topidx tree_left, 575 types<3>::topidx tree_right, 576 int face_left, 577 int face_right, 578 int orientation) = 579 p8est_connectivity_join_faces; 580 581 types<3>::forest *(&functions<3>::new_forest)( 582 MPI_Comm mpicomm, 583 types<3>::connectivity *connectivity, 584 types<3>::locidx min_quadrants, 585 int min_level, 586 int fill_uniform, 587 std::size_t data_size, 588 p8est_init_t init_fn, 589 void * user_pointer) = p8est_new_ext; 590 591 void (&functions<3>::destroy)(types<3>::forest *p8est) = p8est_destroy; 592 593 void (&functions<3>::refine)(types<3>::forest *p8est, 594 int refine_recursive, 595 p8est_refine_t refine_fn, 596 p8est_init_t init_fn) = p8est_refine; 597 598 void (&functions<3>::coarsen)(types<3>::forest *p8est, 599 int coarsen_recursive, 600 p8est_coarsen_t coarsen_fn, 601 p8est_init_t init_fn) = p8est_coarsen; 602 603 void (&functions<3>::balance)(types<3>::forest * p8est, 604 types<3>::balance_type btype, 605 p8est_init_t init_fn) = p8est_balance; 606 607 types<3>::gloidx (&functions<3>::partition)(types<3>::forest *p8est, 608 int partition_for_coarsening, 609 p8est_weight_t weight_fn) = 610 p8est_partition_ext; 611 612 void (&functions<3>::save)(const char * filename, 613 types<3>::forest *p4est, 614 int save_data) = p8est_save; 615 616 types<3>::forest *(&functions<3>::load_ext)( 617 const char * filename, 618 MPI_Comm mpicomm, 619 std::size_t data_size, 620 int load_data, 621 int autopartition, 622 int broadcasthead, 623 void * user_pointer, 624 types<3>::connectivity **p4est) = p8est_load_ext; 625 626 int (&functions<3>::connectivity_save)( 627 const char * filename, 628 types<3>::connectivity *connectivity) = p8est_connectivity_save; 629 630 int (&functions<3>::connectivity_is_valid)( 631 types<3>::connectivity *connectivity) = p8est_connectivity_is_valid; 632 633 types<3>::connectivity *(&functions<3>::connectivity_load)( 634 const char * filename, 635 std::size_t *length) = p8est_connectivity_load; 636 637 unsigned int (&functions<3>::checksum)(types<3>::forest *p8est) = 638 p8est_checksum; 639 640 void (&functions<3>::vtk_write_file)(types<3>::forest *p8est, 641 p8est_geometry_t *, 642 const char *baseName) = 643 p8est_vtk_write_file; 644 645 types<3>::ghost *(&functions<3>::ghost_new)(types<3>::forest * p4est, 646 types<3>::balance_type btype) = 647 p8est_ghost_new; 648 649 void (&functions<3>::ghost_destroy)(types<3>::ghost *ghost) = 650 p8est_ghost_destroy; 651 652 void (&functions<3>::reset_data)(types<3>::forest *p4est, 653 std::size_t data_size, 654 p8est_init_t init_fn, 655 void *user_pointer) = p8est_reset_data; 656 657 std::size_t (&functions<3>::forest_memory_used)(types<3>::forest *p4est) = 658 p8est_memory_used; 659 660 std::size_t (&functions<3>::connectivity_memory_used)( 661 types<3>::connectivity *p4est) = p8est_connectivity_memory_used; 662 663 constexpr unsigned int functions<3>::max_level; 664 665 void (&functions<3>::transfer_fixed)(const types<3>::gloidx *dest_gfq, 666 const types<3>::gloidx *src_gfq, 667 MPI_Comm mpicomm, 668 int tag, 669 void * dest_data, 670 const void * src_data, 671 std::size_t data_size) = 672 p8est_transfer_fixed; 673 674 types<3>::transfer_context *(&functions<3>::transfer_fixed_begin)( 675 const types<3>::gloidx *dest_gfq, 676 const types<3>::gloidx *src_gfq, 677 MPI_Comm mpicomm, 678 int tag, 679 void * dest_data, 680 const void * src_data, 681 std::size_t data_size) = p8est_transfer_fixed_begin; 682 683 void (&functions<3>::transfer_fixed_end)(types<3>::transfer_context *tc) = 684 p8est_transfer_fixed_end; 685 686 void (&functions<3>::transfer_custom)(const types<3>::gloidx *dest_gfq, 687 const types<3>::gloidx *src_gfq, 688 MPI_Comm mpicomm, 689 int tag, 690 void * dest_data, 691 const int * dest_sizes, 692 const void * src_data, 693 const int * src_sizes) = 694 p8est_transfer_custom; 695 696 types<3>::transfer_context *(&functions<3>::transfer_custom_begin)( 697 const types<3>::gloidx *dest_gfq, 698 const types<3>::gloidx *src_gfq, 699 MPI_Comm mpicomm, 700 int tag, 701 void * dest_data, 702 const int * dest_sizes, 703 const void * src_data, 704 const int * src_sizes) = p8est_transfer_custom_begin; 705 706 void (&functions<3>::transfer_custom_end)(types<3>::transfer_context *tc) = 707 p8est_transfer_custom_end; 708 709 710 711 template <int dim> 712 void init_quadrant_children(const typename types<dim>::quadrant & p4est_cell,typename types<dim>::quadrant (& p4est_children)[dealii::GeometryInfo<dim>::max_children_per_cell])713 init_quadrant_children( 714 const typename types<dim>::quadrant &p4est_cell, 715 typename types<dim>::quadrant ( 716 &p4est_children)[dealii::GeometryInfo<dim>::max_children_per_cell]) 717 { 718 for (unsigned int c = 0; 719 c < dealii::GeometryInfo<dim>::max_children_per_cell; 720 ++c) 721 switch (dim) 722 { 723 case 2: 724 P4EST_QUADRANT_INIT(&p4est_children[c]); 725 break; 726 case 3: 727 P8EST_QUADRANT_INIT(&p4est_children[c]); 728 break; 729 default: 730 Assert(false, ExcNotImplemented()); 731 } 732 733 734 functions<dim>::quadrant_childrenv(&p4est_cell, p4est_children); 735 } 736 737 template <int dim> 738 void init_coarse_quadrant(typename types<dim>::quadrant & quad)739 init_coarse_quadrant(typename types<dim>::quadrant &quad) 740 { 741 switch (dim) 742 { 743 case 2: 744 P4EST_QUADRANT_INIT(&quad); 745 break; 746 case 3: 747 P8EST_QUADRANT_INIT(&quad); 748 break; 749 default: 750 Assert(false, ExcNotImplemented()); 751 } 752 functions<dim>::quadrant_set_morton(&quad, 753 /*level=*/0, 754 /*index=*/0); 755 } 756 757 template <int dim> 758 bool quadrant_is_equal(const typename types<dim>::quadrant & q1,const typename types<dim>::quadrant & q2)759 quadrant_is_equal(const typename types<dim>::quadrant &q1, 760 const typename types<dim>::quadrant &q2) 761 { 762 return functions<dim>::quadrant_is_equal(&q1, &q2); 763 } 764 765 766 767 template <int dim> 768 bool quadrant_is_ancestor(const typename types<dim>::quadrant & q1,const typename types<dim>::quadrant & q2)769 quadrant_is_ancestor(const typename types<dim>::quadrant &q1, 770 const typename types<dim>::quadrant &q2) 771 { 772 return functions<dim>::quadrant_is_ancestor(&q1, &q2); 773 } 774 775 template <int dim> 776 bool tree_exists_locally(const typename types<dim>::forest * parallel_forest,const typename types<dim>::topidx coarse_grid_cell)777 tree_exists_locally(const typename types<dim>::forest *parallel_forest, 778 const typename types<dim>::topidx coarse_grid_cell) 779 { 780 Assert(coarse_grid_cell < parallel_forest->connectivity->num_trees, 781 ExcInternalError()); 782 return ((coarse_grid_cell >= parallel_forest->first_local_tree) && 783 (coarse_grid_cell <= parallel_forest->last_local_tree)); 784 } 785 786 787 788 template <> 789 bool quadrant_is_equal(const typename types<1>::quadrant & q1,const typename types<1>::quadrant & q2)790 quadrant_is_equal<1>(const typename types<1>::quadrant &q1, 791 const typename types<1>::quadrant &q2) 792 { 793 return q1 == q2; 794 } 795 796 797 798 template <> quadrant_is_ancestor(types<1>::quadrant const & q1,types<1>::quadrant const & q2)799 bool quadrant_is_ancestor<1>(types<1>::quadrant const &q1, 800 types<1>::quadrant const &q2) 801 { 802 // determine level of quadrants 803 const int level_1 = (q1 << types<1>::max_n_child_indices_bits) >> 804 types<1>::max_n_child_indices_bits; 805 const int level_2 = (q2 << types<1>::max_n_child_indices_bits) >> 806 types<1>::max_n_child_indices_bits; 807 808 // q1 can be an ancestor of q2 if q1's level is smaller 809 if (level_1 >= level_2) 810 return false; 811 812 // extract path of quadrants up to level of possible ancestor q1 813 const int truncated_id_1 = (q1 >> (types<1>::n_bits - 1 - level_1)) 814 << (types<1>::n_bits - 1 - level_1); 815 const int truncated_id_2 = (q2 >> (types<1>::n_bits - 1 - level_1)) 816 << (types<1>::n_bits - 1 - level_1); 817 818 // compare paths 819 return truncated_id_1 == truncated_id_2; 820 } 821 822 823 824 template <> 825 void init_quadrant_children(const typename types<1>::quadrant & q,typename types<1>::quadrant (& p4est_children)[dealii::GeometryInfo<1>::max_children_per_cell])826 init_quadrant_children<1>( 827 const typename types<1>::quadrant &q, 828 typename types<1>::quadrant ( 829 &p4est_children)[dealii::GeometryInfo<1>::max_children_per_cell]) 830 { 831 // determine the current level of quadrant 832 const int level_parent = (q << types<1>::max_n_child_indices_bits) >> 833 types<1>::max_n_child_indices_bits; 834 const int level_child = level_parent + 1; 835 836 // left child: only n_child_indices has to be incremented 837 p4est_children[0] = (q + 1); 838 839 // right child: increment and set a bit to 1 indicating that it is a right 840 // child 841 p4est_children[1] = (q + 1) | (1 << (types<1>::n_bits - 1 - level_child)); 842 } 843 844 845 846 template <> init_coarse_quadrant(typename types<1>::quadrant & quad)847 void init_coarse_quadrant<1>(typename types<1>::quadrant &quad) 848 { 849 quad = 0; 850 } 851 852 } // namespace p4est 853 } // namespace internal 854 855 #endif // DEAL_II_WITH_P4EST 856 857 /*-------------- Explicit Instantiations -------------------------------*/ 858 #include "p4est_wrappers.inst" 859 860 861 DEAL_II_NAMESPACE_CLOSE 862