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 #ifndef dealii_grid_tools_h 17 # define dealii_grid_tools_h 18 19 20 # include <deal.II/base/config.h> 21 22 # include <deal.II/base/bounding_box.h> 23 # include <deal.II/base/geometry_info.h> 24 # include <deal.II/base/std_cxx17/optional.h> 25 26 # include <deal.II/boost_adaptors/bounding_box.h> 27 28 # include <deal.II/dofs/dof_handler.h> 29 30 # include <deal.II/fe/mapping.h> 31 # include <deal.II/fe/mapping_q1.h> 32 33 # include <deal.II/grid/manifold.h> 34 # include <deal.II/grid/tria.h> 35 # include <deal.II/grid/tria_accessor.h> 36 # include <deal.II/grid/tria_iterator.h> 37 38 # include <deal.II/hp/dof_handler.h> 39 40 # include <deal.II/lac/sparsity_tools.h> 41 42 # include <deal.II/numerics/rtree.h> 43 44 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS 45 # include <boost/archive/binary_iarchive.hpp> 46 # include <boost/archive/binary_oarchive.hpp> 47 # include <boost/geometry/index/rtree.hpp> 48 # include <boost/serialization/array.hpp> 49 # include <boost/serialization/vector.hpp> 50 51 # ifdef DEAL_II_WITH_ZLIB 52 # include <boost/iostreams/device/back_inserter.hpp> 53 # include <boost/iostreams/filter/gzip.hpp> 54 # include <boost/iostreams/filtering_stream.hpp> 55 # include <boost/iostreams/stream.hpp> 56 # endif 57 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS 58 59 # include <bitset> 60 # include <list> 61 # include <set> 62 63 DEAL_II_NAMESPACE_OPEN 64 65 // Forward declarations 66 # ifndef DOXYGEN 67 namespace parallel 68 { 69 namespace distributed 70 { 71 template <int, int> 72 class Triangulation; 73 } 74 } // namespace parallel 75 76 namespace hp 77 { 78 template <int, int> 79 class MappingCollection; 80 } 81 82 class SparsityPattern; 83 # endif 84 85 namespace internal 86 { 87 template <int dim, int spacedim, class MeshType> 88 class ActiveCellIterator 89 { 90 public: 91 # ifndef _MSC_VER 92 using type = typename MeshType::active_cell_iterator; 93 # else 94 using type = TriaActiveIterator<dealii::CellAccessor<dim, spacedim>>; 95 # endif 96 }; 97 98 # ifdef _MSC_VER 99 template <int dim, int spacedim> 100 class ActiveCellIterator<dim, spacedim, dealii::DoFHandler<dim, spacedim>> 101 { 102 public: 103 using type = 104 TriaActiveIterator<dealii::DoFCellAccessor<dim, spacedim, false>>; 105 }; 106 # endif 107 108 109 # ifdef _MSC_VER 110 template <int dim, int spacedim> 111 class ActiveCellIterator<dim, spacedim, dealii::hp::DoFHandler<dim, spacedim>> 112 { 113 public: 114 using type = 115 TriaActiveIterator<dealii::DoFCellAccessor<dim, spacedim, false>>; 116 }; 117 # endif 118 } // namespace internal 119 120 /** 121 * This namespace is a collection of algorithms working on triangulations, 122 * such as shifting or rotating triangulations, but also finding a cell that 123 * contains a given point. See the descriptions of the individual functions 124 * for more information. 125 * 126 * @ingroup grid 127 */ 128 namespace GridTools 129 { 130 template <int dim, int spacedim> 131 class Cache; 132 133 /** 134 * @name Information about meshes and cells 135 */ 136 /*@{*/ 137 138 /** 139 * Return the diameter of a triangulation. The diameter is computed using 140 * only the vertices, i.e. if the diameter should be larger than the maximal 141 * distance between boundary vertices due to a higher order mapping, then 142 * this function will not catch this. 143 */ 144 template <int dim, int spacedim> 145 double 146 diameter(const Triangulation<dim, spacedim> &tria); 147 148 /** 149 * Compute the volume (i.e. the dim-dimensional measure) of the 150 * triangulation. We compute the measure using the integral $\sum_K \int_K 1 151 * \; dx$ where $K$ are the cells of the given triangulation. The integral 152 * is approximated via quadrature for which we need the mapping argument. 153 * 154 * If the triangulation is a dim-dimensional one embedded in a higher 155 * dimensional space of dimension spacedim, then the value returned is the 156 * dim-dimensional measure. For example, for a two-dimensional triangulation 157 * in three-dimensional space, the value returned is the area of the surface 158 * so described. (This obviously makes sense since the spacedim-dimensional 159 * measure of a dim-dimensional triangulation would always be zero if dim @< 160 * spacedim. 161 * 162 * This function also works for objects of type 163 * parallel::distributed::Triangulation, in which case the function is a 164 * collective operation. 165 * 166 * @param tria The triangulation. 167 * @param mapping An optional argument used to denote the mapping that 168 * should be used when describing whether cells are bounded by straight or 169 * curved faces. The default is to use a $Q_1$ mapping, which corresponds to 170 * straight lines bounding the cells. 171 * @return The dim-dimensional measure of the domain described by the 172 * triangulation, as discussed above. 173 */ 174 template <int dim, int spacedim> 175 double 176 volume(const Triangulation<dim, spacedim> &tria, 177 const Mapping<dim, spacedim> & mapping = 178 (StaticMappingQ1<dim, spacedim>::mapping)); 179 180 /** 181 * Return an approximation of the diameter of the smallest active cell of a 182 * triangulation. See step-24 for an example of use of this function. 183 * 184 * Notice that, even if you pass a non-trivial mapping, the returned value is 185 * computed only using information on the vertices of the triangulation, 186 * possibly transformed by the mapping. While this is accurate most of the 187 * times, it may fail to give the correct result when the triangulation 188 * contains very distorted cells. 189 */ 190 template <int dim, int spacedim> 191 double 192 minimal_cell_diameter(const Triangulation<dim, spacedim> &triangulation, 193 const Mapping<dim, spacedim> & mapping = 194 (StaticMappingQ1<dim, spacedim>::mapping)); 195 196 /** 197 * Return an approximation of the diameter of the largest active cell of a 198 * triangulation. 199 * 200 * Notice that, even if you pass a non-trivial mapping to this function, the 201 * returned value is computed only using information on the vertices of the 202 * triangulation, possibly transformed by the mapping. While this is accurate 203 * most of the times, it may fail to give the correct result when the 204 * triangulation contains very distorted cells. 205 */ 206 template <int dim, int spacedim> 207 double 208 maximal_cell_diameter(const Triangulation<dim, spacedim> &triangulation, 209 const Mapping<dim, spacedim> & mapping = 210 (StaticMappingQ1<dim, spacedim>::mapping)); 211 212 /** 213 * Given a list of vertices (typically obtained using 214 * Triangulation::get_vertices) as the first, and a list of vertex indices 215 * that characterize a single cell as the second argument, return the 216 * measure (area, volume) of this cell. If this is a real cell, then you can 217 * get the same result using <code>cell-@>measure()</code>, but this 218 * function also works for cells that do not exist except that you make it 219 * up by naming its vertices from the list. 220 */ 221 template <int dim> 222 double 223 cell_measure( 224 const std::vector<Point<dim>> &all_vertices, 225 const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]); 226 227 /** 228 * A variant of cell_measure() accepting an ArrayView instead of a 229 * fixed-sized array for @p vertex_indices. 230 * 231 * The parameter @p vertex_indices is expected to have 232 * GeometryInfo<dim>::vertices_per_cell entries. A std::vector is implicitly 233 * convertible to an ArrayView, so it can be passed directly. See the 234 * ArrayView class for more information. 235 */ 236 template <int dim> 237 double 238 cell_measure(const std::vector<Point<dim>> & all_vertices, 239 const ArrayView<const unsigned int> &vertex_indices); 240 241 /** 242 * A version of the function above that can accept input for nonzero 243 * codimension cases. This function only exists to aid generic programming 244 * and calling it will just raise an exception. 245 */ 246 template <int dim, typename T> 247 double 248 cell_measure(const T &, ...); 249 250 /** 251 * Computes an aspect ratio measure for all locally-owned active cells and 252 * fills a vector with one entry per cell, given a @p triangulation and 253 * @p mapping. The size of the vector that is returned equals the number of 254 * active cells. The vector contains zero for non locally-owned cells. The 255 * aspect ratio of a cell is defined as the ratio of the maximum to minimum 256 * singular value of the Jacobian, taking the maximum over all quadrature 257 * points of a quadrature rule specified via @p quadrature. For example, for 258 * the special case of rectangular elements in 2d with dimensions $a$ and $b$ 259 * ($a \geq b$), this function returns the usual aspect ratio definition 260 * $a/b$. The above definition using singular values is a generalization to 261 * arbitrarily deformed elements. This function is intended to be used for 262 * $d=2,3$ space dimensions, but it can also be used for $d=1$ returning a 263 * value of 1. 264 * 265 * @note Inverted elements do not throw an exception. Instead, a value of inf 266 * is written into the vector in case of inverted elements. 267 * 268 * @note Make sure to use enough quadrature points for a precise calculation 269 * of the aspect ratio in case of deformed elements. 270 * 271 * @note In parallel computations the return value will have the length 272 * n_active_cells but the aspect ratio is only computed for the cells that 273 * are locally owned and placed at index CellAccessor::active_cell_index(), 274 * respectively. All other values are set to 0. 275 */ 276 template <int dim> 277 Vector<double> 278 compute_aspect_ratio_of_cells(const Mapping<dim> & mapping, 279 const Triangulation<dim> &triangulation, 280 const Quadrature<dim> & quadrature); 281 282 /** 283 * Computes the maximum aspect ratio by taking the maximum over all cells. 284 * 285 * @note When running in parallel with a Triangulation that supports MPI, 286 * this is a collective call and the return value is the maximum over all 287 * processors. 288 */ 289 template <int dim> 290 double 291 compute_maximum_aspect_ratio(const Mapping<dim> & mapping, 292 const Triangulation<dim> &triangulation, 293 const Quadrature<dim> & quadrature); 294 295 /** 296 * Compute the smallest box containing the entire triangulation. 297 * 298 * If the input triangulation is a `parallel::distributed::Triangulation`, 299 * then each processor will compute a bounding box enclosing all locally 300 * owned, ghost, and artificial cells. In the case of a domain without curved 301 * boundaries, these bounding boxes will all agree between processors because 302 * the union of the areas occupied by artificial and ghost cells equals the 303 * union of the areas occupied by the cells that other processors own. 304 * However, if the domain has curved boundaries, this is no longer the case. 305 * The bounding box returned may be appropriate for the current processor, 306 * but different from the bounding boxes computed on other processors. 307 */ 308 template <int dim, int spacedim> 309 BoundingBox<spacedim> 310 compute_bounding_box(const Triangulation<dim, spacedim> &triangulation); 311 312 /** 313 * Return the point on the geometrical object @p object closest to the given 314 * point @p trial_point. For example, if @p object is a one-dimensional line 315 * or edge, then the returned point will be a point on the geodesic that 316 * connects the vertices as the manifold associated with the object sees it 317 * (i.e., the geometric line may be curved if it lives in a higher 318 * dimensional space). If the iterator points to a quadrilateral in a higher 319 * dimensional space, then the returned point lies within the convex hull of 320 * the vertices of the quad as seen by the associated manifold. 321 * 322 * @note This projection is usually not well-posed since there may be 323 * multiple points on the object that minimize the distance. The algorithm 324 * used in this function is robust (and the output is guaranteed to be on 325 * the given @p object) but may only provide a few correct digits if the 326 * object has high curvature. If your manifold supports it then the 327 * specialized function Manifold::project_to_manifold() may perform better. 328 */ 329 template <typename Iterator> 330 Point<Iterator::AccessorType::space_dimension> 331 project_to_object( 332 const Iterator & object, 333 const Point<Iterator::AccessorType::space_dimension> &trial_point); 334 335 /** 336 * Return the arrays that define the coarse mesh of a Triangulation. This 337 * function is the inverse of Triangulation::create_triangulation(). 338 * 339 * The return value is a tuple with the vector of vertices, the vector of 340 * cells, and a SubCellData structure. The latter contains additional 341 * information about faces and lines. 342 * 343 * This function is useful in cases where one needs to deconstruct a 344 * Triangulation or manipulate the numbering of the vertices in some way: an 345 * example is GridGenerator::merge_triangulations(). 346 */ 347 template <int dim, int spacedim> 348 std:: 349 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData> 350 get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria); 351 352 /*@}*/ 353 /** 354 * @name Functions supporting the creation of meshes 355 */ 356 /*@{*/ 357 358 /** 359 * Remove vertices that are not referenced by any of the cells. This 360 * function is called by all <tt>GridIn::read_*</tt> functions to eliminate 361 * vertices that are listed in the input files but are not used by the cells 362 * in the input file. While these vertices should not be in the input from 363 * the beginning, they sometimes are, most often when some cells have been 364 * removed by hand without wanting to update the vertex lists, as they might 365 * be lengthy. 366 * 367 * This function is called by all <tt>GridIn::read_*</tt> functions as the 368 * triangulation class requires them to be called with used vertices only. 369 * This is so, since the vertices are copied verbatim by that class, so we 370 * have to eliminate unused vertices beforehand. 371 * 372 * Not implemented for the codimension one case. 373 */ 374 template <int dim, int spacedim> 375 void 376 delete_unused_vertices(std::vector<Point<spacedim>> &vertices, 377 std::vector<CellData<dim>> & cells, 378 SubCellData & subcelldata); 379 380 /** 381 * Remove vertices that are duplicated, due to the input of a structured 382 * grid, for example. If these vertices are not removed, the faces bounded 383 * by these vertices become part of the boundary, even if they are in the 384 * interior of the mesh. 385 * 386 * This function is called by some <tt>GridIn::read_*</tt> functions. Only 387 * the vertices with indices in @p considered_vertices are tested for 388 * equality. This speeds up the algorithm, which is, for worst-case hyper 389 * cube geometries $O(N^{3/2})$ in 2D and $O(N^{5/3})$ in 3D: quite slow. 390 * However, if you wish to consider all vertices, simply pass an empty 391 * vector. In that case, the function fills @p considered_vertices with all 392 * vertices. 393 * 394 * Two vertices are considered equal if their difference in each coordinate 395 * direction is less than @p tol. 396 */ 397 template <int dim, int spacedim> 398 void 399 delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices, 400 std::vector<CellData<dim>> & cells, 401 SubCellData & subcelldata, 402 std::vector<unsigned int> & considered_vertices, 403 const double tol = 1e-12); 404 405 /*@}*/ 406 /** 407 * @name Rotating, stretching and otherwise transforming meshes 408 */ 409 /*@{*/ 410 411 /** 412 * Transform the vertices of the given triangulation by applying the 413 * function object provided as first argument to all its vertices. 414 * 415 * The transformation given as argument is used to transform each vertex. 416 * Its respective type has to offer a function-like syntax, i.e. the 417 * predicate is either an object of a type that has an <tt>operator()</tt>, 418 * or it is a pointer to the function. In either case, argument and return 419 * value have to be of type <tt>Point@<spacedim@></tt>. 420 * 421 * @note The transformations that make sense to use with this function 422 * should have a Jacobian with a positive determinant. For example, 423 * rotation, shearing, stretching, or scaling all satisfy this (though 424 * there is no requirement that the transformation used actually is 425 * linear, as all of these examples are). On the other hand, reflections 426 * or inversions have a negative determinant of the Jacobian. The 427 * current function has no way of asserting a positive determinant 428 * of the Jacobian, but if you happen to use such a transformation, 429 * the result will be a triangulation in which cells have a negative 430 * volume. 431 * 432 * @note If you are using a parallel::distributed::Triangulation you will 433 * have hanging nodes in your local Triangulation even if your "global" mesh 434 * has no hanging nodes. This will cause issues with wrong positioning of 435 * hanging nodes in ghost cells if you call the current functions: The 436 * vertices of all locally owned cells will be correct, but the vertices of 437 * some ghost cells may not. This means that computations like 438 * KellyErrorEstimator may give wrong answers. 439 * 440 * @note This function is in general not compatible with manifolds attached 441 * to the triangulation. For example, in order to refine the grid (using 442 * manifolds) after the grid transformation, you have to make sure that 443 * the original manifold is still valid for the transformed geometry. This 444 * does not hold in general, and it is necessary to clear the manifold and 445 * attach a new one for the transformed geometry in these cases. 446 * If you want to perform refinements according to the original 447 * manifold description attached to the triangulation, you should first do 448 * the refinements, subsequently deactivate all manifolds, and finally call 449 * the transform() function. The result is a triangulation with correctly 450 * transformed vertices, but otherwise straight-sided elements. The 451 * following procedure is recommended 452 * @code 453 * ... 454 * triangulation.refine_global(n_refinements); 455 * triangulation.reset_all_manifolds(); 456 * Transformation<dim> transformation; 457 * GridTools::transform(transformation, triangulation); 458 * ... 459 * @endcode 460 * 461 * This function is used in the "Possibilities for extensions" section of 462 * step-38. It is also used in step-49 and step-53. 463 */ 464 template <int dim, typename Transformation, int spacedim> 465 void 466 transform(const Transformation & transformation, 467 Triangulation<dim, spacedim> &triangulation); 468 469 /** 470 * Shift each vertex of the triangulation by the given shift vector. This 471 * function uses the transform() function above, so the requirements on the 472 * triangulation stated there hold for this function as well. 473 */ 474 template <int dim, int spacedim> 475 void 476 shift(const Tensor<1, spacedim> & shift_vector, 477 Triangulation<dim, spacedim> &triangulation); 478 479 480 /** 481 * Rotate all vertices of the given two-dimensional triangulation in 482 * counter-clockwise sense around the origin of the coordinate system by the 483 * given angle (given in radians, rather than degrees). This function uses 484 * the transform() function above, so the requirements on the triangulation 485 * stated there hold for this function as well. 486 * 487 * @note This function is only supported for dim=2. 488 */ 489 template <int dim> 490 void 491 rotate(const double angle, Triangulation<dim> &triangulation); 492 493 /** 494 * Rotate all vertices of the given @p triangulation in counter-clockwise 495 * direction around the axis with the given index. Otherwise like the 496 * function above. 497 * 498 * @param[in] angle Angle in radians to rotate the Triangulation by. 499 * @param[in] axis Index of the coordinate axis to rotate around, keeping 500 * that coordinate fixed (0=x axis, 1=y axis, 2=z axis). 501 * @param[in,out] triangulation The Triangulation object to rotate. 502 * 503 * @note Implemented for dim=1, 2, and 3. 504 */ 505 template <int dim> 506 void 507 rotate(const double angle, 508 const unsigned int axis, 509 Triangulation<dim, 3> &triangulation); 510 511 /** 512 * Transform the given triangulation smoothly to a different domain where, 513 * typically, each of the vertices at the boundary of the triangulation is 514 * mapped to the corresponding points in the @p new_points map. 515 * 516 * The unknown displacement field $u_d(\mathbf x)$ in direction $d$ is 517 * obtained from the minimization problem \f[ \min\, \int \frac{1}{2} 518 * c(\mathbf x) 519 * \mathbf \nabla u_d(\mathbf x) \cdot 520 * \mathbf \nabla u_d(\mathbf x) 521 * \,\rm d x 522 * \f] 523 * subject to prescribed constraints. The minimizer is obtained by solving the 524 * Laplace equation of the dim components of a displacement field that maps 525 * the current 526 * domain into one described by @p new_points . Linear finite elements with 527 * four Gaussian quadrature points in each direction are used. The difference 528 * between the vertex positions specified in @p new_points and their current 529 * value in @p tria therefore represents the prescribed values of this 530 * displacement field at the boundary of the domain, or more precisely at all 531 * of those locations for which @p new_points provides values (which may be 532 * at part of the boundary, or even in the interior of the domain). The 533 * function then evaluates this displacement field at each unconstrained 534 * vertex and uses it to place the mapped vertex where the displacement 535 * field locates it. Because the solution of the Laplace equation is smooth, 536 * this guarantees a smooth mapping from the old domain to the new one. 537 * 538 * @param[in] new_points The locations where a subset of the existing 539 * vertices are to be placed. Typically, this would be a map from the vertex 540 * indices of all nodes on the boundary to their new locations, thus 541 * completely specifying the geometry of the mapped domain. However, it may 542 * also include interior points if necessary and it does not need to include 543 * all boundary vertices (although you then lose control over the exact 544 * shape of the mapped domain). 545 * 546 * @param[in,out] tria The Triangulation object. This object is changed in- 547 * place, i.e., the previous locations of vertices are overwritten. 548 * 549 * @param[in] coefficient An optional coefficient for the Laplace problem. 550 * Larger values make cells less prone to deformation (effectively 551 * increasing their stiffness). The coefficient is evaluated in the 552 * coordinate system of the old, undeformed configuration of the 553 * triangulation as input, i.e., before the transformation is applied. 554 * Should this function be provided, sensible results can only be expected 555 * if all coefficients are positive. 556 * 557 * @param[in] solve_for_absolute_positions If set to <code>true</code>, the 558 * minimization problem is formulated with respect to the final vertex 559 * positions as opposed to their displacement. The two formulations are 560 * equivalent for 561 * the homogeneous problem (default value of @p coefficient), but they 562 * result in very different mesh motion otherwise. Since in most cases one 563 * will be using a non-constant coefficient in displacement formulation, the 564 * default value of this parameter is <code>false</code>. 565 * 566 * @note This function is not currently implemented for the 1d case. 567 */ 568 template <int dim> 569 void 570 laplace_transform(const std::map<unsigned int, Point<dim>> &new_points, 571 Triangulation<dim> & tria, 572 const Function<dim, double> *coefficient = nullptr, 573 const bool solve_for_absolute_positions = false); 574 575 /** 576 * Return a std::map with all vertices of faces located in the boundary 577 * 578 * @param[in] tria The Triangulation object. 579 */ 580 template <int dim, int spacedim> 581 std::map<unsigned int, Point<spacedim>> 582 get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria); 583 584 /** 585 * Scale the entire triangulation by the given factor. To preserve the 586 * orientation of the triangulation, the factor must be positive. 587 * 588 * This function uses the transform() function above, so the requirements on 589 * the triangulation stated there hold for this function as well. 590 */ 591 template <int dim, int spacedim> 592 void 593 scale(const double scaling_factor, 594 Triangulation<dim, spacedim> &triangulation); 595 596 /** 597 * Distort the given triangulation by randomly moving around all the 598 * vertices of the grid. The direction of movement of each vertex is 599 * random, while the length of the shift vector has a value of @p factor 600 * times the minimal length of the active edges adjacent to this vertex. 601 * Note that @p factor should obviously be well below <tt>0.5</tt>. 602 * 603 * If @p keep_boundary is set to @p true (which is the default), then 604 * boundary vertices are not moved. 605 */ 606 template <int dim, int spacedim> 607 void 608 distort_random(const double factor, 609 Triangulation<dim, spacedim> &triangulation, 610 const bool keep_boundary = true); 611 612 /** 613 * Remove hanging nodes from a grid. If the @p isotropic parameter is set 614 * to @p false (default) this function detects cells with hanging nodes and 615 * refines the neighbours in the direction that removes hanging nodes. 616 * If the @p isotropic parameter is set 617 * to @p true, the neighbours refinement is made in each directions. 618 * In order to remove all hanging nodes this procedure has to be repeated: 619 * this could require a large number of iterations. 620 * To avoid this a max number (@p max_iterations) of iteration is provided. 621 * 622 * Consider the following grid: 623 * @image html remove_hanging_nodes-hanging.png 624 * 625 * @p isotropic == @p false would return: 626 * @image html remove_hanging_nodes-aniso.png 627 * 628 * @p isotropic == @p true would return: 629 * @image html remove_hanging_nodes-isotro.png 630 * 631 * @param[in,out] tria Triangulation to refine. 632 * 633 * @param[in] isotropic If true refine cells in each directions, otherwise 634 * (default value) refine the cell in the direction that removes hanging node. 635 * 636 * @param[in] max_iterations At each step only closest cells to hanging nodes 637 * are refined. The code may require a lot of iterations to remove all 638 * hanging nodes. @p max_iterations is the maximum number of iteration 639 * allowed. If @p max_iterations == numbers::invalid_unsigned_int this 640 * function continues refining until there are no hanging nodes. 641 * 642 * @note In the case of parallel codes, this function should be combined 643 * with GridGenerator::flatten_triangulation. 644 */ 645 template <int dim, int spacedim> 646 void 647 remove_hanging_nodes(Triangulation<dim, spacedim> &tria, 648 const bool isotropic = false, 649 const unsigned int max_iterations = 100); 650 651 /** 652 * Refine a mesh anisotropically such that the resulting mesh is composed by 653 * cells with maximum ratio between dimensions less than @p max_ratio. 654 * This procedure requires an algorithm that may not terminate. Consequently, 655 * it is possible to set a maximum number of iterations through the 656 * @p max_iterations parameter. 657 * 658 * Starting from a cell like this: 659 * @image html remove_anisotropy-coarse.png 660 * 661 * This function would return: 662 * @image html remove_anisotropy-refined.png 663 * 664 * @param[in,out] tria Triangulation to refine. 665 * 666 * @param[in] max_ratio Maximum value allowed among the ratio between 667 * the dimensions of each cell. 668 * 669 * @param[in] max_iterations Maximum number of iterations allowed. 670 * 671 * @note In the case of parallel codes, this function should be combined 672 * with GridGenerator::flatten_triangulation and 673 * GridTools::remove_hanging_nodes. 674 */ 675 template <int dim, int spacedim> 676 void 677 remove_anisotropy(Triangulation<dim, spacedim> &tria, 678 const double max_ratio = 1.6180339887, 679 const unsigned int max_iterations = 5); 680 681 /** 682 * Analyze the boundary cells of a mesh, and if one cell is found at 683 * a corner position (with dim adjacent faces on the boundary), and its 684 * dim-dimensional angle fraction exceeds @p limit_angle_fraction, 685 * refine globally once, and replace the children of such cell 686 * with children where the corner is no longer offending the given angle 687 * fraction. 688 * 689 * If no boundary cells exist with two adjacent faces on the boundary, then 690 * the triangulation is left untouched. If instead we do have cells with dim 691 * adjacent faces on the boundary, then the fraction between the 692 * dim-dimensional 693 * solid angle and dim*pi/2 is checked against the parameter @p limit_angle_fraction. 694 * If it is higher, the grid is refined once, and the children of the 695 * offending cell are replaced with some cells that instead respect the limit. 696 * After this process the triangulation is flattened, and all Manifold objects 697 * are restored as they were in the original triangulation. 698 * 699 * An example is given by the following mesh, obtained by attaching a 700 * SphericalManifold to a mesh generated using GridGenerator::hyper_cube: 701 * 702 * @code 703 * const SphericalManifold<dim> m0; 704 * Triangulation<dim> tria; 705 * GridGenerator::hyper_cube(tria,-1,1); 706 * tria.set_all_manifold_ids_on_boundary(0); 707 * tria.set_manifold(0, m0); 708 * tria.refine_global(4); 709 * @endcode 710 * 711 * <p ALIGN="center"> 712 * @image html regularize_mesh_01.png 713 * </p> 714 * 715 * The four cells that were originally the corners of a square will give you 716 * some troubles during computations, as the jacobian of the transformation 717 * from the reference cell to those cells will go to zero, affecting the error 718 * constants of the finite element estimates. 719 * 720 * Those cells have a corner with an angle that is very close to 180 degrees, 721 * i.e., an angle fraction very close to one. 722 * 723 * The same code, adding a call to regularize_corner_cells: 724 * @code 725 * const SphericalManifold<dim> m0; 726 * Triangulation<dim> tria; 727 * GridGenerator::hyper_cube(tria,-1,1); 728 * tria.set_all_manifold_ids_on_boundary(0); 729 * tria.set_manifold(0, m0); 730 * GridTools::regularize_corner_cells(tria); 731 * tria.refine_global(2); 732 * @endcode 733 * generates a mesh that has a much better behaviour w.r.t. the jacobian of 734 * the Mapping: 735 * 736 * <p ALIGN="center"> 737 * @image html regularize_mesh_02.png 738 * </p> 739 * 740 * This mesh is very similar to the one obtained by GridGenerator::hyper_ball. 741 * However, using GridTools::regularize_corner_cells one has the freedom to 742 * choose when to apply the regularization, i.e., one could in principle first 743 * refine a few times, and then call the regularize_corner_cells function: 744 * 745 * @code 746 * const SphericalManifold<dim> m0; 747 * Triangulation<dim> tria; 748 * GridGenerator::hyper_cube(tria,-1,1); 749 * tria.set_all_manifold_ids_on_boundary(0); 750 * tria.set_manifold(0, m0); 751 * tria.refine_global(2); 752 * GridTools::regularize_corner_cells(tria); 753 * tria.refine_global(1); 754 * @endcode 755 * 756 * This generates the following mesh: 757 * 758 * <p ALIGN="center"> 759 * @image html regularize_mesh_03.png 760 * </p> 761 * 762 * The function is currently implemented only for dim = 2 and 763 * will throw an exception if called with dim = 3. 764 * 765 * @param[in,out] tria Triangulation to regularize. 766 * 767 * @param[in] limit_angle_fraction Maximum ratio of angle or solid 768 * angle that is allowed for a corner element in the mesh. 769 */ 770 template <int dim, int spacedim> 771 void 772 regularize_corner_cells(Triangulation<dim, spacedim> &tria, 773 const double limit_angle_fraction = .75); 774 775 /*@}*/ 776 /** 777 * @name Finding cells and vertices of a triangulation 778 */ 779 /*@{*/ 780 781 /** 782 * Given a Triangulation's @p cache and a list of @p points create the quadrature rules. 783 * 784 * @param[in] cache The triangulation's GridTools::Cache . 785 * @param[in] points The point's vector. 786 * @param[in] cell_hint (optional) A cell iterator for a cell which likely 787 * contains the first point of @p points. 788 * 789 * @return A tuple containing the following information: 790 * - @p cells : A vector of all the cells containing at least one of 791 * the @p points . 792 * - @p qpoints : A vector of vectors of points. @p qpoints[i] contains 793 * the reference positions of all points that fall within the cell @p cells[i] . 794 * - @p indices : A vector of vectors of integers, containing the mapping between 795 * local numbering in @p qpoints , and global index in @p points . 796 * 797 * If @p points[a] and @p points[b] are the only two points that fall in @p cells[c], 798 * then @p qpoints[c][0] and @p qpoints[c][1] are the reference positions of 799 * @p points[a] and @p points[b] in @p cells[c], and @p indices[c][0] = a, 800 * @p indices[c][1] = b. The function 801 * Mapping::transform_unit_to_real(qpoints[c][0]) 802 * returns @p points[a]. 803 * 804 * The algorithm assumes it's easier to look for a point in the cell that was 805 * used previously. For this reason random points are, computationally 806 * speaking, the worst case scenario while points grouped by the cell to which 807 * they belong are the best case. Pre-sorting points, trying to minimize 808 * distances between them, might make the function extremely faster. 809 * 810 * @note If a point is not found inside the mesh, or is lying inside an 811 * artificial cell of a parallel::TriangulationBase, an exception is thrown. 812 * 813 * @note The actual return type of this function, i.e., the type referenced 814 * above as @p return_type, is 815 * @code 816 * std::tuple< 817 * std::vector< 818 * typename Triangulation<dim, spacedim>::active_cell_iterator>, 819 * std::vector<std::vector<Point<dim>>>, 820 * std::vector<std::vector<unsigned int>>> 821 * @endcode 822 * The type is abbreviated in the online documentation to improve readability 823 * of this page. 824 * 825 * @note This function optimizes the search by making use of 826 * GridTools::Cache::get_cell_bounding_boxes_rtree(), which either returns 827 * a cached rtree or builds and stores one. Building an rtree might hinder 828 * the performance if the function is called only once on few points. 829 */ 830 template <int dim, int spacedim> 831 # ifndef DOXYGEN 832 std::tuple< 833 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>, 834 std::vector<std::vector<Point<dim>>>, 835 std::vector<std::vector<unsigned int>>> 836 # else 837 return_type 838 # endif 839 compute_point_locations( 840 const Cache<dim, spacedim> & cache, 841 const std::vector<Point<spacedim>> &points, 842 const typename Triangulation<dim, spacedim>::active_cell_iterator 843 &cell_hint = 844 typename Triangulation<dim, spacedim>::active_cell_iterator()); 845 846 /** 847 * This function is similar to GridTools::compute_point_locations(), 848 * but it tries to find and transform every point of @p points. 849 * 850 * @return A tuple containing four elements; the first three 851 * are documented in GridTools::compute_point_locations(). 852 * The last element of the @p return_type contains the 853 * indices of points which are neither found inside the mesh 854 * nor lie in artificial cells. The @p return_type equals the 855 * following tuple type: 856 * @code 857 * std::tuple< 858 * std::vector< 859 * typename Triangulation<dim,spacedim>::active_cell_iterator>, 860 * std::vector<std::vector<Point<dim>>>, 861 * std::vector<std::vector<unsigned int>>, 862 * std::vector<unsigned int> 863 * > 864 * @endcode 865 * 866 * @note This function optimizes the search by making use of 867 * GridTools::Cache::get_cell_bounding_boxes_rtree(), which either returns 868 * a cached rtree or builds and stores one. Building an rtree might hinder 869 * the performance if the function is called only once on few points. 870 * 871 * For a more detailed documentation see 872 * GridTools::compute_point_locations(). 873 */ 874 template <int dim, int spacedim> 875 # ifndef DOXYGEN 876 std::tuple< 877 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>, 878 std::vector<std::vector<Point<dim>>>, 879 std::vector<std::vector<unsigned int>>, 880 std::vector<unsigned int>> 881 # else 882 return_type 883 # endif 884 compute_point_locations_try_all( 885 const Cache<dim, spacedim> & cache, 886 const std::vector<Point<spacedim>> &points, 887 const typename Triangulation<dim, spacedim>::active_cell_iterator 888 &cell_hint = 889 typename Triangulation<dim, spacedim>::active_cell_iterator()); 890 891 /** 892 * Given a @p cache and a list of 893 * @p local_points for each process, find the points lying on the locally 894 * owned part of the mesh and compute the quadrature rules for them. 895 * Distributed compute point locations is a function similar to 896 * GridTools::compute_point_locations but working for 897 * parallel::TriangulationBase objects and, unlike its serial version, also 898 * for a distributed triangulation (see parallel::distributed::Triangulation). 899 * 900 * @param[in] cache a GridTools::Cache object 901 * @param[in] local_points the array of points owned by the current process. 902 * Every process can have a different array of points which can be empty and 903 * not contained within the locally owned part of the triangulation 904 * @param[in] global_bboxes a vector of vectors of bounding boxes; it 905 * describes the locally owned part of the mesh for each process. The bounding 906 * boxes describing which part of the mesh is locally owned by process with 907 * rank rk are contained in global_bboxes[rk]. The local description can be 908 * obtained from GridTools::compute_mesh_predicate_bounding_box; then the 909 * global one can be obtained using either 910 * GridTools::exchange_local_bounding_boxes or Utilities::MPI::all_gather 911 * @return A tuple containing the quadrature information 912 * 913 * The elements of the output tuple are: 914 * - cells : a vector of all cells containing at least one point. 915 * - qpoints : a vector of vector of points; containing in @p qpoints[i] 916 * the reference positions of all points that fall within the cell @p cells[i] . 917 * - maps : a vector of vector of integers, containing the mapping between 918 * the numbering in qpoints (previous element of the tuple), and the vector 919 * of local points of the process owning the points. 920 * - points : a vector of vector of points. @p points[i][j] is the point in the 921 * real space corresponding. 922 * to @p qpoints[i][j] . Notice @p points are the points lying on the locally 923 * owned part of the mesh; thus these can be either copies of @p local_points 924 * or points received from other processes i.e. local_points for other 925 * processes 926 * - owners : a vector of vectors; @p owners[i][j] contains the rank of 927 * the process owning the point[i][j] (previous element of the tuple). 928 * 929 * The function uses the triangulation's mpi communicator: for this reason it 930 * throws an assert error if the Triangulation is not derived from 931 * parallel::TriangulationBase . 932 * 933 * In a serial execution the first three elements of the tuple are the same 934 * as in GridTools::compute_point_locations . 935 * 936 * Note: this function is a collective operation. 937 * 938 * @note The actual return type of this function, i.e., the type referenced 939 * above as @p return_type, is 940 * @code 941 * std::tuple< 942 * std::vector< 943 * typename Triangulation<dim, spacedim>::active_cell_iterator>, 944 * std::vector<std::vector<Point<dim>>>, 945 * std::vector<std::vector<unsigned int>>, 946 * std::vector<std::vector<Point<spacedim>>>, 947 * std::vector<std::vector<unsigned int>>> 948 * @endcode 949 * The type is abbreviated in the online documentation to improve readability 950 * of this page. 951 */ 952 template <int dim, int spacedim> 953 # ifndef DOXYGEN 954 std::tuple< 955 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>, 956 std::vector<std::vector<Point<dim>>>, 957 std::vector<std::vector<unsigned int>>, 958 std::vector<std::vector<Point<spacedim>>>, 959 std::vector<std::vector<unsigned int>>> 960 # else 961 return_type 962 # endif 963 distributed_compute_point_locations( 964 const GridTools::Cache<dim, spacedim> & cache, 965 const std::vector<Point<spacedim>> & local_points, 966 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes); 967 968 /** 969 * Return a map `vertex index -> Point<spacedim>` containing the used 970 * vertices of the given `container`. The key of the returned map (i.e., 971 * the first element of the pair above) is the global index in the 972 * triangulation, whereas the value of each pair is the physical 973 * location of the corresponding vertex. The used vertices are obtained by 974 * looping over all cells, 975 * and querying for each cell where its vertices are through the (optional) 976 * `mapping` argument. 977 * 978 * In serial Triangulation objects and parallel::shared::Triangulation 979 * objects, the size of the returned map 980 * equals Triangulation::n_used_vertices() (not Triangulation::n_vertices()). 981 * Note that in parallel::distributed::Triangulation objects, only vertices in 982 * locally owned cells and ghost cells are returned, as for all other vertices 983 * their real location might not be known (e.g. for distributed computations 984 * using MappingQEulerian). 985 * 986 * If you use the default `mapping`, the returned map satisfies the following 987 * equality: 988 * 989 * @code 990 * const auto used_vertices = extract_used_vertices(tria); 991 * auto all_vertices = tria.get_vertices(); 992 * 993 * for(const auto &id_and_v : used_vertices) 994 * all_vertices[id_and_v.first] == id_and_v.second; // true 995 * @endcode 996 * 997 * Notice that the above is not satisfied for mappings that change the 998 * location of vertices, like MappingQEulerian. 999 * 1000 * @ref ConceptMeshType "MeshType concept". 1001 * @param container The container to extract vertices from. 1002 * @param mapping The mapping to use to compute the points locations. 1003 */ 1004 template <int dim, int spacedim> 1005 std::map<unsigned int, Point<spacedim>> 1006 extract_used_vertices(const Triangulation<dim, spacedim> &container, 1007 const Mapping<dim, spacedim> & mapping = 1008 StaticMappingQ1<dim, spacedim>::mapping); 1009 1010 /** 1011 * Find and return the index of the closest vertex to a given point in the 1012 * map of vertices passed as the first argument. 1013 * 1014 * @param vertices A map of index->vertex, as returned by 1015 * GridTools::extract_used_vertices(). 1016 * @param p The target point. 1017 * @return The index of the vertex that is closest to the target point `p`. 1018 */ 1019 template <int spacedim> 1020 unsigned int 1021 find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices, 1022 const Point<spacedim> & p); 1023 1024 /** 1025 * Find and return the index of the used vertex (or marked vertex) in a 1026 * given mesh that is located closest to a given point. 1027 * 1028 * This function uses the locations of vertices as stored in the 1029 * triangulation. This is usually sufficient, unless you are using a Mapping 1030 * that moves the vertices around (for example, MappingQEulerian). In this 1031 * case, you should call the function with the same name and with an 1032 * additional Mapping argument. 1033 * 1034 * @param mesh A variable of a type that satisfies the requirements of the 1035 * @ref ConceptMeshType "MeshType concept". 1036 * @param p The point for which we want to find the closest vertex. 1037 * @param marked_vertices An array of bools indicating which 1038 * vertices of @p mesh will be considered within the search 1039 * as the potentially closest vertex. On receiving a non-empty 1040 * @p marked_vertices, the function will 1041 * only search among @p marked_vertices for the closest vertex. 1042 * The size of this array should be equal to the value returned by 1043 * Triangulation::n_vertices() for the triangulation underlying the given mesh 1044 * (as opposed to the value returned by Triangulation::n_used_vertices()). 1045 * @return The index of the closest vertex found. 1046 */ 1047 template <int dim, template <int, int> class MeshType, int spacedim> 1048 unsigned int 1049 find_closest_vertex(const MeshType<dim, spacedim> &mesh, 1050 const Point<spacedim> & p, 1051 const std::vector<bool> & marked_vertices = {}); 1052 1053 /** 1054 * Find and return the index of the used vertex (or marked vertex) in a 1055 * given mesh that is located closest to a given point. Use the given 1056 * mapping to compute the actual location of the vertices. 1057 * 1058 * If the Mapping does not modify the position of the mesh vertices (like, 1059 * for example, MappingQEulerian does), then this function is equivalent to 1060 * the one with the same name, and without the `mapping` argument. 1061 * 1062 * @param mapping A mapping used to compute the vertex locations 1063 * @param mesh A variable of a type that satisfies the requirements of the 1064 * @ref ConceptMeshType "MeshType concept". 1065 * @param p The point for which we want to find the closest vertex. 1066 * @param marked_vertices An array of bools indicating which 1067 * vertices of @p mesh will be considered within the search 1068 * as the potentially closest vertex. On receiving a non-empty 1069 * @p marked_vertices, the function will 1070 * only search among @p marked_vertices for the closest vertex. 1071 * The size of this array should be equal to the value returned by 1072 * Triangulation::n_vertices() for the triangulation underlying the given mesh 1073 * (as opposed to the value returned by Triangulation::n_used_vertices()). 1074 * @return The index of the closest vertex found. 1075 */ 1076 template <int dim, template <int, int> class MeshType, int spacedim> 1077 unsigned int 1078 find_closest_vertex(const Mapping<dim, spacedim> & mapping, 1079 const MeshType<dim, spacedim> &mesh, 1080 const Point<spacedim> & p, 1081 const std::vector<bool> & marked_vertices = {}); 1082 1083 1084 /** 1085 * Find and return a vector of iterators to active cells that surround a 1086 * given vertex with index @p vertex_index. 1087 * 1088 * For locally refined grids, the vertex itself might not be a vertex of all 1089 * adjacent cells that are returned. However, it will always be either a 1090 * vertex of a cell or be a hanging node located on a face or an edge of it. 1091 * 1092 * @param container A variable of a type that satisfies the requirements of 1093 * the 1094 * @ref ConceptMeshType "MeshType concept". 1095 * @param vertex_index The index of the vertex for which we try to find 1096 * adjacent cells. 1097 * @return A vector of cells that lie adjacent to the given vertex. 1098 * 1099 * @note If the point requested does not lie in any of the cells of the mesh 1100 * given, then this function throws an exception of type 1101 * GridTools::ExcPointNotFound. You can catch this exception and decide what 1102 * to do in that case. 1103 * 1104 * @note It isn't entirely clear at this time whether the function does the 1105 * right thing with anisotropically refined meshes. It needs to be checked 1106 * for this case. 1107 */ 1108 template <int dim, template <int, int> class MeshType, int spacedim> 1109 # ifndef _MSC_VER 1110 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator> 1111 # else 1112 std::vector< 1113 typename dealii::internal:: 1114 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type> 1115 # endif 1116 find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container, 1117 const unsigned int vertex_index); 1118 1119 /** 1120 * Find an active cell that surrounds a given point @p p. The return type 1121 * is a pair of an iterator to the active cell along with the unit cell 1122 * coordinates of the point. 1123 * 1124 * The algorithm used in this function proceeds by first looking for the 1125 * vertex located closest to the given point, see 1126 * GridTools::find_closest_vertex(). Secondly, all adjacent cells to this 1127 * vertex are found in the mesh, see 1128 * GridTools::find_cells_adjacent_to_vertex(). Lastly, for each of these 1129 * cells, the function tests whether the point is inside. This check is 1130 * performed using the given @p mapping argument to determine whether cells 1131 * have straight or curved boundaries. 1132 * 1133 * If a point lies on the boundary of two or more cells, then the algorithm 1134 * tries to identify the cell that is of highest refinement level. 1135 * 1136 * If the point requested does not lie in a locally-owned or ghost cell, 1137 * then this function throws an exception of type GridTools::ExcPointNotFound. 1138 * You can catch this exception and decide what to do in that case. Hence, 1139 * for programs that work with partitioned (parallel) triangulations, this 1140 * function should always be called inside a `try`-block unless it is a 1141 * priori clear that the point with which it is called must be inside 1142 * a locally owned or ghost cell (and not close enough to the boundary 1143 * between ghost and artificial cells so that decision which cell it is 1144 * on depends on floating point accuracy). 1145 * 1146 * @param mapping The mapping used to determine whether the given point is 1147 * inside a given cell. 1148 * @param mesh A variable of a type that satisfies the requirements of the 1149 * @ref ConceptMeshType "MeshType concept". 1150 * @param p The point for which we want to find the surrounding cell. 1151 * @param marked_vertices An array of `bool`s indicating whether an 1152 * entry in the vertex array should be considered 1153 * (and the others must be ignored) as the potentially 1154 * closest vertex to the specified point. On specifying a non-default 1155 * @p marked_vertices, find_closest_vertex() would 1156 * only search among @p marked_vertices for the closest vertex. 1157 * The size of this array should be equal to n_vertices() of the 1158 * triangulation (as opposed to n_used_vertices() ). The motivation of using 1159 * @p marked_vertices is to cut down the search space of vertices if one has 1160 * a priori knowledge of a collection of vertices that the point of interest 1161 * may be close to. 1162 * @param tolerance Tolerance in terms of unit cell coordinates. Depending 1163 * on the problem, it might be necessary to adjust the tolerance in order 1164 * to be able to identify a cell. Floating 1165 * point arithmetic implies that a point will, in general, not lie exactly 1166 * on a vertex, edge, or face. In either case, it is not predictable which 1167 * of the cells adjacent to a vertex or an edge/face this function returns. 1168 * Consequently, algorithms that call this function need to take into 1169 * account that the returned cell will only contain the point approximately. 1170 * 1171 * @return A pair of an iterators into the mesh that points to the 1172 * surrounding cell, and of the unit cell coordinates of that point. This 1173 * local position might be located slightly outside an actual unit cell, 1174 * due to numerical roundoff. Therefore, the point returned by this function 1175 * should be projected onto the unit cell, using 1176 * GeometryInfo::project_to_unit_cell(). This is not automatically performed 1177 * by the algorithm. The returned cell can be a locally-owned cell or a 1178 * ghost cell (but not an artificial cell). The returned cell might be a 1179 * ghost cell even if the given point is a vertex of a locally owned cell. 1180 * The reason behind is that this is the only way to guarantee that all 1181 * processors that participate in a parallel triangulation will agree which 1182 * cell contains a point. For example, if two processors come together 1183 * at one vertex and the function is called with this vertex, then one 1184 * processor will return a locally owned cell and the other one a ghost cell. 1185 */ 1186 template <int dim, template <int, int> class MeshType, int spacedim> 1187 # ifndef _MSC_VER 1188 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>> 1189 # else 1190 std::pair<typename dealii::internal:: 1191 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 1192 Point<dim>> 1193 # endif 1194 find_active_cell_around_point(const Mapping<dim, spacedim> & mapping, 1195 const MeshType<dim, spacedim> &mesh, 1196 const Point<spacedim> & p, 1197 const std::vector<bool> &marked_vertices = {}, 1198 const double tolerance = 1.e-10); 1199 1200 /** 1201 * A version of the above function that assumes straight boundaries and 1202 * as a consequence simply calls the above function using MappingQ1 for 1203 * the mapping argument. 1204 * 1205 * @return An iterator into the mesh that points to the surrounding cell. 1206 */ 1207 template <int dim, template <int, int> class MeshType, int spacedim> 1208 # ifndef _MSC_VER 1209 typename MeshType<dim, spacedim>::active_cell_iterator 1210 # else 1211 typename dealii::internal:: 1212 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type 1213 # endif 1214 find_active_cell_around_point(const MeshType<dim, spacedim> &mesh, 1215 const Point<spacedim> & p, 1216 const std::vector<bool> &marked_vertices = {}, 1217 const double tolerance = 1.e-10); 1218 1219 /** 1220 * Another version where we use that mapping on a given 1221 * cell that corresponds to the active finite element index of that cell. 1222 * This is obviously only useful for hp problems, since the active finite 1223 * element index for all other DoF handlers is always zero. 1224 */ 1225 template <int dim, int spacedim> 1226 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator, 1227 Point<dim>> 1228 find_active_cell_around_point( 1229 const hp::MappingCollection<dim, spacedim> &mapping, 1230 const DoFHandler<dim, spacedim> & mesh, 1231 const Point<spacedim> & p, 1232 const double tolerance = 1.e-10); 1233 1234 /** 1235 * Finding an active cell around a point can be very expensive in terms 1236 * of computational costs. This function aims at providing a fast version 1237 * of the above functions by using a space-tree to speed up the geometry 1238 * search. 1239 * 1240 * @param cache Object with information about the space-tree of a triangulation, 1241 * see GridTools::Cache. 1242 * @param p The point for which we want to find the surrounding cell. 1243 * @param cell_hint Gives a hint for the geometry search, which is beneficial 1244 * if a-priori knowledge is available regarding the cell on which the point 1245 * may likely be located. A typical use case would be that this search has 1246 * to be done for an array of points that are close to each other and where 1247 * the adjacent cell of the previous point is a good hint for the next point 1248 * in the array. 1249 * @param marked_vertices See above. 1250 * @param tolerance See above. 1251 * 1252 * 1253 * The following code example shows how to use this function: 1254 * 1255 * @code 1256 * GridTools::Cache<dim, dim> cache(triangulation, mapping); 1257 * auto cell_hint = typename Triangulation<dim, dim>::active_cell_iterator(); 1258 * std::vector<bool> marked_vertices = {}; 1259 * double tolerance = 1.e-10; 1260 * 1261 * std::vector<Point<dim>> points; // a vector of many points 1262 * ... 1263 * 1264 * for(auto p : points) 1265 * { 1266 * try 1267 * { 1268 * auto cell_and_ref_point = GridTools::find_active_cell_around_point( 1269 * cache, p, cell_hint, marked_vertices, tolerance); 1270 * 1271 * // use current cell as hint for the next point 1272 * cell_hint = cell_and_ref_point.first; 1273 * } 1274 * catch(...) 1275 * { 1276 * } 1277 * 1278 * ... 1279 * } 1280 * @endcode 1281 */ 1282 template <int dim, int spacedim> 1283 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator, 1284 Point<dim>> 1285 find_active_cell_around_point( 1286 const Cache<dim, spacedim> &cache, 1287 const Point<spacedim> & p, 1288 const typename Triangulation<dim, spacedim>::active_cell_iterator & 1289 cell_hint = typename Triangulation<dim, spacedim>::active_cell_iterator(), 1290 const std::vector<bool> &marked_vertices = {}, 1291 const double tolerance = 1.e-10); 1292 1293 /** 1294 * A version of the previous function that exploits an already existing 1295 * map between vertices and cells (constructed using the function 1296 * GridTools::vertex_to_cell_map()), a map of vertex_to_cell_centers (obtained 1297 * through GridTools::vertex_to_cell_centers_directions()), and 1298 * optionally an RTree constructed from the used vertices of the 1299 * Triangulation. 1300 * 1301 * @note All of these structures can be queried from a 1302 * GridTools::Cache object. Note, however, that in this case MeshType 1303 * has to be Triangulation, so that it might be more appropriate to directly 1304 * call the function above with argument `cache` in this case. 1305 */ 1306 template <int dim, template <int, int> class MeshType, int spacedim> 1307 # ifndef _MSC_VER 1308 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>> 1309 # else 1310 std::pair<typename dealii::internal:: 1311 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 1312 Point<dim>> 1313 # endif 1314 find_active_cell_around_point( 1315 const Mapping<dim, spacedim> & mapping, 1316 const MeshType<dim, spacedim> &mesh, 1317 const Point<spacedim> & p, 1318 const std::vector< 1319 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>> 1320 & vertex_to_cell_map, 1321 const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers, 1322 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint = 1323 typename MeshType<dim, spacedim>::active_cell_iterator(), 1324 const std::vector<bool> & marked_vertices = {}, 1325 const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree = 1326 RTree<std::pair<Point<spacedim>, unsigned int>>{}, 1327 const double tolerance = 1.e-10); 1328 1329 /** 1330 * As compared to the functions above, this function identifies all cells 1331 * around a point for a given tolerance level `tolerance` in terms of unit 1332 * coordinates. Given a first cell with reference coordinates as parameter 1333 * @p first_cell, e.g. obtained by one of the functions above, all 1334 * corresponding neighboring cells with points in unit coordinates are also 1335 * identified. 1336 * 1337 * This function is useful e.g. for discontinuous function spaces where, for 1338 * the case the given point `p` lies on a vertex, edge or face, several 1339 * cells might hold independent values of the solution that get combined in 1340 * some way in a user code. 1341 * 1342 * This function is used as follows 1343 * @code 1344 * auto first_cell = GridTools::find_active_cell_around_point(...); 1345 * auto all_cells = GridTools::find_all_active_cells_around_point( 1346 * mapping, mesh, p, tolerance, first_cell); 1347 * @endcode 1348 */ 1349 template <int dim, template <int, int> class MeshType, int spacedim> 1350 # ifndef _MSC_VER 1351 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 1352 Point<dim>>> 1353 # else 1354 std::vector<std::pair< 1355 typename dealii::internal:: 1356 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 1357 Point<dim>>> 1358 # endif 1359 find_all_active_cells_around_point( 1360 const Mapping<dim, spacedim> & mapping, 1361 const MeshType<dim, spacedim> &mesh, 1362 const Point<spacedim> & p, 1363 const double tolerance, 1364 const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 1365 Point<dim>> & first_cell); 1366 1367 /** 1368 * A variant of the previous function that internally calls one of the 1369 * functions find_active_cell_around_point() to obtain a first cell, and 1370 * subsequently adds all other cells by calling the function 1371 * find_all_active_cells_around_point() above. 1372 */ 1373 template <int dim, template <int, int> class MeshType, int spacedim> 1374 # ifndef _MSC_VER 1375 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, 1376 Point<dim>>> 1377 # else 1378 std::vector<std::pair< 1379 typename dealii::internal:: 1380 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type, 1381 Point<dim>>> 1382 # endif 1383 find_all_active_cells_around_point( 1384 const Mapping<dim, spacedim> & mapping, 1385 const MeshType<dim, spacedim> &mesh, 1386 const Point<spacedim> & p, 1387 const double tolerance = 1e-10, 1388 const std::vector<bool> & marked_vertices = {}); 1389 1390 /** 1391 * Return a list of all descendants of the given cell that are active. For 1392 * example, if the current cell is once refined but none of its children are 1393 * any further refined, then the returned list will contain all its 1394 * children. 1395 * 1396 * If the current cell is already active, then the returned list is empty 1397 * (because the cell has no children that may be active). 1398 * 1399 * @tparam MeshType A type that satisfies the requirements of the 1400 * @ref ConceptMeshType "MeshType concept". 1401 * @param cell An iterator pointing to a cell of the mesh. 1402 * @return A list of active descendants of the given cell 1403 * 1404 * @note Since in C++ the MeshType template argument can not be deduced from 1405 * a function call, you will have to specify it after the function name, as 1406 * for example in 1407 * @code 1408 * GridTools::get_active_child_cells<DoFHandler<dim> > (cell) 1409 * @endcode 1410 */ 1411 template <class MeshType> 1412 std::vector<typename MeshType::active_cell_iterator> 1413 get_active_child_cells(const typename MeshType::cell_iterator &cell); 1414 1415 /** 1416 * Extract the active cells around a given cell @p cell and return them in 1417 * the vector @p active_neighbors. These neighbors are specifically the 1418 * <i>face</i> neighbors of a cell or, if that neighbor is further 1419 * refined, its active children that border on that face. On the other 1420 * hand, the neighbors returned do not include cells that lie, for 1421 * example, diagonally opposite to a vertex but are not face neighbors 1422 * themselves. (In 3d, it also does not include cells that are 1423 * adjacent to one of the edges of the current cell, but are not 1424 * face neighbors.) 1425 * 1426 * @tparam MeshType A type that satisfies the requirements of the 1427 * @ref ConceptMeshType "MeshType concept". 1428 * @param[in] cell An iterator pointing to a cell of the mesh. 1429 * @param[out] active_neighbors A list of active descendants of the given 1430 * cell 1431 * 1432 * @note Since in C++ the MeshType template argument can not be deduced from 1433 * a function call, you will have to specify it after the function name, as 1434 * for example in 1435 * @code 1436 * GridTools::get_active_neighbors<DoFHandler<dim>>(cell, active_neighbors) 1437 * @endcode 1438 */ 1439 template <class MeshType> 1440 void 1441 get_active_neighbors( 1442 const typename MeshType::active_cell_iterator & cell, 1443 std::vector<typename MeshType::active_cell_iterator> &active_neighbors); 1444 1445 /** 1446 * Extract and return the active cell layer around a subdomain (set of 1447 * active cells) in the @p mesh (i.e. those that share a common set of 1448 * vertices with the subdomain but are not a part of it). Here, the 1449 * "subdomain" consists of exactly all of those cells for which the @p 1450 * predicate returns @p true. 1451 * 1452 * An example of a custom predicate is one that checks for a given material 1453 * id 1454 * @code 1455 * template <int dim> 1456 * bool 1457 * pred_mat_id(const typename Triangulation<dim>::active_cell_iterator & cell) 1458 * { 1459 * return cell->material_id() == 1; 1460 * } 1461 * @endcode 1462 * and we can then extract the layer of cells around this material with the 1463 * following call: 1464 * @code 1465 * GridTools::compute_active_cell_halo_layer(tria, pred_mat_id<dim>); 1466 * @endcode 1467 * 1468 * Predicates that are frequently useful can be found in namespace 1469 * IteratorFilters. For example, it is possible to extract a layer 1470 * of cells around all of those cells with a given material id, 1471 * @code 1472 * GridTools::compute_active_cell_halo_layer( 1473 * tria, IteratorFilters::MaterialIdEqualTo(1, true)); 1474 * @endcode 1475 * or around all cells with one of a set of active FE indices for an 1476 * hp::DoFHandler 1477 * @code 1478 * GridTools::compute_active_cell_halo_layer( 1479 * hp_dof_handler, IteratorFilters::ActiveFEIndexEqualTo({1,2}, true)); 1480 * @endcode 1481 * Note that in the last two examples we ensure that the predicate returns 1482 * true only for locally owned cells. This means that the halo layer will 1483 * not contain any artificial cells. 1484 * 1485 * @tparam MeshType A type that satisfies the requirements of the 1486 * @ref ConceptMeshType "MeshType concept". 1487 * @param[in] mesh A mesh (i.e. objects of type Triangulation, DoFHandler, 1488 * or hp::DoFHandler). 1489 * @param[in] predicate A function (or object of a type with an operator()) 1490 * defining the subdomain around which the halo layer is to be extracted. It 1491 * is a function that takes in an active cell and returns a boolean. 1492 * @return A list of active cells sharing at least one common vertex with 1493 * the predicated subdomain. 1494 */ 1495 template <class MeshType> 1496 std::vector<typename MeshType::active_cell_iterator> 1497 compute_active_cell_halo_layer( 1498 const MeshType &mesh, 1499 const std::function<bool(const typename MeshType::active_cell_iterator &)> 1500 &predicate); 1501 1502 1503 /** 1504 * Extract and return the cell layer around a subdomain (set of 1505 * cells) on a specified level of the @p mesh (i.e. those cells on 1506 * that level that share a common set of vertices with the subdomain 1507 * but are not a part of it). Here, the "subdomain" consists of exactly 1508 * all of those cells for which the @p predicate returns @p true. 1509 */ 1510 template <class MeshType> 1511 std::vector<typename MeshType::cell_iterator> 1512 compute_cell_halo_layer_on_level( 1513 const MeshType &mesh, 1514 const std::function<bool(const typename MeshType::cell_iterator &)> 1515 & predicate, 1516 const unsigned int level); 1517 1518 1519 /** 1520 * Extract and return ghost cells which are the active cell layer around all 1521 * locally owned cells. This is most relevant for 1522 * parallel::shared::Triangulation where it will return a subset of all 1523 * ghost cells on a processor, but for parallel::distributed::Triangulation 1524 * this will return all the ghost cells. 1525 * 1526 * @tparam MeshType A type that satisfies the requirements of the 1527 * @ref ConceptMeshType "MeshType concept". 1528 * @param[in] mesh A mesh (i.e. objects of type Triangulation, DoFHandler, 1529 * or hp::DoFHandler). 1530 * @return A list of ghost cells 1531 */ 1532 template <class MeshType> 1533 std::vector<typename MeshType::active_cell_iterator> 1534 compute_ghost_cell_halo_layer(const MeshType &mesh); 1535 1536 /** 1537 * Extract and return the set of active cells within a geometric distance of 1538 * @p layer_thickness around a subdomain (set of active cells) in the @p mesh. 1539 * Here, the "subdomain" consists of exactly all of 1540 * those cells for which the @p predicate returns @p true. 1541 * 1542 * The function first computes the cells that form the 'surface' of the 1543 * subdomain that consists of all of the active cells for which the predicate 1544 * is true. Using compute_bounding_box(), a bounding box is 1545 * computed for this subdomain and extended by @p layer_thickness. These 1546 * cells are called interior subdomain boundary cells. 1547 * The active cells with all of their vertices outside the extended 1548 * bounding box are ignored. 1549 * The cells that are inside the extended bounding box are then checked for 1550 * their proximity to the interior subdomain boundary cells. This implies 1551 * checking the distance between a pair of arbitrarily oriented cells, 1552 * which is not trivial in general. To simplify this, the algorithm checks 1553 * the distance between the two enclosing spheres of the cells. 1554 * This will definitely result in slightly more cells being marked but 1555 * also greatly simplifies the arithmetic complexity of the algorithm. 1556 * 1557 * @image html active_cell_layer_within_distance.png 1558 * The image shows a mesh generated by subdivided_hyper_rectangle(). The cells 1559 * are marked using three different colors. If the grey colored cells in the 1560 * image are the cells for which the predicate is true, then the function 1561 * compute_active_cell_layer_within_distance() will return a set of cell 1562 * iterators corresponding to the cells colored in red. 1563 * The red colored cells are the active cells that are within a given 1564 * distance to the grey colored cells. 1565 * 1566 * @tparam MeshType A type that satisfies the requirements of the 1567 * @ref ConceptMeshType "MeshType concept". 1568 * @param mesh A mesh (i.e. objects of type Triangulation, DoFHandler, 1569 * or hp::DoFHandler). 1570 * @param predicate A function (or object of a type with an operator()) 1571 * defining the subdomain around which the halo layer is to be extracted. It 1572 * is a function that takes in an active cell and returns a boolean. 1573 * @param layer_thickness specifies the geometric distance within 1574 * which the function searches for active cells from the predicate domain. 1575 * If the minimal distance between the enclosing sphere of the an 1576 * active cell and the enclosing sphere of any of the cells for which 1577 * the @p predicate returns @p true is less than @p layer_thickness, 1578 * then the active cell is an \a active_cell_within_distance. 1579 * @return A list of active cells within a given geometric distance 1580 * @p layer_thickness from the set of active cells for which the @p predicate 1581 * returns @p true. 1582 * 1583 * See compute_active_cell_halo_layer(). 1584 */ 1585 template <class MeshType> 1586 std::vector<typename MeshType::active_cell_iterator> 1587 compute_active_cell_layer_within_distance( 1588 const MeshType &mesh, 1589 const std::function<bool(const typename MeshType::active_cell_iterator &)> 1590 & predicate, 1591 const double layer_thickness); 1592 1593 /** 1594 * Extract and return a set of ghost cells which are within a 1595 * @p layer_thickness around all locally owned cells. 1596 * This is most relevant for parallel::shared::Triangulation 1597 * where it will return a subset of all ghost cells on a process, but for 1598 * parallel::distributed::Triangulation this will return all the ghost cells. 1599 * All the cells for the parallel::shared::Triangulation class that 1600 * are not owned by the current processor can be considered as ghost cells; 1601 * in particular, they do not only form a single layer of cells around the 1602 * locally owned ones. 1603 * 1604 * @tparam MeshType A type that satisfies the requirements of the 1605 * @ref ConceptMeshType "MeshType concept". 1606 * @param mesh A mesh (i.e. objects of type Triangulation, DoFHandler, 1607 * or hp::DoFHandler). 1608 * @param layer_thickness specifies the geometric distance within 1609 * which the function searches for active cells from the locally owned cells. 1610 * @return A subset of ghost cells within a given geometric distance of @p 1611 * layer_thickness from the locally owned cells of a current process. 1612 * 1613 * Also see compute_ghost_cell_halo_layer() and 1614 * compute_active_cell_layer_within_distance(). 1615 */ 1616 template <class MeshType> 1617 std::vector<typename MeshType::active_cell_iterator> 1618 compute_ghost_cell_layer_within_distance(const MeshType &mesh, 1619 const double layer_thickness); 1620 1621 /** 1622 * Compute and return a bounding box, defined through a pair of points 1623 * bottom left and top right, that surrounds a subdomain of the @p mesh. 1624 * Here, the "subdomain" consists of exactly all of those 1625 * active cells for which the @p predicate returns @p true. 1626 * 1627 * For a description of how @p predicate works, 1628 * see compute_active_cell_halo_layer(). 1629 * 1630 * @note This function was written before the BoundingBox class was invented. 1631 * Consequently, it returns a pair of points, rather than a BoundingBox 1632 * object as one may expect. However, BoundingBox has a conversion constructor 1633 * from pairs of points, so the result of this function can still be assigned 1634 * to a BoundingBox object. 1635 */ 1636 template <class MeshType> 1637 std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>> 1638 compute_bounding_box( 1639 const MeshType &mesh, 1640 const std::function<bool(const typename MeshType::active_cell_iterator &)> 1641 &predicate); 1642 1643 /** 1644 * Compute a collection of bounding boxes so that all active cells for which 1645 * the given predicate is true, are completely enclosed in at least one of the 1646 * bounding boxes. Notice the cover is only guaranteed to contain all these 1647 * active cells but it's not necessarily exact i.e. it can include a bigger 1648 * area than their union. 1649 * 1650 * For each cell at a given refinement level containing active cells for which @p predicate is true, 1651 * the function creates a bounding box of its children for which @p predicate is true. 1652 * 1653 * This results in a cover of all active cells for which @p predicate is true; the parameters 1654 * @p allow_merge and @p max_boxes are used to reduce the number of cells at a computational cost and 1655 * covering a bigger n-dimensional volume. 1656 * 1657 * The parameters to control the algorithm are: 1658 * - @p predicate : the property of the cells to enclose e.g. IteratorFilters::LocallyOwnedCell . 1659 * The predicate is tested only on active cells. 1660 * - @p refinement_level : it defines the level at which the initial bounding box are created. The refinement 1661 * should be set to a coarse refinement level. A bounding box is created for 1662 * each active cell at coarser 1663 * level than @p refinement_level; if @p refinement_level is higher than the number of levels of the 1664 * triangulation an exception is thrown. 1665 * - @p allow_merge : This flag allows for box merging and, by default, is false. The algorithm has a cost of 1666 * O(N^2) where N is the number of the bounding boxes created from the 1667 * refinement level; for this reason, if 1668 * the flag is set to true, make sure to choose wisely a coarse enough @p refinement_level. 1669 * - @p max_boxes : the maximum number of bounding boxes to compute. If more are created the smaller ones are 1670 * merged with neighbors. By default after merging the boxes which can be 1671 * expressed as a single one no more boxes are merged. See the 1672 * BoundingBox::get_neighbor_type () function for details. 1673 * Notice only neighboring cells are merged (see the @p get_neighbor_type function in bounding box class): if 1674 * the target number of bounding boxes max_boxes can't be reached by merging 1675 * neighbors an exception is thrown 1676 * 1677 * The following image describes an example of the algorithm with @p refinement_level = 2, @p allow_merge = true 1678 * and @p max_boxes = 1. The cells with the property predicate are in red, the area of a bounding box is 1679 * slightly orange. 1680 * @image html bounding_box_predicate.png 1681 * - 1. In black we can see the cells of the current level. 1682 * - 2. For each cell containing the red area a bounding box is created: by 1683 * default these are returned. 1684 * - 3. Because @p allow_merge = true the number of bounding boxes is reduced while not changing the cover. 1685 * If @p max_boxes was left as default or bigger than 1 these two boxes would be returned. 1686 * - 4. Because @p max_boxes = 1 the smallest bounding box is merged to the bigger one. 1687 * Notice it is important to choose the parameters wisely. For instance, @p allow_merge = false and 1688 * @p refinement_level = 1 returns the very same bounding box but with a 1689 * fraction of the computational cost. 1690 * 1691 * This function does not take into account the curvature of cells and thus it 1692 * is not suited for handling curved geometry: the mapping is assumed to be 1693 * linear. 1694 */ 1695 template <class MeshType> 1696 std::vector<BoundingBox<MeshType::space_dimension>> 1697 compute_mesh_predicate_bounding_box( 1698 const MeshType &mesh, 1699 const std::function<bool(const typename MeshType::active_cell_iterator &)> 1700 & predicate, 1701 const unsigned int refinement_level = 0, 1702 const bool allow_merge = false, 1703 const unsigned int max_boxes = numbers::invalid_unsigned_int); 1704 1705 /** 1706 * Given an array of points, use the global bounding box description obtained 1707 * using GridTools::compute_mesh_predicate_bounding_box to guess, for each of 1708 * them, which process might own it. 1709 * 1710 * @param[in] global_bboxes Vector of bounding boxes describing the portion of 1711 * mesh with a property for each process. 1712 * @param[in] points Array of points to test. 1713 * 1714 * @return A tuple containing the following information: 1715 * - A vector indicized with ranks of processes. For each rank it contains 1716 * a vector of the indices of points it might own. 1717 * - A map from the index <code>unsigned int</code> of the point in @p points 1718 * to the rank of the owner. 1719 * - A map from the index <code>unsigned int</code> of the point in @p points 1720 * to the ranks of the guessed owners. 1721 * 1722 * @note The actual return type of this function, i.e., the type referenced 1723 * above as @p return_type, is 1724 * @code 1725 * std::tuple<std::vector<std::vector<unsigned int>>, 1726 * std::map< unsigned int, unsigned int>, 1727 * std::map< unsigned int, std::vector<unsigned int>>> 1728 * @endcode 1729 * The type is abbreviated in the online documentation to improve readability 1730 * of this page. 1731 */ 1732 template <int spacedim> 1733 # ifndef DOXYGEN 1734 std::tuple<std::vector<std::vector<unsigned int>>, 1735 std::map<unsigned int, unsigned int>, 1736 std::map<unsigned int, std::vector<unsigned int>>> 1737 # else 1738 return_type 1739 # endif 1740 guess_point_owner( 1741 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes, 1742 const std::vector<Point<spacedim>> & points); 1743 1744 1745 /** 1746 * Given a covering rtree (see GridTools::Cache::get_covering_rtree()), and an 1747 * array of points, find a superset of processes which, individually, 1748 * may own the cell containing the points. 1749 * 1750 * For further details see GridTools::guess_point_owner; here only 1751 * different input/output types are reported: 1752 * 1753 * @param[in] covering_rtree RTRee which enables us to identify which 1754 * process(es) in a parallel computation may own the cell that 1755 * surrounds a given point. 1756 * 1757 * @param[in] points A vector of points to consider. 1758 * 1759 * @return A tuple containing the following information: 1760 * - A map indexed by processor ranks. For each rank it contains 1761 * a vector of the indices of points it might own. 1762 * - A map from the index <code>unsigned int</code> of the point in @p points 1763 * to the rank of the owner; these are points for which a single possible 1764 * owner was found. 1765 * - A map from the index <code>unsigned int</code> of the point in @p points 1766 * to the ranks of the guessed owners; these are points for which multiple 1767 * possible owners were found. 1768 * 1769 * @note The actual return type of this function, i.e., the type referenced 1770 * above as @p return_type, is 1771 * @code 1772 * std::tuple<std::map<unsigned int, std::vector<unsigned int>>, 1773 * std::map<unsigned int, unsigned int>, 1774 * std::map<unsigned int, std::vector<unsigned int>>> 1775 * @endcode 1776 * The type is abbreviated in the online documentation to improve readability 1777 * of this page. 1778 */ 1779 template <int spacedim> 1780 # ifndef DOXYGEN 1781 std::tuple<std::map<unsigned int, std::vector<unsigned int>>, 1782 std::map<unsigned int, unsigned int>, 1783 std::map<unsigned int, std::vector<unsigned int>>> 1784 # else 1785 return_type 1786 # endif 1787 guess_point_owner( 1788 const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree, 1789 const std::vector<Point<spacedim>> & points); 1790 1791 1792 /** 1793 * Return the adjacent cells of all the vertices. If a vertex is also a 1794 * hanging node, the associated coarse cell is also returned. The vertices 1795 * are ordered by the vertex index. This is the number returned by the 1796 * function <code>cell-@>vertex_index()</code>. Notice that only the indices 1797 * marked in the array returned by 1798 * Triangulation<dim,spacedim>::get_used_vertices() are used. 1799 */ 1800 template <int dim, int spacedim> 1801 std::vector< 1802 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>> 1803 vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation); 1804 1805 /** 1806 * Return a vector of normalized tensors for each vertex-cell combination of 1807 * the output of GridTools::vertex_to_cell_map() (which is expected as input 1808 * parameter for this function). Each tensor represents a geometric vector 1809 * from the vertex to the respective cell center. 1810 * 1811 * An assertion will be thrown if the size of the input vector is not equal to 1812 * the number of vertices of the triangulation. 1813 * 1814 * result[v][c] is a unit Tensor for vertex index v, indicating the direction 1815 * of the center of the c-th cell with respect to the vertex v. 1816 */ 1817 template <int dim, int spacedim> 1818 std::vector<std::vector<Tensor<1, spacedim>>> 1819 vertex_to_cell_centers_directions( 1820 const Triangulation<dim, spacedim> &mesh, 1821 const std::vector< 1822 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>> 1823 &vertex_to_cells); 1824 1825 1826 /** 1827 * Return the local vertex index of cell @p cell that is closest to 1828 * the given location @p position. The location of the vertices is extracted 1829 * from the (optional) @p mapping argument, to guarantee that the correct 1830 * answer is returned when the underlying mapping modifies the position of the 1831 * vertices. 1832 */ 1833 template <int dim, int spacedim> 1834 unsigned int 1835 find_closest_vertex_of_cell( 1836 const typename Triangulation<dim, spacedim>::active_cell_iterator &cell, 1837 const Point<spacedim> & position, 1838 const Mapping<dim, spacedim> & mapping = 1839 StaticMappingQ1<dim, spacedim>::mapping); 1840 1841 /** 1842 * Compute a globally unique index for each vertex and hanging node 1843 * associated with a locally owned active cell. The vertices of a ghost cell 1844 * that are hanging nodes of a locally owned cells have a global index. 1845 * However, the other vertices of the cells that do not <i>touch</i> an 1846 * active cell do not have a global index on this processor. 1847 * 1848 * The key of the map is the local index of the vertex and the value is the 1849 * global index. The indices need to be recomputed after refinement or 1850 * coarsening and may be different. 1851 */ 1852 template <int dim, int spacedim> 1853 std::map<unsigned int, types::global_vertex_index> 1854 compute_local_to_global_vertex_index_map( 1855 const parallel::distributed::Triangulation<dim, spacedim> &triangulation); 1856 1857 /** 1858 * Return the highest value among ratios between extents in each of the 1859 * coordinate directions of a @p cell. Moreover, return the dimension 1860 * relative to the highest elongation. 1861 * 1862 * @param[in] cell an iterator pointing to the cell. 1863 * 1864 * @return A std::pair<unsigned int, double> such that the @p first value 1865 * is the dimension of the highest elongation and the @p second value is the 1866 * ratio among the dimensions of the @p cell. 1867 */ 1868 template <int dim, int spacedim> 1869 std::pair<unsigned int, double> 1870 get_longest_direction( 1871 typename Triangulation<dim, spacedim>::active_cell_iterator cell); 1872 1873 /*@}*/ 1874 /** 1875 * @name Partitions and subdomains of triangulations 1876 */ 1877 /*@{*/ 1878 1879 /** 1880 * Produce a sparsity pattern in which nonzero entries indicate that two 1881 * cells are connected via a common face. The diagonal entries of the 1882 * sparsity pattern are also set. 1883 * 1884 * The rows and columns refer to the cells as they are traversed in their 1885 * natural order using cell iterators. 1886 */ 1887 template <int dim, int spacedim> 1888 void 1889 get_face_connectivity_of_cells( 1890 const Triangulation<dim, spacedim> &triangulation, 1891 DynamicSparsityPattern & connectivity); 1892 1893 /** 1894 * Produce a sparsity pattern in which nonzero entries indicate that two 1895 * cells are connected via a common vertex. The diagonal entries of the 1896 * sparsity pattern are also set. 1897 * 1898 * The rows and columns refer to the cells as they are traversed in their 1899 * natural order using cell iterators. 1900 */ 1901 template <int dim, int spacedim> 1902 void 1903 get_vertex_connectivity_of_cells( 1904 const Triangulation<dim, spacedim> &triangulation, 1905 DynamicSparsityPattern & connectivity); 1906 1907 /** 1908 * Produce a sparsity pattern for a given level mesh in which nonzero entries 1909 * indicate that two cells are connected via a common vertex. The diagonal 1910 * entries of the sparsity pattern are also set. 1911 * 1912 * The rows and columns refer to the cells as they are traversed in their 1913 * natural order using cell iterators. 1914 */ 1915 template <int dim, int spacedim> 1916 void 1917 get_vertex_connectivity_of_cells_on_level( 1918 const Triangulation<dim, spacedim> &triangulation, 1919 const unsigned int level, 1920 DynamicSparsityPattern & connectivity); 1921 1922 /** 1923 * Use graph partitioner to partition the active cells making up the entire 1924 * domain. After calling this function, the subdomain ids of all active cells 1925 * will have values 1926 * between zero and @p n_partitions-1. You can access the subdomain id of a cell by using 1927 * <tt>cell-@>subdomain_id()</tt>. 1928 * 1929 * Use the third argument to select between partitioning algorithms provided 1930 * by METIS or ZOLTAN. METIS is the default partitioner. 1931 * 1932 * If deal.II was not installed with ZOLTAN or METIS, this function will 1933 * generate an error 1934 * when corresponding partition method is chosen, unless @p n_partitions is one. 1935 * I.e., you can write a program so that it runs in the single-processor 1936 * single-partition case without packages installed, and only requires them 1937 * installed when multiple partitions are required. 1938 * 1939 * @note If the @p cell_weight signal has been attached to the @p triangulation, 1940 * then this will be used and passed to the partitioner. 1941 */ 1942 template <int dim, int spacedim> 1943 void 1944 partition_triangulation(const unsigned int n_partitions, 1945 Triangulation<dim, spacedim> & triangulation, 1946 const SparsityTools::Partitioner partitioner = 1947 SparsityTools::Partitioner::metis); 1948 1949 /** 1950 * This function performs the same operation as the one above, except that 1951 * it takes into consideration a specific set of @p cell_weights, which allow the 1952 * partitioner to balance the graph while taking into consideration the 1953 * computational effort expended on each cell. 1954 * 1955 * @note If the @p cell_weights vector is empty, then no weighting is taken 1956 * into consideration. If not then the size of this vector must equal to the 1957 * number of active cells in the triangulation. 1958 */ 1959 template <int dim, int spacedim> 1960 void 1961 partition_triangulation(const unsigned int n_partitions, 1962 const std::vector<unsigned int> &cell_weights, 1963 Triangulation<dim, spacedim> & triangulation, 1964 const SparsityTools::Partitioner partitioner = 1965 SparsityTools::Partitioner::metis); 1966 1967 /** 1968 * This function does the same as the previous one, i.e. it partitions a 1969 * triangulation using a partitioning algorithm into a number of subdomains 1970 * identified by the <code>cell-@>subdomain_id()</code> flag. 1971 * 1972 * The difference to the previous function is the second argument, a 1973 * sparsity pattern that represents the connectivity pattern between cells. 1974 * 1975 * While the function above builds it directly from the triangulation by 1976 * considering which cells neighbor each other, this function can take a 1977 * more refined connectivity graph. The sparsity pattern needs to be of size 1978 * $N\times N$, where $N$ is the number of active cells in the 1979 * triangulation. If the sparsity pattern contains an entry at position 1980 * $(i,j)$, then this means that cells $i$ and $j$ (in the order in which 1981 * they are traversed by active cell iterators) are to be considered 1982 * connected; partitioning algorithm will then try to partition the domain in 1983 * such a way that (i) the subdomains are of roughly equal size, and (ii) a 1984 * minimal number of connections are broken. 1985 * 1986 * This function is mainly useful in cases where connections between cells 1987 * exist that are not present in the triangulation alone (otherwise the 1988 * previous function would be the simpler one to use). Such connections may 1989 * include that certain parts of the boundary of a domain are coupled 1990 * through symmetric boundary conditions or integrals (e.g. friction contact 1991 * between the two sides of a crack in the domain), or if a numerical scheme 1992 * is used that not only connects immediate neighbors but a larger 1993 * neighborhood of cells (e.g. when solving integral equations). 1994 * 1995 * In addition, this function may be useful in cases where the default 1996 * sparsity pattern is not entirely sufficient. This can happen because the 1997 * default is to just consider face neighbors, not neighboring cells that 1998 * are connected by edges or vertices. While the latter couple when using 1999 * continuous finite elements, they are typically still closely connected in 2000 * the neighborship graph, and partitioning algorithm 2001 * will not usually cut important connections in this case. However, if there 2002 * are vertices in the mesh where many cells (many more than the common 4 or 6 2003 * in 2d and 3d, respectively) come together, then there will be a significant 2004 * number of cells that are connected across a vertex, but several degrees 2005 * removed in the connectivity graph built only using face neighbors. In a 2006 * case like this, partitioning algorithm may sometimes make bad decisions and 2007 * you may want to build your own connectivity graph. 2008 * 2009 * @note If the @p cell_weight signal has been attached to the @p triangulation, 2010 * then this will be used and passed to the partitioner. 2011 */ 2012 template <int dim, int spacedim> 2013 void 2014 partition_triangulation(const unsigned int n_partitions, 2015 const SparsityPattern & cell_connection_graph, 2016 Triangulation<dim, spacedim> &triangulation, 2017 const SparsityTools::Partitioner partitioner = 2018 SparsityTools::Partitioner::metis); 2019 2020 /** 2021 * This function performs the same operation as the one above, except that 2022 * it takes into consideration a specific set of @p cell_weights, which allow the 2023 * partitioner to balance the graph while taking into consideration the 2024 * computational effort expended on each cell. 2025 * 2026 * @note If the @p cell_weights vector is empty, then no weighting is taken 2027 * into consideration. If not then the size of this vector must equal to the 2028 * number of active cells in the triangulation. 2029 */ 2030 template <int dim, int spacedim> 2031 void 2032 partition_triangulation(const unsigned int n_partitions, 2033 const std::vector<unsigned int> &cell_weights, 2034 const SparsityPattern & cell_connection_graph, 2035 Triangulation<dim, spacedim> &triangulation, 2036 const SparsityTools::Partitioner partitioner = 2037 SparsityTools::Partitioner::metis); 2038 2039 /** 2040 * Generates a partitioning of the active cells making up the entire domain 2041 * using the same partitioning scheme as in the p4est library if the flag 2042 * @p group_siblings is set to true (default behavior of this function). 2043 * After calling this function, the subdomain ids of all active cells will 2044 * have values between zero and @p n_partitions-1. You can access the 2045 * subdomain id of a cell by using <tt>cell-@>subdomain_id()</tt>. 2046 * 2047 * @note If the flag @p group_siblings is set to false, children of a 2048 * cell might be placed on different processors even though they are all 2049 * active, which is an assumption made by p4est. By relaxing this, we 2050 * can create partitions owning a single cell (also for refined 2051 * meshes). 2052 */ 2053 template <int dim, int spacedim> 2054 void 2055 partition_triangulation_zorder(const unsigned int n_partitions, 2056 Triangulation<dim, spacedim> &triangulation, 2057 const bool group_siblings = true); 2058 2059 /** 2060 * Partitions the cells of a multigrid hierarchy by assigning level subdomain 2061 * ids using the "youngest child" rule, that is, each cell in the hierarchy is 2062 * owned by the processor who owns its left most child in the forest, and 2063 * active cells have the same subdomain id and level subdomain id. You can 2064 * access the level subdomain id of a cell by using 2065 * <tt>cell-@>level_subdomain_id()</tt>. 2066 * 2067 * Note: This function assumes that the active cells have already been 2068 * partitioned. 2069 */ 2070 template <int dim, int spacedim> 2071 void 2072 partition_multigrid_levels(Triangulation<dim, spacedim> &triangulation); 2073 2074 /** 2075 * For each active cell, return in the output array to which subdomain (as 2076 * given by the <tt>cell->subdomain_id()</tt> function) it belongs. The 2077 * output array is supposed to have the right size already when calling this 2078 * function. 2079 * 2080 * This function returns the association of each cell with one subdomain. If 2081 * you are looking for the association of each @em DoF with a subdomain, use 2082 * the <tt>DoFTools::get_subdomain_association</tt> function. 2083 */ 2084 template <int dim, int spacedim> 2085 void 2086 get_subdomain_association(const Triangulation<dim, spacedim> &triangulation, 2087 std::vector<types::subdomain_id> & subdomain); 2088 2089 /** 2090 * Count how many cells are uniquely associated with the given @p subdomain 2091 * index. 2092 * 2093 * This function may return zero if there are no cells with the given @p 2094 * subdomain index. This can happen, for example, if you try to partition a 2095 * coarse mesh into more partitions (one for each processor) than there are 2096 * cells in the mesh. 2097 * 2098 * This function returns the number of cells associated with one subdomain. 2099 * If you are looking for the association of @em DoFs with this subdomain, 2100 * use the <tt>DoFTools::count_dofs_with_subdomain_association</tt> 2101 * function. 2102 */ 2103 template <int dim, int spacedim> 2104 unsigned int 2105 count_cells_with_subdomain_association( 2106 const Triangulation<dim, spacedim> &triangulation, 2107 const types::subdomain_id subdomain); 2108 2109 2110 /** 2111 * For a triangulation, return a mask that represents which of its vertices 2112 * are "owned" by the current process in the same way as we talk about 2113 * locally owned cells or degrees of freedom (see 2114 * @ref GlossLocallyOwnedCell 2115 * and 2116 * @ref GlossLocallyOwnedDof). 2117 * For the purpose of this function, we define a locally owned vertex as 2118 * follows: a vertex is owned by that processor with the smallest subdomain 2119 * id (which equals the MPI rank of that processor) among all owners of 2120 * cells adjacent to this vertex. In other words, vertices that are in the 2121 * interior of a partition of the triangulation are owned by the owner of 2122 * this partition; for vertices that lie on the boundary between two or more 2123 * partitions, the owner is the processor with the least subdomain_id among 2124 * all adjacent subdomains. 2125 * 2126 * For sequential triangulations (as opposed to, for example, 2127 * parallel::distributed::Triangulation), every user vertex is of course 2128 * owned by the current processor, i.e., the function returns 2129 * Triangulation::get_used_vertices(). For parallel triangulations, the 2130 * returned mask is a subset of what Triangulation::get_used_vertices() 2131 * returns. 2132 * 2133 * @param triangulation The triangulation of which the function evaluates 2134 * which vertices are locally owned. 2135 * @return The subset of vertices, as described above. The length of the 2136 * returned array equals Triangulation.n_vertices() and may, consequently, 2137 * be larger than Triangulation::n_used_vertices(). 2138 */ 2139 template <int dim, int spacedim> 2140 std::vector<bool> 2141 get_locally_owned_vertices(const Triangulation<dim, spacedim> &triangulation); 2142 2143 /*@}*/ 2144 /** 2145 * @name Comparing different meshes 2146 */ 2147 /*@{*/ 2148 2149 /** 2150 * Given two meshes (i.e. objects of type Triangulation, DoFHandler, or 2151 * hp::DoFHandler) that are based on the same coarse mesh, this function 2152 * figures out a set of cells that are matched between the two meshes and 2153 * where at most one of the meshes is more refined on this cell. In other 2154 * words, it finds the smallest cells that are common to both meshes, and 2155 * that together completely cover the domain. 2156 * 2157 * This function is useful, for example, in time-dependent or nonlinear 2158 * application, where one has to integrate a solution defined on one mesh 2159 * (e.g., the one from the previous time step or nonlinear iteration) 2160 * against the shape functions of another mesh (the next time step, the next 2161 * nonlinear iteration). If, for example, the new mesh is finer, then one 2162 * has to obtain the solution on the coarse mesh (mesh_1) and interpolate it 2163 * to the children of the corresponding cell of mesh_2. Conversely, if the 2164 * new mesh is coarser, one has to express the coarse cell shape function by 2165 * a linear combination of fine cell shape functions. In either case, one 2166 * needs to loop over the finest cells that are common to both 2167 * triangulations. This function returns a list of pairs of matching 2168 * iterators to cells in the two meshes that can be used to this end. 2169 * 2170 * Note that the list of these iterators is not necessarily ordered, and 2171 * does also not necessarily coincide with the order in which cells are 2172 * traversed in one, or both, of the meshes given as arguments. 2173 * 2174 * @tparam MeshType A type that satisfies the requirements of the 2175 * @ref ConceptMeshType "MeshType concept". 2176 * 2177 * @note This function can only be used with 2178 * parallel::distributed::Triangulation when both meshes use the same 2179 * Triangulation since, with a distributed Triangulation, not all cells are 2180 * stored locally, so the resulting list may not cover the entire domain. 2181 */ 2182 template <typename MeshType> 2183 std::list<std::pair<typename MeshType::cell_iterator, 2184 typename MeshType::cell_iterator>> 2185 get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2); 2186 2187 /** 2188 * Return true if the two triangulations are based on the same coarse mesh. 2189 * This is determined by checking whether they have the same number of cells 2190 * on the coarsest level, and then checking that they have the same 2191 * vertices. 2192 * 2193 * The two meshes may have different refinement histories beyond the coarse 2194 * mesh. 2195 */ 2196 template <int dim, int spacedim> 2197 bool 2198 have_same_coarse_mesh(const Triangulation<dim, spacedim> &mesh_1, 2199 const Triangulation<dim, spacedim> &mesh_2); 2200 2201 /** 2202 * The same function as above, but working on arguments of type DoFHandler, 2203 * or hp::DoFHandler. This function is provided to allow calling 2204 * have_same_coarse_mesh for all types of containers representing 2205 * triangulations or the classes built on triangulations. 2206 * 2207 * @tparam MeshType A type that satisfies the requirements of the 2208 * @ref ConceptMeshType "MeshType concept". 2209 */ 2210 template <typename MeshType> 2211 bool 2212 have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2); 2213 2214 /*@}*/ 2215 /** 2216 * @name Dealing with distorted cells 2217 */ 2218 /*@{*/ 2219 2220 /** 2221 * Given a triangulation and a list of cells whose children have become 2222 * distorted as a result of mesh refinement, try to fix these cells up by 2223 * moving the center node around. 2224 * 2225 * The function returns a list of cells with distorted children that 2226 * couldn't be fixed up for whatever reason. The returned list is therefore 2227 * a subset of the input argument. 2228 * 2229 * For a definition of the concept of distorted cells, see the 2230 * @ref GlossDistorted "glossary entry". 2231 * The first argument passed to the current function is typically the 2232 * exception thrown by the Triangulation::execute_coarsening_and_refinement 2233 * function. 2234 */ 2235 template <int dim, int spacedim> 2236 typename Triangulation<dim, spacedim>::DistortedCellList 2237 fix_up_distorted_child_cells( 2238 const typename Triangulation<dim, spacedim>::DistortedCellList 2239 & distorted_cells, 2240 Triangulation<dim, spacedim> &triangulation); 2241 2242 2243 2244 /*@}*/ 2245 /** 2246 * @name Extracting and creating patches of cells 2247 * 2248 * These functions extract and create patches of cells surrounding a single 2249 * cell, and creating triangulation out of them. 2250 */ 2251 /*@{*/ 2252 2253 2254 /** 2255 * This function returns a list of all the active neighbor cells of the 2256 * given, active cell. Here, a neighbor is defined as one having at least 2257 * part of a face in common with the given cell, but not edge (in 3d) or 2258 * vertex neighbors (in 2d and 3d). 2259 * 2260 * The first element of the returned list is the cell provided as argument. 2261 * The remaining ones are neighbors: The function loops over all faces of 2262 * that given cell and checks if that face is not on the boundary of the 2263 * domain. Then, if the neighbor cell does not have any children (i.e., it 2264 * is either at the same refinement level as the current cell, or coarser) 2265 * then this neighbor cell is added to the list of cells. Otherwise, if the 2266 * neighbor cell is refined and therefore has children, then this function 2267 * loops over all subfaces of current face adds the neighbors behind these 2268 * sub-faces to the list to be returned. 2269 * 2270 * @tparam MeshType A type that satisfies the requirements of the 2271 * @ref ConceptMeshType "MeshType concept". 2272 * In C++, the compiler can not determine <code>MeshType</code> from the 2273 * function call. You need to specify it as an explicit template argument 2274 * following the function name. 2275 * @param[in] cell An iterator pointing to a cell of the mesh. 2276 * @return A list of active cells that form the patch around the given cell 2277 * 2278 * @note Patches are often used in defining error estimators that require 2279 * the solution of a local problem on the patch surrounding each of the 2280 * cells of the mesh. This also requires manipulating the degrees of freedom 2281 * associated with the cells of a patch. To this end, there are further 2282 * functions working on patches in namespace DoFTools. 2283 * 2284 * @note In the context of a parallel distributed computation, it only makes 2285 * sense to call this function on locally owned cells. This is because the 2286 * neighbors of locally owned cells are either locally owned themselves, or 2287 * ghost cells. For both, we know that these are in fact the real cells of 2288 * the complete, parallel triangulation. We can also query the degrees of 2289 * freedom on these. 2290 */ 2291 template <class MeshType> 2292 std::vector<typename MeshType::active_cell_iterator> 2293 get_patch_around_cell(const typename MeshType::active_cell_iterator &cell); 2294 2295 2296 /** 2297 * This function takes a vector of active cells (hereafter named @p 2298 * patch_cells) as input argument, and returns a vector of their parent 2299 * cells with the coarsest common level of refinement. In other words, find 2300 * that set of cells living at the same refinement level so that all cells 2301 * in the input vector are children of the cells in the set, or are in the 2302 * set itself. 2303 * 2304 * @tparam Container In C++, the compiler can not determine the type of 2305 * <code>Container</code> from the function call. You need to specify it as 2306 * an explicit template argument following the function name. This type has 2307 * to satisfy the requirements of a mesh container (see 2308 * @ref ConceptMeshType). 2309 * 2310 * @param[in] patch_cells A vector of active cells for which this function 2311 * finds the parents at the coarsest common level. This vector of cells 2312 * typically results from calling the function 2313 * GridTools::get_patch_around_cell(). 2314 * @return A list of cells with the coarsest common level of refinement of 2315 * the input cells. 2316 */ 2317 template <class Container> 2318 std::vector<typename Container::cell_iterator> 2319 get_cells_at_coarsest_common_level( 2320 const std::vector<typename Container::active_cell_iterator> &patch_cells); 2321 2322 /** 2323 * This function constructs a Triangulation (named @p local_triangulation) 2324 * from a given vector of active cells. This vector (which we think of the 2325 * cells corresponding to a "patch") contains active cells that are part of 2326 * an existing global Triangulation. The goal of this function is to build a 2327 * local Triangulation that contains only the active cells given in @p patch 2328 * (and potentially a minimum number of additional cells required to form a 2329 * valid Triangulation). The function also returns a map that allows to 2330 * identify the cells in the output Triangulation and corresponding cells in 2331 * the input list. 2332 * 2333 * The function copies the location of vertices of cells from the cells of the 2334 * source triangulation to the triangulation that is built from the list of 2335 * patch cells. This adds support for triangulations which have been 2336 * perturbed or smoothed in some manner which makes the triangulation 2337 * deviate from the standard deal.II refinement strategy of placing new 2338 * vertices at midpoints of faces or edges. 2339 * 2340 * The operation implemented by this function is frequently used in the 2341 * definition of error estimators that need to solve "local" problems on 2342 * each cell and its neighbors. A similar construction is necessary in the 2343 * definition of the Clement interpolation operator in which one needs to 2344 * solve a local problem on all cells within the support of a shape 2345 * function. This function then builds a complete Triangulation from a list 2346 * of cells that make up such a patch; one can then later attach a 2347 * DoFHandler to such a Triangulation. 2348 * 2349 * If the list of input cells contains only cells at the same refinement 2350 * level, then the output Triangulation simply consists of a Triangulation 2351 * containing only exactly these patch cells. On the other hand, if the 2352 * input cells live on different refinement levels, i.e., the Triangulation 2353 * of which they are part is adaptively refined, then the construction of 2354 * the output Triangulation is not so simple because the coarsest level of a 2355 * Triangulation can not contain hanging nodes. Rather, we first have to 2356 * find the common refinement level of all input cells, along with their 2357 * common parents (see GridTools::get_cells_at_coarsest_common_level()), 2358 * build a Triangulation from those, and then adaptively refine it so that 2359 * the input cells all also exist in the output Triangulation. 2360 * 2361 * A consequence of this procedure is that the output Triangulation may 2362 * contain more active cells than the ones that exist in the input vector. 2363 * On the other hand, one typically wants to solve the local problem not on 2364 * the entire output Triangulation, but only on those cells of it that 2365 * correspond to cells in the input list. In this case, a user typically 2366 * wants to assign degrees of freedom only on cells that are part of the 2367 * "patch", and somehow ignore those excessive cells. The current function 2368 * supports this common requirement by setting the user flag for the cells 2369 * in the output Triangulation that match with cells in the input list. 2370 * Cells which are not part of the original patch will not have their @p 2371 * user_flag set; we can then avoid assigning degrees of freedom using the 2372 * FE_Nothing<dim> element. 2373 * 2374 * @tparam Container In C++, the compiler can not determine the type of 2375 * <code>Container</code> from the function call. You need to specify it as 2376 * an explicit template argument following the function name. This type that 2377 * satisfies the requirements of a mesh container (see 2378 * @ref ConceptMeshType). 2379 * 2380 * @param[in] patch A vector of active cells from a common triangulation. 2381 * These cells may or may not all be at the same refinement level. 2382 * @param[out] local_triangulation A triangulation whose active cells 2383 * correspond to the given vector of active cells in @p patch. 2384 * @param[out] patch_to_global_tria_map A map between the local 2385 * triangulation which is built as explained above, and the cell iterators 2386 * in the input list. 2387 */ 2388 template <class Container> 2389 void 2390 build_triangulation_from_patch( 2391 const std::vector<typename Container::active_cell_iterator> &patch, 2392 Triangulation<Container::dimension, Container::space_dimension> 2393 &local_triangulation, 2394 std::map< 2395 typename Triangulation<Container::dimension, 2396 Container::space_dimension>::active_cell_iterator, 2397 typename Container::active_cell_iterator> &patch_to_global_tria_map); 2398 2399 /** 2400 * This function runs through the degrees of freedom defined by the 2401 * DoFHandler and for each dof constructs a vector of 2402 * active_cell_iterators representing the cells of support of the associated 2403 * basis element at that degree of freedom. This function was originally 2404 * designed for the implementation of local projections, for instance the 2405 * Clement interpolant, in conjunction with other local patch functions like 2406 * GridTools::build_triangulation_from_patch. 2407 * 2408 * DoFHandler's built on top of Triangulation or 2409 * parallel:distributed::Triangulation are supported and handled 2410 * appropriately. 2411 * 2412 * The result is the patch of cells representing the support of the basis 2413 * element associated to the degree of freedom. For instance using an FE_Q 2414 * finite element, we obtain the standard patch of cells touching the degree 2415 * of freedom and then add other cells that take care of possible hanging node 2416 * constraints. Using a FE_DGQ finite element, the degrees of freedom are 2417 * logically considered to be "interior" to the cells so the patch would 2418 * consist exclusively of the single cell on which the degree of freedom is 2419 * located. 2420 * 2421 * @param[in] dof_handler The DoFHandler which could be built on a 2422 * Triangulation or a parallel::distributed::Triangulation with a finite 2423 * element that has degrees of freedom that are logically associated to a 2424 * vertex, line, quad, or hex. 2425 * @return A map from the global_dof_index of 2426 * degrees of freedom on locally relevant cells to vectors containing 2427 * DoFHandler::active_cell_iterators of cells in the support of the basis 2428 * function at that degree of freedom. 2429 */ 2430 template <int dim, int spacedim> 2431 std::map< 2432 types::global_dof_index, 2433 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>> 2434 get_dof_to_support_patch_map(DoFHandler<dim, spacedim> &dof_handler); 2435 2436 2437 /*@}*/ 2438 2439 /** 2440 * @name Dealing with periodic domains 2441 */ 2442 /*@{*/ 2443 2444 /** 2445 * Data type that provides all information necessary to create periodicity 2446 * constraints and a periodic p4est forest with respect to two 'periodic' 2447 * cell faces. 2448 */ 2449 template <typename CellIterator> 2450 struct PeriodicFacePair 2451 { 2452 /** 2453 * The cells associated with the two 'periodic' faces. 2454 */ 2455 CellIterator cell[2]; 2456 2457 /** 2458 * The local face indices (with respect to the specified cells) of the two 2459 * 'periodic' faces. 2460 */ 2461 unsigned int face_idx[2]; 2462 2463 /** 2464 * The relative orientation of the first face with respect to the second 2465 * face as described in orthogonal_equality() and 2466 * DoFTools::make_periodicity_constraints() (and stored as a bitset). 2467 */ 2468 std::bitset<3> orientation; 2469 2470 /** 2471 * A @p dim $\times$ @p dim rotation matrix that describes how vector 2472 * valued DoFs of the first face should be modified prior to constraining 2473 * to the DoFs of the second face. 2474 * 2475 * The rotation matrix is used in DoFTools::make_periodicity_constraints() 2476 * by applying the rotation to all vector valued blocks listed in the 2477 * parameter @p first_vector_components of the finite element space. For 2478 * more details see DoFTools::make_periodicity_constraints() and the 2479 * glossary 2480 * @ref GlossPeriodicConstraints "glossary entry on periodic conditions". 2481 */ 2482 FullMatrix<double> matrix; 2483 }; 2484 2485 2486 /** 2487 * An orthogonal equality test for faces. 2488 * 2489 * @p face1 and @p face2 are considered equal, if a one to one matching 2490 * between its vertices can be achieved via an orthogonal equality relation. 2491 * 2492 * Here, two vertices <tt>v_1</tt> and <tt>v_2</tt> are considered equal, if 2493 * $M\cdot v_1 + offset - v_2$ is parallel to the unit vector in unit 2494 * direction @p direction. If the parameter @p matrix is a reference to a 2495 * spacedim x spacedim matrix, $M$ is set to @p matrix, otherwise $M$ is the 2496 * identity matrix. 2497 * 2498 * If the matching was successful, the _relative_ orientation of @p face1 2499 * with respect to @p face2 is returned in the bitset @p orientation, where 2500 * @code 2501 * orientation[0] -> face_orientation 2502 * orientation[1] -> face_flip 2503 * orientation[2] -> face_rotation 2504 * @endcode 2505 * 2506 * In 2D <tt>face_orientation</tt> is always <tt>true</tt>, 2507 * <tt>face_rotation</tt> is always <tt>false</tt>, and face_flip has the 2508 * meaning of <tt>line_flip</tt>. More precisely in 3d: 2509 * 2510 * <tt>face_orientation</tt>: <tt>true</tt> if @p face1 and @p face2 have 2511 * the same orientation. Otherwise, the vertex indices of @p face1 match the 2512 * vertex indices of @p face2 in the following manner: 2513 * 2514 * @code 2515 * face1: face2: 2516 * 2517 * 1 - 3 2 - 3 2518 * | | <--> | | 2519 * 0 - 2 0 - 1 2520 * @endcode 2521 * 2522 * <tt>face_flip</tt>: <tt>true</tt> if the matched vertices are rotated by 2523 * 180 degrees: 2524 * 2525 * @code 2526 * face1: face2: 2527 * 2528 * 1 - 0 2 - 3 2529 * | | <--> | | 2530 * 3 - 2 0 - 1 2531 * @endcode 2532 * 2533 * <tt>face_rotation</tt>: <tt>true</tt> if the matched vertices are rotated 2534 * by 90 degrees counterclockwise: 2535 * 2536 * @code 2537 * face1: face2: 2538 * 2539 * 0 - 2 2 - 3 2540 * | | <--> | | 2541 * 1 - 3 0 - 1 2542 * @endcode 2543 * 2544 * and any combination of that... More information on the topic can be found 2545 * in the 2546 * @ref GlossFaceOrientation "glossary" 2547 * article. 2548 */ 2549 template <typename FaceIterator> 2550 bool orthogonal_equality( 2551 std::bitset<3> & orientation, 2552 const FaceIterator & face1, 2553 const FaceIterator & face2, 2554 const int direction, 2555 const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset = 2556 Tensor<1, FaceIterator::AccessorType::space_dimension>(), 2557 const FullMatrix<double> &matrix = FullMatrix<double>()); 2558 2559 2560 /** 2561 * Same function as above, but doesn't return the actual orientation 2562 */ 2563 template <typename FaceIterator> 2564 bool 2565 orthogonal_equality( 2566 const FaceIterator & face1, 2567 const FaceIterator & face2, 2568 const int direction, 2569 const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset = 2570 Tensor<1, FaceIterator::AccessorType::space_dimension>(), 2571 const FullMatrix<double> &matrix = FullMatrix<double>()); 2572 2573 2574 /** 2575 * This function will collect periodic face pairs on the coarsest mesh level 2576 * of the given @p mesh (a Triangulation or DoFHandler) and add them to the 2577 * vector @p matched_pairs leaving the original contents intact. 2578 * 2579 * Define a 'first' boundary as all boundary faces having boundary_id @p 2580 * b_id1 and a 'second' boundary consisting of all faces belonging to @p 2581 * b_id2. 2582 * 2583 * This function tries to match all faces belonging to the first boundary 2584 * with faces belonging to the second boundary with the help of 2585 * orthogonal_equality(). 2586 * 2587 * The bitset that is returned inside of PeriodicFacePair encodes the 2588 * _relative_ orientation of the first face with respect to the second face, 2589 * see the documentation of orthogonal_equality() for further details. 2590 * 2591 * The @p direction refers to the space direction in which periodicity is 2592 * enforced. When matching periodic faces this vector component is ignored. 2593 * 2594 * The @p offset is a vector tangential to the faces that is added to the 2595 * location of vertices of the 'first' boundary when attempting to match 2596 * them to the corresponding vertices of the 'second' boundary. This can be 2597 * used to implement conditions such as $u(0,y)=u(1,y+1)$. 2598 * 2599 * Optionally, a $dim\times dim$ rotation @p matrix can be specified that 2600 * describes how vector valued DoFs of the first face should be modified 2601 * prior to constraining to the DoFs of the second face. The @p matrix is 2602 * used in two places. First, @p matrix will be supplied to 2603 * orthogonal_equality() and used for matching faces: Two vertices $v_1$ and 2604 * $v_2$ match if $\text{matrix}\cdot v_1 + \text{offset} - v_2$ is parallel 2605 * to the unit vector in unit direction @p direction. (For more details see 2606 * DoFTools::make_periodicity_constraints(), the glossary 2607 * @ref GlossPeriodicConstraints "glossary entry on periodic conditions" 2608 * and step-45). Second, @p matrix will be stored in the PeriodicFacePair 2609 * collection @p matched_pairs for further use. 2610 * 2611 * @tparam MeshType A type that satisfies the requirements of the 2612 * @ref ConceptMeshType "MeshType concept". 2613 * 2614 * @note The created std::vector can be used in 2615 * DoFTools::make_periodicity_constraints() and in 2616 * parallel::distributed::Triangulation::add_periodicity() to enforce 2617 * periodicity algebraically. 2618 * 2619 * @note Because elements will be added to @p matched_pairs (and existing 2620 * entries will be preserved), it is possible to call this function several 2621 * times with different boundary ids to generate a vector with all periodic 2622 * pairs. 2623 * 2624 * @note Since the periodic face pairs are found on the coarsest mesh level, 2625 * it is necessary to ensure that the coarsest level faces have the correct 2626 * boundary indicators set. In general, this means that one must first set 2627 * all boundary indicators on the coarse grid before performing any global 2628 * or local grid refinement. 2629 */ 2630 template <typename MeshType> 2631 void 2632 collect_periodic_faces( 2633 const MeshType & mesh, 2634 const types::boundary_id b_id1, 2635 const types::boundary_id b_id2, 2636 const int direction, 2637 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>> 2638 & matched_pairs, 2639 const Tensor<1, MeshType::space_dimension> &offset = 2640 dealii::Tensor<1, MeshType::space_dimension>(), 2641 const FullMatrix<double> &matrix = FullMatrix<double>()); 2642 2643 2644 /** 2645 * This compatibility version of collect_periodic_faces() only works on 2646 * grids with cells in 2647 * @ref GlossFaceOrientation "standard orientation". 2648 * 2649 * Instead of defining a 'first' and 'second' boundary with the help of two 2650 * boundary_ids this function defines a 'left' boundary as all faces with 2651 * local face index <code>2*dimension</code> and boundary indicator @p b_id 2652 * and, similarly, a 'right' boundary consisting of all face with local face 2653 * index <code>2*dimension+1</code> and boundary indicator @p b_id. Faces 2654 * with coordinates only differing in the @p direction component are 2655 * identified. 2656 * 2657 * This function will collect periodic face pairs on the coarsest mesh level 2658 * and add them to @p matched_pairs leaving the original contents intact. 2659 * 2660 * See above function for further details. 2661 * 2662 * @note This version of collect_periodic_faces() will not work on 2663 * meshes with cells not in 2664 * @ref GlossFaceOrientation "standard orientation". 2665 */ 2666 template <typename MeshType> 2667 void 2668 collect_periodic_faces( 2669 const MeshType & mesh, 2670 const types::boundary_id b_id, 2671 const int direction, 2672 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>> 2673 & matched_pairs, 2674 const dealii::Tensor<1, MeshType::space_dimension> &offset = 2675 dealii::Tensor<1, MeshType::space_dimension>(), 2676 const FullMatrix<double> &matrix = FullMatrix<double>()); 2677 2678 /*@}*/ 2679 /** 2680 * @name Dealing with boundary and manifold ids 2681 */ 2682 /*@{*/ 2683 2684 /** 2685 * Copy boundary ids to manifold ids on faces and edges at the boundary. The 2686 * default manifold_id for new Triangulation objects is 2687 * numbers::flat_manifold_id. This function copies the boundary_ids of 2688 * the boundary faces and edges to the manifold_ids of the same faces and 2689 * edges, allowing the user to change the boundary_ids and use them for 2690 * boundary conditions regardless of the geometry, which will use 2691 * manifold_ids to create new points. Only active cells will be iterated 2692 * over. This is a function you'd typically call when there is only one 2693 * active level on your Triangulation. Mesh refinement will then inherit 2694 * these indicators to child cells, faces, and edges. 2695 * 2696 * The optional parameter @p reset_boundary_ids, indicates whether this 2697 * function should reset the boundary_ids of boundary faces and edges to its 2698 * default value 0 after copying its value to the manifold_id. By default, 2699 * boundary_ids are left untouched. 2700 * 2701 * @ingroup manifold 2702 * @relatesalso boundary 2703 */ 2704 template <int dim, int spacedim> 2705 void 2706 copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria, 2707 const bool reset_boundary_ids = false); 2708 2709 /** 2710 * Map the given boundary ids to the given manifold ids on faces and 2711 * edges at the boundary. 2712 * 2713 * This function copies the boundary ids of the boundary faces and 2714 * edges that are present in the parameter @p src_boundary_ids to 2715 * the corresponding manifold id in @p dst_manifold_ids, of the same 2716 * faces and edges. 2717 * 2718 * If the optional parameter @p reset_boundary_ids is non empty, 2719 * each boundary id in @p src_boundary_ids, is replaced with the 2720 * corresponding boundary id in @p reset_boundary_ids. 2721 * 2722 * An exception is thrown if the size of the input vectors do not 2723 * match. If a boundary id indicated in @p src_boundary_ids is not 2724 * present in the triangulation, it is simply ignored during the 2725 * process. 2726 * 2727 * @ingroup manifold 2728 * @relatesalso boundary 2729 */ 2730 template <int dim, int spacedim> 2731 void 2732 map_boundary_to_manifold_ids( 2733 const std::vector<types::boundary_id> &src_boundary_ids, 2734 const std::vector<types::manifold_id> &dst_manifold_ids, 2735 Triangulation<dim, spacedim> & tria, 2736 const std::vector<types::boundary_id> &reset_boundary_ids = {}); 2737 2738 /** 2739 * Copy material ids to manifold ids. The default manifold_id for new 2740 * Triangulation objects is numbers::flat_manifold_id. When refinements 2741 * occurs, the Triangulation asks where to locate new points to the 2742 * underlying manifold. 2743 * 2744 * When reading a Triangulation from a supported input format, typical 2745 * information that can be stored in a file are boundary conditions for 2746 * boundary faces (which we store in the boundary_id of the faces), material 2747 * types for cells (which we store in the material_id of the cells) and in 2748 * some cases subdomain ids for cells (which we store in the subdomain_id of 2749 * the cell). 2750 * 2751 * If you read one of these grids into a Triangulation, you might still want 2752 * to use the material_id specified in the input file as a manifold_id 2753 * description. In this case you can associate a Manifold object to internal 2754 * cells, and this object will be used by the Triangulation to query 2755 * Manifold objects for new points. This function iterates over active cells 2756 * and copies the material_ids to the manifold_ids. 2757 * 2758 * The optional parameter @p compute_face_ids, indicates whether this 2759 * function should also set the manifold_ids of the faces (both for internal 2760 * faces and for faces on the boundary). If set to true, then each face will 2761 * get a manifold_id equal to the minimum of the surrounding manifold_ids, 2762 * ensuring that a unique manifold id is selected for each face of the 2763 * Triangulation. By default, face manifold_ids are not computed. 2764 * 2765 * @ingroup manifold 2766 */ 2767 template <int dim, int spacedim> 2768 void 2769 copy_material_to_manifold_id(Triangulation<dim, spacedim> &tria, 2770 const bool compute_face_ids = false); 2771 2772 /** 2773 * Propagate manifold indicators associated with the cells of the 2774 * Triangulation @p tria to their co-dimension one and two objects. 2775 * 2776 * This function sets the @p manifold_id of faces and edges (both on the 2777 * interior and on the boundary) to the value returned by the 2778 * @p disambiguation_function method, called with the set of 2779 * manifold indicators of the cells that share the same face or edge. 2780 * 2781 * By default, the @p disambiguation_function returns 2782 * numbers::flat_manifold_id when the set has size greater than one (i.e., 2783 * when it is not possible to decide what manifold indicator a face or edge 2784 * should have according to the manifold indicators of the adjacent cells) 2785 * and it returns the manifold indicator contained in the set when it has 2786 * dimension one (i.e., when all adjacent cells and faces have the same 2787 * manifold indicator). 2788 * 2789 * The parameter @p overwrite_only_flat_manifold_ids allows you to specify 2790 * what to do when a face or an edge already has a manifold indicator 2791 * different from numbers::flat_manifold_id. If the flag is @p true, the edge 2792 * or face will maintain its original manifold indicator. 2793 * If it is @p false, then also the manifold indicator of these faces and edges 2794 * is set according to the return value of the @p disambiguation_function. 2795 */ 2796 template <int dim, int spacedim> 2797 void 2798 assign_co_dimensional_manifold_indicators( 2799 Triangulation<dim, spacedim> & tria, 2800 const std::function<types::manifold_id( 2801 const std::set<types::manifold_id> &)> &disambiguation_function = 2802 [](const std::set<types::manifold_id> &manifold_ids) { 2803 if (manifold_ids.size() == 1) 2804 return *manifold_ids.begin(); 2805 else 2806 return numbers::flat_manifold_id; 2807 }, 2808 bool overwrite_only_flat_manifold_ids = true); 2809 /*@}*/ 2810 2811 /** 2812 * Exchange arbitrary data of type @p DataType provided by the function 2813 * objects from locally owned cells to ghost cells on other processors. 2814 * 2815 * After this call, you typically will have received data from @p unpack on 2816 * every ghost cell as it was given by @p pack on the owning processor. 2817 * Whether you do or do not receive information to @p unpack on a given 2818 * ghost cell depends on whether the @p pack function decided that 2819 * something needs to be sent. It does so using the std_cxx17::optional 2820 * mechanism: if the std_cxx17::optional return object of the @p pack 2821 * function is empty, then this implies that no data has to be sent for 2822 * the locally owned cell it was called on. In that case, @p unpack will 2823 * also not be called on the ghost cell that corresponds to it on the 2824 * receiving side. On the other hand, if the std_cxx17::optional object is 2825 * not empty, then the data stored within it will be sent to the received 2826 * and the @p unpack function called with it. 2827 * 2828 * @tparam DataType The type of the data to be communicated. It is assumed 2829 * to be serializable by boost::serialization. In many cases, this 2830 * data type can not be deduced by the compiler, e.g., if you provide 2831 * lambda functions for the second and third argument 2832 * to this function. In this case, you have to explicitly specify 2833 * the @p DataType as a template argument to the function call. 2834 * @tparam MeshType The type of @p mesh. 2835 * 2836 * @param mesh A variable of a type that satisfies the requirements of the 2837 * @ref ConceptMeshType "MeshType concept". 2838 * @param pack The function that will be called on each locally owned cell 2839 * that is a ghost cell somewhere else. As mentioned above, the function 2840 * may return a regular data object of type @p DataType to indicate 2841 * that data should be sent, or an empty 2842 * <code>std_cxx17::optional@<DataType@></code> to indicate that nothing has 2843 * to be sent for this cell. 2844 * @param unpack The function that will be called for each ghost cell 2845 * for which data was sent, i.e., for which the @p pack function 2846 * on the sending side returned a non-empty std_cxx17::optional object. 2847 * The @p unpack function is then called with the data sent by the 2848 * processor that owns that cell. 2849 * @param cell_filter Only cells are communicated where this filter function returns 2850 * the value `true`. In the default case, the function returns true on all 2851 * cells and thus, all relevant cells are communicated. 2852 * 2853 * <h4> An example </h4> 2854 * 2855 * Here is an example that shows how this function is to be used 2856 * in a concrete context. It is taken from the code that makes 2857 * sure that the @p active_fe_index (a single unsigned integer) is 2858 * transported from locally owned cells where one can set it in 2859 * hp::DoFHandler objects, to the corresponding ghost cells on 2860 * other processors to ensure that one can query the right value 2861 * also on those processors: 2862 * @code 2863 * using active_cell_iterator = 2864 * typename dealii::hp::DoFHandler<dim,spacedim>::active_cell_iterator; 2865 * auto pack = [] (const active_cell_iterator &cell) -> unsigned int 2866 * { 2867 * return cell->active_fe_index(); 2868 * }; 2869 * 2870 * auto unpack = [] (const active_cell_iterator &cell, 2871 * const unsigned int active_fe_index) -> void 2872 * { 2873 * cell->set_active_fe_index(active_fe_index); 2874 * }; 2875 * 2876 * GridTools::exchange_cell_data_to_ghosts< 2877 * unsigned int, dealii::hp::DoFHandler<dim,spacedim>> (dof_handler, 2878 * pack, 2879 * unpack); 2880 * @endcode 2881 * 2882 * You will notice that the @p pack lambda function returns an `unsigned int`, 2883 * not a `std_cxx17::optional<unsigned int>`. The former converts 2884 * automatically to the latter, implying that data will always be transported 2885 * to the other processor. 2886 * 2887 * (In reality, the @p unpack function needs to be a bit more 2888 * complicated because it is not allowed to call 2889 * DoFAccessor::set_active_fe_index() on ghost cells. Rather, the 2890 * @p unpack function directly accesses internal data structures. But 2891 * you get the idea -- the code could, just as well, have exchanged 2892 * material ids, user indices, boundary indicators, or any kind of other 2893 * data with similar calls as the ones above.) 2894 */ 2895 template <typename DataType, typename MeshType> 2896 void 2897 exchange_cell_data_to_ghosts( 2898 const MeshType & mesh, 2899 const std::function<std_cxx17::optional<DataType>( 2900 const typename MeshType::active_cell_iterator &)> &pack, 2901 const std::function<void(const typename MeshType::active_cell_iterator &, 2902 const DataType &)> & unpack, 2903 const std::function<bool(const typename MeshType::active_cell_iterator &)> 2904 &cell_filter = 2905 [](const typename MeshType::active_cell_iterator &) { return true; }); 2906 2907 /** 2908 * Exchange arbitrary data of type @p DataType provided by the function 2909 * objects from locally owned level cells to ghost level cells on other 2910 * processes. 2911 * 2912 * In addition to the parameters of exchange_cell_data_to_ghosts(), this 2913 * function allows to provide a @p cell_filter function, which can be used to only 2914 * communicate marked cells. In the default case, all relevant cells are 2915 * communicated. 2916 */ 2917 template <typename DataType, typename MeshType> 2918 void 2919 exchange_cell_data_to_level_ghosts( 2920 const MeshType & mesh, 2921 const std::function<std_cxx17::optional<DataType>( 2922 const typename MeshType::level_cell_iterator &)> &pack, 2923 const std::function<void(const typename MeshType::level_cell_iterator &, 2924 const DataType &)> & unpack, 2925 const std::function<bool(const typename MeshType::level_cell_iterator &)> 2926 &cell_filter = 2927 [](const typename MeshType::level_cell_iterator &) { return true; }); 2928 2929 /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding 2930 * boxes @p local_bboxes. 2931 * 2932 * This function is meant to exchange bounding boxes describing the locally 2933 * owned cells in a distributed triangulation obtained with the function 2934 * GridTools::compute_mesh_predicate_bounding_box . 2935 * 2936 * The output vector's size is the number of processes of the MPI 2937 * communicator: 2938 * its i-th entry contains the vector @p local_bboxes of the i-th process. 2939 */ 2940 template <int spacedim> 2941 std::vector<std::vector<BoundingBox<spacedim>>> 2942 exchange_local_bounding_boxes( 2943 const std::vector<BoundingBox<spacedim>> &local_bboxes, 2944 MPI_Comm mpi_communicator); 2945 2946 /** 2947 * In this collective operation each process provides a vector 2948 * of bounding boxes and a communicator. 2949 * All these vectors are gathered on each of the processes, 2950 * organized in a search tree which, and then returned. 2951 * 2952 * The idea is that the vector of bounding boxes describes a 2953 * relevant property of the computations on each process 2954 * individually, which could also be of use to other processes. An 2955 * example would be if the input vector of bounding boxes 2956 * corresponded to a covering of the locally owned partition of a 2957 * mesh (see 2958 * @ref GlossLocallyOwnedCell) 2959 * of a 2960 * parallel::distributed::Triangulation object. While these may 2961 * overlap the bounding boxes of other processes, finding which 2962 * process owns the cell that encloses a given point is vastly 2963 * easier if the process trying to figure this out has a list of 2964 * bounding boxes for each of the other processes at hand. 2965 * 2966 * The returned search tree object is an r-tree with packing 2967 * algorithm, which is provided by boost library. See 2968 * https://www.boost.org/doc/libs/1_67_0/libs/geometry/doc/html/geometry/spatial_indexes/introduction.html 2969 * for more information. 2970 * 2971 * In the returned tree, each node contains a pair of elements: 2972 * the first being a bounding box, 2973 * the second being the rank of the process whose local description 2974 * contains the bounding box. 2975 * 2976 * @note This function is a collective operation. 2977 */ 2978 template <int spacedim> 2979 RTree<std::pair<BoundingBox<spacedim>, unsigned int>> 2980 build_global_description_tree( 2981 const std::vector<BoundingBox<spacedim>> &local_description, 2982 MPI_Comm mpi_communicator); 2983 2984 /** 2985 * Collect for a given triangulation all locally relevant vertices that 2986 * coincide due to periodicity. 2987 * 2988 * Coinciding vertices are put into a group, e.g.: [1, 25, 51], which is 2989 * labeled by an arbitrary element from it, e.g.: "1". All coinciding vertices 2990 * store the label to its group, so that they can quickly access all the 2991 * coinciding vertices in that group: e.g.: 51 -> "1" -> [1, 25, 51] 2992 * 2993 * @param[in] tria Triangulation. 2994 * @param[out] coinciding_vertex_groups A map of equivalence classes (of 2995 * coinciding vertices) labeled by an arbitrary element from them. 2996 * Vertices not coinciding are ignored. 2997 * @param[out] vertex_to_coinciding_vertex_group Map of a vertex to the label 2998 * of a group of coinciding vertices. Vertices not contained in 2999 * this vector are not coinciding with any other vertex. 3000 */ 3001 template <int dim, int spacedim> 3002 void 3003 collect_coinciding_vertices( 3004 const Triangulation<dim, spacedim> & tria, 3005 std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups, 3006 std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group); 3007 3008 /** 3009 * Return a map that, for each vertex, lists all the processes whose 3010 * subdomains are adjacent to that vertex. 3011 * 3012 * @param[in] tria Triangulation. 3013 */ 3014 template <int dim, int spacedim> 3015 std::map<unsigned int, std::set<dealii::types::subdomain_id>> 3016 compute_vertices_with_ghost_neighbors( 3017 const Triangulation<dim, spacedim> &tria); 3018 3019 /** 3020 * A structure that allows the transfer of cell data of type @p T from one processor 3021 * to another. It corresponds to a packed buffer that stores a vector of 3022 * CellId and a vector of type @p T. 3023 * 3024 * This class facilitates the transfer by providing the save/load functions 3025 * that are able to pack up the vector of CellId's and the associated 3026 * data of type @p T into a stream. 3027 * 3028 * Type @p T is assumed to be serializable by <code>boost::serialization</code> (for 3029 * example <code>unsigned int</code> or <code>std::vector@<double@></code>). 3030 */ 3031 template <int dim, typename T> 3032 struct CellDataTransferBuffer 3033 { 3034 /** 3035 * A vector to store IDs of cells to be transferred. 3036 */ 3037 std::vector<CellId> cell_ids; 3038 3039 /** 3040 * A vector of cell data to be transferred. 3041 */ 3042 std::vector<T> data; 3043 3044 /** 3045 * Write the data of this object to a stream for the purpose of 3046 * serialization. 3047 * 3048 * @pre The user is responsible to keep the size of @p data 3049 * equal to the size as @p cell_ids . 3050 */ 3051 template <class Archive> 3052 void 3053 save(Archive &ar, const unsigned int version) const; 3054 3055 /** 3056 * Read the data of this object from a stream for the purpose of 3057 * serialization. Throw away the previous content. 3058 */ 3059 template <class Archive> 3060 void 3061 load(Archive &ar, const unsigned int version); 3062 3063 # ifdef DOXYGEN 3064 /** 3065 * Write and read the data of this object from a stream for the purpose 3066 * of serialization. 3067 */ 3068 template <class Archive> 3069 void 3070 serialize(Archive &archive, const unsigned int version); 3071 # else 3072 // This macro defines the serialize() method that is compatible with 3073 // the templated save() and load() method that have been implemented. 3074 BOOST_SERIALIZATION_SPLIT_MEMBER() 3075 # endif 3076 }; 3077 3078 /** 3079 * @name Exceptions 3080 */ 3081 /*@{*/ 3082 3083 /** 3084 * Exception 3085 */ 3086 DeclException1(ExcInvalidNumberOfPartitions, 3087 int, 3088 << "The number of partitions you gave is " << arg1 3089 << ", but must be greater than zero."); 3090 /** 3091 * Exception 3092 */ 3093 DeclException1(ExcNonExistentSubdomain, 3094 int, 3095 << "The subdomain id " << arg1 3096 << " has no cells associated with it."); 3097 /** 3098 * Exception 3099 */ 3100 DeclException0(ExcTriangulationHasBeenRefined); 3101 3102 /** 3103 * Exception 3104 */ 3105 DeclException1(ExcScalingFactorNotPositive, 3106 double, 3107 << "The scaling factor must be positive, but it is " << arg1 3108 << "."); 3109 /** 3110 * Exception 3111 */ 3112 template <int N> 3113 DeclException1(ExcPointNotFoundInCoarseGrid, 3114 Point<N>, 3115 << "The point <" << arg1 3116 << "> could not be found inside any of the " 3117 << "coarse grid cells."); 3118 /** 3119 * Exception 3120 */ 3121 template <int N> 3122 DeclException1(ExcPointNotFound, 3123 Point<N>, 3124 << "The point <" << arg1 3125 << "> could not be found inside any of the " 3126 << "subcells of a coarse grid cell."); 3127 3128 /** 3129 * Exception 3130 */ 3131 DeclException1(ExcVertexNotUsed, 3132 unsigned int, 3133 << "The given vertex with index " << arg1 3134 << " is not used in the given triangulation."); 3135 3136 3137 /*@}*/ 3138 3139 } /*namespace GridTools*/ 3140 3141 3142 3143 /* ----------------- Template function --------------- */ 3144 3145 # ifndef DOXYGEN 3146 3147 namespace GridTools 3148 { 3149 template <int dim> 3150 double cell_measure(const std::vector<Point<dim>> & all_vertices,const unsigned int (& indices)[GeometryInfo<dim>::vertices_per_cell])3151 cell_measure( 3152 const std::vector<Point<dim>> &all_vertices, 3153 const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell]) 3154 { 3155 // We forward call to the ArrayView version: 3156 const ArrayView<const unsigned int> view( 3157 indices, GeometryInfo<dim>::vertices_per_cell); 3158 return cell_measure(all_vertices, view); 3159 } 3160 3161 template <int dim, typename T> 3162 double cell_measure(const T &,...)3163 cell_measure(const T &, ...) 3164 { 3165 Assert(false, ExcNotImplemented()); 3166 return std::numeric_limits<double>::quiet_NaN(); 3167 } 3168 3169 3170 3171 // This specialization is defined here so that the general template in the 3172 // source file doesn't need to have further 1D overloads for the internal 3173 // functions it calls. 3174 template <> 3175 inline Triangulation<1, 1>::DistortedCellList fix_up_distorted_child_cells(const Triangulation<1,1>::DistortedCellList &,Triangulation<1,1> &)3176 fix_up_distorted_child_cells(const Triangulation<1, 1>::DistortedCellList &, 3177 Triangulation<1, 1> &) 3178 { 3179 return {}; 3180 } 3181 3182 3183 3184 template <int dim, typename Predicate, int spacedim> 3185 void transform(const Predicate & predicate,Triangulation<dim,spacedim> & triangulation)3186 transform(const Predicate & predicate, 3187 Triangulation<dim, spacedim> &triangulation) 3188 { 3189 std::vector<bool> treated_vertices(triangulation.n_vertices(), false); 3190 3191 // loop over all active cells, and 3192 // transform those vertices that 3193 // have not yet been touched. note 3194 // that we get to all vertices in 3195 // the triangulation by only 3196 // visiting the active cells. 3197 typename Triangulation<dim, spacedim>::active_cell_iterator 3198 cell = triangulation.begin_active(), 3199 endc = triangulation.end(); 3200 for (; cell != endc; ++cell) 3201 for (const unsigned int v : cell->vertex_indices()) 3202 if (treated_vertices[cell->vertex_index(v)] == false) 3203 { 3204 // transform this vertex 3205 cell->vertex(v) = predicate(cell->vertex(v)); 3206 // and mark it as treated 3207 treated_vertices[cell->vertex_index(v)] = true; 3208 }; 3209 3210 3211 // now fix any vertices on hanging nodes so that we don't create any holes 3212 if (dim == 2) 3213 { 3214 typename Triangulation<dim, spacedim>::active_cell_iterator 3215 cell = triangulation.begin_active(), 3216 endc = triangulation.end(); 3217 for (; cell != endc; ++cell) 3218 for (const unsigned int face : cell->face_indices()) 3219 if (cell->face(face)->has_children() && 3220 !cell->face(face)->at_boundary()) 3221 { 3222 Assert(cell->reference_cell_type() == 3223 ReferenceCell::get_hypercube(dim), 3224 ExcNotImplemented()); 3225 3226 // this line has children 3227 cell->face(face)->child(0)->vertex(1) = 3228 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) / 3229 2; 3230 } 3231 } 3232 else if (dim == 3) 3233 { 3234 typename Triangulation<dim, spacedim>::active_cell_iterator 3235 cell = triangulation.begin_active(), 3236 endc = triangulation.end(); 3237 for (; cell != endc; ++cell) 3238 for (const unsigned int face : cell->face_indices()) 3239 if (cell->face(face)->has_children() && 3240 !cell->face(face)->at_boundary()) 3241 { 3242 Assert(cell->reference_cell_type() == 3243 ReferenceCell::get_hypercube(dim), 3244 ExcNotImplemented()); 3245 3246 // this face has hanging nodes 3247 cell->face(face)->child(0)->vertex(1) = 3248 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) / 3249 2.0; 3250 cell->face(face)->child(0)->vertex(2) = 3251 (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) / 3252 2.0; 3253 cell->face(face)->child(1)->vertex(3) = 3254 (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) / 3255 2.0; 3256 cell->face(face)->child(2)->vertex(3) = 3257 (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 3258 2.0; 3259 3260 // center of the face 3261 cell->face(face)->child(0)->vertex(3) = 3262 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) + 3263 cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 3264 4.0; 3265 } 3266 } 3267 3268 // Make sure FEValues notices that the mesh has changed 3269 triangulation.signals.mesh_movement(); 3270 } 3271 3272 3273 3274 template <class MeshType> 3275 std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(const typename MeshType::cell_iterator & cell)3276 get_active_child_cells(const typename MeshType::cell_iterator &cell) 3277 { 3278 std::vector<typename MeshType::active_cell_iterator> child_cells; 3279 3280 if (cell->has_children()) 3281 { 3282 for (unsigned int child = 0; child < cell->n_children(); ++child) 3283 if (cell->child(child)->has_children()) 3284 { 3285 const std::vector<typename MeshType::active_cell_iterator> 3286 children = get_active_child_cells<MeshType>(cell->child(child)); 3287 child_cells.insert(child_cells.end(), 3288 children.begin(), 3289 children.end()); 3290 } 3291 else 3292 child_cells.push_back(cell->child(child)); 3293 } 3294 3295 return child_cells; 3296 } 3297 3298 3299 3300 template <class MeshType> 3301 void get_active_neighbors(const typename MeshType::active_cell_iterator & cell,std::vector<typename MeshType::active_cell_iterator> & active_neighbors)3302 get_active_neighbors( 3303 const typename MeshType::active_cell_iterator & cell, 3304 std::vector<typename MeshType::active_cell_iterator> &active_neighbors) 3305 { 3306 active_neighbors.clear(); 3307 for (const unsigned int n : cell->face_indices()) 3308 if (!cell->at_boundary(n)) 3309 { 3310 if (MeshType::dimension == 1) 3311 { 3312 // check children of neighbor. note 3313 // that in 1d children of the neighbor 3314 // may be further refined. In 1d the 3315 // case is simple since we know what 3316 // children bound to the present cell 3317 typename MeshType::cell_iterator neighbor_child = 3318 cell->neighbor(n); 3319 if (!neighbor_child->is_active()) 3320 { 3321 while (neighbor_child->has_children()) 3322 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0); 3323 3324 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell, 3325 ExcInternalError()); 3326 } 3327 active_neighbors.push_back(neighbor_child); 3328 } 3329 else 3330 { 3331 if (cell->face(n)->has_children()) 3332 // this neighbor has children. find 3333 // out which border to the present 3334 // cell 3335 for (unsigned int c = 0; 3336 c < cell->face(n)->number_of_children(); 3337 ++c) 3338 active_neighbors.push_back( 3339 cell->neighbor_child_on_subface(n, c)); 3340 else 3341 { 3342 // the neighbor must be active 3343 // himself 3344 Assert(cell->neighbor(n)->is_active(), ExcInternalError()); 3345 active_neighbors.push_back(cell->neighbor(n)); 3346 } 3347 } 3348 } 3349 } 3350 3351 3352 3353 namespace internal 3354 { 3355 namespace ProjectToObject 3356 { 3357 /** 3358 * The method GridTools::project_to_object requires taking derivatives 3359 * along the surface of a simplex. In general these cannot be 3360 * approximated with finite differences but special differences of the 3361 * form 3362 * 3363 * df/dx_i - df/dx_j 3364 * 3365 * <em>can</em> be approximated. This <code>struct</code> just stores 3366 * the two derivatives approximated by the stencil (in the case of the 3367 * example above <code>i</code> and <code>j</code>). 3368 */ 3369 struct CrossDerivative 3370 { 3371 const unsigned int direction_0; 3372 const unsigned int direction_1; 3373 3374 CrossDerivative(const unsigned int d0, const unsigned int d1); 3375 }; 3376 CrossDerivative(const unsigned int d0,const unsigned int d1)3377 inline CrossDerivative::CrossDerivative(const unsigned int d0, 3378 const unsigned int d1) 3379 : direction_0(d0) 3380 , direction_1(d1) 3381 {} 3382 3383 3384 3385 /** 3386 * Standard second-order approximation to the first derivative with a 3387 * two-point centered scheme. This is used below in a 1D Newton method. 3388 */ 3389 template <typename F> 3390 inline auto 3391 centered_first_difference(const double center, 3392 const double step, 3393 const F &f) -> decltype(f(center) - f(center)) 3394 { 3395 return (f(center + step) - f(center - step)) / (2.0 * step); 3396 } 3397 3398 3399 3400 /** 3401 * Standard second-order approximation to the second derivative with a 3402 * three-point centered scheme. This is used below in a 1D Newton method. 3403 */ 3404 template <typename F> 3405 inline auto 3406 centered_second_difference(const double center, 3407 const double step, 3408 const F &f) -> decltype(f(center) - f(center)) 3409 { 3410 return (f(center + step) - 2.0 * f(center) + f(center - step)) / 3411 (step * step); 3412 } 3413 3414 3415 3416 /** 3417 * Fourth order approximation of the derivative 3418 * 3419 * df/dx_i - df/dx_j 3420 * 3421 * where <code>i</code> and <code>j</code> are specified by @p 3422 * cross_derivative. The derivative approximation is at @p center with a 3423 * step size of @p step and function @p f. 3424 */ 3425 template <int structdim, typename F> 3426 inline auto 3427 cross_stencil( 3428 const CrossDerivative cross_derivative, 3429 const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> ¢er, 3430 const double step, 3431 const F &f) -> decltype(f(center) - f(center)) 3432 { 3433 Tensor<1, GeometryInfo<structdim>::vertices_per_cell> simplex_vector; 3434 simplex_vector[cross_derivative.direction_0] = 0.5 * step; 3435 simplex_vector[cross_derivative.direction_1] = -0.5 * step; 3436 return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) - 3437 1.0 / 3.0 * f(center - simplex_vector) + 3438 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) / 3439 step; 3440 } 3441 3442 3443 3444 /** 3445 * The optimization algorithm used in GridTools::project_to_object is 3446 * essentially a gradient descent method. This function computes entries 3447 * in the gradient of the objective function; see the description in the 3448 * comments inside GridTools::project_to_object for more information. 3449 */ 3450 template <int spacedim, int structdim, typename F> 3451 inline double gradient_entry(const unsigned int row_n,const unsigned int dependent_direction,const Point<spacedim> & p0,const Tensor<1,GeometryInfo<structdim>::vertices_per_cell> & center,const double step,const F & f)3452 gradient_entry( 3453 const unsigned int row_n, 3454 const unsigned int dependent_direction, 3455 const Point<spacedim> &p0, 3456 const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> ¢er, 3457 const double step, 3458 const F & f) 3459 { 3460 Assert(row_n < GeometryInfo<structdim>::vertices_per_cell && 3461 dependent_direction < 3462 GeometryInfo<structdim>::vertices_per_cell, 3463 ExcMessage("This function assumes that the last weight is a " 3464 "dependent variable (and hence we cannot take its " 3465 "derivative directly).")); 3466 Assert(row_n != dependent_direction, 3467 ExcMessage( 3468 "We cannot differentiate with respect to the variable " 3469 "that is assumed to be dependent.")); 3470 3471 const Point<spacedim> manifold_point = f(center); 3472 const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>( 3473 {row_n, dependent_direction}, center, step, f); 3474 double entry = 0.0; 3475 for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n) 3476 entry += 3477 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n]; 3478 return entry; 3479 } 3480 3481 /** 3482 * Project onto a d-linear object. This is more accurate than the 3483 * general algorithm in project_to_object but only works for geometries 3484 * described by linear, bilinear, or trilinear mappings. 3485 */ 3486 template <typename Iterator, int spacedim, int structdim> 3487 Point<spacedim> project_to_d_linear_object(const Iterator & object,const Point<spacedim> & trial_point)3488 project_to_d_linear_object(const Iterator & object, 3489 const Point<spacedim> &trial_point) 3490 { 3491 // let's look at this for simplicity for a quad (structdim==2) in a 3492 // space with spacedim>2 (notate trial_point by y): all points on the 3493 // surface are given by 3494 // x(\xi) = sum_i v_i phi_x(\xi) 3495 // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the 3496 // reference coordinates of the quad. so what we are trying to do is 3497 // find a point x on the surface that is closest to the point y. there 3498 // are different ways to solve this problem, but in the end it's a 3499 // nonlinear problem and we have to find reference coordinates \xi so 3500 // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function 3501 // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree 3502 // 2*structdim that we'd like to minimize. unless structdim==1, we'll 3503 // have to use a Newton method to find the answer. This leads to the 3504 // following formulation of Newton steps: 3505 // 3506 // Given \xi_k, find \delta\xi_k so that 3507 // H_k \delta\xi_k = - F_k 3508 // where H_k is an approximation to the second derivatives of J at 3509 // \xi_k, and F_k is the first derivative of J. We'll iterate this a 3510 // number of times until the right hand side is small enough. As a 3511 // stopping criterion, we terminate if ||\delta\xi||<eps. 3512 // 3513 // As for the Hessian, the best choice would be 3514 // H_k = J''(\xi_k) 3515 // but we'll opt for the simpler Gauss-Newton form 3516 // H_k = A^T A 3517 // i.e. 3518 // (H_k)_{nm} = \sum_{i,j} v_i*v_j * 3519 // \partial_n phi_i * 3520 // \partial_m phi_j 3521 // we start at xi=(0.5, 0.5). 3522 Point<structdim> xi; 3523 for (unsigned int d = 0; d < structdim; ++d) 3524 xi[d] = 0.5; 3525 3526 Point<spacedim> x_k; 3527 for (const unsigned int i : GeometryInfo<structdim>::vertex_indices()) 3528 x_k += object->vertex(i) * 3529 GeometryInfo<structdim>::d_linear_shape_function(xi, i); 3530 3531 do 3532 { 3533 Tensor<1, structdim> F_k; 3534 for (const unsigned int i : 3535 GeometryInfo<structdim>::vertex_indices()) 3536 F_k += 3537 (x_k - trial_point) * object->vertex(i) * 3538 GeometryInfo<structdim>::d_linear_shape_function_gradient(xi, 3539 i); 3540 3541 Tensor<2, structdim> H_k; 3542 for (const unsigned int i : 3543 GeometryInfo<structdim>::vertex_indices()) 3544 for (const unsigned int j : 3545 GeometryInfo<structdim>::vertex_indices()) 3546 { 3547 Tensor<2, structdim> tmp = outer_product( 3548 GeometryInfo<structdim>::d_linear_shape_function_gradient( 3549 xi, i), 3550 GeometryInfo<structdim>::d_linear_shape_function_gradient( 3551 xi, j)); 3552 H_k += (object->vertex(i) * object->vertex(j)) * tmp; 3553 } 3554 3555 const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k; 3556 xi += delta_xi; 3557 3558 x_k = Point<spacedim>(); 3559 for (const unsigned int i : 3560 GeometryInfo<structdim>::vertex_indices()) 3561 x_k += object->vertex(i) * 3562 GeometryInfo<structdim>::d_linear_shape_function(xi, i); 3563 3564 if (delta_xi.norm() < 1e-7) 3565 break; 3566 } 3567 while (true); 3568 3569 return x_k; 3570 } 3571 } // namespace ProjectToObject 3572 } // namespace internal 3573 3574 3575 3576 namespace internal 3577 { 3578 // We hit an internal compiler error in ICC 15 if we define this as a lambda 3579 // inside the project_to_object function below. 3580 template <int structdim> 3581 inline bool weights_are_ok(const Tensor<1,GeometryInfo<structdim>::vertices_per_cell> & v)3582 weights_are_ok( 3583 const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v) 3584 { 3585 // clang has trouble figuring out structdim here, so define it 3586 // again: 3587 static const std::size_t n_vertices_per_cell = 3588 Tensor<1, GeometryInfo<structdim>::vertices_per_cell>:: 3589 n_independent_components; 3590 std::array<double, n_vertices_per_cell> copied_weights; 3591 for (unsigned int i = 0; i < n_vertices_per_cell; ++i) 3592 { 3593 copied_weights[i] = v[i]; 3594 if (v[i] < 0.0 || v[i] > 1.0) 3595 return false; 3596 } 3597 3598 // check the sum: try to avoid some roundoff errors by summing in order 3599 std::sort(copied_weights.begin(), copied_weights.end()); 3600 const double sum = 3601 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0); 3602 return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc 3603 } 3604 } // namespace internal 3605 3606 template <typename Iterator> 3607 Point<Iterator::AccessorType::space_dimension> project_to_object(const Iterator & object,const Point<Iterator::AccessorType::space_dimension> & trial_point)3608 project_to_object( 3609 const Iterator & object, 3610 const Point<Iterator::AccessorType::space_dimension> &trial_point) 3611 { 3612 const int spacedim = Iterator::AccessorType::space_dimension; 3613 const int structdim = Iterator::AccessorType::structure_dimension; 3614 3615 Point<spacedim> projected_point = trial_point; 3616 3617 if (structdim >= spacedim) 3618 return projected_point; 3619 else if (structdim == 1 || structdim == 2) 3620 { 3621 using namespace internal::ProjectToObject; 3622 // Try to use the special flat algorithm for quads (this is better 3623 // than the general algorithm in 3D). This does not take into account 3624 // whether projected_point is outside the quad, but we optimize along 3625 // lines below anyway: 3626 const int dim = Iterator::AccessorType::dimension; 3627 const Manifold<dim, spacedim> &manifold = object->get_manifold(); 3628 if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>( 3629 &manifold) != nullptr) 3630 { 3631 projected_point = 3632 project_to_d_linear_object<Iterator, spacedim, structdim>( 3633 object, trial_point); 3634 } 3635 else 3636 { 3637 // We want to find a point on the convex hull (defined by the 3638 // vertices of the object and the manifold description) that is 3639 // relatively close to the trial point. This has a few issues: 3640 // 3641 // 1. For a general convex hull we are not guaranteed that a unique 3642 // minimum exists. 3643 // 2. The independent variables in the optimization process are the 3644 // weights given to Manifold::get_new_point, which must sum to 1, 3645 // so we cannot use standard finite differences to approximate a 3646 // gradient. 3647 // 3648 // There is not much we can do about 1., but for 2. we can derive 3649 // finite difference stencils that work on a structdim-dimensional 3650 // simplex and rewrite the optimization problem to use those 3651 // instead. Consider the structdim 2 case and let 3652 // 3653 // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1, 3654 // c2, c3}) 3655 // 3656 // where {c0, c1, c2, c3} are the weights for the four vertices on 3657 // the quadrilateral. We seek to minimize the Euclidean distance 3658 // between F(...) and trial_point. We can solve for c3 in terms of 3659 // the other weights and get, for one coordinate direction 3660 // 3661 // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2) 3662 // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...)) 3663 // 3664 // where we substitute back in for c3 after taking the 3665 // derivative. We can compute a stencil for the cross derivative 3666 // d/dc0 - d/dc3: this is exactly what cross_stencil approximates 3667 // (and gradient_entry computes the sum over the independent 3668 // variables). Below, we somewhat arbitrarily pick the last 3669 // component as the dependent one. 3670 // 3671 // Since we can now calculate derivatives of the objective 3672 // function we can use gradient descent to minimize it. 3673 // 3674 // Of course, this is much simpler in the structdim = 1 case (we 3675 // could rewrite the projection as a 1D optimization problem), but 3676 // to reduce the potential for bugs we use the same code in both 3677 // cases. 3678 const double step_size = object->diameter() / 64.0; 3679 3680 constexpr unsigned int n_vertices_per_cell = 3681 GeometryInfo<structdim>::vertices_per_cell; 3682 3683 std::array<Point<spacedim>, n_vertices_per_cell> vertices; 3684 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell; 3685 ++vertex_n) 3686 vertices[vertex_n] = object->vertex(vertex_n); 3687 3688 auto get_point_from_weights = 3689 [&](const Tensor<1, n_vertices_per_cell> &weights) 3690 -> Point<spacedim> { 3691 return object->get_manifold().get_new_point( 3692 make_array_view(vertices.begin(), vertices.end()), 3693 make_array_view(weights.begin_raw(), weights.end_raw())); 3694 }; 3695 3696 // pick the initial weights as (normalized) inverse distances from 3697 // the trial point: 3698 Tensor<1, n_vertices_per_cell> guess_weights; 3699 double guess_weights_sum = 0.0; 3700 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell; 3701 ++vertex_n) 3702 { 3703 const double distance = 3704 vertices[vertex_n].distance(trial_point); 3705 if (distance == 0.0) 3706 { 3707 guess_weights = 0.0; 3708 guess_weights[vertex_n] = 1.0; 3709 guess_weights_sum = 1.0; 3710 break; 3711 } 3712 else 3713 { 3714 guess_weights[vertex_n] = 1.0 / distance; 3715 guess_weights_sum += guess_weights[vertex_n]; 3716 } 3717 } 3718 guess_weights /= guess_weights_sum; 3719 Assert(internal::weights_are_ok<structdim>(guess_weights), 3720 ExcInternalError()); 3721 3722 // The optimization algorithm consists of two parts: 3723 // 3724 // 1. An outer loop where we apply the gradient descent algorithm. 3725 // 2. An inner loop where we do a line search to find the optimal 3726 // length of the step one should take in the gradient direction. 3727 // 3728 for (unsigned int outer_n = 0; outer_n < 40; ++outer_n) 3729 { 3730 const unsigned int dependent_direction = 3731 n_vertices_per_cell - 1; 3732 Tensor<1, n_vertices_per_cell> current_gradient; 3733 for (unsigned int row_n = 0; row_n < n_vertices_per_cell; 3734 ++row_n) 3735 { 3736 if (row_n != dependent_direction) 3737 { 3738 current_gradient[row_n] = 3739 gradient_entry<spacedim, structdim>( 3740 row_n, 3741 dependent_direction, 3742 trial_point, 3743 guess_weights, 3744 step_size, 3745 get_point_from_weights); 3746 3747 current_gradient[dependent_direction] -= 3748 current_gradient[row_n]; 3749 } 3750 } 3751 3752 // We need to travel in the -gradient direction, as noted 3753 // above, but we may not want to take a full step in that 3754 // direction; instead, guess that we will go -0.5*gradient and 3755 // do quasi-Newton iteration to pick the best multiplier. The 3756 // goal is to find a scalar alpha such that 3757 // 3758 // F(x - alpha g) 3759 // 3760 // is minimized, where g is the gradient and F is the 3761 // objective function. To find the optimal value we find roots 3762 // of the derivative of the objective function with respect to 3763 // alpha by Newton iteration, where we approximate the first 3764 // and second derivatives of F(x - alpha g) with centered 3765 // finite differences. 3766 double gradient_weight = -0.5; 3767 auto gradient_weight_objective_function = 3768 [&](const double gradient_weight_guess) -> double { 3769 return (trial_point - 3770 get_point_from_weights(guess_weights + 3771 gradient_weight_guess * 3772 current_gradient)) 3773 .norm_square(); 3774 }; 3775 3776 for (unsigned int inner_n = 0; inner_n < 10; ++inner_n) 3777 { 3778 const double update_numerator = centered_first_difference( 3779 gradient_weight, 3780 step_size, 3781 gradient_weight_objective_function); 3782 const double update_denominator = 3783 centered_second_difference( 3784 gradient_weight, 3785 step_size, 3786 gradient_weight_objective_function); 3787 3788 // avoid division by zero. Note that we limit the gradient 3789 // weight below 3790 if (std::abs(update_denominator) == 0.0) 3791 break; 3792 gradient_weight = 3793 gradient_weight - update_numerator / update_denominator; 3794 3795 // Put a fairly lenient bound on the largest possible 3796 // gradient (things tend to be locally flat, so the gradient 3797 // itself is usually small) 3798 if (std::abs(gradient_weight) > 10) 3799 { 3800 gradient_weight = -10.0; 3801 break; 3802 } 3803 } 3804 3805 // It only makes sense to take convex combinations with weights 3806 // between zero and one. If the update takes us outside of this 3807 // region then rescale the update to stay within the region and 3808 // try again 3809 Tensor<1, n_vertices_per_cell> tentative_weights = 3810 guess_weights + gradient_weight * current_gradient; 3811 3812 double new_gradient_weight = gradient_weight; 3813 for (unsigned int iteration_count = 0; iteration_count < 40; 3814 ++iteration_count) 3815 { 3816 if (internal::weights_are_ok<structdim>(tentative_weights)) 3817 break; 3818 3819 for (unsigned int i = 0; i < n_vertices_per_cell; ++i) 3820 { 3821 if (tentative_weights[i] < 0.0) 3822 { 3823 tentative_weights -= 3824 (tentative_weights[i] / current_gradient[i]) * 3825 current_gradient; 3826 } 3827 if (tentative_weights[i] < 0.0 || 3828 1.0 < tentative_weights[i]) 3829 { 3830 new_gradient_weight /= 2.0; 3831 tentative_weights = 3832 guess_weights + 3833 new_gradient_weight * current_gradient; 3834 } 3835 } 3836 } 3837 3838 // the update might still send us outside the valid region, so 3839 // check again and quit if the update is still not valid 3840 if (!internal::weights_are_ok<structdim>(tentative_weights)) 3841 break; 3842 3843 // if we cannot get closer by traveling in the gradient 3844 // direction then quit 3845 if (get_point_from_weights(tentative_weights) 3846 .distance(trial_point) < 3847 get_point_from_weights(guess_weights).distance(trial_point)) 3848 guess_weights = tentative_weights; 3849 else 3850 break; 3851 Assert(internal::weights_are_ok<structdim>(guess_weights), 3852 ExcInternalError()); 3853 } 3854 Assert(internal::weights_are_ok<structdim>(guess_weights), 3855 ExcInternalError()); 3856 projected_point = get_point_from_weights(guess_weights); 3857 } 3858 3859 // if structdim == 2 and the optimal point is not on the interior then 3860 // we may be able to get a more accurate result by projecting onto the 3861 // lines. 3862 if (structdim == 2) 3863 { 3864 std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell> 3865 line_projections; 3866 for (unsigned int line_n = 0; 3867 line_n < GeometryInfo<structdim>::lines_per_cell; 3868 ++line_n) 3869 { 3870 line_projections[line_n] = 3871 project_to_object(object->line(line_n), trial_point); 3872 } 3873 std::sort(line_projections.begin(), 3874 line_projections.end(), 3875 [&](const Point<spacedim> &a, const Point<spacedim> &b) { 3876 return a.distance(trial_point) < 3877 b.distance(trial_point); 3878 }); 3879 if (line_projections[0].distance(trial_point) < 3880 projected_point.distance(trial_point)) 3881 projected_point = line_projections[0]; 3882 } 3883 } 3884 else 3885 { 3886 Assert(false, ExcNotImplemented()); 3887 return projected_point; 3888 } 3889 3890 return projected_point; 3891 } 3892 3893 3894 3895 template <int dim, typename T> 3896 template <class Archive> 3897 void save(Archive & ar,const unsigned int)3898 CellDataTransferBuffer<dim, T>::save(Archive &ar, 3899 const unsigned int /*version*/) const 3900 { 3901 Assert(cell_ids.size() == data.size(), 3902 ExcDimensionMismatch(cell_ids.size(), data.size())); 3903 // archive the cellids in an efficient binary format 3904 const std::size_t n_cells = cell_ids.size(); 3905 ar & n_cells; 3906 for (const auto &id : cell_ids) 3907 { 3908 CellId::binary_type binary_cell_id = id.template to_binary<dim>(); 3909 ar & binary_cell_id; 3910 } 3911 3912 ar &data; 3913 } 3914 3915 3916 3917 template <int dim, typename T> 3918 template <class Archive> 3919 void load(Archive & ar,const unsigned int)3920 CellDataTransferBuffer<dim, T>::load(Archive &ar, 3921 const unsigned int /*version*/) 3922 { 3923 std::size_t n_cells; 3924 ar & n_cells; 3925 cell_ids.clear(); 3926 cell_ids.reserve(n_cells); 3927 for (unsigned int c = 0; c < n_cells; ++c) 3928 { 3929 CellId::binary_type value; 3930 ar & value; 3931 cell_ids.emplace_back(value); 3932 } 3933 ar &data; 3934 } 3935 3936 3937 namespace internal 3938 { 3939 template <typename DataType, 3940 typename MeshType, 3941 typename MeshCellIteratorType> 3942 void exchange_cell_data(const MeshType & mesh,const std::function<std_cxx17::optional<DataType> (const MeshCellIteratorType &)> & pack,const std::function<void (const MeshCellIteratorType &,const DataType &)> & unpack,const std::function<bool (const MeshCellIteratorType &)> & cell_filter,const std::function<void (const std::function<void (const MeshCellIteratorType &,const types::subdomain_id)> &)> & process_cells,const std::function<std::set<types::subdomain_id> (const parallel::TriangulationBase<MeshType::dimension,MeshType::space_dimension> &)> & compute_ghost_owners)3943 exchange_cell_data( 3944 const MeshType &mesh, 3945 const std::function< 3946 std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack, 3947 const std::function<void(const MeshCellIteratorType &, const DataType &)> 3948 & unpack, 3949 const std::function<bool(const MeshCellIteratorType &)> & cell_filter, 3950 const std::function<void( 3951 const std::function<void(const MeshCellIteratorType &, 3952 const types::subdomain_id)> &)> &process_cells, 3953 const std::function<std::set<types::subdomain_id>( 3954 const parallel::TriangulationBase<MeshType::dimension, 3955 MeshType::space_dimension> &)> 3956 &compute_ghost_owners) 3957 { 3958 # ifndef DEAL_II_WITH_MPI 3959 (void)mesh; 3960 (void)pack; 3961 (void)unpack; 3962 (void)cell_filter; 3963 (void)process_cells; 3964 (void)compute_ghost_owners; 3965 Assert(false, 3966 ExcMessage("GridTools::exchange_cell_data() requires MPI.")); 3967 # else 3968 constexpr int dim = MeshType::dimension; 3969 constexpr int spacedim = MeshType::space_dimension; 3970 auto tria = 3971 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>( 3972 &mesh.get_triangulation()); 3973 Assert( 3974 tria != nullptr, 3975 ExcMessage( 3976 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations.")); 3977 3978 // build list of cells to request for each neighbor 3979 std::set<dealii::types::subdomain_id> ghost_owners = 3980 compute_ghost_owners(*tria); 3981 std::map<dealii::types::subdomain_id, 3982 std::vector<typename CellId::binary_type>> 3983 neighbor_cell_list; 3984 3985 for (const auto ghost_owner : ghost_owners) 3986 neighbor_cell_list[ghost_owner] = {}; 3987 3988 process_cells([&](const auto &cell, const auto key) { 3989 if (cell_filter(cell)) 3990 neighbor_cell_list[key].emplace_back( 3991 cell->id().template to_binary<spacedim>()); 3992 }); 3993 3994 Assert(ghost_owners.size() == neighbor_cell_list.size(), 3995 ExcInternalError()); 3996 3997 3998 // Before sending & receiving, make sure we protect this section with 3999 // a mutex: 4000 static Utilities::MPI::CollectiveMutex mutex; 4001 Utilities::MPI::CollectiveMutex::ScopedLock lock( 4002 mutex, tria->get_communicator()); 4003 4004 const int mpi_tag = 4005 Utilities::MPI::internal::Tags::exchange_cell_data_request; 4006 const int mpi_tag_reply = 4007 Utilities::MPI::internal::Tags::exchange_cell_data_reply; 4008 4009 // send our requests: 4010 std::vector<MPI_Request> requests(ghost_owners.size()); 4011 { 4012 unsigned int idx = 0; 4013 for (const auto &it : neighbor_cell_list) 4014 { 4015 // send the data about the relevant cells 4016 const int ierr = MPI_Isend(it.second.data(), 4017 it.second.size() * sizeof(it.second[0]), 4018 MPI_BYTE, 4019 it.first, 4020 mpi_tag, 4021 tria->get_communicator(), 4022 &requests[idx]); 4023 AssertThrowMPI(ierr); 4024 ++idx; 4025 } 4026 } 4027 4028 using DestinationToBufferMap = 4029 std::map<dealii::types::subdomain_id, 4030 GridTools::CellDataTransferBuffer<dim, DataType>>; 4031 DestinationToBufferMap destination_to_data_buffer_map; 4032 4033 // receive requests and reply with the ghost indices 4034 std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send( 4035 ghost_owners.size()); 4036 std::vector<std::vector<types::global_dof_index>> 4037 send_dof_numbers_and_indices(ghost_owners.size()); 4038 std::vector<MPI_Request> reply_requests(ghost_owners.size()); 4039 std::vector<std::vector<char>> sendbuffers(ghost_owners.size()); 4040 4041 for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx) 4042 { 4043 MPI_Status status; 4044 int ierr = MPI_Probe(MPI_ANY_SOURCE, 4045 mpi_tag, 4046 tria->get_communicator(), 4047 &status); 4048 AssertThrowMPI(ierr); 4049 4050 int len; 4051 ierr = MPI_Get_count(&status, MPI_BYTE, &len); 4052 AssertThrowMPI(ierr); 4053 Assert(len % sizeof(cell_data_to_send[idx][0]) == 0, 4054 ExcInternalError()); 4055 4056 const unsigned int n_cells = 4057 len / sizeof(typename CellId::binary_type); 4058 cell_data_to_send[idx].resize(n_cells); 4059 4060 ierr = MPI_Recv(cell_data_to_send[idx].data(), 4061 len, 4062 MPI_BYTE, 4063 status.MPI_SOURCE, 4064 status.MPI_TAG, 4065 tria->get_communicator(), 4066 &status); 4067 AssertThrowMPI(ierr); 4068 4069 // store data for each cell 4070 for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c) 4071 { 4072 const auto cell = 4073 CellId(cell_data_to_send[idx][c]).to_cell(*tria); 4074 4075 MeshCellIteratorType mesh_it(tria, 4076 cell->level(), 4077 cell->index(), 4078 &mesh); 4079 const std_cxx17::optional<DataType> data = pack(mesh_it); 4080 4081 if (data) 4082 { 4083 typename DestinationToBufferMap::iterator p = 4084 destination_to_data_buffer_map 4085 .insert(std::make_pair( 4086 idx, 4087 GridTools::CellDataTransferBuffer<dim, DataType>())) 4088 .first; 4089 4090 p->second.cell_ids.emplace_back(cell->id()); 4091 p->second.data.emplace_back(*data); 4092 } 4093 } 4094 4095 // send reply 4096 GridTools::CellDataTransferBuffer<dim, DataType> &data = 4097 destination_to_data_buffer_map[idx]; 4098 4099 sendbuffers[idx] = 4100 Utilities::pack(data, /*enable_compression*/ false); 4101 ierr = MPI_Isend(sendbuffers[idx].data(), 4102 sendbuffers[idx].size(), 4103 MPI_BYTE, 4104 status.MPI_SOURCE, 4105 mpi_tag_reply, 4106 tria->get_communicator(), 4107 &reply_requests[idx]); 4108 AssertThrowMPI(ierr); 4109 } 4110 4111 // finally receive the replies 4112 std::vector<char> receive; 4113 for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx) 4114 { 4115 MPI_Status status; 4116 int ierr = MPI_Probe(MPI_ANY_SOURCE, 4117 mpi_tag_reply, 4118 tria->get_communicator(), 4119 &status); 4120 AssertThrowMPI(ierr); 4121 4122 int len; 4123 ierr = MPI_Get_count(&status, MPI_BYTE, &len); 4124 AssertThrowMPI(ierr); 4125 4126 receive.resize(len); 4127 4128 char *ptr = receive.data(); 4129 ierr = MPI_Recv(ptr, 4130 len, 4131 MPI_BYTE, 4132 status.MPI_SOURCE, 4133 status.MPI_TAG, 4134 tria->get_communicator(), 4135 &status); 4136 AssertThrowMPI(ierr); 4137 4138 auto cellinfo = 4139 Utilities::unpack<CellDataTransferBuffer<dim, DataType>>( 4140 receive, /*enable_compression*/ false); 4141 4142 DataType *data = cellinfo.data.data(); 4143 for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data) 4144 { 4145 const typename Triangulation<dim, spacedim>::cell_iterator 4146 tria_cell = cellinfo.cell_ids[c].to_cell(*tria); 4147 4148 MeshCellIteratorType cell(tria, 4149 tria_cell->level(), 4150 tria_cell->index(), 4151 &mesh); 4152 4153 unpack(cell, *data); 4154 } 4155 } 4156 4157 // make sure that all communication is finished 4158 // when we leave this function. 4159 if (requests.size() > 0) 4160 { 4161 const int ierr = 4162 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); 4163 AssertThrowMPI(ierr); 4164 } 4165 if (reply_requests.size() > 0) 4166 { 4167 const int ierr = MPI_Waitall(reply_requests.size(), 4168 reply_requests.data(), 4169 MPI_STATUSES_IGNORE); 4170 AssertThrowMPI(ierr); 4171 } 4172 4173 4174 # endif // DEAL_II_WITH_MPI 4175 } 4176 4177 } // namespace internal 4178 4179 template <typename DataType, typename MeshType> 4180 void exchange_cell_data_to_ghosts(const MeshType & mesh,const std::function<std_cxx17::optional<DataType> (const typename MeshType::active_cell_iterator &)> & pack,const std::function<void (const typename MeshType::active_cell_iterator &,const DataType &)> & unpack,const std::function<bool (const typename MeshType::active_cell_iterator &)> & cell_filter)4181 exchange_cell_data_to_ghosts( 4182 const MeshType & mesh, 4183 const std::function<std_cxx17::optional<DataType>( 4184 const typename MeshType::active_cell_iterator &)> &pack, 4185 const std::function<void(const typename MeshType::active_cell_iterator &, 4186 const DataType &)> & unpack, 4187 const std::function<bool(const typename MeshType::active_cell_iterator &)> 4188 &cell_filter) 4189 { 4190 # ifndef DEAL_II_WITH_MPI 4191 (void)mesh; 4192 (void)pack; 4193 (void)unpack; 4194 (void)cell_filter; 4195 Assert(false, 4196 ExcMessage( 4197 "GridTools::exchange_cell_data_to_ghosts() requires MPI.")); 4198 # else 4199 internal::exchange_cell_data<DataType, 4200 MeshType, 4201 typename MeshType::active_cell_iterator>( 4202 mesh, 4203 pack, 4204 unpack, 4205 cell_filter, 4206 [&](const auto &process) { 4207 for (const auto &cell : mesh.active_cell_iterators()) 4208 if (cell->is_ghost()) 4209 process(cell, cell->subdomain_id()); 4210 }, 4211 [](const auto &tria) { return tria.ghost_owners(); }); 4212 # endif 4213 } 4214 4215 4216 4217 template <typename DataType, typename MeshType> 4218 void exchange_cell_data_to_level_ghosts(const MeshType & mesh,const std::function<std_cxx17::optional<DataType> (const typename MeshType::level_cell_iterator &)> & pack,const std::function<void (const typename MeshType::level_cell_iterator &,const DataType &)> & unpack,const std::function<bool (const typename MeshType::level_cell_iterator &)> & cell_filter)4219 exchange_cell_data_to_level_ghosts( 4220 const MeshType & mesh, 4221 const std::function<std_cxx17::optional<DataType>( 4222 const typename MeshType::level_cell_iterator &)> &pack, 4223 const std::function<void(const typename MeshType::level_cell_iterator &, 4224 const DataType &)> & unpack, 4225 const std::function<bool(const typename MeshType::level_cell_iterator &)> 4226 &cell_filter) 4227 { 4228 # ifndef DEAL_II_WITH_MPI 4229 (void)mesh; 4230 (void)pack; 4231 (void)unpack; 4232 (void)cell_filter; 4233 Assert(false, 4234 ExcMessage( 4235 "GridTools::exchange_cell_data_to_level_ghosts() requires MPI.")); 4236 # else 4237 internal::exchange_cell_data<DataType, 4238 MeshType, 4239 typename MeshType::level_cell_iterator>( 4240 mesh, 4241 pack, 4242 unpack, 4243 cell_filter, 4244 [&](const auto &process) { 4245 for (const auto &cell : mesh.cell_iterators()) 4246 if (cell->level_subdomain_id() != 4247 dealii::numbers::artificial_subdomain_id && 4248 !cell->is_locally_owned_on_level()) 4249 process(cell, cell->level_subdomain_id()); 4250 }, 4251 [](const auto &tria) { return tria.level_ghost_owners(); }); 4252 # endif 4253 } 4254 } // namespace GridTools 4255 4256 # endif 4257 4258 DEAL_II_NAMESPACE_CLOSE 4259 4260 /*---------------------------- grid_tools.h ---------------------------*/ 4261 /* end of #ifndef dealii_grid_tools_h */ 4262 #endif 4263 /*---------------------------- grid_tools.h ---------------------------*/ 4264