1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2001 - 2020 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 #include <deal.II/base/geometry_info.h> 17 #include <deal.II/base/point.h> 18 #include <deal.II/base/tensor.h> 19 20 #include <deal.II/distributed/fully_distributed_tria.h> 21 #include <deal.II/distributed/shared_tria.h> 22 #include <deal.II/distributed/tria.h> 23 24 #include <deal.II/dofs/dof_accessor.h> 25 #include <deal.II/dofs/dof_handler.h> 26 27 #include <deal.II/fe/mapping_q1.h> 28 #include <deal.II/fe/mapping_q_generic.h> 29 30 #include <deal.II/grid/filtered_iterator.h> 31 #include <deal.II/grid/grid_tools.h> 32 #include <deal.II/grid/tria.h> 33 #include <deal.II/grid/tria_accessor.h> 34 #include <deal.II/grid/tria_iterator.h> 35 36 #include <deal.II/hp/dof_handler.h> 37 #include <deal.II/hp/mapping_collection.h> 38 39 #include <deal.II/lac/full_matrix.h> 40 41 #include <algorithm> 42 #include <array> 43 #include <cmath> 44 #include <list> 45 #include <map> 46 #include <numeric> 47 #include <set> 48 #include <vector> 49 50 51 DEAL_II_NAMESPACE_OPEN 52 53 namespace GridTools 54 { 55 template <int dim, template <int, int> class MeshType, int spacedim> 56 unsigned int find_closest_vertex(const MeshType<dim,spacedim> & mesh,const Point<spacedim> & p,const std::vector<bool> & marked_vertices)57 find_closest_vertex(const MeshType<dim, spacedim> &mesh, 58 const Point<spacedim> & p, 59 const std::vector<bool> & marked_vertices) 60 { 61 // first get the underlying 62 // triangulation from the 63 // mesh and determine vertices 64 // and used vertices 65 const Triangulation<dim, spacedim> &tria = mesh.get_triangulation(); 66 67 const std::vector<Point<spacedim>> &vertices = tria.get_vertices(); 68 69 Assert(tria.get_vertices().size() == marked_vertices.size() || 70 marked_vertices.size() == 0, 71 ExcDimensionMismatch(tria.get_vertices().size(), 72 marked_vertices.size())); 73 74 // If p is an element of marked_vertices, 75 // and q is that of used_Vertices, 76 // the vector marked_vertices does NOT 77 // contain unused vertices if p implies q. 78 // I.e., if p is true q must be true 79 // (if p is false, q could be false or true). 80 // p implies q logic is encapsulated in ~p|q. 81 Assert( 82 marked_vertices.size() == 0 || 83 std::equal(marked_vertices.begin(), 84 marked_vertices.end(), 85 tria.get_used_vertices().begin(), 86 [](bool p, bool q) { return !p || q; }), 87 ExcMessage( 88 "marked_vertices should be a subset of used vertices in the triangulation " 89 "but marked_vertices contains one or more vertices that are not used vertices!")); 90 91 // In addition, if a vector bools 92 // is specified (marked_vertices) 93 // marking all the vertices which 94 // could be the potentially closest 95 // vertex to the point, use it instead 96 // of used vertices 97 const std::vector<bool> &used = (marked_vertices.size() == 0) ? 98 tria.get_used_vertices() : 99 marked_vertices; 100 101 // At the beginning, the first 102 // used vertex is the closest one 103 std::vector<bool>::const_iterator first = 104 std::find(used.begin(), used.end(), true); 105 106 // Assert that at least one vertex 107 // is actually used 108 Assert(first != used.end(), ExcInternalError()); 109 110 unsigned int best_vertex = std::distance(used.begin(), first); 111 double best_dist = (p - vertices[best_vertex]).norm_square(); 112 113 // For all remaining vertices, test 114 // whether they are any closer 115 for (unsigned int j = best_vertex + 1; j < vertices.size(); j++) 116 if (used[j]) 117 { 118 double dist = (p - vertices[j]).norm_square(); 119 if (dist < best_dist) 120 { 121 best_vertex = j; 122 best_dist = dist; 123 } 124 } 125 126 return best_vertex; 127 } 128 129 130 131 template <int dim, template <int, int> class MeshType, int spacedim> 132 unsigned int find_closest_vertex(const Mapping<dim,spacedim> & mapping,const MeshType<dim,spacedim> & mesh,const Point<spacedim> & p,const std::vector<bool> & marked_vertices)133 find_closest_vertex(const Mapping<dim, spacedim> & mapping, 134 const MeshType<dim, spacedim> &mesh, 135 const Point<spacedim> & p, 136 const std::vector<bool> & marked_vertices) 137 { 138 // Take a shortcut in the simple case. 139 if (mapping.preserves_vertex_locations() == true) 140 return find_closest_vertex(mesh, p, marked_vertices); 141 142 // first get the underlying 143 // triangulation from the 144 // mesh and determine vertices 145 // and used vertices 146 const Triangulation<dim, spacedim> &tria = mesh.get_triangulation(); 147 148 auto vertices = extract_used_vertices(tria, mapping); 149 150 Assert(tria.get_vertices().size() == marked_vertices.size() || 151 marked_vertices.size() == 0, 152 ExcDimensionMismatch(tria.get_vertices().size(), 153 marked_vertices.size())); 154 155 // If p is an element of marked_vertices, 156 // and q is that of used_Vertices, 157 // the vector marked_vertices does NOT 158 // contain unused vertices if p implies q. 159 // I.e., if p is true q must be true 160 // (if p is false, q could be false or true). 161 // p implies q logic is encapsulated in ~p|q. 162 Assert( 163 marked_vertices.size() == 0 || 164 std::equal(marked_vertices.begin(), 165 marked_vertices.end(), 166 tria.get_used_vertices().begin(), 167 [](bool p, bool q) { return !p || q; }), 168 ExcMessage( 169 "marked_vertices should be a subset of used vertices in the triangulation " 170 "but marked_vertices contains one or more vertices that are not used vertices!")); 171 172 // Remove from the map unwanted elements. 173 if (marked_vertices.size()) 174 for (auto it = vertices.begin(); it != vertices.end();) 175 { 176 if (marked_vertices[it->first] == false) 177 { 178 vertices.erase(it++); 179 } 180 else 181 { 182 ++it; 183 } 184 } 185 186 return find_closest_vertex(vertices, p); 187 } 188 189 190 191 template <int dim, template <int, int> class MeshType, int spacedim> 192 #ifndef _MSC_VER 193 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator> 194 #else 195 std::vector< 196 typename dealii::internal:: 197 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type> 198 #endif find_cells_adjacent_to_vertex(const MeshType<dim,spacedim> & mesh,const unsigned int vertex)199 find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh, 200 const unsigned int vertex) 201 { 202 // make sure that the given vertex is 203 // an active vertex of the underlying 204 // triangulation 205 AssertIndexRange(vertex, mesh.get_triangulation().n_vertices()); 206 Assert(mesh.get_triangulation().get_used_vertices()[vertex], 207 ExcVertexNotUsed(vertex)); 208 209 // use a set instead of a vector 210 // to ensure that cells are inserted only 211 // once 212 std::set<typename dealii::internal:: 213 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type> 214 adjacent_cells; 215 216 typename dealii::internal:: 217 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type 218 cell = mesh.begin_active(), 219 endc = mesh.end(); 220 221 // go through all active cells and look if the vertex is part of that cell 222 // 223 // in 1d, this is all we need to care about. in 2d/3d we also need to worry 224 // that the vertex might be a hanging node on a face or edge of a cell; in 225 // this case, we would want to add those cells as well on whose faces the 226 // vertex is located but for which it is not a vertex itself. 227 // 228 // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging 229 // node can only be in the middle of a face and we can query the neighboring 230 // cell from the current cell. on the other hand, in 3d a hanging node 231 // vertex can also be on an edge but there can be many other cells on 232 // this edge and we can not access them from the cell we are currently 233 // on. 234 // 235 // so, in the 3d case, if we run the algorithm as in 2d, we catch all 236 // those cells for which the vertex we seek is on a *subface*, but we 237 // miss the case of cells for which the vertex we seek is on a 238 // sub-edge for which there is no corresponding sub-face (because the 239 // immediate neighbor behind this face is not refined), see for example 240 // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we 241 // haven't yet found the vertex for the current cell we also need to 242 // look at the mid-points of edges 243 // 244 // as a final note, deciding whether a neighbor is actually coarser is 245 // simple in the case of isotropic refinement (we just need to look at 246 // the level of the current and the neighboring cell). however, this 247 // isn't so simple if we have used anisotropic refinement since then 248 // the level of a cell is not indicative of whether it is coarser or 249 // not than the current cell. ultimately, we want to add all cells on 250 // which the vertex is, independent of whether they are coarser or 251 // finer and so in the 2d case below we simply add *any* *active* neighbor. 252 // in the worst case, we add cells multiple times to the adjacent_cells 253 // list, but std::set throws out those cells already entered 254 for (; cell != endc; ++cell) 255 { 256 for (const unsigned int v : cell->vertex_indices()) 257 if (cell->vertex_index(v) == vertex) 258 { 259 // OK, we found a cell that contains 260 // the given vertex. We add it 261 // to the list. 262 adjacent_cells.insert(cell); 263 264 // as explained above, in 2+d we need to check whether 265 // this vertex is on a face behind which there is a 266 // (possibly) coarser neighbor. if this is the case, 267 // then we need to also add this neighbor 268 if (dim >= 2) 269 for (unsigned int vface = 0; vface < dim; vface++) 270 { 271 const unsigned int face = 272 GeometryInfo<dim>::vertex_to_face[v][vface]; 273 274 if (!cell->at_boundary(face) && 275 cell->neighbor(face)->is_active()) 276 { 277 // there is a (possibly) coarser cell behind a 278 // face to which the vertex belongs. the 279 // vertex we are looking at is then either a 280 // vertex of that coarser neighbor, or it is a 281 // hanging node on one of the faces of that 282 // cell. in either case, it is adjacent to the 283 // vertex, so add it to the list as well (if 284 // the cell was already in the list then the 285 // std::set makes sure that we get it only 286 // once) 287 adjacent_cells.insert(cell->neighbor(face)); 288 } 289 } 290 291 // in any case, we have found a cell, so go to the next cell 292 goto next_cell; 293 } 294 295 // in 3d also loop over the edges 296 if (dim >= 3) 297 { 298 for (unsigned int e = 0; e < cell->n_lines(); ++e) 299 if (cell->line(e)->has_children()) 300 // the only place where this vertex could have been 301 // hiding is on the mid-edge point of the edge we 302 // are looking at 303 if (cell->line(e)->child(0)->vertex_index(1) == vertex) 304 { 305 adjacent_cells.insert(cell); 306 307 // jump out of this tangle of nested loops 308 goto next_cell; 309 } 310 } 311 312 // in more than 3d we would probably have to do the same as 313 // above also for even lower-dimensional objects 314 Assert(dim <= 3, ExcNotImplemented()); 315 316 // move on to the next cell if we have found the 317 // vertex on the current one 318 next_cell:; 319 } 320 321 // if this was an active vertex then there needs to have been 322 // at least one cell to which it is adjacent! 323 Assert(adjacent_cells.size() > 0, ExcInternalError()); 324 325 // return the result as a vector, rather than the set we built above 326 return std::vector< 327 typename dealii::internal:: 328 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>( 329 adjacent_cells.begin(), adjacent_cells.end()); 330 } 331 332 333 334 namespace 335 { 336 template <int dim, template <int, int> class MeshType, int spacedim> 337 void find_active_cell_around_point_internal(const MeshType<dim,spacedim> & mesh,std::set<typename MeshType<dim,spacedim>::active_cell_iterator> & searched_cells,std::set<typename MeshType<dim,spacedim>::active_cell_iterator> & adjacent_cells)338 find_active_cell_around_point_internal( 339 const MeshType<dim, spacedim> &mesh, 340 #ifndef _MSC_VER 341 std::set<typename MeshType<dim, spacedim>::active_cell_iterator> 342 &searched_cells, 343 std::set<typename MeshType<dim, spacedim>::active_cell_iterator> 344 &adjacent_cells) 345 #else 346 std::set< 347 typename dealii::internal:: 348 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type> 349 &searched_cells, 350 std::set< 351 typename dealii::internal:: 352 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type> 353 &adjacent_cells) 354 #endif 355 { 356 #ifndef _MSC_VER 357 using cell_iterator = 358 typename MeshType<dim, spacedim>::active_cell_iterator; 359 #else 360 using cell_iterator = typename dealii::internal:: 361 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type; 362 #endif 363 364 // update the searched cells 365 searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end()); 366 // now we to collect all neighbors 367 // of the cells in adjacent_cells we 368 // have not yet searched. 369 std::set<cell_iterator> adjacent_cells_new; 370 371 typename std::set<cell_iterator>::const_iterator cell = 372 adjacent_cells.begin(), 373 endc = 374 adjacent_cells.end(); 375 for (; cell != endc; ++cell) 376 { 377 std::vector<cell_iterator> active_neighbors; 378 get_active_neighbors<MeshType<dim, spacedim>>(*cell, 379 active_neighbors); 380 for (unsigned int i = 0; i < active_neighbors.size(); ++i) 381 if (searched_cells.find(active_neighbors[i]) == 382 searched_cells.end()) 383 adjacent_cells_new.insert(active_neighbors[i]); 384 } 385 adjacent_cells.clear(); 386 adjacent_cells.insert(adjacent_cells_new.begin(), 387 adjacent_cells_new.end()); 388 if (adjacent_cells.size() == 0) 389 { 390 // we haven't found any other cell that would be a 391 // neighbor of a previously found cell, but we know 392 // that we haven't checked all cells yet. that means 393 // that the domain is disconnected. in that case, 394 // choose the first previously untouched cell we 395 // can find 396 cell_iterator it = mesh.begin_active(); 397 for (; it != mesh.end(); ++it) 398 if (searched_cells.find(it) == searched_cells.end()) 399 { 400 adjacent_cells.insert(it); 401 break; 402 } 403 } 404 } 405 } // namespace 406 407 408 409 template <int dim, template <int, int> class MeshType, int spacedim> 410 #ifndef _MSC_VER 411 typename MeshType<dim, spacedim>::active_cell_iterator 412 #else 413 typename dealii::internal:: 414 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type 415 #endif find_active_cell_around_point(const MeshType<dim,spacedim> & mesh,const Point<spacedim> & p,const std::vector<bool> & marked_vertices,const double tolerance)416 find_active_cell_around_point(const MeshType<dim, spacedim> &mesh, 417 const Point<spacedim> & p, 418 const std::vector<bool> & marked_vertices, 419 const double tolerance) 420 { 421 return find_active_cell_around_point<dim, MeshType, spacedim>( 422 StaticMappingQ1<dim, spacedim>::mapping, 423 mesh, 424 p, 425 marked_vertices, 426 tolerance) 427 .first; 428 } 429 430 431 432 template <int dim, template <int, int> class MeshType, int spacedim> 433 #ifndef _MSC_VER 434 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>> 435 #else 436 std::pair<typename dealii::internal:: 437 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 438 Point<dim>> 439 #endif find_active_cell_around_point(const Mapping<dim,spacedim> & mapping,const MeshType<dim,spacedim> & mesh,const Point<spacedim> & p,const std::vector<bool> & marked_vertices,const double tolerance)440 find_active_cell_around_point(const Mapping<dim, spacedim> & mapping, 441 const MeshType<dim, spacedim> &mesh, 442 const Point<spacedim> & p, 443 const std::vector<bool> & marked_vertices, 444 const double tolerance) 445 { 446 using active_cell_iterator = typename dealii::internal:: 447 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type; 448 449 // The best distance is set to the 450 // maximum allowable distance from 451 // the unit cell; we assume a 452 // max. deviation of the given tolerance 453 double best_distance = tolerance; 454 int best_level = -1; 455 std::pair<active_cell_iterator, Point<dim>> best_cell; 456 457 // Find closest vertex and determine 458 // all adjacent cells 459 std::vector<active_cell_iterator> adjacent_cells_tmp = 460 find_cells_adjacent_to_vertex( 461 mesh, find_closest_vertex(mapping, mesh, p, marked_vertices)); 462 463 // Make sure that we have found 464 // at least one cell adjacent to vertex. 465 Assert(adjacent_cells_tmp.size() > 0, ExcInternalError()); 466 467 // Copy all the cells into a std::set 468 std::set<active_cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), 469 adjacent_cells_tmp.end()); 470 std::set<active_cell_iterator> searched_cells; 471 472 // Determine the maximal number of cells 473 // in the grid. 474 // As long as we have not found 475 // the cell and have not searched 476 // every cell in the triangulation, 477 // we keep on looking. 478 const unsigned int n_active_cells = 479 mesh.get_triangulation().n_active_cells(); 480 bool found = false; 481 unsigned int cells_searched = 0; 482 while (!found && cells_searched < n_active_cells) 483 { 484 typename std::set<active_cell_iterator>::const_iterator 485 cell = adjacent_cells.begin(), 486 endc = adjacent_cells.end(); 487 for (; cell != endc; ++cell) 488 { 489 if ((*cell)->is_artificial() == false) 490 { 491 try 492 { 493 const Point<dim> p_cell = 494 mapping.transform_real_to_unit_cell(*cell, p); 495 496 // calculate the infinity norm of 497 // the distance vector to the unit cell. 498 const double dist = 499 GeometryInfo<dim>::distance_to_unit_cell(p_cell); 500 501 // We compare if the point is inside the 502 // unit cell (or at least not too far 503 // outside). If it is, it is also checked 504 // that the cell has a more refined state 505 if ((dist < best_distance) || 506 ((dist == best_distance) && 507 ((*cell)->level() > best_level))) 508 { 509 found = true; 510 best_distance = dist; 511 best_level = (*cell)->level(); 512 best_cell = std::make_pair(*cell, p_cell); 513 } 514 } 515 catch ( 516 typename MappingQGeneric<dim, 517 spacedim>::ExcTransformationFailed &) 518 { 519 // ok, the transformation 520 // failed presumably 521 // because the point we 522 // are looking for lies 523 // outside the current 524 // cell. this means that 525 // the current cell can't 526 // be the cell around the 527 // point, so just ignore 528 // this cell and move on 529 // to the next 530 } 531 } 532 } 533 534 // update the number of cells searched 535 cells_searched += adjacent_cells.size(); 536 537 // if the user provided a custom mask for vertices, 538 // terminate the search without trying to expand the search 539 // to all cells of the triangulation, as done below. 540 if (marked_vertices.size() > 0) 541 cells_searched = n_active_cells; 542 543 // if we have not found the cell in 544 // question and have not yet searched every 545 // cell, we expand our search to 546 // all the not already searched neighbors of 547 // the cells in adjacent_cells. This is 548 // what find_active_cell_around_point_internal 549 // is for. 550 if (!found && cells_searched < n_active_cells) 551 { 552 find_active_cell_around_point_internal<dim, MeshType, spacedim>( 553 mesh, searched_cells, adjacent_cells); 554 } 555 } 556 557 AssertThrow(best_cell.first.state() == IteratorState::valid, 558 ExcPointNotFound<spacedim>(p)); 559 560 return best_cell; 561 } 562 563 564 565 template <int dim, template <int, int> class MeshType, int spacedim> 566 #ifndef _MSC_VER 567 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 568 Point<dim>>> 569 #else 570 std::vector<std::pair< 571 typename dealii::internal:: 572 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 573 Point<dim>>> 574 #endif find_all_active_cells_around_point(const Mapping<dim,spacedim> & mapping,const MeshType<dim,spacedim> & mesh,const Point<spacedim> & p,const double tolerance,const std::vector<bool> & marked_vertices)575 find_all_active_cells_around_point(const Mapping<dim, spacedim> & mapping, 576 const MeshType<dim, spacedim> &mesh, 577 const Point<spacedim> & p, 578 const double tolerance, 579 const std::vector<bool> &marked_vertices) 580 { 581 try 582 { 583 const auto cell_and_point = find_active_cell_around_point( 584 mapping, mesh, p, marked_vertices, tolerance); 585 586 return find_all_active_cells_around_point( 587 mapping, mesh, p, tolerance, cell_and_point); 588 } 589 catch (ExcPointNotFound<spacedim> &) 590 {} 591 592 return {}; 593 } 594 595 596 597 template <int dim, template <int, int> class MeshType, int spacedim> 598 #ifndef _MSC_VER 599 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 600 Point<dim>>> 601 #else 602 std::vector<std::pair< 603 typename dealii::internal:: 604 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 605 Point<dim>>> 606 #endif find_all_active_cells_around_point(const Mapping<dim,spacedim> & mapping,const MeshType<dim,spacedim> & mesh,const Point<spacedim> & p,const double tolerance,const std::pair<typename MeshType<dim,spacedim>::active_cell_iterator,Point<dim>> & first_cell)607 find_all_active_cells_around_point( 608 const Mapping<dim, spacedim> & mapping, 609 const MeshType<dim, spacedim> &mesh, 610 const Point<spacedim> & p, 611 const double tolerance, 612 const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 613 Point<dim>> & first_cell) 614 { 615 std::vector< 616 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 617 Point<dim>>> 618 cells_and_points; 619 620 // insert the fist cell and point into the vector 621 cells_and_points.push_back(first_cell); 622 623 // check if the given point is on the surface of the unit cell. if yes, 624 // need to find all neighbors 625 const Point<dim> unit_point = cells_and_points.front().second; 626 const auto my_cell = cells_and_points.front().first; 627 Tensor<1, dim> distance_to_center; 628 unsigned int n_dirs_at_threshold = 0; 629 unsigned int last_point_at_threshold = numbers::invalid_unsigned_int; 630 for (unsigned int d = 0; d < dim; ++d) 631 { 632 distance_to_center[d] = std::abs(unit_point[d] - 0.5); 633 if (distance_to_center[d] > 0.5 - tolerance) 634 { 635 ++n_dirs_at_threshold; 636 last_point_at_threshold = d; 637 } 638 } 639 640 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator> 641 cells_to_add; 642 // point is within face -> only need neighbor 643 if (n_dirs_at_threshold == 1) 644 { 645 unsigned int neighbor_index = 646 2 * last_point_at_threshold + 647 (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0); 648 if (!my_cell->at_boundary(neighbor_index)) 649 cells_to_add.push_back(my_cell->neighbor(neighbor_index)); 650 } 651 // corner point -> use all neighbors 652 else if (n_dirs_at_threshold == dim) 653 { 654 unsigned int local_vertex_index = 0; 655 for (unsigned int d = 0; d < dim; ++d) 656 local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d; 657 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator> 658 cells = find_cells_adjacent_to_vertex( 659 mesh, my_cell->vertex_index(local_vertex_index)); 660 for (const auto &cell : cells) 661 if (cell != my_cell) 662 cells_to_add.push_back(cell); 663 } 664 // point on line in 3D: We cannot simply take the intersection between 665 // the two vertices of cells because of hanging nodes. So instead we 666 // list the vertices around both points and then select the 667 // appropriate cells according to the result of read_to_unit_cell 668 // below. 669 else if (n_dirs_at_threshold == 2) 670 { 671 std::pair<unsigned int, unsigned int> vertex_indices[3]; 672 unsigned int count_vertex_indices = 0; 673 unsigned int free_direction = numbers::invalid_unsigned_int; 674 for (unsigned int d = 0; d < dim; ++d) 675 { 676 if (distance_to_center[d] > 0.5 - tolerance) 677 { 678 vertex_indices[count_vertex_indices].first = d; 679 vertex_indices[count_vertex_indices].second = 680 unit_point[d] > 0.5 ? 1 : 0; 681 count_vertex_indices++; 682 } 683 else 684 free_direction = d; 685 } 686 687 AssertDimension(count_vertex_indices, 2); 688 Assert(free_direction != numbers::invalid_unsigned_int, 689 ExcInternalError()); 690 691 const unsigned int first_vertex = 692 (vertex_indices[0].second << vertex_indices[0].first) + 693 (vertex_indices[1].second << vertex_indices[1].first); 694 for (unsigned int d = 0; d < 2; ++d) 695 { 696 auto tentative_cells = find_cells_adjacent_to_vertex( 697 mesh, 698 my_cell->vertex_index(first_vertex + (d << free_direction))); 699 for (const auto &cell : tentative_cells) 700 { 701 bool cell_not_yet_present = true; 702 for (const auto &other_cell : cells_to_add) 703 if (cell == other_cell) 704 { 705 cell_not_yet_present = false; 706 break; 707 } 708 if (cell_not_yet_present) 709 cells_to_add.push_back(cell); 710 } 711 } 712 } 713 714 const double original_distance_to_unit_cell = 715 GeometryInfo<dim>::distance_to_unit_cell(unit_point); 716 for (const auto &cell : cells_to_add) 717 { 718 if (cell != my_cell) 719 try 720 { 721 const Point<dim> p_unit = 722 mapping.transform_real_to_unit_cell(cell, p); 723 if (GeometryInfo<dim>::distance_to_unit_cell(p_unit) < 724 original_distance_to_unit_cell + tolerance) 725 cells_and_points.emplace_back(cell, p_unit); 726 } 727 catch (typename Mapping<dim>::ExcTransformationFailed &) 728 {} 729 } 730 731 std::sort( 732 cells_and_points.begin(), 733 cells_and_points.end(), 734 [](const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 735 Point<dim>> &a, 736 const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 737 Point<dim>> &b) { return a.first < b.first; }); 738 739 return cells_and_points; 740 } 741 742 743 744 template <class MeshType> 745 std::vector<typename MeshType::active_cell_iterator> compute_active_cell_halo_layer(const MeshType & mesh,const std::function<bool (const typename MeshType::active_cell_iterator &)> & predicate)746 compute_active_cell_halo_layer( 747 const MeshType &mesh, 748 const std::function<bool(const typename MeshType::active_cell_iterator &)> 749 &predicate) 750 { 751 std::vector<typename MeshType::active_cell_iterator> active_halo_layer; 752 std::vector<bool> locally_active_vertices_on_subdomain( 753 mesh.get_triangulation().n_vertices(), false); 754 755 // Find the cells for which the predicate is true 756 // These are the cells around which we wish to construct 757 // the halo layer 758 for (const auto &cell : mesh.active_cell_iterators()) 759 if (predicate(cell)) // True predicate --> Part of subdomain 760 for (const auto v : cell->vertex_indices()) 761 locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true; 762 763 // Find the cells that do not conform to the predicate 764 // but share a vertex with the selected subdomain 765 // These comprise the halo layer 766 for (const auto &cell : mesh.active_cell_iterators()) 767 if (!predicate(cell)) // False predicate --> Potential halo cell 768 for (const auto v : cell->vertex_indices()) 769 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == 770 true) 771 { 772 active_halo_layer.push_back(cell); 773 break; 774 } 775 776 return active_halo_layer; 777 } 778 779 780 781 template <class MeshType> 782 std::vector<typename MeshType::cell_iterator> compute_cell_halo_layer_on_level(const MeshType & mesh,const std::function<bool (const typename MeshType::cell_iterator &)> & predicate,const unsigned int level)783 compute_cell_halo_layer_on_level( 784 const MeshType &mesh, 785 const std::function<bool(const typename MeshType::cell_iterator &)> 786 & predicate, 787 const unsigned int level) 788 { 789 std::vector<typename MeshType::cell_iterator> level_halo_layer; 790 std::vector<bool> locally_active_vertices_on_level_subdomain( 791 mesh.get_triangulation().n_vertices(), false); 792 793 // Find the cells for which the predicate is true 794 // These are the cells around which we wish to construct 795 // the halo layer 796 for (typename MeshType::cell_iterator cell = mesh.begin(level); 797 cell != mesh.end(level); 798 ++cell) 799 if (predicate(cell)) // True predicate --> Part of subdomain 800 for (const unsigned int v : cell->vertex_indices()) 801 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] = 802 true; 803 804 // Find the cells that do not conform to the predicate 805 // but share a vertex with the selected subdomain on that level 806 // These comprise the halo layer 807 for (typename MeshType::cell_iterator cell = mesh.begin(level); 808 cell != mesh.end(level); 809 ++cell) 810 if (!predicate(cell)) // False predicate --> Potential halo cell 811 for (const unsigned int v : cell->vertex_indices()) 812 if (locally_active_vertices_on_level_subdomain[cell->vertex_index( 813 v)] == true) 814 { 815 level_halo_layer.push_back(cell); 816 break; 817 } 818 819 return level_halo_layer; 820 } 821 822 823 namespace 824 { 825 template <class MeshType> 826 bool contains_locally_owned_cells(const std::vector<typename MeshType::active_cell_iterator> & cells)827 contains_locally_owned_cells( 828 const std::vector<typename MeshType::active_cell_iterator> &cells) 829 { 830 for (typename std::vector< 831 typename MeshType::active_cell_iterator>::const_iterator it = 832 cells.begin(); 833 it != cells.end(); 834 ++it) 835 { 836 if ((*it)->is_locally_owned()) 837 return true; 838 } 839 return false; 840 } 841 842 template <class MeshType> 843 bool contains_artificial_cells(const std::vector<typename MeshType::active_cell_iterator> & cells)844 contains_artificial_cells( 845 const std::vector<typename MeshType::active_cell_iterator> &cells) 846 { 847 for (typename std::vector< 848 typename MeshType::active_cell_iterator>::const_iterator it = 849 cells.begin(); 850 it != cells.end(); 851 ++it) 852 { 853 if ((*it)->is_artificial()) 854 return true; 855 } 856 return false; 857 } 858 } // namespace 859 860 861 862 template <class MeshType> 863 std::vector<typename MeshType::active_cell_iterator> compute_ghost_cell_halo_layer(const MeshType & mesh)864 compute_ghost_cell_halo_layer(const MeshType &mesh) 865 { 866 std::function<bool(const typename MeshType::active_cell_iterator &)> 867 predicate = IteratorFilters::LocallyOwnedCell(); 868 869 const std::vector<typename MeshType::active_cell_iterator> 870 active_halo_layer = compute_active_cell_halo_layer(mesh, predicate); 871 872 // Check that we never return locally owned or artificial cells 873 // What is left should only be the ghost cells 874 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) == false, 875 ExcMessage("Halo layer contains locally owned cells")); 876 Assert(contains_artificial_cells<MeshType>(active_halo_layer) == false, 877 ExcMessage("Halo layer contains artificial cells")); 878 879 return active_halo_layer; 880 } 881 882 883 884 template <class MeshType> 885 std::vector<typename MeshType::active_cell_iterator> compute_active_cell_layer_within_distance(const MeshType & mesh,const std::function<bool (const typename MeshType::active_cell_iterator &)> & predicate,const double layer_thickness)886 compute_active_cell_layer_within_distance( 887 const MeshType &mesh, 888 const std::function<bool(const typename MeshType::active_cell_iterator &)> 889 & predicate, 890 const double layer_thickness) 891 { 892 std::vector<typename MeshType::active_cell_iterator> 893 subdomain_boundary_cells, active_cell_layer_within_distance; 894 std::vector<bool> vertices_outside_subdomain( 895 mesh.get_triangulation().n_vertices(), false); 896 897 const unsigned int spacedim = MeshType::space_dimension; 898 899 unsigned int n_non_predicate_cells = 0; // Number of non predicate cells 900 901 // Find the layer of cells for which predicate is true and that 902 // are on the boundary with other cells. These are 903 // subdomain boundary cells. 904 905 // Find the cells for which the predicate is false 906 // These are the cells which are around the predicate subdomain 907 for (const auto &cell : mesh.active_cell_iterators()) 908 if (!predicate(cell)) // Negation of predicate --> Not Part of subdomain 909 { 910 for (const unsigned int v : cell->vertex_indices()) 911 vertices_outside_subdomain[cell->vertex_index(v)] = true; 912 n_non_predicate_cells++; 913 } 914 915 // If all the active cells conform to the predicate 916 // or if none of the active cells conform to the predicate 917 // there is no active cell layer around the predicate 918 // subdomain (within any distance) 919 if (n_non_predicate_cells == 0 || 920 n_non_predicate_cells == mesh.get_triangulation().n_active_cells()) 921 return std::vector<typename MeshType::active_cell_iterator>(); 922 923 // Find the cells that conform to the predicate 924 // but share a vertex with the cell not in the predicate subdomain 925 for (const auto &cell : mesh.active_cell_iterators()) 926 if (predicate(cell)) // True predicate --> Potential boundary cell of the 927 // subdomain 928 for (const unsigned int v : cell->vertex_indices()) 929 if (vertices_outside_subdomain[cell->vertex_index(v)] == true) 930 { 931 subdomain_boundary_cells.push_back(cell); 932 break; // No need to go through remaining vertices 933 } 934 935 // To cheaply filter out some cells located far away from the predicate 936 // subdomain, get the bounding box of the predicate subdomain. 937 std::pair<Point<spacedim>, Point<spacedim>> bounding_box = 938 compute_bounding_box(mesh, predicate); 939 940 // DOUBLE_EPSILON to compare really close double values 941 const double DOUBLE_EPSILON = 100. * std::numeric_limits<double>::epsilon(); 942 943 // Add layer_thickness to the bounding box 944 for (unsigned int d = 0; d < spacedim; ++d) 945 { 946 bounding_box.first[d] -= (layer_thickness + DOUBLE_EPSILON); // minp 947 bounding_box.second[d] += (layer_thickness + DOUBLE_EPSILON); // maxp 948 } 949 950 std::vector<Point<spacedim>> 951 subdomain_boundary_cells_centers; // cache all the subdomain boundary 952 // cells centers here 953 std::vector<double> 954 subdomain_boundary_cells_radii; // cache all the subdomain boundary cells 955 // radii 956 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size()); 957 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size()); 958 // compute cell radius for each boundary cell of the predicate subdomain 959 for (typename std::vector<typename MeshType::active_cell_iterator>:: 960 const_iterator subdomain_boundary_cell_iterator = 961 subdomain_boundary_cells.begin(); 962 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end(); 963 ++subdomain_boundary_cell_iterator) 964 { 965 const std::pair<Point<spacedim>, double> 966 &subdomain_boundary_cell_enclosing_ball = 967 (*subdomain_boundary_cell_iterator)->enclosing_ball(); 968 969 subdomain_boundary_cells_centers.push_back( 970 subdomain_boundary_cell_enclosing_ball.first); 971 subdomain_boundary_cells_radii.push_back( 972 subdomain_boundary_cell_enclosing_ball.second); 973 } 974 AssertThrow(subdomain_boundary_cells_radii.size() == 975 subdomain_boundary_cells_centers.size(), 976 ExcInternalError()); 977 978 // Find the cells that are within layer_thickness of predicate subdomain 979 // boundary distance but are inside the extended bounding box. Most cells 980 // might be outside the extended bounding box, so we could skip them. Those 981 // cells that are inside the extended bounding box but are not part of the 982 // predicate subdomain are possible candidates to be within the distance to 983 // the boundary cells of the predicate subdomain. 984 for (const auto &cell : mesh.active_cell_iterators()) 985 { 986 // Ignore all the cells that are in the predicate subdomain 987 if (predicate(cell)) 988 continue; 989 990 const std::pair<Point<spacedim>, double> &cell_enclosing_ball = 991 cell->enclosing_ball(); 992 993 const Point<spacedim> cell_enclosing_ball_center = 994 cell_enclosing_ball.first; 995 const double cell_enclosing_ball_radius = cell_enclosing_ball.second; 996 997 bool cell_inside = true; // reset for each cell 998 999 for (unsigned int d = 0; d < spacedim; ++d) 1000 cell_inside &= 1001 (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius > 1002 bounding_box.first[d]) && 1003 (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius < 1004 bounding_box.second[d]); 1005 // cell_inside is true if its enclosing ball intersects the extended 1006 // bounding box 1007 1008 // Ignore all the cells that are outside the extended bounding box 1009 if (cell_inside) 1010 for (unsigned int i = 0; i < subdomain_boundary_cells_radii.size(); 1011 ++i) 1012 if (cell_enclosing_ball_center.distance_square( 1013 subdomain_boundary_cells_centers[i]) < 1014 Utilities::fixed_power<2>(cell_enclosing_ball_radius + 1015 subdomain_boundary_cells_radii[i] + 1016 layer_thickness + DOUBLE_EPSILON)) 1017 { 1018 active_cell_layer_within_distance.push_back(cell); 1019 break; // Exit the loop checking all the remaining subdomain 1020 // boundary cells 1021 } 1022 } 1023 return active_cell_layer_within_distance; 1024 } 1025 1026 1027 1028 template <class MeshType> 1029 std::vector<typename MeshType::active_cell_iterator> compute_ghost_cell_layer_within_distance(const MeshType & mesh,const double layer_thickness)1030 compute_ghost_cell_layer_within_distance(const MeshType &mesh, 1031 const double layer_thickness) 1032 { 1033 IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate; 1034 std::function<bool(const typename MeshType::active_cell_iterator &)> 1035 predicate(locally_owned_cell_predicate); 1036 1037 const std::vector<typename MeshType::active_cell_iterator> 1038 ghost_cell_layer_within_distance = 1039 compute_active_cell_layer_within_distance(mesh, 1040 predicate, 1041 layer_thickness); 1042 1043 // Check that we never return locally owned or artificial cells 1044 // What is left should only be the ghost cells 1045 Assert( 1046 contains_locally_owned_cells<MeshType>( 1047 ghost_cell_layer_within_distance) == false, 1048 ExcMessage( 1049 "Ghost cells within layer_thickness contains locally owned cells.")); 1050 Assert( 1051 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) == 1052 false, 1053 ExcMessage( 1054 "Ghost cells within layer_thickness contains artificial cells. " 1055 "The function compute_ghost_cell_layer_within_distance " 1056 "is probably called while using parallel::distributed::Triangulation. " 1057 "In such case please refer to the description of this function.")); 1058 1059 return ghost_cell_layer_within_distance; 1060 } 1061 1062 1063 1064 template <class MeshType> 1065 std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>> compute_bounding_box(const MeshType & mesh,const std::function<bool (const typename MeshType::active_cell_iterator &)> & predicate)1066 compute_bounding_box( 1067 const MeshType &mesh, 1068 const std::function<bool(const typename MeshType::active_cell_iterator &)> 1069 &predicate) 1070 { 1071 std::vector<bool> locally_active_vertices_on_subdomain( 1072 mesh.get_triangulation().n_vertices(), false); 1073 1074 const unsigned int spacedim = MeshType::space_dimension; 1075 1076 // Two extreme points can define the bounding box 1077 // around the active cells that conform to the given predicate. 1078 Point<MeshType::space_dimension> maxp, minp; 1079 1080 // initialize minp and maxp with the first predicate cell center 1081 for (const auto &cell : mesh.active_cell_iterators()) 1082 if (predicate(cell)) 1083 { 1084 minp = cell->center(); 1085 maxp = cell->center(); 1086 break; 1087 } 1088 1089 // Run through all the cells to check if it belongs to predicate domain, 1090 // if it belongs to the predicate domain, extend the bounding box. 1091 for (const auto &cell : mesh.active_cell_iterators()) 1092 if (predicate(cell)) // True predicate --> Part of subdomain 1093 for (const unsigned int v : cell->vertex_indices()) 1094 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == 1095 false) 1096 { 1097 locally_active_vertices_on_subdomain[cell->vertex_index(v)] = 1098 true; 1099 for (unsigned int d = 0; d < spacedim; ++d) 1100 { 1101 minp[d] = std::min(minp[d], cell->vertex(v)[d]); 1102 maxp[d] = std::max(maxp[d], cell->vertex(v)[d]); 1103 } 1104 } 1105 1106 return std::make_pair(minp, maxp); 1107 } 1108 1109 1110 1111 template <typename MeshType> 1112 std::list<std::pair<typename MeshType::cell_iterator, 1113 typename MeshType::cell_iterator>> get_finest_common_cells(const MeshType & mesh_1,const MeshType & mesh_2)1114 get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2) 1115 { 1116 Assert(have_same_coarse_mesh(mesh_1, mesh_2), 1117 ExcMessage("The two meshes must be represent triangulations that " 1118 "have the same coarse meshes")); 1119 // We will allow the output to contain ghost cells when we have shared 1120 // Triangulations (i.e., so that each processor will get exactly the same 1121 // list of cell pairs), but not when we have two distributed 1122 // Triangulations (so that all active cells are partitioned by processor). 1123 // Non-parallel Triangulations have no ghost or artificial cells, so they 1124 // work the same way as shared Triangulations here. 1125 bool remove_ghost_cells = false; 1126 #ifdef DEAL_II_WITH_MPI 1127 { 1128 constexpr int dim = MeshType::dimension; 1129 constexpr int spacedim = MeshType::space_dimension; 1130 if (dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> 1131 *>(&mesh_1.get_triangulation()) != nullptr || 1132 dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim> 1133 *>(&mesh_2.get_triangulation()) != nullptr) 1134 { 1135 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(), 1136 ExcMessage("This function can only be used with meshes " 1137 "corresponding to distributed Triangulations when " 1138 "both Triangulations are equal.")); 1139 remove_ghost_cells = true; 1140 } 1141 } 1142 #endif 1143 1144 // the algorithm goes as follows: first, we fill a list with pairs of 1145 // iterators common to the two meshes on the coarsest level. then we 1146 // traverse the list; each time, we find a pair of iterators for which 1147 // both correspond to non-active cells, we delete this item and push the 1148 // pairs of iterators to their children to the back. if these again both 1149 // correspond to non-active cells, we will get to the later on for further 1150 // consideration 1151 using CellList = std::list<std::pair<typename MeshType::cell_iterator, 1152 typename MeshType::cell_iterator>>; 1153 CellList cell_list; 1154 1155 // first push the coarse level cells 1156 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0), 1157 cell_2 = mesh_2.begin(0); 1158 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2) 1159 cell_list.emplace_back(cell_1, cell_2); 1160 1161 // then traverse list as described above 1162 typename CellList::iterator cell_pair = cell_list.begin(); 1163 while (cell_pair != cell_list.end()) 1164 { 1165 // if both cells in this pair have children, then erase this element 1166 // and push their children instead 1167 if (cell_pair->first->has_children() && 1168 cell_pair->second->has_children()) 1169 { 1170 Assert(cell_pair->first->refinement_case() == 1171 cell_pair->second->refinement_case(), 1172 ExcNotImplemented()); 1173 for (unsigned int c = 0; c < cell_pair->first->n_children(); ++c) 1174 cell_list.emplace_back(cell_pair->first->child(c), 1175 cell_pair->second->child(c)); 1176 1177 // erasing an iterator keeps other iterators valid, so already 1178 // advance the present iterator by one and then delete the element 1179 // we've visited before 1180 const auto previous_cell_pair = cell_pair; 1181 ++cell_pair; 1182 cell_list.erase(previous_cell_pair); 1183 } 1184 else 1185 { 1186 // at least one cell is active 1187 if (remove_ghost_cells && 1188 ((cell_pair->first->is_active() && 1189 !cell_pair->first->is_locally_owned()) || 1190 (cell_pair->second->is_active() && 1191 !cell_pair->second->is_locally_owned()))) 1192 { 1193 // we only exclude ghost cells for distributed Triangulations 1194 const auto previous_cell_pair = cell_pair; 1195 ++cell_pair; 1196 cell_list.erase(previous_cell_pair); 1197 } 1198 else 1199 ++cell_pair; 1200 } 1201 } 1202 1203 // just to make sure everything is ok, validate that all pairs have at 1204 // least one active iterator or have different refinement_cases 1205 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end(); 1206 ++cell_pair) 1207 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() || 1208 (cell_pair->first->refinement_case() != 1209 cell_pair->second->refinement_case()), 1210 ExcInternalError()); 1211 1212 return cell_list; 1213 } 1214 1215 1216 1217 template <int dim, int spacedim> 1218 bool have_same_coarse_mesh(const Triangulation<dim,spacedim> & mesh_1,const Triangulation<dim,spacedim> & mesh_2)1219 have_same_coarse_mesh(const Triangulation<dim, spacedim> &mesh_1, 1220 const Triangulation<dim, spacedim> &mesh_2) 1221 { 1222 // make sure the two meshes have 1223 // the same number of coarse cells 1224 if (mesh_1.n_cells(0) != mesh_2.n_cells(0)) 1225 return false; 1226 1227 // if so, also make sure they have 1228 // the same vertices on the cells 1229 // of the coarse mesh 1230 typename Triangulation<dim, spacedim>::cell_iterator cell_1 = 1231 mesh_1.begin(0), 1232 cell_2 = 1233 mesh_2.begin(0), 1234 endc = mesh_1.end(0); 1235 for (; cell_1 != endc; ++cell_1, ++cell_2) 1236 { 1237 if (cell_1->n_vertices() != cell_2->n_vertices()) 1238 return false; 1239 for (const unsigned int v : cell_1->vertex_indices()) 1240 if (cell_1->vertex(v) != cell_2->vertex(v)) 1241 return false; 1242 } 1243 1244 // if we've gotten through all 1245 // this, then the meshes really 1246 // seem to have a common coarse 1247 // mesh 1248 return true; 1249 } 1250 1251 1252 1253 template <typename MeshType> 1254 bool have_same_coarse_mesh(const MeshType & mesh_1,const MeshType & mesh_2)1255 have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2) 1256 { 1257 return have_same_coarse_mesh(mesh_1.get_triangulation(), 1258 mesh_2.get_triangulation()); 1259 } 1260 1261 1262 1263 template <int dim, int spacedim> 1264 std::pair<typename dealii::DoFHandler<dim, spacedim>::active_cell_iterator, 1265 Point<dim>> find_active_cell_around_point(const hp::MappingCollection<dim,spacedim> & mapping,const DoFHandler<dim,spacedim> & mesh,const Point<spacedim> & p,const double tolerance)1266 find_active_cell_around_point( 1267 const hp::MappingCollection<dim, spacedim> &mapping, 1268 const DoFHandler<dim, spacedim> & mesh, 1269 const Point<spacedim> & p, 1270 const double tolerance) 1271 { 1272 Assert((mapping.size() == 1) || 1273 (mapping.size() == mesh.get_fe_collection().size()), 1274 ExcMessage("Mapping collection needs to have either size 1 " 1275 "or size equal to the number of elements in " 1276 "the FECollection.")); 1277 1278 using cell_iterator = 1279 typename dealii::DoFHandler<dim, spacedim>::active_cell_iterator; 1280 1281 std::pair<cell_iterator, Point<dim>> best_cell; 1282 // If we have only one element in the MappingCollection, 1283 // we use find_active_cell_around_point using only one 1284 // mapping. 1285 if (mapping.size() == 1) 1286 { 1287 const std::vector<bool> marked_vertices = {}; 1288 best_cell = find_active_cell_around_point( 1289 mapping[0], mesh, p, marked_vertices, tolerance); 1290 } 1291 else 1292 { 1293 // The best distance is set to the 1294 // maximum allowable distance from 1295 // the unit cell 1296 double best_distance = tolerance; 1297 int best_level = -1; 1298 1299 1300 // Find closest vertex and determine 1301 // all adjacent cells 1302 unsigned int vertex = find_closest_vertex(mesh, p); 1303 1304 std::vector<cell_iterator> adjacent_cells_tmp = 1305 find_cells_adjacent_to_vertex(mesh, vertex); 1306 1307 // Make sure that we have found 1308 // at least one cell adjacent to vertex. 1309 Assert(adjacent_cells_tmp.size() > 0, ExcInternalError()); 1310 1311 // Copy all the cells into a std::set 1312 std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), 1313 adjacent_cells_tmp.end()); 1314 std::set<cell_iterator> searched_cells; 1315 1316 // Determine the maximal number of cells 1317 // in the grid. 1318 // As long as we have not found 1319 // the cell and have not searched 1320 // every cell in the triangulation, 1321 // we keep on looking. 1322 const unsigned int n_cells = mesh.get_triangulation().n_cells(); 1323 bool found = false; 1324 unsigned int cells_searched = 0; 1325 while (!found && cells_searched < n_cells) 1326 { 1327 typename std::set<cell_iterator>::const_iterator 1328 cell = adjacent_cells.begin(), 1329 endc = adjacent_cells.end(); 1330 for (; cell != endc; ++cell) 1331 { 1332 try 1333 { 1334 const Point<dim> p_cell = 1335 mapping[(*cell)->active_fe_index()] 1336 .transform_real_to_unit_cell(*cell, p); 1337 1338 1339 // calculate the infinity norm of 1340 // the distance vector to the unit cell. 1341 const double dist = 1342 GeometryInfo<dim>::distance_to_unit_cell(p_cell); 1343 1344 // We compare if the point is inside the 1345 // unit cell (or at least not too far 1346 // outside). If it is, it is also checked 1347 // that the cell has a more refined state 1348 if (dist < best_distance || (dist == best_distance && 1349 (*cell)->level() > best_level)) 1350 { 1351 found = true; 1352 best_distance = dist; 1353 best_level = (*cell)->level(); 1354 best_cell = std::make_pair(*cell, p_cell); 1355 } 1356 } 1357 catch ( 1358 typename MappingQGeneric<dim, 1359 spacedim>::ExcTransformationFailed &) 1360 { 1361 // ok, the transformation 1362 // failed presumably 1363 // because the point we 1364 // are looking for lies 1365 // outside the current 1366 // cell. this means that 1367 // the current cell can't 1368 // be the cell around the 1369 // point, so just ignore 1370 // this cell and move on 1371 // to the next 1372 } 1373 } 1374 // update the number of cells searched 1375 cells_searched += adjacent_cells.size(); 1376 // if we have not found the cell in 1377 // question and have not yet searched every 1378 // cell, we expand our search to 1379 // all the not already searched neighbors of 1380 // the cells in adjacent_cells. 1381 if (!found && cells_searched < n_cells) 1382 { 1383 find_active_cell_around_point_internal<dim, 1384 dealii::DoFHandler, 1385 spacedim>( 1386 mesh, searched_cells, adjacent_cells); 1387 } 1388 } 1389 } 1390 1391 AssertThrow(best_cell.first.state() == IteratorState::valid, 1392 ExcPointNotFound<spacedim>(p)); 1393 1394 return best_cell; 1395 } 1396 1397 1398 template <class MeshType> 1399 std::vector<typename MeshType::active_cell_iterator> get_patch_around_cell(const typename MeshType::active_cell_iterator & cell)1400 get_patch_around_cell(const typename MeshType::active_cell_iterator &cell) 1401 { 1402 Assert(cell->is_locally_owned(), 1403 ExcMessage("This function only makes sense if the cell for " 1404 "which you are asking for a patch, is locally " 1405 "owned.")); 1406 1407 std::vector<typename MeshType::active_cell_iterator> patch; 1408 patch.push_back(cell); 1409 for (const unsigned int face_number : cell->face_indices()) 1410 if (cell->face(face_number)->at_boundary() == false) 1411 { 1412 if (cell->neighbor(face_number)->has_children() == false) 1413 patch.push_back(cell->neighbor(face_number)); 1414 else 1415 // the neighbor is refined. in 2d/3d, we can simply ask for the 1416 // children of the neighbor because they can not be further refined 1417 // and, consequently, the children is active 1418 if (MeshType::dimension > 1) 1419 { 1420 for (unsigned int subface = 0; 1421 subface < cell->face(face_number)->n_children(); 1422 ++subface) 1423 patch.push_back( 1424 cell->neighbor_child_on_subface(face_number, subface)); 1425 } 1426 else 1427 { 1428 // in 1d, we need to work a bit harder: iterate until we find 1429 // the child by going from cell to child to child etc 1430 typename MeshType::cell_iterator neighbor = 1431 cell->neighbor(face_number); 1432 while (neighbor->has_children()) 1433 neighbor = neighbor->child(1 - face_number); 1434 1435 Assert(neighbor->neighbor(1 - face_number) == cell, 1436 ExcInternalError()); 1437 patch.push_back(neighbor); 1438 } 1439 } 1440 return patch; 1441 } 1442 1443 1444 1445 template <class Container> 1446 std::vector<typename Container::cell_iterator> get_cells_at_coarsest_common_level(const std::vector<typename Container::active_cell_iterator> & patch)1447 get_cells_at_coarsest_common_level( 1448 const std::vector<typename Container::active_cell_iterator> &patch) 1449 { 1450 Assert(patch.size() > 0, 1451 ExcMessage( 1452 "Vector containing patch cells should not be an empty vector!")); 1453 // In order to extract the set of cells with the coarsest common level from 1454 // the give vector of cells: First it finds the number associated with the 1455 // minimum level of refinement, namely "min_level" 1456 int min_level = patch[0]->level(); 1457 1458 for (unsigned int i = 0; i < patch.size(); ++i) 1459 min_level = std::min(min_level, patch[i]->level()); 1460 std::set<typename Container::cell_iterator> uniform_cells; 1461 typename std::vector< 1462 typename Container::active_cell_iterator>::const_iterator patch_cell; 1463 // it loops through all cells of the input vector 1464 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell) 1465 { 1466 // If the refinement level of each cell i the loop be equal to the 1467 // min_level, so that that cell inserted into the set of uniform_cells, 1468 // as the set of cells with the coarsest common refinement level 1469 if ((*patch_cell)->level() == min_level) 1470 uniform_cells.insert(*patch_cell); 1471 else 1472 // If not, it asks for the parent of the cell, until it finds the 1473 // parent cell with the refinement level equal to the min_level and 1474 // inserts that parent cell into the the set of uniform_cells, as the 1475 // set of cells with the coarsest common refinement level. 1476 { 1477 typename Container::cell_iterator parent = *patch_cell; 1478 1479 while (parent->level() > min_level) 1480 parent = parent->parent(); 1481 uniform_cells.insert(parent); 1482 } 1483 } 1484 1485 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(), 1486 uniform_cells.end()); 1487 } 1488 1489 1490 1491 template <class Container> 1492 void build_triangulation_from_patch(const std::vector<typename Container::active_cell_iterator> & patch,Triangulation<Container::dimension,Container::space_dimension> & local_triangulation,std::map<typename Triangulation<Container::dimension,Container::space_dimension>::active_cell_iterator,typename Container::active_cell_iterator> & patch_to_global_tria_map)1493 build_triangulation_from_patch( 1494 const std::vector<typename Container::active_cell_iterator> &patch, 1495 Triangulation<Container::dimension, Container::space_dimension> 1496 &local_triangulation, 1497 std::map< 1498 typename Triangulation<Container::dimension, 1499 Container::space_dimension>::active_cell_iterator, 1500 typename Container::active_cell_iterator> &patch_to_global_tria_map) 1501 1502 { 1503 const std::vector<typename Container::cell_iterator> uniform_cells = 1504 get_cells_at_coarsest_common_level<Container>(patch); 1505 // First it creates triangulation from the vector of "uniform_cells" 1506 local_triangulation.clear(); 1507 std::vector<Point<Container::space_dimension>> vertices; 1508 const unsigned int n_uniform_cells = uniform_cells.size(); 1509 std::vector<CellData<Container::dimension>> cells(n_uniform_cells); 1510 unsigned int k = 0; // for enumerating cells 1511 unsigned int i = 0; // for enumerating vertices 1512 typename std::vector<typename Container::cell_iterator>::const_iterator 1513 uniform_cell; 1514 for (uniform_cell = uniform_cells.begin(); 1515 uniform_cell != uniform_cells.end(); 1516 ++uniform_cell) 1517 { 1518 for (const unsigned int v : (*uniform_cell)->vertex_indices()) 1519 { 1520 Point<Container::space_dimension> position = 1521 (*uniform_cell)->vertex(v); 1522 bool repeat_vertex = false; 1523 1524 for (unsigned int m = 0; m < i; ++m) 1525 { 1526 if (position == vertices[m]) 1527 { 1528 repeat_vertex = true; 1529 cells[k].vertices[v] = m; 1530 break; 1531 } 1532 } 1533 if (repeat_vertex == false) 1534 { 1535 vertices.push_back(position); 1536 cells[k].vertices[v] = i; 1537 i = i + 1; 1538 } 1539 1540 } // for vertices_per_cell 1541 k = k + 1; 1542 } 1543 local_triangulation.create_triangulation(vertices, cells, SubCellData()); 1544 Assert(local_triangulation.n_active_cells() == uniform_cells.size(), 1545 ExcInternalError()); 1546 local_triangulation.clear_user_flags(); 1547 unsigned int index = 0; 1548 // Create a map between cells of class DoFHandler into class Triangulation 1549 std::map<typename Triangulation<Container::dimension, 1550 Container::space_dimension>::cell_iterator, 1551 typename Container::cell_iterator> 1552 patch_to_global_tria_map_tmp; 1553 for (typename Triangulation<Container::dimension, 1554 Container::space_dimension>::cell_iterator 1555 coarse_cell = local_triangulation.begin(); 1556 coarse_cell != local_triangulation.end(); 1557 ++coarse_cell, ++index) 1558 { 1559 patch_to_global_tria_map_tmp.insert( 1560 std::make_pair(coarse_cell, uniform_cells[index])); 1561 // To ensure that the cells with the same coordinates (here, we compare 1562 // their centers) are mapped into each other. 1563 1564 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <= 1565 1e-15 * coarse_cell->diameter(), 1566 ExcInternalError()); 1567 } 1568 bool refinement_necessary; 1569 // In this loop we start to do refinement on the above coarse triangulation 1570 // to reach to the same level of refinement as the patch cells are really on 1571 do 1572 { 1573 refinement_necessary = false; 1574 for (const auto &active_tria_cell : 1575 local_triangulation.active_cell_iterators()) 1576 { 1577 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children()) 1578 { 1579 active_tria_cell->set_refine_flag(); 1580 refinement_necessary = true; 1581 } 1582 else 1583 for (unsigned int i = 0; i < patch.size(); ++i) 1584 { 1585 // Even though vertices may not be exactly the same, the 1586 // appropriate cells will match since == for TriAccessors 1587 // checks only cell level and index. 1588 if (patch_to_global_tria_map_tmp[active_tria_cell] == 1589 patch[i]) 1590 { 1591 // adjust the cell vertices of the local_triangulation to 1592 // match cell vertices of the global triangulation 1593 for (const unsigned int v : 1594 active_tria_cell->vertex_indices()) 1595 active_tria_cell->vertex(v) = patch[i]->vertex(v); 1596 1597 Assert(active_tria_cell->center().distance( 1598 patch_to_global_tria_map_tmp[active_tria_cell] 1599 ->center()) <= 1600 1e-15 * active_tria_cell->diameter(), 1601 ExcInternalError()); 1602 1603 active_tria_cell->set_user_flag(); 1604 break; 1605 } 1606 } 1607 } 1608 1609 if (refinement_necessary) 1610 { 1611 local_triangulation.execute_coarsening_and_refinement(); 1612 1613 for (typename Triangulation< 1614 Container::dimension, 1615 Container::space_dimension>::cell_iterator cell = 1616 local_triangulation.begin(); 1617 cell != local_triangulation.end(); 1618 ++cell) 1619 { 1620 if (patch_to_global_tria_map_tmp.find(cell) != 1621 patch_to_global_tria_map_tmp.end()) 1622 { 1623 if (cell->has_children()) 1624 { 1625 // Note: Since the cell got children, then it should not 1626 // be in the map anymore children may be added into the 1627 // map, instead 1628 1629 // these children may not yet be in the map 1630 for (unsigned int c = 0; c < cell->n_children(); ++c) 1631 { 1632 if (patch_to_global_tria_map_tmp.find(cell->child( 1633 c)) == patch_to_global_tria_map_tmp.end()) 1634 { 1635 patch_to_global_tria_map_tmp.insert( 1636 std::make_pair( 1637 cell->child(c), 1638 patch_to_global_tria_map_tmp[cell]->child( 1639 c))); 1640 1641 // One might be tempted to assert that the cell 1642 // being added here has the same center as the 1643 // equivalent cell in the global triangulation, 1644 // but it may not be the case. For 1645 // triangulations that have been perturbed or 1646 // smoothed, the cell indices and levels may be 1647 // the same, but the vertex locations may not. 1648 // We adjust the vertices of the cells that have 1649 // no children (ie the active cells) to be 1650 // consistent with the global triangulation 1651 // later on and add assertions at that time 1652 // to guarantee the cells in the 1653 // local_triangulation are physically at the 1654 // same locations of the cells in the patch of 1655 // the global triangulation. 1656 } 1657 } 1658 // The parent cell whose children were added 1659 // into the map should be deleted from the map 1660 patch_to_global_tria_map_tmp.erase(cell); 1661 } 1662 } 1663 } 1664 } 1665 } 1666 while (refinement_necessary); 1667 1668 1669 // Last assertion check to make sure we have the right cells and centers 1670 // in the map, and hence the correct vertices of the triangulation 1671 for (typename Triangulation<Container::dimension, 1672 Container::space_dimension>::cell_iterator 1673 cell = local_triangulation.begin(); 1674 cell != local_triangulation.end(); 1675 ++cell) 1676 { 1677 if (cell->user_flag_set()) 1678 { 1679 Assert(patch_to_global_tria_map_tmp.find(cell) != 1680 patch_to_global_tria_map_tmp.end(), 1681 ExcInternalError()); 1682 1683 Assert(cell->center().distance( 1684 patch_to_global_tria_map_tmp[cell]->center()) <= 1685 1e-15 * cell->diameter(), 1686 ExcInternalError()); 1687 } 1688 } 1689 1690 1691 typename std::map< 1692 typename Triangulation<Container::dimension, 1693 Container::space_dimension>::cell_iterator, 1694 typename Container::cell_iterator>::iterator 1695 map_tmp_it = patch_to_global_tria_map_tmp.begin(), 1696 map_tmp_end = patch_to_global_tria_map_tmp.end(); 1697 // Now we just need to take the temporary map of pairs of type cell_iterator 1698 // "patch_to_global_tria_map_tmp" making pair of active_cell_iterators so 1699 // that filling out the final map "patch_to_global_tria_map" 1700 for (; map_tmp_it != map_tmp_end; ++map_tmp_it) 1701 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second; 1702 } 1703 1704 1705 1706 template <int dim, int spacedim> 1707 std::map< 1708 types::global_dof_index, 1709 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>> get_dof_to_support_patch_map(DoFHandler<dim,spacedim> & dof_handler)1710 get_dof_to_support_patch_map(DoFHandler<dim, spacedim> &dof_handler) 1711 { 1712 // This is the map from global_dof_index to 1713 // a set of cells on patch. We first map into 1714 // a set because it is very likely that we 1715 // will attempt to add a cell more than once 1716 // to a particular patch and we want to preserve 1717 // uniqueness of cell iterators. std::set does this 1718 // automatically for us. Later after it is all 1719 // constructed, we will copy to a map of vectors 1720 // since that is the preferred output for other 1721 // functions. 1722 std::map<types::global_dof_index, 1723 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>> 1724 dof_to_set_of_cells_map; 1725 1726 std::vector<types::global_dof_index> local_dof_indices; 1727 std::vector<types::global_dof_index> local_face_dof_indices; 1728 std::vector<types::global_dof_index> local_line_dof_indices; 1729 1730 // a place to save the dof_handler user flags and restore them later 1731 // to maintain const of dof_handler. 1732 std::vector<bool> user_flags; 1733 1734 1735 // in 3d, we need pointers from active lines to the 1736 // active parent lines, so we construct it as needed. 1737 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator, 1738 typename DoFHandler<dim, spacedim>::line_iterator> 1739 lines_to_parent_lines_map; 1740 if (dim == 3) 1741 { 1742 // save user flags as they will be modified and then later restored 1743 dof_handler.get_triangulation().save_user_flags(user_flags); 1744 const_cast<dealii::Triangulation<dim, spacedim> &>( 1745 dof_handler.get_triangulation()) 1746 .clear_user_flags(); 1747 1748 1749 typename DoFHandler<dim, spacedim>::active_cell_iterator 1750 cell = dof_handler.begin_active(), 1751 endc = dof_handler.end(); 1752 for (; cell != endc; ++cell) 1753 { 1754 // We only want lines that are locally_relevant 1755 // although it doesn't hurt to have lines that 1756 // are children of ghost cells since there are 1757 // few and we don't have to use them. 1758 if (cell->is_artificial() == false) 1759 { 1760 for (unsigned int l = 0; l < cell->n_lines(); ++l) 1761 if (cell->line(l)->has_children()) 1762 for (unsigned int c = 0; c < cell->line(l)->n_children(); 1763 ++c) 1764 { 1765 lines_to_parent_lines_map[cell->line(l)->child(c)] = 1766 cell->line(l); 1767 // set flags to know that child 1768 // line has an active parent. 1769 cell->line(l)->child(c)->set_user_flag(); 1770 } 1771 } 1772 } 1773 } 1774 1775 1776 // We loop through all cells and add cell to the 1777 // map for the dofs that it immediately touches 1778 // and then account for all the other dofs of 1779 // which it is a part, mainly the ones that must 1780 // be added on account of adaptivity hanging node 1781 // constraints. 1782 typename DoFHandler<dim, spacedim>::active_cell_iterator 1783 cell = dof_handler.begin_active(), 1784 endc = dof_handler.end(); 1785 for (; cell != endc; ++cell) 1786 { 1787 // Need to loop through all cells that could 1788 // be in the patch of dofs on locally_owned 1789 // cells including ghost cells 1790 if (cell->is_artificial() == false) 1791 { 1792 const unsigned int n_dofs_per_cell = 1793 cell->get_fe().n_dofs_per_cell(); 1794 local_dof_indices.resize(n_dofs_per_cell); 1795 1796 // Take care of adding cell pointer to each 1797 // dofs that exists on cell. 1798 cell->get_dof_indices(local_dof_indices); 1799 for (unsigned int i = 0; i < n_dofs_per_cell; ++i) 1800 dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell); 1801 1802 // In the case of the adjacent cell (over 1803 // faces or edges) being more refined, we 1804 // want to add all of the children to the 1805 // patch since the support function at that 1806 // dof could be non-zero along that entire 1807 // face (or line). 1808 1809 // Take care of dofs on neighbor faces 1810 for (const unsigned int f : cell->face_indices()) 1811 { 1812 if (cell->face(f)->has_children()) 1813 { 1814 for (unsigned int c = 0; c < cell->face(f)->n_children(); 1815 ++c) 1816 { 1817 // Add cell to dofs of all subfaces 1818 // 1819 // *-------------------*----------*---------* 1820 // | | add cell | | 1821 // | |<- to dofs| | 1822 // | |of subface| | 1823 // | cell *----------*---------* 1824 // | | add cell | | 1825 // | |<- to dofs| | 1826 // | |of subface| | 1827 // *-------------------*----------*---------* 1828 // 1829 Assert(cell->face(f)->child(c)->has_children() == false, 1830 ExcInternalError()); 1831 1832 const unsigned int n_dofs_per_face = 1833 cell->get_fe().n_dofs_per_face(); 1834 local_face_dof_indices.resize(n_dofs_per_face); 1835 1836 cell->face(f)->child(c)->get_dof_indices( 1837 local_face_dof_indices); 1838 for (unsigned int i = 0; i < n_dofs_per_face; ++i) 1839 dof_to_set_of_cells_map[local_face_dof_indices[i]] 1840 .insert(cell); 1841 } 1842 } 1843 else if ((cell->face(f)->at_boundary() == false) && 1844 (cell->neighbor_is_coarser(f))) 1845 { 1846 // Add cell to dofs of parent face and all 1847 // child faces of parent face 1848 // 1849 // *-------------------*----------*---------* 1850 // | | | | 1851 // | | cell | | 1852 // | add cell | | | 1853 // | to dofs -> *----------*---------* 1854 // | of parent | add cell | | 1855 // | face |<- to dofs| | 1856 // | |of subface| | 1857 // *-------------------*----------*---------* 1858 // 1859 1860 // Add cell to all dofs of parent face 1861 std::pair<unsigned int, unsigned int> 1862 neighbor_face_no_subface_no = 1863 cell->neighbor_of_coarser_neighbor(f); 1864 unsigned int face_no = neighbor_face_no_subface_no.first; 1865 unsigned int subface = neighbor_face_no_subface_no.second; 1866 1867 const unsigned int n_dofs_per_face = 1868 cell->get_fe().n_dofs_per_face(); 1869 local_face_dof_indices.resize(n_dofs_per_face); 1870 1871 cell->neighbor(f)->face(face_no)->get_dof_indices( 1872 local_face_dof_indices); 1873 for (unsigned int i = 0; i < n_dofs_per_face; ++i) 1874 dof_to_set_of_cells_map[local_face_dof_indices[i]].insert( 1875 cell); 1876 1877 // Add cell to all dofs of children of 1878 // parent face 1879 for (unsigned int c = 0; 1880 c < cell->neighbor(f)->face(face_no)->n_children(); 1881 ++c) 1882 { 1883 if (c != subface) // don't repeat work on dofs of 1884 // original cell 1885 { 1886 const unsigned int n_dofs_per_face = 1887 cell->get_fe().n_dofs_per_face(); 1888 local_face_dof_indices.resize(n_dofs_per_face); 1889 1890 Assert(cell->neighbor(f) 1891 ->face(face_no) 1892 ->child(c) 1893 ->has_children() == false, 1894 ExcInternalError()); 1895 cell->neighbor(f) 1896 ->face(face_no) 1897 ->child(c) 1898 ->get_dof_indices(local_face_dof_indices); 1899 for (unsigned int i = 0; i < n_dofs_per_face; ++i) 1900 dof_to_set_of_cells_map[local_face_dof_indices[i]] 1901 .insert(cell); 1902 } 1903 } 1904 } 1905 } 1906 1907 1908 // If 3d, take care of dofs on lines in the 1909 // same pattern as faces above. That is, if 1910 // a cell's line has children, distribute 1911 // cell to dofs of children of line, and 1912 // if cell's line has an active parent, then 1913 // distribute cell to dofs on parent line 1914 // and dofs on all children of parent line. 1915 if (dim == 3) 1916 { 1917 for (unsigned int l = 0; l < cell->n_lines(); ++l) 1918 { 1919 if (cell->line(l)->has_children()) 1920 { 1921 for (unsigned int c = 0; 1922 c < cell->line(l)->n_children(); 1923 ++c) 1924 { 1925 Assert(cell->line(l)->child(c)->has_children() == 1926 false, 1927 ExcInternalError()); 1928 1929 // dofs_per_line returns number of dofs 1930 // on line not including the vertices of the line. 1931 const unsigned int n_dofs_per_line = 1932 2 * cell->get_fe().n_dofs_per_vertex() + 1933 cell->get_fe().n_dofs_per_line(); 1934 local_line_dof_indices.resize(n_dofs_per_line); 1935 1936 cell->line(l)->child(c)->get_dof_indices( 1937 local_line_dof_indices); 1938 for (unsigned int i = 0; i < n_dofs_per_line; ++i) 1939 dof_to_set_of_cells_map[local_line_dof_indices[i]] 1940 .insert(cell); 1941 } 1942 } 1943 // user flag was set above to denote that 1944 // an active parent line exists so add 1945 // cell to dofs of parent and all it's 1946 // children 1947 else if (cell->line(l)->user_flag_set() == true) 1948 { 1949 typename DoFHandler<dim, spacedim>::line_iterator 1950 parent_line = 1951 lines_to_parent_lines_map[cell->line(l)]; 1952 Assert(parent_line->has_children(), ExcInternalError()); 1953 1954 // dofs_per_line returns number of dofs 1955 // on line not including the vertices of the line. 1956 const unsigned int n_dofs_per_line = 1957 2 * cell->get_fe().n_dofs_per_vertex() + 1958 cell->get_fe().n_dofs_per_line(); 1959 local_line_dof_indices.resize(n_dofs_per_line); 1960 1961 parent_line->get_dof_indices(local_line_dof_indices); 1962 for (unsigned int i = 0; i < n_dofs_per_line; ++i) 1963 dof_to_set_of_cells_map[local_line_dof_indices[i]] 1964 .insert(cell); 1965 1966 for (unsigned int c = 0; c < parent_line->n_children(); 1967 ++c) 1968 { 1969 Assert(parent_line->child(c)->has_children() == 1970 false, 1971 ExcInternalError()); 1972 1973 const unsigned int n_dofs_per_line = 1974 2 * cell->get_fe().n_dofs_per_vertex() + 1975 cell->get_fe().n_dofs_per_line(); 1976 local_line_dof_indices.resize(n_dofs_per_line); 1977 1978 parent_line->child(c)->get_dof_indices( 1979 local_line_dof_indices); 1980 for (unsigned int i = 0; i < n_dofs_per_line; ++i) 1981 dof_to_set_of_cells_map[local_line_dof_indices[i]] 1982 .insert(cell); 1983 } 1984 } 1985 } // for lines l 1986 } // if dim == 3 1987 } // if cell->is_locally_owned() 1988 } // for cells 1989 1990 1991 if (dim == 3) 1992 { 1993 // finally, restore user flags that were changed above 1994 // to when we constructed the pointers to parent of lines 1995 // Since dof_handler is const, we must leave it unchanged. 1996 const_cast<dealii::Triangulation<dim, spacedim> &>( 1997 dof_handler.get_triangulation()) 1998 .load_user_flags(user_flags); 1999 } 2000 2001 // Finally, we copy map of sets to 2002 // map of vectors using the std::vector::assign() function 2003 std::map< 2004 types::global_dof_index, 2005 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>> 2006 dof_to_cell_patches; 2007 2008 typename std::map< 2009 types::global_dof_index, 2010 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>:: 2011 iterator it = dof_to_set_of_cells_map.begin(), 2012 it_end = dof_to_set_of_cells_map.end(); 2013 for (; it != it_end; ++it) 2014 dof_to_cell_patches[it->first].assign(it->second.begin(), 2015 it->second.end()); 2016 2017 return dof_to_cell_patches; 2018 } 2019 2020 /* 2021 * Internally used in collect_periodic_faces 2022 */ 2023 template <typename CellIterator> 2024 void match_periodic_face_pairs(std::set<std::pair<CellIterator,unsigned int>> & pairs1,std::set<std::pair<typename identity<CellIterator>::type,unsigned int>> & pairs2,const int direction,std::vector<PeriodicFacePair<CellIterator>> & matched_pairs,const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> & offset,const FullMatrix<double> & matrix)2025 match_periodic_face_pairs( 2026 std::set<std::pair<CellIterator, unsigned int>> &pairs1, 2027 std::set<std::pair<typename identity<CellIterator>::type, unsigned int>> 2028 & pairs2, 2029 const int direction, 2030 std::vector<PeriodicFacePair<CellIterator>> &matched_pairs, 2031 const dealii::Tensor<1, CellIterator::AccessorType::space_dimension> 2032 & offset, 2033 const FullMatrix<double> &matrix) 2034 { 2035 static const int space_dim = CellIterator::AccessorType::space_dimension; 2036 (void)space_dim; 2037 AssertIndexRange(direction, space_dim); 2038 2039 #ifdef DEBUG 2040 { 2041 constexpr int dim = CellIterator::AccessorType::dimension; 2042 constexpr int spacedim = CellIterator::AccessorType::space_dimension; 2043 // For parallel::fullydistributed::Triangulation there might be unmatched 2044 // faces on periodic boundaries on the coarse grid, which results that 2045 // this assert is not fulfilled (not a bug!). See also the discussion in 2046 // the method collect_periodic_faces. 2047 if (!(((pairs1.size() > 0) && 2048 (dynamic_cast<const parallel::fullydistributed:: 2049 Triangulation<dim, spacedim> *>( 2050 &pairs1.begin()->first->get_triangulation()) != nullptr)) || 2051 ((pairs2.size() > 0) && 2052 (dynamic_cast< 2053 const parallel::fullydistributed::Triangulation<dim, spacedim> 2054 *>(&pairs2.begin()->first->get_triangulation()) != nullptr)))) 2055 Assert(pairs1.size() == pairs2.size(), 2056 ExcMessage("Unmatched faces on periodic boundaries")); 2057 } 2058 #endif 2059 2060 unsigned int n_matches = 0; 2061 2062 // Match with a complexity of O(n^2). This could be improved... 2063 std::bitset<3> orientation; 2064 using PairIterator = 2065 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator; 2066 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1) 2067 { 2068 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2) 2069 { 2070 const CellIterator cell1 = it1->first; 2071 const CellIterator cell2 = it2->first; 2072 const unsigned int face_idx1 = it1->second; 2073 const unsigned int face_idx2 = it2->second; 2074 if (GridTools::orthogonal_equality(orientation, 2075 cell1->face(face_idx1), 2076 cell2->face(face_idx2), 2077 direction, 2078 offset, 2079 matrix)) 2080 { 2081 // We have a match, so insert the matching pairs and 2082 // remove the matched cell in pairs2 to speed up the 2083 // matching: 2084 const PeriodicFacePair<CellIterator> matched_face = { 2085 {cell1, cell2}, {face_idx1, face_idx2}, orientation, matrix}; 2086 matched_pairs.push_back(matched_face); 2087 pairs2.erase(it2); 2088 ++n_matches; 2089 break; 2090 } 2091 } 2092 } 2093 2094 // Assure that all faces are matched if not 2095 // parallel::fullydistributed::Triangulation is used. This is related to the 2096 // fact that the faces might not be successfully matched on the coarse grid 2097 // (not a bug!). See also the comment above and in the 2098 // method collect_periodic_faces. 2099 { 2100 constexpr int dim = CellIterator::AccessorType::dimension; 2101 constexpr int spacedim = CellIterator::AccessorType::space_dimension; 2102 if (!(((pairs1.size() > 0) && 2103 (dynamic_cast<const parallel::fullydistributed:: 2104 Triangulation<dim, spacedim> *>( 2105 &pairs1.begin()->first->get_triangulation()) != nullptr)) || 2106 ((pairs2.size() > 0) && 2107 (dynamic_cast< 2108 const parallel::fullydistributed::Triangulation<dim, spacedim> 2109 *>(&pairs2.begin()->first->get_triangulation()) != nullptr)))) 2110 AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0, 2111 ExcMessage("Unmatched faces on periodic boundaries")); 2112 } 2113 } 2114 2115 2116 2117 template <typename MeshType> 2118 void collect_periodic_faces(const MeshType & mesh,const types::boundary_id b_id,const int direction,std::vector<PeriodicFacePair<typename MeshType::cell_iterator>> & matched_pairs,const Tensor<1,MeshType::space_dimension> & offset,const FullMatrix<double> & matrix)2119 collect_periodic_faces( 2120 const MeshType & mesh, 2121 const types::boundary_id b_id, 2122 const int direction, 2123 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>> 2124 & matched_pairs, 2125 const Tensor<1, MeshType::space_dimension> &offset, 2126 const FullMatrix<double> & matrix) 2127 { 2128 static const int dim = MeshType::dimension; 2129 static const int space_dim = MeshType::space_dimension; 2130 (void)dim; 2131 (void)space_dim; 2132 AssertIndexRange(direction, space_dim); 2133 2134 Assert(dim == space_dim, ExcNotImplemented()); 2135 2136 // Loop over all cells on the highest level and collect all boundary 2137 // faces 2*direction and 2*direction*1: 2138 2139 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1; 2140 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2; 2141 2142 for (typename MeshType::cell_iterator cell = mesh.begin(0); 2143 cell != mesh.end(0); 2144 ++cell) 2145 { 2146 const typename MeshType::face_iterator face_1 = 2147 cell->face(2 * direction); 2148 const typename MeshType::face_iterator face_2 = 2149 cell->face(2 * direction + 1); 2150 2151 if (face_1->at_boundary() && face_1->boundary_id() == b_id) 2152 { 2153 const std::pair<typename MeshType::cell_iterator, unsigned int> 2154 pair1 = std::make_pair(cell, 2 * direction); 2155 pairs1.insert(pair1); 2156 } 2157 2158 if (face_2->at_boundary() && face_2->boundary_id() == b_id) 2159 { 2160 const std::pair<typename MeshType::cell_iterator, unsigned int> 2161 pair2 = std::make_pair(cell, 2 * direction + 1); 2162 pairs2.insert(pair2); 2163 } 2164 } 2165 2166 Assert(pairs1.size() == pairs2.size(), 2167 ExcMessage("Unmatched faces on periodic boundaries")); 2168 2169 Assert(pairs1.size() > 0, 2170 ExcMessage("No new periodic face pairs have been found. " 2171 "Are you sure that you've selected the correct boundary " 2172 "id's and that the coarsest level mesh is colorized?")); 2173 2174 #ifdef DEBUG 2175 const unsigned int size_old = matched_pairs.size(); 2176 #endif 2177 2178 // and call match_periodic_face_pairs that does the actual matching: 2179 match_periodic_face_pairs( 2180 pairs1, pairs2, direction, matched_pairs, offset, matrix); 2181 2182 #ifdef DEBUG 2183 // check for standard orientation 2184 const unsigned int size_new = matched_pairs.size(); 2185 for (unsigned int i = size_old; i < size_new; ++i) 2186 { 2187 Assert(matched_pairs[i].orientation == 1, 2188 ExcMessage( 2189 "Found a face match with non standard orientation. " 2190 "This function is only suitable for meshes with cells " 2191 "in default orientation")); 2192 } 2193 #endif 2194 } 2195 2196 2197 2198 template <typename MeshType> 2199 void collect_periodic_faces(const MeshType & mesh,const types::boundary_id b_id1,const types::boundary_id b_id2,const int direction,std::vector<PeriodicFacePair<typename MeshType::cell_iterator>> & matched_pairs,const Tensor<1,MeshType::space_dimension> & offset,const FullMatrix<double> & matrix)2200 collect_periodic_faces( 2201 const MeshType & mesh, 2202 const types::boundary_id b_id1, 2203 const types::boundary_id b_id2, 2204 const int direction, 2205 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>> 2206 & matched_pairs, 2207 const Tensor<1, MeshType::space_dimension> &offset, 2208 const FullMatrix<double> & matrix) 2209 { 2210 static const int dim = MeshType::dimension; 2211 static const int space_dim = MeshType::space_dimension; 2212 (void)dim; 2213 (void)space_dim; 2214 AssertIndexRange(direction, space_dim); 2215 2216 // Loop over all cells on the highest level and collect all boundary 2217 // faces belonging to b_id1 and b_id2: 2218 2219 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1; 2220 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2; 2221 2222 for (typename MeshType::cell_iterator cell = mesh.begin(0); 2223 cell != mesh.end(0); 2224 ++cell) 2225 { 2226 for (unsigned int i : cell->face_indices()) 2227 { 2228 const typename MeshType::face_iterator face = cell->face(i); 2229 if (face->at_boundary() && face->boundary_id() == b_id1) 2230 { 2231 const std::pair<typename MeshType::cell_iterator, unsigned int> 2232 pair1 = std::make_pair(cell, i); 2233 pairs1.insert(pair1); 2234 } 2235 2236 if (face->at_boundary() && face->boundary_id() == b_id2) 2237 { 2238 const std::pair<typename MeshType::cell_iterator, unsigned int> 2239 pair2 = std::make_pair(cell, i); 2240 pairs2.insert(pair2); 2241 } 2242 } 2243 } 2244 2245 // Assure that all faces are matched on the coare grid. This requirement 2246 // can only fulfilled if a process owns the complete coarse grid. This is 2247 // not the case for a parallel::fullydistributed::Triangulation, i.e. this 2248 // requirement has not to be met neither (consider faces on the outside of a 2249 // ghost cell that are periodic but for which the ghost neighbor doesn't 2250 // exist). 2251 if (!(((pairs1.size() > 0) && 2252 (dynamic_cast< 2253 const parallel::fullydistributed::Triangulation<dim, space_dim> 2254 *>(&pairs1.begin()->first->get_triangulation()) != nullptr)) || 2255 ((pairs2.size() > 0) && 2256 (dynamic_cast< 2257 const parallel::fullydistributed::Triangulation<dim, space_dim> 2258 *>(&pairs2.begin()->first->get_triangulation()) != nullptr)))) 2259 Assert(pairs1.size() == pairs2.size(), 2260 ExcMessage("Unmatched faces on periodic boundaries")); 2261 2262 Assert( 2263 (pairs1.size() > 0 || 2264 (dynamic_cast< 2265 const parallel::fullydistributed::Triangulation<dim, space_dim> *>( 2266 &mesh.begin()->get_triangulation()) != nullptr)), 2267 ExcMessage("No new periodic face pairs have been found. " 2268 "Are you sure that you've selected the correct boundary " 2269 "id's and that the coarsest level mesh is colorized?")); 2270 2271 // and call match_periodic_face_pairs that does the actual matching: 2272 match_periodic_face_pairs( 2273 pairs1, pairs2, direction, matched_pairs, offset, matrix); 2274 } 2275 2276 2277 2278 /* 2279 * Internally used in orthogonal_equality 2280 * 2281 * An orthogonal equality test for points: 2282 * 2283 * point1 and point2 are considered equal, if 2284 * matrix.point1 + offset - point2 2285 * is parallel to the unit vector in <direction> 2286 */ 2287 template <int spacedim> 2288 inline bool orthogonal_equality(const Point<spacedim> & point1,const Point<spacedim> & point2,const int direction,const Tensor<1,spacedim> & offset,const FullMatrix<double> & matrix)2289 orthogonal_equality(const Point<spacedim> & point1, 2290 const Point<spacedim> & point2, 2291 const int direction, 2292 const Tensor<1, spacedim> &offset, 2293 const FullMatrix<double> & matrix) 2294 { 2295 AssertIndexRange(direction, spacedim); 2296 2297 Assert(matrix.m() == matrix.n(), ExcInternalError()); 2298 2299 Point<spacedim> distance; 2300 2301 if (matrix.m() == spacedim) 2302 for (int i = 0; i < spacedim; ++i) 2303 for (int j = 0; j < spacedim; ++j) 2304 distance(i) += matrix(i, j) * point1(j); 2305 else 2306 distance = point1; 2307 2308 distance += offset - point2; 2309 2310 for (int i = 0; i < spacedim; ++i) 2311 { 2312 // Only compare coordinate-components != direction: 2313 if (i == direction) 2314 continue; 2315 2316 if (std::abs(distance(i)) > 1.e-10) 2317 return false; 2318 } 2319 2320 return true; 2321 } 2322 2323 2324 /* 2325 * Internally used in orthogonal_equality 2326 * 2327 * A lookup table to transform vertex matchings to orientation flags of 2328 * the form (face_orientation, face_flip, face_rotation) 2329 * 2330 * See the comment on the next function as well as the detailed 2331 * documentation of make_periodicity_constraints and 2332 * collect_periodic_faces for details 2333 */ 2334 template <int dim> 2335 struct OrientationLookupTable 2336 {}; 2337 2338 template <> 2339 struct OrientationLookupTable<1> 2340 { 2341 using MATCH_T = 2342 std::array<unsigned int, GeometryInfo<1>::vertices_per_face>; 2343 static inline std::bitset<3> lookupGridTools::OrientationLookupTable2344 lookup(const MATCH_T &) 2345 { 2346 // The 1D case is trivial 2347 return 1; // [true ,false,false] 2348 } 2349 }; 2350 2351 template <> 2352 struct OrientationLookupTable<2> 2353 { 2354 using MATCH_T = 2355 std::array<unsigned int, GeometryInfo<2>::vertices_per_face>; 2356 static inline std::bitset<3> lookupGridTools::OrientationLookupTable2357 lookup(const MATCH_T &matching) 2358 { 2359 // In 2D matching faces (=lines) results in two cases: Either 2360 // they are aligned or flipped. We store this "line_flip" 2361 // property somewhat sloppy as "face_flip" 2362 // (always: face_orientation = true, face_rotation = false) 2363 2364 static const MATCH_T m_tff = {{0, 1}}; 2365 if (matching == m_tff) 2366 return 1; // [true ,false,false] 2367 static const MATCH_T m_ttf = {{1, 0}}; 2368 if (matching == m_ttf) 2369 return 3; // [true ,true ,false] 2370 Assert(false, ExcInternalError()); 2371 // what follows is dead code, but it avoids warnings about the lack 2372 // of a return value 2373 return 0; 2374 } 2375 }; 2376 2377 template <> 2378 struct OrientationLookupTable<3> 2379 { 2380 using MATCH_T = 2381 std::array<unsigned int, GeometryInfo<3>::vertices_per_face>; 2382 static inline std::bitset<3> lookupGridTools::OrientationLookupTable2383 lookup(const MATCH_T &matching) 2384 { 2385 // The full fledged 3D case. *Yay* 2386 // See the documentation in include/deal.II/base/geometry_info.h 2387 // as well as the actual implementation in source/grid/tria.cc 2388 // for more details... 2389 2390 static const MATCH_T m_tff = {{0, 1, 2, 3}}; 2391 if (matching == m_tff) 2392 return 1; // [true ,false,false] 2393 static const MATCH_T m_tft = {{1, 3, 0, 2}}; 2394 if (matching == m_tft) 2395 return 5; // [true ,false,true ] 2396 static const MATCH_T m_ttf = {{3, 2, 1, 0}}; 2397 if (matching == m_ttf) 2398 return 3; // [true ,true ,false] 2399 static const MATCH_T m_ttt = {{2, 0, 3, 1}}; 2400 if (matching == m_ttt) 2401 return 7; // [true ,true ,true ] 2402 static const MATCH_T m_fff = {{0, 2, 1, 3}}; 2403 if (matching == m_fff) 2404 return 0; // [false,false,false] 2405 static const MATCH_T m_fft = {{2, 3, 0, 1}}; 2406 if (matching == m_fft) 2407 return 4; // [false,false,true ] 2408 static const MATCH_T m_ftf = {{3, 1, 2, 0}}; 2409 if (matching == m_ftf) 2410 return 2; // [false,true ,false] 2411 static const MATCH_T m_ftt = {{1, 0, 3, 2}}; 2412 if (matching == m_ftt) 2413 return 6; // [false,true ,true ] 2414 Assert(false, ExcInternalError()); 2415 // what follows is dead code, but it avoids warnings about the lack 2416 // of a return value 2417 return 0; 2418 } 2419 }; 2420 2421 2422 2423 template <typename FaceIterator> orthogonal_equality(std::bitset<3> & orientation,const FaceIterator & face1,const FaceIterator & face2,const int direction,const Tensor<1,FaceIterator::AccessorType::space_dimension> & offset,const FullMatrix<double> & matrix)2424 inline bool orthogonal_equality( 2425 std::bitset<3> & orientation, 2426 const FaceIterator & face1, 2427 const FaceIterator & face2, 2428 const int direction, 2429 const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset, 2430 const FullMatrix<double> & matrix) 2431 { 2432 Assert(matrix.m() == matrix.n(), 2433 ExcMessage("The supplied matrix must be a square matrix")); 2434 2435 static const int dim = FaceIterator::AccessorType::dimension; 2436 2437 // Do a full matching of the face vertices: 2438 2439 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching; 2440 2441 AssertDimension(face1->n_vertices(), face2->n_vertices()); 2442 2443 std::set<unsigned int> face2_vertices; 2444 for (unsigned int i = 0; i < face1->n_vertices(); ++i) 2445 face2_vertices.insert(i); 2446 2447 for (unsigned int i = 0; i < face1->n_vertices(); ++i) 2448 { 2449 for (std::set<unsigned int>::iterator it = face2_vertices.begin(); 2450 it != face2_vertices.end(); 2451 ++it) 2452 { 2453 if (orthogonal_equality(face1->vertex(i), 2454 face2->vertex(*it), 2455 direction, 2456 offset, 2457 matrix)) 2458 { 2459 matching[i] = *it; 2460 face2_vertices.erase(it); 2461 break; // jump out of the innermost loop 2462 } 2463 } 2464 } 2465 2466 // And finally, a lookup to determine the ordering bitmask: 2467 if (face2_vertices.empty()) 2468 orientation = OrientationLookupTable<dim>::lookup(matching); 2469 2470 return face2_vertices.empty(); 2471 } 2472 2473 2474 2475 template <typename FaceIterator> 2476 inline bool orthogonal_equality(const FaceIterator & face1,const FaceIterator & face2,const int direction,const Tensor<1,FaceIterator::AccessorType::space_dimension> & offset,const FullMatrix<double> & matrix)2477 orthogonal_equality( 2478 const FaceIterator & face1, 2479 const FaceIterator & face2, 2480 const int direction, 2481 const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset, 2482 const FullMatrix<double> & matrix) 2483 { 2484 // Call the function above with a dummy orientation array 2485 std::bitset<3> dummy; 2486 return orthogonal_equality(dummy, face1, face2, direction, offset, matrix); 2487 } 2488 2489 2490 2491 } // namespace GridTools 2492 2493 2494 #include "grid_tools_dof_handlers.inst" 2495 2496 2497 DEAL_II_NAMESPACE_CLOSE 2498