Reference documentation for deal.II version Git 5b8b897cb2 2021-04-22 22:24:19 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.h
Go to the documentation of this file.
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 
25 
27 
29 
31 
32 # include <deal.II/fe/mapping.h>
33 # include <deal.II/fe/mapping_q1.h>
34 
35 # include <deal.II/grid/manifold.h>
36 # include <deal.II/grid/tria.h>
39 
40 # include <deal.II/hp/dof_handler.h>
41 
43 
44 # include <deal.II/numerics/rtree.h>
45 
47 # include <boost/archive/binary_iarchive.hpp>
48 # include <boost/archive/binary_oarchive.hpp>
49 # include <boost/geometry/index/rtree.hpp>
50 # include <boost/random/mersenne_twister.hpp>
51 # include <boost/serialization/array.hpp>
52 # include <boost/serialization/vector.hpp>
53 
54 # ifdef DEAL_II_WITH_ZLIB
55 # include <boost/iostreams/device/back_inserter.hpp>
56 # include <boost/iostreams/filter/gzip.hpp>
57 # include <boost/iostreams/filtering_stream.hpp>
58 # include <boost/iostreams/stream.hpp>
59 # endif
61 
62 # include <bitset>
63 # include <list>
64 # include <set>
65 
67 
68 // Forward declarations
69 # ifndef DOXYGEN
70 namespace parallel
71 {
72  namespace distributed
73  {
74  template <int, int>
75  class Triangulation;
76  }
77 } // namespace parallel
78 
79 namespace hp
80 {
81  template <int, int>
82  class MappingCollection;
83 }
84 
85 class SparsityPattern;
86 # endif
87 
88 namespace internal
89 {
90  template <int dim, int spacedim, class MeshType>
92  {
93  public:
94 # ifndef _MSC_VER
95  using type = typename MeshType::active_cell_iterator;
96 # else
98 # endif
99  };
100 
101 # ifdef _MSC_VER
102  template <int dim, int spacedim>
103  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
104  {
105  public:
106  using type =
108  };
109 # endif
110 
111 
112 # ifdef _MSC_VER
113  template <int dim, int spacedim>
114  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim>>
115  {
116  public:
117  using type =
119  };
120 # endif
121 } // namespace internal
122 
131 namespace GridTools
132 {
133  template <int dim, int spacedim>
134  class Cache;
135 
140 
147  template <int dim, int spacedim>
148  double
150 
177  template <int dim, int spacedim>
178  double
180  const Mapping<dim, spacedim> & mapping =
181  (ReferenceCells::get_hypercube<dim>()
182  .template get_default_linear_mapping<dim, spacedim>()));
183 
194  template <int dim, int spacedim>
195  double
198  const Mapping<dim, spacedim> & mapping =
199  (ReferenceCells::get_hypercube<dim>()
200  .template get_default_linear_mapping<dim, spacedim>()));
201 
212  template <int dim, int spacedim>
213  double
215  const Triangulation<dim, spacedim> &triangulation,
216  const Mapping<dim, spacedim> & mapping =
217  (ReferenceCells::get_hypercube<dim>()
218  .template get_default_linear_mapping<dim, spacedim>()));
219 
229  template <int dim>
230  double
231  cell_measure(
232  const std::vector<Point<dim>> &all_vertices,
234 
244  template <int dim>
245  double
246  cell_measure(const std::vector<Point<dim>> & all_vertices,
248 
254  template <int dim, typename T>
255  double
256  cell_measure(const T &, ...);
257 
280  template <int dim, int spacedim>
281  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
283 
310  template <int dim>
311  Vector<double>
313  const Triangulation<dim> &triangulation,
314  const Quadrature<dim> & quadrature);
315 
323  template <int dim>
324  double
326  const Triangulation<dim> &triangulation,
327  const Quadrature<dim> & quadrature);
328 
342  template <int dim, int spacedim>
345 
363  template <typename Iterator>
366  const Iterator & object,
368 
381  template <int dim, int spacedim>
382  std::
383  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
385 
391 
408  template <int dim, int spacedim>
409  void
411  std::vector<CellData<dim>> & cells,
412  SubCellData & subcelldata);
413 
431  template <int dim, int spacedim>
432  void
433  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
434  std::vector<CellData<dim>> & cells,
435  SubCellData & subcelldata,
436  std::vector<unsigned int> & considered_vertices,
437  const double tol = 1e-12);
438 
457  template <int dim, int spacedim>
458  void
460  const std::vector<Point<spacedim>> &all_vertices,
461  std::vector<CellData<dim>> & cells);
462 
472  template <int dim>
473  void
474  consistently_order_cells(std::vector<CellData<dim>> &cells);
475 
481 
535  template <int dim, typename Transformation, int spacedim>
536  void
537  transform(const Transformation & transformation,
538  Triangulation<dim, spacedim> &triangulation);
539 
545  template <int dim, int spacedim>
546  void
547  shift(const Tensor<1, spacedim> & shift_vector,
548  Triangulation<dim, spacedim> &triangulation);
549 
550 
560  template <int dim>
561  void
562  rotate(const double angle, Triangulation<dim> &triangulation);
563 
576  template <int dim>
577  void
578  rotate(const double angle,
579  const unsigned int axis,
580  Triangulation<dim, 3> &triangulation);
581 
639  template <int dim>
640  void
641  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
642  Triangulation<dim> & tria,
643  const Function<dim, double> *coefficient = nullptr,
644  const bool solve_for_absolute_positions = false);
645 
651  template <int dim, int spacedim>
652  std::map<unsigned int, Point<spacedim>>
654 
662  template <int dim, int spacedim>
663  void
664  scale(const double scaling_factor,
665  Triangulation<dim, spacedim> &triangulation);
666 
681  template <int dim, int spacedim>
682  void
684  const double factor,
685  Triangulation<dim, spacedim> &triangulation,
686  const bool keep_boundary = true,
687  const unsigned int seed = boost::random::mt19937::default_seed);
688 
722  template <int dim, int spacedim>
723  void
725  const bool isotropic = false,
726  const unsigned int max_iterations = 100);
727 
752  template <int dim, int spacedim>
753  void
755  const double max_ratio = 1.6180339887,
756  const unsigned int max_iterations = 5);
757 
847  template <int dim, int spacedim>
848  void
850  const double limit_angle_fraction = .75);
851 
857 
907  template <int dim, int spacedim>
908 # ifndef DOXYGEN
909  std::tuple<
910  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
911  std::vector<std::vector<Point<dim>>>,
912  std::vector<std::vector<unsigned int>>>
913 # else
914  return_type
915 # endif
917  const Cache<dim, spacedim> & cache,
918  const std::vector<Point<spacedim>> &points,
920  &cell_hint =
922 
951  template <int dim, int spacedim>
952 # ifndef DOXYGEN
953  std::tuple<
954  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
955  std::vector<std::vector<Point<dim>>>,
956  std::vector<std::vector<unsigned int>>,
957  std::vector<unsigned int>>
958 # else
959  return_type
960 # endif
962  const Cache<dim, spacedim> & cache,
963  const std::vector<Point<spacedim>> &points,
965  &cell_hint =
967 
1037  template <int dim, int spacedim>
1038 # ifndef DOXYGEN
1039  std::tuple<
1040  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1041  std::vector<std::vector<Point<dim>>>,
1042  std::vector<std::vector<unsigned int>>,
1043  std::vector<std::vector<Point<spacedim>>>,
1044  std::vector<std::vector<unsigned int>>>
1045 # else
1046  return_type
1047 # endif
1049  const GridTools::Cache<dim, spacedim> & cache,
1050  const std::vector<Point<spacedim>> & local_points,
1051  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1052  const double tolerance = 1e-10);
1053 
1054  namespace internal
1055  {
1070  template <int dim, int spacedim>
1072  {
1079  std::vector<std::tuple<std::pair<int, int>,
1080  unsigned int,
1081  unsigned int,
1082  Point<dim>,
1084  unsigned int>>
1086 
1090  std::vector<unsigned int> send_ranks;
1091 
1097  std::vector<unsigned int> send_ptrs;
1098 
1109  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1111 
1115  std::vector<unsigned int> recv_ranks;
1116 
1122  std::vector<unsigned int> recv_ptrs;
1123  };
1124 
1134  template <int dim, int spacedim>
1137  const GridTools::Cache<dim, spacedim> & cache,
1138  const std::vector<Point<spacedim>> & points,
1139  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1140  const double tolerance,
1141  const bool perform_handshake);
1142 
1143  } // namespace internal
1144 
1181  template <int dim, int spacedim>
1182  std::map<unsigned int, Point<spacedim>>
1184  const Triangulation<dim, spacedim> &container,
1185  const Mapping<dim, spacedim> & mapping =
1186  (ReferenceCells::get_hypercube<dim>()
1187  .template get_default_linear_mapping<dim, spacedim>()));
1188 
1198  template <int spacedim>
1199  unsigned int
1200  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1201  const Point<spacedim> & p);
1202 
1226  template <int dim, template <int, int> class MeshType, int spacedim>
1227  unsigned int
1228  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1229  const Point<spacedim> & p,
1230  const std::vector<bool> & marked_vertices = {});
1231 
1255  template <int dim, template <int, int> class MeshType, int spacedim>
1256  unsigned int
1258  const MeshType<dim, spacedim> &mesh,
1259  const Point<spacedim> & p,
1260  const std::vector<bool> & marked_vertices = {});
1261 
1262 
1287  template <int dim, template <int, int> class MeshType, int spacedim>
1288 # ifndef _MSC_VER
1289  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1290 # else
1291  std::vector<
1292  typename ::internal::
1293  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1294 # endif
1295  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1296  const unsigned int vertex_index);
1297 
1365  template <int dim, template <int, int> class MeshType, int spacedim>
1366 # ifndef _MSC_VER
1367  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1368 # else
1369  std::pair<typename ::internal::
1370  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1371  Point<dim>>
1372 # endif
1374  const MeshType<dim, spacedim> &mesh,
1375  const Point<spacedim> & p,
1376  const std::vector<bool> &marked_vertices = {},
1377  const double tolerance = 1.e-10);
1378 
1386  template <int dim, template <int, int> class MeshType, int spacedim>
1387 # ifndef _MSC_VER
1388  typename MeshType<dim, spacedim>::active_cell_iterator
1389 # else
1390  typename ::internal::
1391  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1392 # endif
1393  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1394  const Point<spacedim> & p,
1395  const std::vector<bool> &marked_vertices = {},
1396  const double tolerance = 1.e-10);
1397 
1404  template <int dim, int spacedim>
1405  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1406  Point<dim>>
1408  const hp::MappingCollection<dim, spacedim> &mapping,
1409  const DoFHandler<dim, spacedim> & mesh,
1410  const Point<spacedim> & p,
1411  const double tolerance = 1.e-10);
1412 
1461  template <int dim, int spacedim>
1462  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1463  Point<dim>>
1465  const Cache<dim, spacedim> &cache,
1466  const Point<spacedim> & p,
1469  const std::vector<bool> &marked_vertices = {},
1470  const double tolerance = 1.e-10);
1471 
1485  template <int dim, template <int, int> class MeshType, int spacedim>
1486 # ifndef _MSC_VER
1487  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1488 # else
1489  std::pair<typename ::internal::
1490  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1491  Point<dim>>
1492 # endif
1494  const Mapping<dim, spacedim> & mapping,
1495  const MeshType<dim, spacedim> &mesh,
1496  const Point<spacedim> & p,
1497  const std::vector<
1498  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1500  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1501  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1502  typename MeshType<dim, spacedim>::active_cell_iterator(),
1503  const std::vector<bool> & marked_vertices = {},
1504  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1505  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1506  const double tolerance = 1.e-10);
1507 
1528  template <int dim, template <int, int> class MeshType, int spacedim>
1529 # ifndef _MSC_VER
1530  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1531  Point<dim>>>
1532 # else
1533  std::vector<std::pair<
1534  typename ::internal::
1535  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1536  Point<dim>>>
1537 # endif
1539  const Mapping<dim, spacedim> & mapping,
1540  const MeshType<dim, spacedim> &mesh,
1541  const Point<spacedim> & p,
1542  const double tolerance,
1543  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1544  Point<dim>> & first_cell);
1545 
1552  template <int dim, template <int, int> class MeshType, int spacedim>
1553 # ifndef _MSC_VER
1554  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1555  Point<dim>>>
1556 # else
1557  std::vector<std::pair<
1558  typename ::internal::
1559  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1560  Point<dim>>>
1561 # endif
1563  const Mapping<dim, spacedim> & mapping,
1564  const MeshType<dim, spacedim> &mesh,
1565  const Point<spacedim> & p,
1566  const double tolerance = 1e-10,
1567  const std::vector<bool> & marked_vertices = {});
1568 
1590  template <class MeshType>
1591  std::vector<typename MeshType::active_cell_iterator>
1592  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1593 
1618  template <class MeshType>
1619  void
1621  const typename MeshType::active_cell_iterator & cell,
1622  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1623 
1673  template <class MeshType>
1674  std::vector<typename MeshType::active_cell_iterator>
1676  const MeshType &mesh,
1677  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1678  &predicate);
1679 
1680 
1688  template <class MeshType>
1689  std::vector<typename MeshType::cell_iterator>
1691  const MeshType &mesh,
1692  const std::function<bool(const typename MeshType::cell_iterator &)>
1693  & predicate,
1694  const unsigned int level);
1695 
1696 
1709  template <class MeshType>
1710  std::vector<typename MeshType::active_cell_iterator>
1711  compute_ghost_cell_halo_layer(const MeshType &mesh);
1712 
1761  template <class MeshType>
1762  std::vector<typename MeshType::active_cell_iterator>
1764  const MeshType &mesh,
1765  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1766  & predicate,
1767  const double layer_thickness);
1768 
1791  template <class MeshType>
1792  std::vector<typename MeshType::active_cell_iterator>
1793  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1794  const double layer_thickness);
1795 
1811  template <class MeshType>
1812  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1814  const MeshType &mesh,
1815  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1816  &predicate);
1817 
1870  template <class MeshType>
1871  std::vector<BoundingBox<MeshType::space_dimension>>
1873  const MeshType &mesh,
1874  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1875  & predicate,
1876  const unsigned int refinement_level = 0,
1877  const bool allow_merge = false,
1878  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1879 
1907  template <int spacedim>
1908 # ifndef DOXYGEN
1909  std::tuple<std::vector<std::vector<unsigned int>>,
1910  std::map<unsigned int, unsigned int>,
1911  std::map<unsigned int, std::vector<unsigned int>>>
1912 # else
1913  return_type
1914 # endif
1916  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1917  const std::vector<Point<spacedim>> & points);
1918 
1919 
1954  template <int spacedim>
1955 # ifndef DOXYGEN
1956  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1957  std::map<unsigned int, unsigned int>,
1958  std::map<unsigned int, std::vector<unsigned int>>>
1959 # else
1960  return_type
1961 # endif
1963  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1964  const std::vector<Point<spacedim>> & points);
1965 
1966 
1975  template <int dim, int spacedim>
1976  std::vector<
1977  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1978  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1979 
1992  template <int dim, int spacedim>
1993  std::vector<std::vector<Tensor<1, spacedim>>>
1995  const Triangulation<dim, spacedim> &mesh,
1996  const std::vector<
1998  &vertex_to_cells);
1999 
2000 
2008  template <int dim, int spacedim>
2009  unsigned int
2012  const Point<spacedim> & position,
2013  const Mapping<dim, spacedim> & mapping =
2014  (ReferenceCells::get_hypercube<dim>()
2015  .template get_default_linear_mapping<dim, spacedim>()));
2016 
2028  template <int dim, int spacedim>
2029  std::map<unsigned int, types::global_vertex_index>
2032 
2044  template <int dim, int spacedim>
2045  std::pair<unsigned int, double>
2048 
2054 
2063  template <int dim, int spacedim>
2064  void
2066  const Triangulation<dim, spacedim> &triangulation,
2067  DynamicSparsityPattern & connectivity);
2068 
2077  template <int dim, int spacedim>
2078  void
2080  const Triangulation<dim, spacedim> &triangulation,
2081  DynamicSparsityPattern & connectivity);
2082 
2091  template <int dim, int spacedim>
2092  void
2094  const Triangulation<dim, spacedim> &triangulation,
2095  const unsigned int level,
2096  DynamicSparsityPattern & connectivity);
2097 
2118  template <int dim, int spacedim>
2119  void
2120  partition_triangulation(const unsigned int n_partitions,
2121  Triangulation<dim, spacedim> & triangulation,
2122  const SparsityTools::Partitioner partitioner =
2124 
2135  template <int dim, int spacedim>
2136  void
2137  partition_triangulation(const unsigned int n_partitions,
2138  const std::vector<unsigned int> &cell_weights,
2139  Triangulation<dim, spacedim> & triangulation,
2140  const SparsityTools::Partitioner partitioner =
2142 
2188  template <int dim, int spacedim>
2189  void
2190  partition_triangulation(const unsigned int n_partitions,
2191  const SparsityPattern & cell_connection_graph,
2192  Triangulation<dim, spacedim> &triangulation,
2193  const SparsityTools::Partitioner partitioner =
2195 
2206  template <int dim, int spacedim>
2207  void
2208  partition_triangulation(const unsigned int n_partitions,
2209  const std::vector<unsigned int> &cell_weights,
2210  const SparsityPattern & cell_connection_graph,
2211  Triangulation<dim, spacedim> &triangulation,
2212  const SparsityTools::Partitioner partitioner =
2214 
2229  template <int dim, int spacedim>
2230  void
2231  partition_triangulation_zorder(const unsigned int n_partitions,
2232  Triangulation<dim, spacedim> &triangulation,
2233  const bool group_siblings = true);
2234 
2246  template <int dim, int spacedim>
2247  void
2249 
2257  template <int dim, int spacedim>
2258  std::vector<types::subdomain_id>
2260  const std::vector<CellId> & cell_ids);
2261 
2272  template <int dim, int spacedim>
2273  void
2275  std::vector<types::subdomain_id> & subdomain);
2276 
2291  template <int dim, int spacedim>
2292  unsigned int
2294  const Triangulation<dim, spacedim> &triangulation,
2295  const types::subdomain_id subdomain);
2296 
2326  template <int dim, int spacedim>
2327  std::vector<bool>
2329 
2335 
2369  template <typename MeshType>
2370  std::list<std::pair<typename MeshType::cell_iterator,
2371  typename MeshType::cell_iterator>>
2372  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2373 
2383  template <int dim, int spacedim>
2384  bool
2386  const Triangulation<dim, spacedim> &mesh_2);
2387 
2397  template <typename MeshType>
2398  bool
2399  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2400 
2406 
2422  template <int dim, int spacedim>
2426  & distorted_cells,
2427  Triangulation<dim, spacedim> &triangulation);
2428 
2429 
2430 
2439 
2440 
2478  template <class MeshType>
2479  std::vector<typename MeshType::active_cell_iterator>
2480  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2481 
2482 
2504  template <class Container>
2505  std::vector<typename Container::cell_iterator>
2507  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2508 
2575  template <class Container>
2576  void
2578  const std::vector<typename Container::active_cell_iterator> &patch,
2580  &local_triangulation,
2581  std::map<
2582  typename Triangulation<Container::dimension,
2583  Container::space_dimension>::active_cell_iterator,
2584  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2585 
2617  template <int dim, int spacedim>
2618  std::map<
2620  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2622 
2623 
2630 
2636  template <typename CellIterator>
2638  {
2642  CellIterator cell[2];
2643 
2648  unsigned int face_idx[2];
2649 
2655  std::bitset<3> orientation;
2656 
2670  };
2671 
2672 
2736  template <typename FaceIterator>
2737  bool orthogonal_equality(
2738  std::bitset<3> & orientation,
2739  const FaceIterator & face1,
2740  const FaceIterator & face2,
2741  const int direction,
2745 
2746 
2750  template <typename FaceIterator>
2751  bool
2753  const FaceIterator & face1,
2754  const FaceIterator & face2,
2755  const int direction,
2759 
2760 
2817  template <typename MeshType>
2818  void
2820  const MeshType & mesh,
2821  const types::boundary_id b_id1,
2822  const types::boundary_id b_id2,
2823  const int direction,
2825  & matched_pairs,
2826  const Tensor<1, MeshType::space_dimension> &offset =
2829 
2830 
2853  template <typename MeshType>
2854  void
2856  const MeshType & mesh,
2857  const types::boundary_id b_id,
2858  const int direction,
2860  & matched_pairs,
2861  const ::Tensor<1, MeshType::space_dimension> &offset =
2864 
2870 
2891  template <int dim, int spacedim>
2892  void
2894  const bool reset_boundary_ids = false);
2895 
2917  template <int dim, int spacedim>
2918  void
2920  const std::vector<types::boundary_id> &src_boundary_ids,
2921  const std::vector<types::manifold_id> &dst_manifold_ids,
2923  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2924 
2954  template <int dim, int spacedim>
2955  void
2957  const bool compute_face_ids = false);
2958 
2983  template <int dim, int spacedim>
2984  void
2987  const std::function<types::manifold_id(
2988  const std::set<types::manifold_id> &)> &disambiguation_function =
2989  [](const std::set<types::manifold_id> &manifold_ids) {
2990  if (manifold_ids.size() == 1)
2991  return *manifold_ids.begin();
2992  else
2994  },
2995  bool overwrite_only_flat_manifold_ids = true);
3082  template <typename DataType, typename MeshType>
3083  void
3085  const MeshType & mesh,
3086  const std::function<std_cxx17::optional<DataType>(
3087  const typename MeshType::active_cell_iterator &)> &pack,
3088  const std::function<void(const typename MeshType::active_cell_iterator &,
3089  const DataType &)> & unpack,
3090  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3091  &cell_filter =
3093 
3104  template <typename DataType, typename MeshType>
3105  void
3107  const MeshType & mesh,
3108  const std::function<std_cxx17::optional<DataType>(
3109  const typename MeshType::level_cell_iterator &)> &pack,
3110  const std::function<void(const typename MeshType::level_cell_iterator &,
3111  const DataType &)> & unpack,
3112  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3114  true});
3115 
3116  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3117  * boxes @p local_bboxes.
3118  *
3119  * This function is meant to exchange bounding boxes describing the locally
3120  * owned cells in a distributed triangulation obtained with the function
3121  * GridTools::compute_mesh_predicate_bounding_box .
3122  *
3123  * The output vector's size is the number of processes of the MPI
3124  * communicator:
3125  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3126  */
3127  template <int spacedim>
3128  std::vector<std::vector<BoundingBox<spacedim>>>
3130  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3131  const MPI_Comm & mpi_communicator);
3132 
3165  template <int spacedim>
3168  const std::vector<BoundingBox<spacedim>> &local_description,
3169  const MPI_Comm & mpi_communicator);
3170 
3188  template <int dim, int spacedim>
3189  void
3191  const Triangulation<dim, spacedim> & tria,
3192  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3193  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3194 
3201  template <int dim, int spacedim>
3202  std::map<unsigned int, std::set<::types::subdomain_id>>
3204  const Triangulation<dim, spacedim> &tria);
3205 
3218  template <int dim, typename T>
3220  {
3224  std::vector<CellId> cell_ids;
3225 
3229  std::vector<T> data;
3230 
3239  template <class Archive>
3240  void
3241  save(Archive &ar, const unsigned int version) const;
3242 
3249  template <class Archive>
3250  void
3251  load(Archive &ar, const unsigned int version);
3252 
3253 # ifdef DOXYGEN
3254 
3259  template <class Archive>
3260  void
3261  serialize(Archive &archive, const unsigned int version);
3262 # else
3263  // This macro defines the serialize() method that is compatible with
3264  // the templated save() and load() method that have been implemented.
3265  BOOST_SERIALIZATION_SPLIT_MEMBER()
3266 # endif
3267  };
3268 
3273 
3278  int,
3279  << "The number of partitions you gave is " << arg1
3280  << ", but must be greater than zero.");
3285  int,
3286  << "The subdomain id " << arg1
3287  << " has no cells associated with it.");
3292 
3297  double,
3298  << "The scaling factor must be positive, but it is " << arg1
3299  << ".");
3303  template <int N>
3305  Point<N>,
3306  << "The point <" << arg1
3307  << "> could not be found inside any of the "
3308  << "coarse grid cells.");
3312  template <int N>
3314  Point<N>,
3315  << "The point <" << arg1
3316  << "> could not be found inside any of the "
3317  << "subcells of a coarse grid cell.");
3318 
3323  unsigned int,
3324  << "The given vertex with index " << arg1
3325  << " is not used in the given triangulation.");
3326 
3329 } /*namespace GridTools*/
3330 
3331 
3342  "The edges of the mesh are not consistently orientable.");
3343 
3344 
3345 
3346 /* ----------------- Template function --------------- */
3347 
3348 # ifndef DOXYGEN
3349 
3350 namespace GridTools
3351 {
3352  template <int dim>
3353  double
3354  cell_measure(
3355  const std::vector<Point<dim>> &all_vertices,
3356  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3357  {
3358  // We forward call to the ArrayView version:
3359  const ArrayView<const unsigned int> view(
3360  indices, GeometryInfo<dim>::vertices_per_cell);
3361  return cell_measure(all_vertices, view);
3362  }
3363 
3364  template <int dim, typename T>
3365  double
3366  cell_measure(const T &, ...)
3367  {
3368  Assert(false, ExcNotImplemented());
3369  return std::numeric_limits<double>::quiet_NaN();
3370  }
3371 
3372 
3373 
3374  // This specialization is defined here so that the general template in the
3375  // source file doesn't need to have further 1D overloads for the internal
3376  // functions it calls.
3377  template <>
3381  {
3382  return {};
3383  }
3384 
3385 
3386 
3387  template <int dim, typename Predicate, int spacedim>
3388  void
3389  transform(const Predicate & predicate,
3390  Triangulation<dim, spacedim> &triangulation)
3391  {
3392  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3393 
3394  // loop over all active cells, and
3395  // transform those vertices that
3396  // have not yet been touched. note
3397  // that we get to all vertices in
3398  // the triangulation by only
3399  // visiting the active cells.
3401  cell = triangulation.begin_active(),
3402  endc = triangulation.end();
3403  for (; cell != endc; ++cell)
3404  for (const unsigned int v : cell->vertex_indices())
3405  if (treated_vertices[cell->vertex_index(v)] == false)
3406  {
3407  // transform this vertex
3408  cell->vertex(v) = predicate(cell->vertex(v));
3409  // and mark it as treated
3410  treated_vertices[cell->vertex_index(v)] = true;
3411  };
3412 
3413 
3414  // now fix any vertices on hanging nodes so that we don't create any holes
3415  if (dim == 2)
3416  {
3418  cell = triangulation.begin_active(),
3419  endc = triangulation.end();
3420  for (; cell != endc; ++cell)
3421  for (const unsigned int face : cell->face_indices())
3422  if (cell->face(face)->has_children() &&
3423  !cell->face(face)->at_boundary())
3424  {
3425  Assert(cell->reference_cell() ==
3426  ReferenceCells::get_hypercube<dim>(),
3427  ExcNotImplemented());
3428 
3429  // this line has children
3430  cell->face(face)->child(0)->vertex(1) =
3431  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3432  2;
3433  }
3434  }
3435  else if (dim == 3)
3436  {
3438  cell = triangulation.begin_active(),
3439  endc = triangulation.end();
3440  for (; cell != endc; ++cell)
3441  for (const unsigned int face : cell->face_indices())
3442  if (cell->face(face)->has_children() &&
3443  !cell->face(face)->at_boundary())
3444  {
3445  Assert(cell->reference_cell() ==
3446  ReferenceCells::get_hypercube<dim>(),
3447  ExcNotImplemented());
3448 
3449  // this face has hanging nodes
3450  cell->face(face)->child(0)->vertex(1) =
3451  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3452  2.0;
3453  cell->face(face)->child(0)->vertex(2) =
3454  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3455  2.0;
3456  cell->face(face)->child(1)->vertex(3) =
3457  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3458  2.0;
3459  cell->face(face)->child(2)->vertex(3) =
3460  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3461  2.0;
3462 
3463  // center of the face
3464  cell->face(face)->child(0)->vertex(3) =
3465  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3466  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3467  4.0;
3468  }
3469  }
3470 
3471  // Make sure FEValues notices that the mesh has changed
3472  triangulation.signals.mesh_movement();
3473  }
3474 
3475 
3476 
3477  template <class MeshType>
3478  std::vector<typename MeshType::active_cell_iterator>
3479  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3480  {
3481  std::vector<typename MeshType::active_cell_iterator> child_cells;
3482 
3483  if (cell->has_children())
3484  {
3485  for (unsigned int child = 0; child < cell->n_children(); ++child)
3486  if (cell->child(child)->has_children())
3487  {
3488  const std::vector<typename MeshType::active_cell_iterator>
3489  children = get_active_child_cells<MeshType>(cell->child(child));
3490  child_cells.insert(child_cells.end(),
3491  children.begin(),
3492  children.end());
3493  }
3494  else
3495  child_cells.push_back(cell->child(child));
3496  }
3497 
3498  return child_cells;
3499  }
3500 
3501 
3502 
3503  template <class MeshType>
3504  void
3506  const typename MeshType::active_cell_iterator & cell,
3507  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3508  {
3509  active_neighbors.clear();
3510  for (const unsigned int n : cell->face_indices())
3511  if (!cell->at_boundary(n))
3512  {
3513  if (MeshType::dimension == 1)
3514  {
3515  // check children of neighbor. note
3516  // that in 1d children of the neighbor
3517  // may be further refined. In 1d the
3518  // case is simple since we know what
3519  // children bound to the present cell
3520  typename MeshType::cell_iterator neighbor_child =
3521  cell->neighbor(n);
3522  if (!neighbor_child->is_active())
3523  {
3524  while (neighbor_child->has_children())
3525  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3526 
3527  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3528  ExcInternalError());
3529  }
3530  active_neighbors.push_back(neighbor_child);
3531  }
3532  else
3533  {
3534  if (cell->face(n)->has_children())
3535  // this neighbor has children. find
3536  // out which border to the present
3537  // cell
3538  for (unsigned int c = 0;
3539  c < cell->face(n)->n_active_descendants();
3540  ++c)
3541  active_neighbors.push_back(
3542  cell->neighbor_child_on_subface(n, c));
3543  else
3544  {
3545  // the neighbor must be active
3546  // himself
3547  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3548  active_neighbors.push_back(cell->neighbor(n));
3549  }
3550  }
3551  }
3552  }
3553 
3554 
3555 
3556  namespace internal
3557  {
3558  namespace ProjectToObject
3559  {
3572  struct CrossDerivative
3573  {
3574  const unsigned int direction_0;
3575  const unsigned int direction_1;
3576 
3577  CrossDerivative(const unsigned int d0, const unsigned int d1);
3578  };
3579 
3580  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3581  const unsigned int d1)
3582  : direction_0(d0)
3583  , direction_1(d1)
3584  {}
3585 
3586 
3587 
3592  template <typename F>
3593  inline auto
3594  centered_first_difference(const double center,
3595  const double step,
3596  const F &f) -> decltype(f(center) - f(center))
3597  {
3598  return (f(center + step) - f(center - step)) / (2.0 * step);
3599  }
3600 
3601 
3602 
3607  template <typename F>
3608  inline auto
3609  centered_second_difference(const double center,
3610  const double step,
3611  const F &f) -> decltype(f(center) - f(center))
3612  {
3613  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3614  (step * step);
3615  }
3616 
3617 
3618 
3628  template <int structdim, typename F>
3629  inline auto
3630  cross_stencil(
3631  const CrossDerivative cross_derivative,
3633  const double step,
3634  const F &f) -> decltype(f(center) - f(center))
3635  {
3637  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3638  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3639  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3640  1.0 / 3.0 * f(center - simplex_vector) +
3641  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3642  step;
3643  }
3644 
3645 
3646 
3653  template <int spacedim, int structdim, typename F>
3654  inline double
3655  gradient_entry(
3656  const unsigned int row_n,
3657  const unsigned int dependent_direction,
3658  const Point<spacedim> &p0,
3660  const double step,
3661  const F & f)
3662  {
3664  dependent_direction <
3666  ExcMessage("This function assumes that the last weight is a "
3667  "dependent variable (and hence we cannot take its "
3668  "derivative directly)."));
3669  Assert(row_n != dependent_direction,
3670  ExcMessage(
3671  "We cannot differentiate with respect to the variable "
3672  "that is assumed to be dependent."));
3673 
3674  const Point<spacedim> manifold_point = f(center);
3675  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3676  {row_n, dependent_direction}, center, step, f);
3677  double entry = 0.0;
3678  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3679  entry +=
3680  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3681  return entry;
3682  }
3683 
3689  template <typename Iterator, int spacedim, int structdim>
3691  project_to_d_linear_object(const Iterator & object,
3692  const Point<spacedim> &trial_point)
3693  {
3694  // let's look at this for simplicity for a quad (structdim==2) in a
3695  // space with spacedim>2 (notate trial_point by y): all points on the
3696  // surface are given by
3697  // x(\xi) = sum_i v_i phi_x(\xi)
3698  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3699  // reference coordinates of the quad. so what we are trying to do is
3700  // find a point x on the surface that is closest to the point y. there
3701  // are different ways to solve this problem, but in the end it's a
3702  // nonlinear problem and we have to find reference coordinates \xi so
3703  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3704  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3705  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3706  // have to use a Newton method to find the answer. This leads to the
3707  // following formulation of Newton steps:
3708  //
3709  // Given \xi_k, find \delta\xi_k so that
3710  // H_k \delta\xi_k = - F_k
3711  // where H_k is an approximation to the second derivatives of J at
3712  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3713  // number of times until the right hand side is small enough. As a
3714  // stopping criterion, we terminate if ||\delta\xi||<eps.
3715  //
3716  // As for the Hessian, the best choice would be
3717  // H_k = J''(\xi_k)
3718  // but we'll opt for the simpler Gauss-Newton form
3719  // H_k = A^T A
3720  // i.e.
3721  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3722  // \partial_n phi_i *
3723  // \partial_m phi_j
3724  // we start at xi=(0.5, 0.5).
3725  Point<structdim> xi;
3726  for (unsigned int d = 0; d < structdim; ++d)
3727  xi[d] = 0.5;
3728 
3729  Point<spacedim> x_k;
3730  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3731  x_k += object->vertex(i) *
3733 
3734  do
3735  {
3737  for (const unsigned int i :
3739  F_k +=
3740  (x_k - trial_point) * object->vertex(i) *
3742  i);
3743 
3745  for (const unsigned int i :
3747  for (const unsigned int j :
3749  {
3752  xi, i),
3754  xi, j));
3755  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3756  }
3757 
3758  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3759  xi += delta_xi;
3760 
3761  x_k = Point<spacedim>();
3762  for (const unsigned int i :
3764  x_k += object->vertex(i) *
3766 
3767  if (delta_xi.norm() < 1e-7)
3768  break;
3769  }
3770  while (true);
3771 
3772  return x_k;
3773  }
3774  } // namespace ProjectToObject
3775  } // namespace internal
3776 
3777 
3778 
3779  namespace internal
3780  {
3781  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3782  // inside the project_to_object function below.
3783  template <int structdim>
3784  inline bool
3785  weights_are_ok(
3787  {
3788  // clang has trouble figuring out structdim here, so define it
3789  // again:
3790  static const std::size_t n_vertices_per_cell =
3792  n_independent_components;
3793  std::array<double, n_vertices_per_cell> copied_weights;
3794  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3795  {
3796  copied_weights[i] = v[i];
3797  if (v[i] < 0.0 || v[i] > 1.0)
3798  return false;
3799  }
3800 
3801  // check the sum: try to avoid some roundoff errors by summing in order
3802  std::sort(copied_weights.begin(), copied_weights.end());
3803  const double sum =
3804  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3805  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3806  }
3807  } // namespace internal
3808 
3809  template <typename Iterator>
3812  const Iterator & object,
3814  {
3815  const int spacedim = Iterator::AccessorType::space_dimension;
3816  const int structdim = Iterator::AccessorType::structure_dimension;
3817 
3818  Point<spacedim> projected_point = trial_point;
3819 
3820  if (structdim >= spacedim)
3821  return projected_point;
3822  else if (structdim == 1 || structdim == 2)
3823  {
3824  using namespace internal::ProjectToObject;
3825  // Try to use the special flat algorithm for quads (this is better
3826  // than the general algorithm in 3D). This does not take into account
3827  // whether projected_point is outside the quad, but we optimize along
3828  // lines below anyway:
3829  const int dim = Iterator::AccessorType::dimension;
3830  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3831  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3832  &manifold) != nullptr)
3833  {
3834  projected_point =
3835  project_to_d_linear_object<Iterator, spacedim, structdim>(
3836  object, trial_point);
3837  }
3838  else
3839  {
3840  // We want to find a point on the convex hull (defined by the
3841  // vertices of the object and the manifold description) that is
3842  // relatively close to the trial point. This has a few issues:
3843  //
3844  // 1. For a general convex hull we are not guaranteed that a unique
3845  // minimum exists.
3846  // 2. The independent variables in the optimization process are the
3847  // weights given to Manifold::get_new_point, which must sum to 1,
3848  // so we cannot use standard finite differences to approximate a
3849  // gradient.
3850  //
3851  // There is not much we can do about 1., but for 2. we can derive
3852  // finite difference stencils that work on a structdim-dimensional
3853  // simplex and rewrite the optimization problem to use those
3854  // instead. Consider the structdim 2 case and let
3855  //
3856  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3857  // c2, c3})
3858  //
3859  // where {c0, c1, c2, c3} are the weights for the four vertices on
3860  // the quadrilateral. We seek to minimize the Euclidean distance
3861  // between F(...) and trial_point. We can solve for c3 in terms of
3862  // the other weights and get, for one coordinate direction
3863  //
3864  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3865  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3866  //
3867  // where we substitute back in for c3 after taking the
3868  // derivative. We can compute a stencil for the cross derivative
3869  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3870  // (and gradient_entry computes the sum over the independent
3871  // variables). Below, we somewhat arbitrarily pick the last
3872  // component as the dependent one.
3873  //
3874  // Since we can now calculate derivatives of the objective
3875  // function we can use gradient descent to minimize it.
3876  //
3877  // Of course, this is much simpler in the structdim = 1 case (we
3878  // could rewrite the projection as a 1D optimization problem), but
3879  // to reduce the potential for bugs we use the same code in both
3880  // cases.
3881  const double step_size = object->diameter() / 64.0;
3882 
3883  constexpr unsigned int n_vertices_per_cell =
3885 
3886  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3887  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3888  ++vertex_n)
3889  vertices[vertex_n] = object->vertex(vertex_n);
3890 
3891  auto get_point_from_weights =
3892  [&](const Tensor<1, n_vertices_per_cell> &weights)
3893  -> Point<spacedim> {
3894  return object->get_manifold().get_new_point(
3895  make_array_view(vertices.begin(), vertices.end()),
3896  make_array_view(weights.begin_raw(), weights.end_raw()));
3897  };
3898 
3899  // pick the initial weights as (normalized) inverse distances from
3900  // the trial point:
3901  Tensor<1, n_vertices_per_cell> guess_weights;
3902  double guess_weights_sum = 0.0;
3903  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3904  ++vertex_n)
3905  {
3906  const double distance =
3907  vertices[vertex_n].distance(trial_point);
3908  if (distance == 0.0)
3909  {
3910  guess_weights = 0.0;
3911  guess_weights[vertex_n] = 1.0;
3912  guess_weights_sum = 1.0;
3913  break;
3914  }
3915  else
3916  {
3917  guess_weights[vertex_n] = 1.0 / distance;
3918  guess_weights_sum += guess_weights[vertex_n];
3919  }
3920  }
3921  guess_weights /= guess_weights_sum;
3922  Assert(internal::weights_are_ok<structdim>(guess_weights),
3923  ExcInternalError());
3924 
3925  // The optimization algorithm consists of two parts:
3926  //
3927  // 1. An outer loop where we apply the gradient descent algorithm.
3928  // 2. An inner loop where we do a line search to find the optimal
3929  // length of the step one should take in the gradient direction.
3930  //
3931  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3932  {
3933  const unsigned int dependent_direction =
3934  n_vertices_per_cell - 1;
3935  Tensor<1, n_vertices_per_cell> current_gradient;
3936  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3937  ++row_n)
3938  {
3939  if (row_n != dependent_direction)
3940  {
3941  current_gradient[row_n] =
3942  gradient_entry<spacedim, structdim>(
3943  row_n,
3944  dependent_direction,
3945  trial_point,
3946  guess_weights,
3947  step_size,
3948  get_point_from_weights);
3949 
3950  current_gradient[dependent_direction] -=
3951  current_gradient[row_n];
3952  }
3953  }
3954 
3955  // We need to travel in the -gradient direction, as noted
3956  // above, but we may not want to take a full step in that
3957  // direction; instead, guess that we will go -0.5*gradient and
3958  // do quasi-Newton iteration to pick the best multiplier. The
3959  // goal is to find a scalar alpha such that
3960  //
3961  // F(x - alpha g)
3962  //
3963  // is minimized, where g is the gradient and F is the
3964  // objective function. To find the optimal value we find roots
3965  // of the derivative of the objective function with respect to
3966  // alpha by Newton iteration, where we approximate the first
3967  // and second derivatives of F(x - alpha g) with centered
3968  // finite differences.
3969  double gradient_weight = -0.5;
3970  auto gradient_weight_objective_function =
3971  [&](const double gradient_weight_guess) -> double {
3972  return (trial_point -
3973  get_point_from_weights(guess_weights +
3974  gradient_weight_guess *
3975  current_gradient))
3976  .norm_square();
3977  };
3978 
3979  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3980  {
3981  const double update_numerator = centered_first_difference(
3982  gradient_weight,
3983  step_size,
3984  gradient_weight_objective_function);
3985  const double update_denominator =
3986  centered_second_difference(
3987  gradient_weight,
3988  step_size,
3989  gradient_weight_objective_function);
3990 
3991  // avoid division by zero. Note that we limit the gradient
3992  // weight below
3993  if (std::abs(update_denominator) == 0.0)
3994  break;
3995  gradient_weight =
3996  gradient_weight - update_numerator / update_denominator;
3997 
3998  // Put a fairly lenient bound on the largest possible
3999  // gradient (things tend to be locally flat, so the gradient
4000  // itself is usually small)
4001  if (std::abs(gradient_weight) > 10)
4002  {
4003  gradient_weight = -10.0;
4004  break;
4005  }
4006  }
4007 
4008  // It only makes sense to take convex combinations with weights
4009  // between zero and one. If the update takes us outside of this
4010  // region then rescale the update to stay within the region and
4011  // try again
4012  Tensor<1, n_vertices_per_cell> tentative_weights =
4013  guess_weights + gradient_weight * current_gradient;
4014 
4015  double new_gradient_weight = gradient_weight;
4016  for (unsigned int iteration_count = 0; iteration_count < 40;
4017  ++iteration_count)
4018  {
4019  if (internal::weights_are_ok<structdim>(tentative_weights))
4020  break;
4021 
4022  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4023  {
4024  if (tentative_weights[i] < 0.0)
4025  {
4026  tentative_weights -=
4027  (tentative_weights[i] / current_gradient[i]) *
4028  current_gradient;
4029  }
4030  if (tentative_weights[i] < 0.0 ||
4031  1.0 < tentative_weights[i])
4032  {
4033  new_gradient_weight /= 2.0;
4034  tentative_weights =
4035  guess_weights +
4036  new_gradient_weight * current_gradient;
4037  }
4038  }
4039  }
4040 
4041  // the update might still send us outside the valid region, so
4042  // check again and quit if the update is still not valid
4043  if (!internal::weights_are_ok<structdim>(tentative_weights))
4044  break;
4045 
4046  // if we cannot get closer by traveling in the gradient
4047  // direction then quit
4048  if (get_point_from_weights(tentative_weights)
4049  .distance(trial_point) <
4050  get_point_from_weights(guess_weights).distance(trial_point))
4051  guess_weights = tentative_weights;
4052  else
4053  break;
4054  Assert(internal::weights_are_ok<structdim>(guess_weights),
4055  ExcInternalError());
4056  }
4057  Assert(internal::weights_are_ok<structdim>(guess_weights),
4058  ExcInternalError());
4059  projected_point = get_point_from_weights(guess_weights);
4060  }
4061 
4062  // if structdim == 2 and the optimal point is not on the interior then
4063  // we may be able to get a more accurate result by projecting onto the
4064  // lines.
4065  if (structdim == 2)
4066  {
4067  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4068  line_projections;
4069  for (unsigned int line_n = 0;
4070  line_n < GeometryInfo<structdim>::lines_per_cell;
4071  ++line_n)
4072  {
4073  line_projections[line_n] =
4074  project_to_object(object->line(line_n), trial_point);
4075  }
4076  std::sort(line_projections.begin(),
4077  line_projections.end(),
4078  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4079  return a.distance(trial_point) <
4080  b.distance(trial_point);
4081  });
4082  if (line_projections[0].distance(trial_point) <
4083  projected_point.distance(trial_point))
4084  projected_point = line_projections[0];
4085  }
4086  }
4087  else
4088  {
4089  Assert(false, ExcNotImplemented());
4090  return projected_point;
4091  }
4092 
4093  return projected_point;
4094  }
4095 
4096 
4097 
4098  template <int dim, typename T>
4099  template <class Archive>
4100  void
4102  const unsigned int /*version*/) const
4103  {
4104  Assert(cell_ids.size() == data.size(),
4105  ExcDimensionMismatch(cell_ids.size(), data.size()));
4106  // archive the cellids in an efficient binary format
4107  const std::size_t n_cells = cell_ids.size();
4108  ar & n_cells;
4109  for (const auto &id : cell_ids)
4110  {
4111  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
4112  ar & binary_cell_id;
4113  }
4114 
4115  ar &data;
4116  }
4117 
4118 
4119 
4120  template <int dim, typename T>
4121  template <class Archive>
4122  void
4124  const unsigned int /*version*/)
4125  {
4126  std::size_t n_cells;
4127  ar & n_cells;
4128  cell_ids.clear();
4129  cell_ids.reserve(n_cells);
4130  for (unsigned int c = 0; c < n_cells; ++c)
4131  {
4132  CellId::binary_type value;
4133  ar & value;
4134  cell_ids.emplace_back(value);
4135  }
4136  ar &data;
4137  }
4138 
4139 
4140  namespace internal
4141  {
4142  template <typename DataType,
4143  typename MeshType,
4144  typename MeshCellIteratorType>
4145  inline void
4146  exchange_cell_data(
4147  const MeshType &mesh,
4148  const std::function<
4149  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4150  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4151  & unpack,
4152  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4153  const std::function<void(
4154  const std::function<void(const MeshCellIteratorType &,
4155  const types::subdomain_id)> &)> &process_cells,
4156  const std::function<std::set<types::subdomain_id>(
4157  const parallel::TriangulationBase<MeshType::dimension,
4158  MeshType::space_dimension> &)>
4159  &compute_ghost_owners)
4160  {
4161 # ifndef DEAL_II_WITH_MPI
4162  (void)mesh;
4163  (void)pack;
4164  (void)unpack;
4165  (void)cell_filter;
4166  (void)process_cells;
4167  (void)compute_ghost_owners;
4168  Assert(false, ExcNeedsMPI());
4169 # else
4170  constexpr int dim = MeshType::dimension;
4171  constexpr int spacedim = MeshType::space_dimension;
4172  auto tria =
4173  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4174  &mesh.get_triangulation());
4175  Assert(
4176  tria != nullptr,
4177  ExcMessage(
4178  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4179 
4180  if (const auto tria = dynamic_cast<
4182  &mesh.get_triangulation()))
4183  {
4184  Assert(
4185  tria->with_artificial_cells(),
4186  ExcMessage(
4187  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4188  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4189  "operate on a single layer ghost cells. However, you have "
4190  "given a Triangulation object of type "
4191  "parallel::shared::Triangulation without artificial cells "
4192  "resulting in arbitrary numbers of ghost layers."));
4193  }
4194 
4195  // build list of cells to request for each neighbor
4196  std::set<::types::subdomain_id> ghost_owners =
4197  compute_ghost_owners(*tria);
4198  std::map<::types::subdomain_id,
4199  std::vector<typename CellId::binary_type>>
4200  neighbor_cell_list;
4201 
4202  for (const auto ghost_owner : ghost_owners)
4203  neighbor_cell_list[ghost_owner] = {};
4204 
4205  process_cells([&](const auto &cell, const auto key) {
4206  if (cell_filter(cell))
4207  neighbor_cell_list[key].emplace_back(
4208  cell->id().template to_binary<spacedim>());
4209  });
4210 
4211  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4212  ExcInternalError());
4213 
4214 
4215  // Before sending & receiving, make sure we protect this section with
4216  // a mutex:
4217  static Utilities::MPI::CollectiveMutex mutex;
4219  mutex, tria->get_communicator());
4220 
4221  const int mpi_tag =
4223  const int mpi_tag_reply =
4225 
4226  // send our requests:
4227  std::vector<MPI_Request> requests(ghost_owners.size());
4228  {
4229  unsigned int idx = 0;
4230  for (const auto &it : neighbor_cell_list)
4231  {
4232  // send the data about the relevant cells
4233  const int ierr = MPI_Isend(it.second.data(),
4234  it.second.size() * sizeof(it.second[0]),
4235  MPI_BYTE,
4236  it.first,
4237  mpi_tag,
4238  tria->get_communicator(),
4239  &requests[idx]);
4240  AssertThrowMPI(ierr);
4241  ++idx;
4242  }
4243  }
4244 
4245  using DestinationToBufferMap =
4246  std::map<::types::subdomain_id,
4248  DestinationToBufferMap destination_to_data_buffer_map;
4249 
4250  // receive requests and reply with the ghost indices
4251  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4252  ghost_owners.size());
4253  std::vector<std::vector<types::global_dof_index>>
4254  send_dof_numbers_and_indices(ghost_owners.size());
4255  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4256  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4257 
4258  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4259  {
4260  MPI_Status status;
4261  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4262  mpi_tag,
4263  tria->get_communicator(),
4264  &status);
4265  AssertThrowMPI(ierr);
4266 
4267  int len;
4268  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4269  AssertThrowMPI(ierr);
4270  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4271  ExcInternalError());
4272 
4273  const unsigned int n_cells =
4274  len / sizeof(typename CellId::binary_type);
4275  cell_data_to_send[idx].resize(n_cells);
4276 
4277  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4278  len,
4279  MPI_BYTE,
4280  status.MPI_SOURCE,
4281  status.MPI_TAG,
4282  tria->get_communicator(),
4283  &status);
4284  AssertThrowMPI(ierr);
4285 
4286  // store data for each cell
4287  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4288  {
4289  const auto cell =
4290  tria->create_cell_iterator(CellId(cell_data_to_send[idx][c]));
4291 
4292  MeshCellIteratorType mesh_it(tria,
4293  cell->level(),
4294  cell->index(),
4295  &mesh);
4296  const std_cxx17::optional<DataType> data = pack(mesh_it);
4297 
4298  if (data)
4299  {
4300  typename DestinationToBufferMap::iterator p =
4301  destination_to_data_buffer_map
4302  .insert(std::make_pair(
4303  idx,
4304  GridTools::CellDataTransferBuffer<dim, DataType>()))
4305  .first;
4306 
4307  p->second.cell_ids.emplace_back(cell->id());
4308  p->second.data.emplace_back(*data);
4309  }
4310  }
4311 
4312  // send reply
4313  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4314  destination_to_data_buffer_map[idx];
4315 
4316  sendbuffers[idx] =
4317  Utilities::pack(data, /*enable_compression*/ false);
4318  ierr = MPI_Isend(sendbuffers[idx].data(),
4319  sendbuffers[idx].size(),
4320  MPI_BYTE,
4321  status.MPI_SOURCE,
4322  mpi_tag_reply,
4323  tria->get_communicator(),
4324  &reply_requests[idx]);
4325  AssertThrowMPI(ierr);
4326  }
4327 
4328  // finally receive the replies
4329  std::vector<char> receive;
4330  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4331  {
4332  MPI_Status status;
4333  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4334  mpi_tag_reply,
4335  tria->get_communicator(),
4336  &status);
4337  AssertThrowMPI(ierr);
4338 
4339  int len;
4340  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4341  AssertThrowMPI(ierr);
4342 
4343  receive.resize(len);
4344 
4345  char *ptr = receive.data();
4346  ierr = MPI_Recv(ptr,
4347  len,
4348  MPI_BYTE,
4349  status.MPI_SOURCE,
4350  status.MPI_TAG,
4351  tria->get_communicator(),
4352  &status);
4353  AssertThrowMPI(ierr);
4354 
4355  auto cellinfo =
4356  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4357  receive, /*enable_compression*/ false);
4358 
4359  DataType *data = cellinfo.data.data();
4360  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4361  {
4363  tria_cell = tria->create_cell_iterator(cellinfo.cell_ids[c]);
4364 
4365  MeshCellIteratorType cell(tria,
4366  tria_cell->level(),
4367  tria_cell->index(),
4368  &mesh);
4369 
4370  unpack(cell, *data);
4371  }
4372  }
4373 
4374  // make sure that all communication is finished
4375  // when we leave this function.
4376  if (requests.size() > 0)
4377  {
4378  const int ierr =
4379  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4380  AssertThrowMPI(ierr);
4381  }
4382  if (reply_requests.size() > 0)
4383  {
4384  const int ierr = MPI_Waitall(reply_requests.size(),
4385  reply_requests.data(),
4386  MPI_STATUSES_IGNORE);
4387  AssertThrowMPI(ierr);
4388  }
4389 
4390 
4391 # endif // DEAL_II_WITH_MPI
4392  }
4393 
4394  } // namespace internal
4395 
4396  template <typename DataType, typename MeshType>
4397  inline void
4399  const MeshType & mesh,
4400  const std::function<std_cxx17::optional<DataType>(
4401  const typename MeshType::active_cell_iterator &)> &pack,
4402  const std::function<void(const typename MeshType::active_cell_iterator &,
4403  const DataType &)> & unpack,
4404  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4405  &cell_filter)
4406  {
4407 # ifndef DEAL_II_WITH_MPI
4408  (void)mesh;
4409  (void)pack;
4410  (void)unpack;
4411  (void)cell_filter;
4412  Assert(false, ExcNeedsMPI());
4413 # else
4414  internal::exchange_cell_data<DataType,
4415  MeshType,
4416  typename MeshType::active_cell_iterator>(
4417  mesh,
4418  pack,
4419  unpack,
4420  cell_filter,
4421  [&](const auto &process) {
4422  for (const auto &cell : mesh.active_cell_iterators())
4423  if (cell->is_ghost())
4424  process(cell, cell->subdomain_id());
4425  },
4426  [](const auto &tria) { return tria.ghost_owners(); });
4427 # endif
4428  }
4429 
4430 
4431 
4432  template <typename DataType, typename MeshType>
4433  inline void
4435  const MeshType & mesh,
4436  const std::function<std_cxx17::optional<DataType>(
4437  const typename MeshType::level_cell_iterator &)> &pack,
4438  const std::function<void(const typename MeshType::level_cell_iterator &,
4439  const DataType &)> & unpack,
4440  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4441  &cell_filter)
4442  {
4443 # ifndef DEAL_II_WITH_MPI
4444  (void)mesh;
4445  (void)pack;
4446  (void)unpack;
4447  (void)cell_filter;
4448  Assert(false, ExcNeedsMPI());
4449 # else
4450  internal::exchange_cell_data<DataType,
4451  MeshType,
4452  typename MeshType::level_cell_iterator>(
4453  mesh,
4454  pack,
4455  unpack,
4456  cell_filter,
4457  [&](const auto &process) {
4458  for (const auto &cell : mesh.cell_iterators())
4459  if (cell->level_subdomain_id() !=
4461  !cell->is_locally_owned_on_level())
4462  process(cell, cell->level_subdomain_id());
4463  },
4464  [](const auto &tria) { return tria.level_ghost_owners(); });
4465 # endif
4466  }
4467 } // namespace GridTools
4468 
4469 # endif
4470 
4472 
4473 /*---------------------------- grid_tools.h ---------------------------*/
4474 /* end of #ifndef dealii_grid_tools_h */
4475 #endif
4476 /*---------------------------- grid_tools.h ---------------------------*/
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5049
void map_boundary_to_manifold_ids(const std::vector< types::boundary_id > &src_boundary_ids, const std::vector< types::manifold_id > &dst_manifold_ids, Triangulation< dim, spacedim > &tria, const std::vector< types::boundary_id > &reset_boundary_ids={})
Definition: grid_tools.cc:4834
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:6358
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:213
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=always_return< typename MeshType::level_cell_iterator, bool >{ true})
unsigned int n_vertices() const
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static ::ExceptionBase & ExcNeedsMPI()
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2186
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:4809
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1922
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:81
virtual MPI_Comm get_communicator() const
Definition: tria.cc:10195
static ::ExceptionBase & ExcPointNotFoundInCoarseGrid(Point< N > arg1)
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:312
Contents is actually a matrix.
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5017
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:2935
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< std::set< typename MeshType< dim, spacedim >::active_cell_iterator >> &vertex_to_cell_map, const std::vector< std::vector< Tensor< 1, spacedim >>> &vertex_to_cell_centers, const typename MeshType< dim, spacedim >::active_cell_iterator &cell_hint=typename MeshType< dim, spacedim >::active_cell_iterator(), const std::vector< bool > &marked_vertices={}, const RTree< std::pair< Point< spacedim >, unsigned int >> &used_vertices_rtree=RTree< std::pair< Point< spacedim >, unsigned int >>{}, const double tolerance=1.e-10)
Definition: grid_tools.cc:2760
void load(Archive &ar, const unsigned int version)
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=always_return< typename MeshType::active_cell_iterator, bool >{true})
return_type guess_point_owner(const RTree< std::pair< BoundingBox< spacedim >, unsigned int >> &covering_rtree, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3212
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const double tolerance, const bool perform_handshake)
Definition: grid_tools.cc:5960
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:5111
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2043
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3262
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12025
std::map< unsigned int, types::global_vertex_index > compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3311
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
std::pair< Point< MeshType::space_dimension >, Point< MeshType::space_dimension > > compute_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
std::bitset< 3 > orientation
Definition: grid_tools.h:2655
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:6198
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4095
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:145
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
const double angle
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:409
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:6262
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:626
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3698
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
void partition_triangulation(const unsigned int n_partitions, const std::vector< unsigned int > &cell_weights, const SparsityPattern &cell_connection_graph, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:3909
cell_iterator end() const
Definition: tria.cc:12116
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:534
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:4274
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1085
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:95
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
std::map< types::global_dof_index, std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandler< dim, spacedim > &dof_handler)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
static const char T
T sum(const T &t, const MPI_Comm &mpi_communicator)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3727
bool orthogonal_equality(const FaceIterator &face1, const FaceIterator &face2, const int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
#define Assert(cond, exc)
Definition: exceptions.h:1473
Signals signals
Definition: tria.h:2295
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3224
Abstract base class for mapping classes.
Definition: mapping.h:303
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:4301
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1257
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:697
#define DeclException0(Exception0)
Definition: exceptions.h:470
std::vector< BoundingBox< MeshType::space_dimension > > compute_mesh_predicate_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const unsigned int refinement_level=0, const bool allow_merge=false, const unsigned int max_boxes=numbers::invalid_unsigned_int)
Definition: grid_tools.cc:3017
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
unsigned int level
Definition: grid_out.cc:4590
void save(Archive &ar, const unsigned int version) const
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5082
Point< 3 > vertices[4]
void build_triangulation_from_patch(const std::vector< typename Container::active_cell_iterator > &patch, Triangulation< Container::dimension, Container::space_dimension > &local_triangulation, std::map< typename Triangulation< Container::dimension, Container::space_dimension >::active_cell_iterator, typename Container::active_cell_iterator > &patch_to_global_tria_map)
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:6421
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:288
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:4901
Triangulation< dim, spacedim >::DistortedCellList fix_up_distorted_child_cells(const typename Triangulation< dim, spacedim >::DistortedCellList &distorted_cells, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4767
cell_iterator create_cell_iterator(const CellId &cell_id) const
Definition: tria.cc:12091
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1218
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: hp.h:117
std::vector< std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_active_cells_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const double tolerance=1e-10, const std::vector< bool > &marked_vertices={})
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12655
Definition: cell_id.h:70
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator >> &vertex_to_cells)
Definition: grid_tools.cc:2697
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
unsigned int global_dof_index
Definition: types.h:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
__global__ void set(Number *val, const Number s, const size_type N)
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:5462
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1754
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:446
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition: grid_tools.cc:2218
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:137
std::vector< typename MeshType< dim, spacedim >::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType< dim, spacedim > &container, const unsigned int vertex_index)
Definition: grid_tools.cc:2569
double cell_measure(const T &,...)
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:378
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:4217
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3664
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1321
Point< 3 > center
FullMatrix< double > matrix
Definition: grid_tools.h:2669
void assign_co_dimensional_manifold_indicators(Triangulation< dim, spacedim > &tria, const std::function< types::manifold_id(const std::set< types::manifold_id > &)> &disambiguation_function=[](const std::set< types::manifold_id > &manifold_ids) { if(manifold_ids.size()==1) return *manifold_ids.begin();else return numbers::flat_manifold_id;}, bool overwrite_only_flat_manifold_ids=true)
Definition: grid_tools.cc:4929
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:2032
static ::ExceptionBase & ExcNotImplemented()
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:5426
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4247
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1110
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:3990
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
numbers::NumberTraits< Number >::real_type norm() const
unsigned int find_closest_vertex(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={})
Definition: grid_tools.cc:2507
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
void invert_all_negative_measure_cells(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:863
static ::ExceptionBase & ExcMeshNotOrientable()
static ::ExceptionBase & ExcPointNotFound(Point< N > arg1)
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id, const int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator >> &matched_pairs, const ::Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2023
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_layer_within_distance(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const double layer_thickness)
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4231
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
static ::ExceptionBase & ExcInternalError()
void delete_duplicated_vertices(std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
Definition: grid_tools.cc:731
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6544