Reference documentation for deal.II version Git 1dc1051882 2021-04-22 23:57:03 +0200
\(\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,
233  const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]);
234 
244  template <int dim>
245  double
246  cell_measure(const std::vector<Point<dim>> & all_vertices,
247  const ArrayView<const unsigned int> &vertex_indices);
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 
468 
522  template <int dim, typename Transformation, int spacedim>
523  void
524  transform(const Transformation & transformation,
525  Triangulation<dim, spacedim> &triangulation);
526 
532  template <int dim, int spacedim>
533  void
534  shift(const Tensor<1, spacedim> & shift_vector,
535  Triangulation<dim, spacedim> &triangulation);
536 
537 
547  template <int dim>
548  void
549  rotate(const double angle, Triangulation<dim> &triangulation);
550 
563  template <int dim>
564  void
565  rotate(const double angle,
566  const unsigned int axis,
567  Triangulation<dim, 3> &triangulation);
568 
626  template <int dim>
627  void
628  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
629  Triangulation<dim> & tria,
630  const Function<dim, double> *coefficient = nullptr,
631  const bool solve_for_absolute_positions = false);
632 
638  template <int dim, int spacedim>
639  std::map<unsigned int, Point<spacedim>>
641 
649  template <int dim, int spacedim>
650  void
651  scale(const double scaling_factor,
652  Triangulation<dim, spacedim> &triangulation);
653 
668  template <int dim, int spacedim>
669  void
671  const double factor,
672  Triangulation<dim, spacedim> &triangulation,
673  const bool keep_boundary = true,
674  const unsigned int seed = boost::random::mt19937::default_seed);
675 
709  template <int dim, int spacedim>
710  void
712  const bool isotropic = false,
713  const unsigned int max_iterations = 100);
714 
739  template <int dim, int spacedim>
740  void
742  const double max_ratio = 1.6180339887,
743  const unsigned int max_iterations = 5);
744 
834  template <int dim, int spacedim>
835  void
837  const double limit_angle_fraction = .75);
838 
844 
894  template <int dim, int spacedim>
895 # ifndef DOXYGEN
896  std::tuple<
897  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
898  std::vector<std::vector<Point<dim>>>,
899  std::vector<std::vector<unsigned int>>>
900 # else
901  return_type
902 # endif
904  const Cache<dim, spacedim> & cache,
905  const std::vector<Point<spacedim>> &points,
907  &cell_hint =
909 
938  template <int dim, int spacedim>
939 # ifndef DOXYGEN
940  std::tuple<
941  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
942  std::vector<std::vector<Point<dim>>>,
943  std::vector<std::vector<unsigned int>>,
944  std::vector<unsigned int>>
945 # else
946  return_type
947 # endif
949  const Cache<dim, spacedim> & cache,
950  const std::vector<Point<spacedim>> &points,
952  &cell_hint =
954 
1024  template <int dim, int spacedim>
1025 # ifndef DOXYGEN
1026  std::tuple<
1027  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1028  std::vector<std::vector<Point<dim>>>,
1029  std::vector<std::vector<unsigned int>>,
1030  std::vector<std::vector<Point<spacedim>>>,
1031  std::vector<std::vector<unsigned int>>>
1032 # else
1033  return_type
1034 # endif
1036  const GridTools::Cache<dim, spacedim> & cache,
1037  const std::vector<Point<spacedim>> & local_points,
1038  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1039  const double tolerance = 1e-10);
1040 
1041  namespace internal
1042  {
1057  template <int dim, int spacedim>
1059  {
1066  std::vector<std::tuple<std::pair<int, int>,
1067  unsigned int,
1068  unsigned int,
1069  Point<dim>,
1071  unsigned int>>
1073 
1077  std::vector<unsigned int> send_ranks;
1078 
1084  std::vector<unsigned int> send_ptrs;
1085 
1096  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1098 
1102  std::vector<unsigned int> recv_ranks;
1103 
1109  std::vector<unsigned int> recv_ptrs;
1110  };
1111 
1121  template <int dim, int spacedim>
1124  const GridTools::Cache<dim, spacedim> & cache,
1125  const std::vector<Point<spacedim>> & points,
1126  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1127  const double tolerance,
1128  const bool perform_handshake);
1129 
1130  } // namespace internal
1131 
1168  template <int dim, int spacedim>
1169  std::map<unsigned int, Point<spacedim>>
1171  const Triangulation<dim, spacedim> &container,
1172  const Mapping<dim, spacedim> & mapping =
1173  (ReferenceCells::get_hypercube<dim>()
1174  .template get_default_linear_mapping<dim, spacedim>()));
1175 
1185  template <int spacedim>
1186  unsigned int
1187  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1188  const Point<spacedim> & p);
1189 
1213  template <int dim, template <int, int> class MeshType, int spacedim>
1214  unsigned int
1215  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1216  const Point<spacedim> & p,
1217  const std::vector<bool> & marked_vertices = {});
1218 
1242  template <int dim, template <int, int> class MeshType, int spacedim>
1243  unsigned int
1245  const MeshType<dim, spacedim> &mesh,
1246  const Point<spacedim> & p,
1247  const std::vector<bool> & marked_vertices = {});
1248 
1249 
1274  template <int dim, template <int, int> class MeshType, int spacedim>
1275 # ifndef _MSC_VER
1276  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1277 # else
1278  std::vector<
1279  typename ::internal::
1280  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1281 # endif
1282  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1283  const unsigned int vertex_index);
1284 
1352  template <int dim, template <int, int> class MeshType, int spacedim>
1353 # ifndef _MSC_VER
1354  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1355 # else
1356  std::pair<typename ::internal::
1357  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1358  Point<dim>>
1359 # endif
1361  const MeshType<dim, spacedim> &mesh,
1362  const Point<spacedim> & p,
1363  const std::vector<bool> &marked_vertices = {},
1364  const double tolerance = 1.e-10);
1365 
1373  template <int dim, template <int, int> class MeshType, int spacedim>
1374 # ifndef _MSC_VER
1375  typename MeshType<dim, spacedim>::active_cell_iterator
1376 # else
1377  typename ::internal::
1378  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1379 # endif
1380  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1381  const Point<spacedim> & p,
1382  const std::vector<bool> &marked_vertices = {},
1383  const double tolerance = 1.e-10);
1384 
1391  template <int dim, int spacedim>
1392  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1393  Point<dim>>
1395  const hp::MappingCollection<dim, spacedim> &mapping,
1396  const DoFHandler<dim, spacedim> & mesh,
1397  const Point<spacedim> & p,
1398  const double tolerance = 1.e-10);
1399 
1448  template <int dim, int spacedim>
1449  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1450  Point<dim>>
1452  const Cache<dim, spacedim> &cache,
1453  const Point<spacedim> & p,
1456  const std::vector<bool> &marked_vertices = {},
1457  const double tolerance = 1.e-10);
1458 
1472  template <int dim, template <int, int> class MeshType, int spacedim>
1473 # ifndef _MSC_VER
1474  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1475 # else
1476  std::pair<typename ::internal::
1477  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1478  Point<dim>>
1479 # endif
1481  const Mapping<dim, spacedim> & mapping,
1482  const MeshType<dim, spacedim> &mesh,
1483  const Point<spacedim> & p,
1484  const std::vector<
1485  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1487  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1488  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1489  typename MeshType<dim, spacedim>::active_cell_iterator(),
1490  const std::vector<bool> & marked_vertices = {},
1491  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1492  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1493  const double tolerance = 1.e-10);
1494 
1515  template <int dim, template <int, int> class MeshType, int spacedim>
1516 # ifndef _MSC_VER
1517  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1518  Point<dim>>>
1519 # else
1520  std::vector<std::pair<
1521  typename ::internal::
1522  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1523  Point<dim>>>
1524 # endif
1526  const Mapping<dim, spacedim> & mapping,
1527  const MeshType<dim, spacedim> &mesh,
1528  const Point<spacedim> & p,
1529  const double tolerance,
1530  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1531  Point<dim>> & first_cell);
1532 
1539  template <int dim, template <int, int> class MeshType, int spacedim>
1540 # ifndef _MSC_VER
1541  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1542  Point<dim>>>
1543 # else
1544  std::vector<std::pair<
1545  typename ::internal::
1546  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1547  Point<dim>>>
1548 # endif
1550  const Mapping<dim, spacedim> & mapping,
1551  const MeshType<dim, spacedim> &mesh,
1552  const Point<spacedim> & p,
1553  const double tolerance = 1e-10,
1554  const std::vector<bool> & marked_vertices = {});
1555 
1577  template <class MeshType>
1578  std::vector<typename MeshType::active_cell_iterator>
1579  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1580 
1605  template <class MeshType>
1606  void
1608  const typename MeshType::active_cell_iterator & cell,
1609  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1610 
1660  template <class MeshType>
1661  std::vector<typename MeshType::active_cell_iterator>
1663  const MeshType &mesh,
1664  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1665  &predicate);
1666 
1667 
1675  template <class MeshType>
1676  std::vector<typename MeshType::cell_iterator>
1678  const MeshType &mesh,
1679  const std::function<bool(const typename MeshType::cell_iterator &)>
1680  & predicate,
1681  const unsigned int level);
1682 
1683 
1696  template <class MeshType>
1697  std::vector<typename MeshType::active_cell_iterator>
1698  compute_ghost_cell_halo_layer(const MeshType &mesh);
1699 
1748  template <class MeshType>
1749  std::vector<typename MeshType::active_cell_iterator>
1751  const MeshType &mesh,
1752  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1753  & predicate,
1754  const double layer_thickness);
1755 
1778  template <class MeshType>
1779  std::vector<typename MeshType::active_cell_iterator>
1780  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1781  const double layer_thickness);
1782 
1798  template <class MeshType>
1799  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1801  const MeshType &mesh,
1802  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1803  &predicate);
1804 
1857  template <class MeshType>
1858  std::vector<BoundingBox<MeshType::space_dimension>>
1860  const MeshType &mesh,
1861  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1862  & predicate,
1863  const unsigned int refinement_level = 0,
1864  const bool allow_merge = false,
1865  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1866 
1894  template <int spacedim>
1895 # ifndef DOXYGEN
1896  std::tuple<std::vector<std::vector<unsigned int>>,
1897  std::map<unsigned int, unsigned int>,
1898  std::map<unsigned int, std::vector<unsigned int>>>
1899 # else
1900  return_type
1901 # endif
1903  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1904  const std::vector<Point<spacedim>> & points);
1905 
1906 
1941  template <int spacedim>
1942 # ifndef DOXYGEN
1943  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1944  std::map<unsigned int, unsigned int>,
1945  std::map<unsigned int, std::vector<unsigned int>>>
1946 # else
1947  return_type
1948 # endif
1950  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1951  const std::vector<Point<spacedim>> & points);
1952 
1953 
1962  template <int dim, int spacedim>
1963  std::vector<
1964  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1965  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1966 
1979  template <int dim, int spacedim>
1980  std::vector<std::vector<Tensor<1, spacedim>>>
1982  const Triangulation<dim, spacedim> &mesh,
1983  const std::vector<
1985  &vertex_to_cells);
1986 
1987 
1995  template <int dim, int spacedim>
1996  unsigned int
1999  const Point<spacedim> & position,
2000  const Mapping<dim, spacedim> & mapping =
2001  (ReferenceCells::get_hypercube<dim>()
2002  .template get_default_linear_mapping<dim, spacedim>()));
2003 
2015  template <int dim, int spacedim>
2016  std::map<unsigned int, types::global_vertex_index>
2019 
2031  template <int dim, int spacedim>
2032  std::pair<unsigned int, double>
2035 
2041 
2050  template <int dim, int spacedim>
2051  void
2053  const Triangulation<dim, spacedim> &triangulation,
2054  DynamicSparsityPattern & connectivity);
2055 
2064  template <int dim, int spacedim>
2065  void
2067  const Triangulation<dim, spacedim> &triangulation,
2068  DynamicSparsityPattern & connectivity);
2069 
2078  template <int dim, int spacedim>
2079  void
2081  const Triangulation<dim, spacedim> &triangulation,
2082  const unsigned int level,
2083  DynamicSparsityPattern & connectivity);
2084 
2105  template <int dim, int spacedim>
2106  void
2107  partition_triangulation(const unsigned int n_partitions,
2108  Triangulation<dim, spacedim> & triangulation,
2109  const SparsityTools::Partitioner partitioner =
2111 
2122  template <int dim, int spacedim>
2123  void
2124  partition_triangulation(const unsigned int n_partitions,
2125  const std::vector<unsigned int> &cell_weights,
2126  Triangulation<dim, spacedim> & triangulation,
2127  const SparsityTools::Partitioner partitioner =
2129 
2175  template <int dim, int spacedim>
2176  void
2177  partition_triangulation(const unsigned int n_partitions,
2178  const SparsityPattern & cell_connection_graph,
2179  Triangulation<dim, spacedim> &triangulation,
2180  const SparsityTools::Partitioner partitioner =
2182 
2193  template <int dim, int spacedim>
2194  void
2195  partition_triangulation(const unsigned int n_partitions,
2196  const std::vector<unsigned int> &cell_weights,
2197  const SparsityPattern & cell_connection_graph,
2198  Triangulation<dim, spacedim> &triangulation,
2199  const SparsityTools::Partitioner partitioner =
2201 
2216  template <int dim, int spacedim>
2217  void
2218  partition_triangulation_zorder(const unsigned int n_partitions,
2219  Triangulation<dim, spacedim> &triangulation,
2220  const bool group_siblings = true);
2221 
2233  template <int dim, int spacedim>
2234  void
2236 
2244  template <int dim, int spacedim>
2245  std::vector<types::subdomain_id>
2247  const std::vector<CellId> & cell_ids);
2248 
2259  template <int dim, int spacedim>
2260  void
2262  std::vector<types::subdomain_id> & subdomain);
2263 
2278  template <int dim, int spacedim>
2279  unsigned int
2281  const Triangulation<dim, spacedim> &triangulation,
2282  const types::subdomain_id subdomain);
2283 
2313  template <int dim, int spacedim>
2314  std::vector<bool>
2316 
2322 
2356  template <typename MeshType>
2357  std::list<std::pair<typename MeshType::cell_iterator,
2358  typename MeshType::cell_iterator>>
2359  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2360 
2370  template <int dim, int spacedim>
2371  bool
2373  const Triangulation<dim, spacedim> &mesh_2);
2374 
2384  template <typename MeshType>
2385  bool
2386  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2387 
2393 
2409  template <int dim, int spacedim>
2413  & distorted_cells,
2414  Triangulation<dim, spacedim> &triangulation);
2415 
2416 
2417 
2426 
2427 
2465  template <class MeshType>
2466  std::vector<typename MeshType::active_cell_iterator>
2467  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2468 
2469 
2491  template <class Container>
2492  std::vector<typename Container::cell_iterator>
2494  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2495 
2562  template <class Container>
2563  void
2565  const std::vector<typename Container::active_cell_iterator> &patch,
2567  &local_triangulation,
2568  std::map<
2569  typename Triangulation<Container::dimension,
2570  Container::space_dimension>::active_cell_iterator,
2571  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2572 
2604  template <int dim, int spacedim>
2605  std::map<
2607  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2609 
2610 
2617 
2623  template <typename CellIterator>
2625  {
2629  CellIterator cell[2];
2630 
2635  unsigned int face_idx[2];
2636 
2642  std::bitset<3> orientation;
2643 
2657  };
2658 
2659 
2723  template <typename FaceIterator>
2724  bool orthogonal_equality(
2725  std::bitset<3> & orientation,
2726  const FaceIterator & face1,
2727  const FaceIterator & face2,
2728  const int direction,
2732 
2733 
2737  template <typename FaceIterator>
2738  bool
2740  const FaceIterator & face1,
2741  const FaceIterator & face2,
2742  const int direction,
2746 
2747 
2804  template <typename MeshType>
2805  void
2807  const MeshType & mesh,
2808  const types::boundary_id b_id1,
2809  const types::boundary_id b_id2,
2810  const int direction,
2812  & matched_pairs,
2813  const Tensor<1, MeshType::space_dimension> &offset =
2816 
2817 
2840  template <typename MeshType>
2841  void
2843  const MeshType & mesh,
2844  const types::boundary_id b_id,
2845  const int direction,
2847  & matched_pairs,
2848  const ::Tensor<1, MeshType::space_dimension> &offset =
2851 
2857 
2878  template <int dim, int spacedim>
2879  void
2881  const bool reset_boundary_ids = false);
2882 
2904  template <int dim, int spacedim>
2905  void
2907  const std::vector<types::boundary_id> &src_boundary_ids,
2908  const std::vector<types::manifold_id> &dst_manifold_ids,
2910  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2911 
2941  template <int dim, int spacedim>
2942  void
2944  const bool compute_face_ids = false);
2945 
2970  template <int dim, int spacedim>
2971  void
2974  const std::function<types::manifold_id(
2975  const std::set<types::manifold_id> &)> &disambiguation_function =
2976  [](const std::set<types::manifold_id> &manifold_ids) {
2977  if (manifold_ids.size() == 1)
2978  return *manifold_ids.begin();
2979  else
2981  },
2982  bool overwrite_only_flat_manifold_ids = true);
3069  template <typename DataType, typename MeshType>
3070  void
3072  const MeshType & mesh,
3073  const std::function<std_cxx17::optional<DataType>(
3074  const typename MeshType::active_cell_iterator &)> &pack,
3075  const std::function<void(const typename MeshType::active_cell_iterator &,
3076  const DataType &)> & unpack,
3077  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3078  &cell_filter =
3080 
3091  template <typename DataType, typename MeshType>
3092  void
3094  const MeshType & mesh,
3095  const std::function<std_cxx17::optional<DataType>(
3096  const typename MeshType::level_cell_iterator &)> &pack,
3097  const std::function<void(const typename MeshType::level_cell_iterator &,
3098  const DataType &)> & unpack,
3099  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3101  true});
3102 
3103  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3104  * boxes @p local_bboxes.
3105  *
3106  * This function is meant to exchange bounding boxes describing the locally
3107  * owned cells in a distributed triangulation obtained with the function
3108  * GridTools::compute_mesh_predicate_bounding_box .
3109  *
3110  * The output vector's size is the number of processes of the MPI
3111  * communicator:
3112  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3113  */
3114  template <int spacedim>
3115  std::vector<std::vector<BoundingBox<spacedim>>>
3117  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3118  const MPI_Comm & mpi_communicator);
3119 
3152  template <int spacedim>
3155  const std::vector<BoundingBox<spacedim>> &local_description,
3156  const MPI_Comm & mpi_communicator);
3157 
3175  template <int dim, int spacedim>
3176  void
3178  const Triangulation<dim, spacedim> & tria,
3179  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3180  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3181 
3188  template <int dim, int spacedim>
3189  std::map<unsigned int, std::set<::types::subdomain_id>>
3191  const Triangulation<dim, spacedim> &tria);
3192 
3205  template <int dim, typename T>
3207  {
3211  std::vector<CellId> cell_ids;
3212 
3216  std::vector<T> data;
3217 
3226  template <class Archive>
3227  void
3228  save(Archive &ar, const unsigned int version) const;
3229 
3236  template <class Archive>
3237  void
3238  load(Archive &ar, const unsigned int version);
3239 
3240 # ifdef DOXYGEN
3241 
3246  template <class Archive>
3247  void
3248  serialize(Archive &archive, const unsigned int version);
3249 # else
3250  // This macro defines the serialize() method that is compatible with
3251  // the templated save() and load() method that have been implemented.
3252  BOOST_SERIALIZATION_SPLIT_MEMBER()
3253 # endif
3254  };
3255 
3260 
3265  int,
3266  << "The number of partitions you gave is " << arg1
3267  << ", but must be greater than zero.");
3272  int,
3273  << "The subdomain id " << arg1
3274  << " has no cells associated with it.");
3279 
3284  double,
3285  << "The scaling factor must be positive, but it is " << arg1
3286  << ".");
3290  template <int N>
3292  Point<N>,
3293  << "The point <" << arg1
3294  << "> could not be found inside any of the "
3295  << "coarse grid cells.");
3299  template <int N>
3301  Point<N>,
3302  << "The point <" << arg1
3303  << "> could not be found inside any of the "
3304  << "subcells of a coarse grid cell.");
3305 
3310  unsigned int,
3311  << "The given vertex with index " << arg1
3312  << " is not used in the given triangulation.");
3313 
3314 
3317 } /*namespace GridTools*/
3318 
3319 
3320 
3321 /* ----------------- Template function --------------- */
3322 
3323 # ifndef DOXYGEN
3324 
3325 namespace GridTools
3326 {
3327  template <int dim>
3328  double
3329  cell_measure(
3330  const std::vector<Point<dim>> &all_vertices,
3331  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3332  {
3333  // We forward call to the ArrayView version:
3334  const ArrayView<const unsigned int> view(
3335  indices, GeometryInfo<dim>::vertices_per_cell);
3336  return cell_measure(all_vertices, view);
3337  }
3338 
3339  template <int dim, typename T>
3340  double
3341  cell_measure(const T &, ...)
3342  {
3343  Assert(false, ExcNotImplemented());
3344  return std::numeric_limits<double>::quiet_NaN();
3345  }
3346 
3347 
3348 
3349  // This specialization is defined here so that the general template in the
3350  // source file doesn't need to have further 1D overloads for the internal
3351  // functions it calls.
3352  template <>
3356  {
3357  return {};
3358  }
3359 
3360 
3361 
3362  template <int dim, typename Predicate, int spacedim>
3363  void
3364  transform(const Predicate & predicate,
3365  Triangulation<dim, spacedim> &triangulation)
3366  {
3367  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3368 
3369  // loop over all active cells, and
3370  // transform those vertices that
3371  // have not yet been touched. note
3372  // that we get to all vertices in
3373  // the triangulation by only
3374  // visiting the active cells.
3376  cell = triangulation.begin_active(),
3377  endc = triangulation.end();
3378  for (; cell != endc; ++cell)
3379  for (const unsigned int v : cell->vertex_indices())
3380  if (treated_vertices[cell->vertex_index(v)] == false)
3381  {
3382  // transform this vertex
3383  cell->vertex(v) = predicate(cell->vertex(v));
3384  // and mark it as treated
3385  treated_vertices[cell->vertex_index(v)] = true;
3386  };
3387 
3388 
3389  // now fix any vertices on hanging nodes so that we don't create any holes
3390  if (dim == 2)
3391  {
3393  cell = triangulation.begin_active(),
3394  endc = triangulation.end();
3395  for (; cell != endc; ++cell)
3396  for (const unsigned int face : cell->face_indices())
3397  if (cell->face(face)->has_children() &&
3398  !cell->face(face)->at_boundary())
3399  {
3400  Assert(cell->reference_cell() ==
3401  ReferenceCells::get_hypercube<dim>(),
3402  ExcNotImplemented());
3403 
3404  // this line has children
3405  cell->face(face)->child(0)->vertex(1) =
3406  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3407  2;
3408  }
3409  }
3410  else if (dim == 3)
3411  {
3413  cell = triangulation.begin_active(),
3414  endc = triangulation.end();
3415  for (; cell != endc; ++cell)
3416  for (const unsigned int face : cell->face_indices())
3417  if (cell->face(face)->has_children() &&
3418  !cell->face(face)->at_boundary())
3419  {
3420  Assert(cell->reference_cell() ==
3421  ReferenceCells::get_hypercube<dim>(),
3422  ExcNotImplemented());
3423 
3424  // this face has hanging nodes
3425  cell->face(face)->child(0)->vertex(1) =
3426  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3427  2.0;
3428  cell->face(face)->child(0)->vertex(2) =
3429  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3430  2.0;
3431  cell->face(face)->child(1)->vertex(3) =
3432  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3433  2.0;
3434  cell->face(face)->child(2)->vertex(3) =
3435  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3436  2.0;
3437 
3438  // center of the face
3439  cell->face(face)->child(0)->vertex(3) =
3440  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3441  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3442  4.0;
3443  }
3444  }
3445 
3446  // Make sure FEValues notices that the mesh has changed
3447  triangulation.signals.mesh_movement();
3448  }
3449 
3450 
3451 
3452  template <class MeshType>
3453  std::vector<typename MeshType::active_cell_iterator>
3454  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3455  {
3456  std::vector<typename MeshType::active_cell_iterator> child_cells;
3457 
3458  if (cell->has_children())
3459  {
3460  for (unsigned int child = 0; child < cell->n_children(); ++child)
3461  if (cell->child(child)->has_children())
3462  {
3463  const std::vector<typename MeshType::active_cell_iterator>
3464  children = get_active_child_cells<MeshType>(cell->child(child));
3465  child_cells.insert(child_cells.end(),
3466  children.begin(),
3467  children.end());
3468  }
3469  else
3470  child_cells.push_back(cell->child(child));
3471  }
3472 
3473  return child_cells;
3474  }
3475 
3476 
3477 
3478  template <class MeshType>
3479  void
3481  const typename MeshType::active_cell_iterator & cell,
3482  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3483  {
3484  active_neighbors.clear();
3485  for (const unsigned int n : cell->face_indices())
3486  if (!cell->at_boundary(n))
3487  {
3488  if (MeshType::dimension == 1)
3489  {
3490  // check children of neighbor. note
3491  // that in 1d children of the neighbor
3492  // may be further refined. In 1d the
3493  // case is simple since we know what
3494  // children bound to the present cell
3495  typename MeshType::cell_iterator neighbor_child =
3496  cell->neighbor(n);
3497  if (!neighbor_child->is_active())
3498  {
3499  while (neighbor_child->has_children())
3500  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3501 
3502  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3503  ExcInternalError());
3504  }
3505  active_neighbors.push_back(neighbor_child);
3506  }
3507  else
3508  {
3509  if (cell->face(n)->has_children())
3510  // this neighbor has children. find
3511  // out which border to the present
3512  // cell
3513  for (unsigned int c = 0;
3514  c < cell->face(n)->n_active_descendants();
3515  ++c)
3516  active_neighbors.push_back(
3517  cell->neighbor_child_on_subface(n, c));
3518  else
3519  {
3520  // the neighbor must be active
3521  // himself
3522  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3523  active_neighbors.push_back(cell->neighbor(n));
3524  }
3525  }
3526  }
3527  }
3528 
3529 
3530 
3531  namespace internal
3532  {
3533  namespace ProjectToObject
3534  {
3547  struct CrossDerivative
3548  {
3549  const unsigned int direction_0;
3550  const unsigned int direction_1;
3551 
3552  CrossDerivative(const unsigned int d0, const unsigned int d1);
3553  };
3554 
3555  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3556  const unsigned int d1)
3557  : direction_0(d0)
3558  , direction_1(d1)
3559  {}
3560 
3561 
3562 
3567  template <typename F>
3568  inline auto
3569  centered_first_difference(const double center,
3570  const double step,
3571  const F &f) -> decltype(f(center) - f(center))
3572  {
3573  return (f(center + step) - f(center - step)) / (2.0 * step);
3574  }
3575 
3576 
3577 
3582  template <typename F>
3583  inline auto
3584  centered_second_difference(const double center,
3585  const double step,
3586  const F &f) -> decltype(f(center) - f(center))
3587  {
3588  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3589  (step * step);
3590  }
3591 
3592 
3593 
3603  template <int structdim, typename F>
3604  inline auto
3605  cross_stencil(
3606  const CrossDerivative cross_derivative,
3608  const double step,
3609  const F &f) -> decltype(f(center) - f(center))
3610  {
3612  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3613  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3614  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3615  1.0 / 3.0 * f(center - simplex_vector) +
3616  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3617  step;
3618  }
3619 
3620 
3621 
3628  template <int spacedim, int structdim, typename F>
3629  inline double
3630  gradient_entry(
3631  const unsigned int row_n,
3632  const unsigned int dependent_direction,
3633  const Point<spacedim> &p0,
3635  const double step,
3636  const F & f)
3637  {
3639  dependent_direction <
3641  ExcMessage("This function assumes that the last weight is a "
3642  "dependent variable (and hence we cannot take its "
3643  "derivative directly)."));
3644  Assert(row_n != dependent_direction,
3645  ExcMessage(
3646  "We cannot differentiate with respect to the variable "
3647  "that is assumed to be dependent."));
3648 
3649  const Point<spacedim> manifold_point = f(center);
3650  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3651  {row_n, dependent_direction}, center, step, f);
3652  double entry = 0.0;
3653  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3654  entry +=
3655  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3656  return entry;
3657  }
3658 
3664  template <typename Iterator, int spacedim, int structdim>
3666  project_to_d_linear_object(const Iterator & object,
3667  const Point<spacedim> &trial_point)
3668  {
3669  // let's look at this for simplicity for a quad (structdim==2) in a
3670  // space with spacedim>2 (notate trial_point by y): all points on the
3671  // surface are given by
3672  // x(\xi) = sum_i v_i phi_x(\xi)
3673  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3674  // reference coordinates of the quad. so what we are trying to do is
3675  // find a point x on the surface that is closest to the point y. there
3676  // are different ways to solve this problem, but in the end it's a
3677  // nonlinear problem and we have to find reference coordinates \xi so
3678  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3679  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3680  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3681  // have to use a Newton method to find the answer. This leads to the
3682  // following formulation of Newton steps:
3683  //
3684  // Given \xi_k, find \delta\xi_k so that
3685  // H_k \delta\xi_k = - F_k
3686  // where H_k is an approximation to the second derivatives of J at
3687  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3688  // number of times until the right hand side is small enough. As a
3689  // stopping criterion, we terminate if ||\delta\xi||<eps.
3690  //
3691  // As for the Hessian, the best choice would be
3692  // H_k = J''(\xi_k)
3693  // but we'll opt for the simpler Gauss-Newton form
3694  // H_k = A^T A
3695  // i.e.
3696  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3697  // \partial_n phi_i *
3698  // \partial_m phi_j
3699  // we start at xi=(0.5, 0.5).
3700  Point<structdim> xi;
3701  for (unsigned int d = 0; d < structdim; ++d)
3702  xi[d] = 0.5;
3703 
3704  Point<spacedim> x_k;
3705  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3706  x_k += object->vertex(i) *
3708 
3709  do
3710  {
3712  for (const unsigned int i :
3714  F_k +=
3715  (x_k - trial_point) * object->vertex(i) *
3717  i);
3718 
3720  for (const unsigned int i :
3722  for (const unsigned int j :
3724  {
3727  xi, i),
3729  xi, j));
3730  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3731  }
3732 
3733  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3734  xi += delta_xi;
3735 
3736  x_k = Point<spacedim>();
3737  for (const unsigned int i :
3739  x_k += object->vertex(i) *
3741 
3742  if (delta_xi.norm() < 1e-7)
3743  break;
3744  }
3745  while (true);
3746 
3747  return x_k;
3748  }
3749  } // namespace ProjectToObject
3750  } // namespace internal
3751 
3752 
3753 
3754  namespace internal
3755  {
3756  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3757  // inside the project_to_object function below.
3758  template <int structdim>
3759  inline bool
3760  weights_are_ok(
3762  {
3763  // clang has trouble figuring out structdim here, so define it
3764  // again:
3765  static const std::size_t n_vertices_per_cell =
3767  n_independent_components;
3768  std::array<double, n_vertices_per_cell> copied_weights;
3769  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3770  {
3771  copied_weights[i] = v[i];
3772  if (v[i] < 0.0 || v[i] > 1.0)
3773  return false;
3774  }
3775 
3776  // check the sum: try to avoid some roundoff errors by summing in order
3777  std::sort(copied_weights.begin(), copied_weights.end());
3778  const double sum =
3779  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3780  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3781  }
3782  } // namespace internal
3783 
3784  template <typename Iterator>
3787  const Iterator & object,
3789  {
3790  const int spacedim = Iterator::AccessorType::space_dimension;
3791  const int structdim = Iterator::AccessorType::structure_dimension;
3792 
3793  Point<spacedim> projected_point = trial_point;
3794 
3795  if (structdim >= spacedim)
3796  return projected_point;
3797  else if (structdim == 1 || structdim == 2)
3798  {
3799  using namespace internal::ProjectToObject;
3800  // Try to use the special flat algorithm for quads (this is better
3801  // than the general algorithm in 3D). This does not take into account
3802  // whether projected_point is outside the quad, but we optimize along
3803  // lines below anyway:
3804  const int dim = Iterator::AccessorType::dimension;
3805  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3806  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3807  &manifold) != nullptr)
3808  {
3809  projected_point =
3810  project_to_d_linear_object<Iterator, spacedim, structdim>(
3811  object, trial_point);
3812  }
3813  else
3814  {
3815  // We want to find a point on the convex hull (defined by the
3816  // vertices of the object and the manifold description) that is
3817  // relatively close to the trial point. This has a few issues:
3818  //
3819  // 1. For a general convex hull we are not guaranteed that a unique
3820  // minimum exists.
3821  // 2. The independent variables in the optimization process are the
3822  // weights given to Manifold::get_new_point, which must sum to 1,
3823  // so we cannot use standard finite differences to approximate a
3824  // gradient.
3825  //
3826  // There is not much we can do about 1., but for 2. we can derive
3827  // finite difference stencils that work on a structdim-dimensional
3828  // simplex and rewrite the optimization problem to use those
3829  // instead. Consider the structdim 2 case and let
3830  //
3831  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3832  // c2, c3})
3833  //
3834  // where {c0, c1, c2, c3} are the weights for the four vertices on
3835  // the quadrilateral. We seek to minimize the Euclidean distance
3836  // between F(...) and trial_point. We can solve for c3 in terms of
3837  // the other weights and get, for one coordinate direction
3838  //
3839  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3840  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3841  //
3842  // where we substitute back in for c3 after taking the
3843  // derivative. We can compute a stencil for the cross derivative
3844  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3845  // (and gradient_entry computes the sum over the independent
3846  // variables). Below, we somewhat arbitrarily pick the last
3847  // component as the dependent one.
3848  //
3849  // Since we can now calculate derivatives of the objective
3850  // function we can use gradient descent to minimize it.
3851  //
3852  // Of course, this is much simpler in the structdim = 1 case (we
3853  // could rewrite the projection as a 1D optimization problem), but
3854  // to reduce the potential for bugs we use the same code in both
3855  // cases.
3856  const double step_size = object->diameter() / 64.0;
3857 
3858  constexpr unsigned int n_vertices_per_cell =
3860 
3861  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3862  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3863  ++vertex_n)
3864  vertices[vertex_n] = object->vertex(vertex_n);
3865 
3866  auto get_point_from_weights =
3867  [&](const Tensor<1, n_vertices_per_cell> &weights)
3868  -> Point<spacedim> {
3869  return object->get_manifold().get_new_point(
3870  make_array_view(vertices.begin(), vertices.end()),
3871  make_array_view(weights.begin_raw(), weights.end_raw()));
3872  };
3873 
3874  // pick the initial weights as (normalized) inverse distances from
3875  // the trial point:
3876  Tensor<1, n_vertices_per_cell> guess_weights;
3877  double guess_weights_sum = 0.0;
3878  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3879  ++vertex_n)
3880  {
3881  const double distance =
3882  vertices[vertex_n].distance(trial_point);
3883  if (distance == 0.0)
3884  {
3885  guess_weights = 0.0;
3886  guess_weights[vertex_n] = 1.0;
3887  guess_weights_sum = 1.0;
3888  break;
3889  }
3890  else
3891  {
3892  guess_weights[vertex_n] = 1.0 / distance;
3893  guess_weights_sum += guess_weights[vertex_n];
3894  }
3895  }
3896  guess_weights /= guess_weights_sum;
3897  Assert(internal::weights_are_ok<structdim>(guess_weights),
3898  ExcInternalError());
3899 
3900  // The optimization algorithm consists of two parts:
3901  //
3902  // 1. An outer loop where we apply the gradient descent algorithm.
3903  // 2. An inner loop where we do a line search to find the optimal
3904  // length of the step one should take in the gradient direction.
3905  //
3906  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3907  {
3908  const unsigned int dependent_direction =
3909  n_vertices_per_cell - 1;
3910  Tensor<1, n_vertices_per_cell> current_gradient;
3911  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3912  ++row_n)
3913  {
3914  if (row_n != dependent_direction)
3915  {
3916  current_gradient[row_n] =
3917  gradient_entry<spacedim, structdim>(
3918  row_n,
3919  dependent_direction,
3920  trial_point,
3921  guess_weights,
3922  step_size,
3923  get_point_from_weights);
3924 
3925  current_gradient[dependent_direction] -=
3926  current_gradient[row_n];
3927  }
3928  }
3929 
3930  // We need to travel in the -gradient direction, as noted
3931  // above, but we may not want to take a full step in that
3932  // direction; instead, guess that we will go -0.5*gradient and
3933  // do quasi-Newton iteration to pick the best multiplier. The
3934  // goal is to find a scalar alpha such that
3935  //
3936  // F(x - alpha g)
3937  //
3938  // is minimized, where g is the gradient and F is the
3939  // objective function. To find the optimal value we find roots
3940  // of the derivative of the objective function with respect to
3941  // alpha by Newton iteration, where we approximate the first
3942  // and second derivatives of F(x - alpha g) with centered
3943  // finite differences.
3944  double gradient_weight = -0.5;
3945  auto gradient_weight_objective_function =
3946  [&](const double gradient_weight_guess) -> double {
3947  return (trial_point -
3948  get_point_from_weights(guess_weights +
3949  gradient_weight_guess *
3950  current_gradient))
3951  .norm_square();
3952  };
3953 
3954  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3955  {
3956  const double update_numerator = centered_first_difference(
3957  gradient_weight,
3958  step_size,
3959  gradient_weight_objective_function);
3960  const double update_denominator =
3961  centered_second_difference(
3962  gradient_weight,
3963  step_size,
3964  gradient_weight_objective_function);
3965 
3966  // avoid division by zero. Note that we limit the gradient
3967  // weight below
3968  if (std::abs(update_denominator) == 0.0)
3969  break;
3970  gradient_weight =
3971  gradient_weight - update_numerator / update_denominator;
3972 
3973  // Put a fairly lenient bound on the largest possible
3974  // gradient (things tend to be locally flat, so the gradient
3975  // itself is usually small)
3976  if (std::abs(gradient_weight) > 10)
3977  {
3978  gradient_weight = -10.0;
3979  break;
3980  }
3981  }
3982 
3983  // It only makes sense to take convex combinations with weights
3984  // between zero and one. If the update takes us outside of this
3985  // region then rescale the update to stay within the region and
3986  // try again
3987  Tensor<1, n_vertices_per_cell> tentative_weights =
3988  guess_weights + gradient_weight * current_gradient;
3989 
3990  double new_gradient_weight = gradient_weight;
3991  for (unsigned int iteration_count = 0; iteration_count < 40;
3992  ++iteration_count)
3993  {
3994  if (internal::weights_are_ok<structdim>(tentative_weights))
3995  break;
3996 
3997  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3998  {
3999  if (tentative_weights[i] < 0.0)
4000  {
4001  tentative_weights -=
4002  (tentative_weights[i] / current_gradient[i]) *
4003  current_gradient;
4004  }
4005  if (tentative_weights[i] < 0.0 ||
4006  1.0 < tentative_weights[i])
4007  {
4008  new_gradient_weight /= 2.0;
4009  tentative_weights =
4010  guess_weights +
4011  new_gradient_weight * current_gradient;
4012  }
4013  }
4014  }
4015 
4016  // the update might still send us outside the valid region, so
4017  // check again and quit if the update is still not valid
4018  if (!internal::weights_are_ok<structdim>(tentative_weights))
4019  break;
4020 
4021  // if we cannot get closer by traveling in the gradient
4022  // direction then quit
4023  if (get_point_from_weights(tentative_weights)
4024  .distance(trial_point) <
4025  get_point_from_weights(guess_weights).distance(trial_point))
4026  guess_weights = tentative_weights;
4027  else
4028  break;
4029  Assert(internal::weights_are_ok<structdim>(guess_weights),
4030  ExcInternalError());
4031  }
4032  Assert(internal::weights_are_ok<structdim>(guess_weights),
4033  ExcInternalError());
4034  projected_point = get_point_from_weights(guess_weights);
4035  }
4036 
4037  // if structdim == 2 and the optimal point is not on the interior then
4038  // we may be able to get a more accurate result by projecting onto the
4039  // lines.
4040  if (structdim == 2)
4041  {
4042  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4043  line_projections;
4044  for (unsigned int line_n = 0;
4045  line_n < GeometryInfo<structdim>::lines_per_cell;
4046  ++line_n)
4047  {
4048  line_projections[line_n] =
4049  project_to_object(object->line(line_n), trial_point);
4050  }
4051  std::sort(line_projections.begin(),
4052  line_projections.end(),
4053  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4054  return a.distance(trial_point) <
4055  b.distance(trial_point);
4056  });
4057  if (line_projections[0].distance(trial_point) <
4058  projected_point.distance(trial_point))
4059  projected_point = line_projections[0];
4060  }
4061  }
4062  else
4063  {
4064  Assert(false, ExcNotImplemented());
4065  return projected_point;
4066  }
4067 
4068  return projected_point;
4069  }
4070 
4071 
4072 
4073  template <int dim, typename T>
4074  template <class Archive>
4075  void
4077  const unsigned int /*version*/) const
4078  {
4079  Assert(cell_ids.size() == data.size(),
4080  ExcDimensionMismatch(cell_ids.size(), data.size()));
4081  // archive the cellids in an efficient binary format
4082  const std::size_t n_cells = cell_ids.size();
4083  ar & n_cells;
4084  for (const auto &id : cell_ids)
4085  {
4086  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
4087  ar & binary_cell_id;
4088  }
4089 
4090  ar &data;
4091  }
4092 
4093 
4094 
4095  template <int dim, typename T>
4096  template <class Archive>
4097  void
4099  const unsigned int /*version*/)
4100  {
4101  std::size_t n_cells;
4102  ar & n_cells;
4103  cell_ids.clear();
4104  cell_ids.reserve(n_cells);
4105  for (unsigned int c = 0; c < n_cells; ++c)
4106  {
4107  CellId::binary_type value;
4108  ar & value;
4109  cell_ids.emplace_back(value);
4110  }
4111  ar &data;
4112  }
4113 
4114 
4115  namespace internal
4116  {
4117  template <typename DataType,
4118  typename MeshType,
4119  typename MeshCellIteratorType>
4120  inline void
4121  exchange_cell_data(
4122  const MeshType &mesh,
4123  const std::function<
4124  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4125  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4126  & unpack,
4127  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4128  const std::function<void(
4129  const std::function<void(const MeshCellIteratorType &,
4130  const types::subdomain_id)> &)> &process_cells,
4131  const std::function<std::set<types::subdomain_id>(
4132  const parallel::TriangulationBase<MeshType::dimension,
4133  MeshType::space_dimension> &)>
4134  &compute_ghost_owners)
4135  {
4136 # ifndef DEAL_II_WITH_MPI
4137  (void)mesh;
4138  (void)pack;
4139  (void)unpack;
4140  (void)cell_filter;
4141  (void)process_cells;
4142  (void)compute_ghost_owners;
4143  Assert(false, ExcNeedsMPI());
4144 # else
4145  constexpr int dim = MeshType::dimension;
4146  constexpr int spacedim = MeshType::space_dimension;
4147  auto tria =
4148  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4149  &mesh.get_triangulation());
4150  Assert(
4151  tria != nullptr,
4152  ExcMessage(
4153  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4154 
4155  if (const auto tria = dynamic_cast<
4157  &mesh.get_triangulation()))
4158  {
4159  Assert(
4160  tria->with_artificial_cells(),
4161  ExcMessage(
4162  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4163  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4164  "operate on a single layer ghost cells. However, you have "
4165  "given a Triangulation object of type "
4166  "parallel::shared::Triangulation without artificial cells "
4167  "resulting in arbitrary numbers of ghost layers."));
4168  }
4169 
4170  // build list of cells to request for each neighbor
4171  std::set<::types::subdomain_id> ghost_owners =
4172  compute_ghost_owners(*tria);
4173  std::map<::types::subdomain_id,
4174  std::vector<typename CellId::binary_type>>
4175  neighbor_cell_list;
4176 
4177  for (const auto ghost_owner : ghost_owners)
4178  neighbor_cell_list[ghost_owner] = {};
4179 
4180  process_cells([&](const auto &cell, const auto key) {
4181  if (cell_filter(cell))
4182  neighbor_cell_list[key].emplace_back(
4183  cell->id().template to_binary<spacedim>());
4184  });
4185 
4186  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4187  ExcInternalError());
4188 
4189 
4190  // Before sending & receiving, make sure we protect this section with
4191  // a mutex:
4192  static Utilities::MPI::CollectiveMutex mutex;
4194  mutex, tria->get_communicator());
4195 
4196  const int mpi_tag =
4198  const int mpi_tag_reply =
4200 
4201  // send our requests:
4202  std::vector<MPI_Request> requests(ghost_owners.size());
4203  {
4204  unsigned int idx = 0;
4205  for (const auto &it : neighbor_cell_list)
4206  {
4207  // send the data about the relevant cells
4208  const int ierr = MPI_Isend(it.second.data(),
4209  it.second.size() * sizeof(it.second[0]),
4210  MPI_BYTE,
4211  it.first,
4212  mpi_tag,
4213  tria->get_communicator(),
4214  &requests[idx]);
4215  AssertThrowMPI(ierr);
4216  ++idx;
4217  }
4218  }
4219 
4220  using DestinationToBufferMap =
4221  std::map<::types::subdomain_id,
4223  DestinationToBufferMap destination_to_data_buffer_map;
4224 
4225  // receive requests and reply with the ghost indices
4226  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4227  ghost_owners.size());
4228  std::vector<std::vector<types::global_dof_index>>
4229  send_dof_numbers_and_indices(ghost_owners.size());
4230  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4231  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4232 
4233  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4234  {
4235  MPI_Status status;
4236  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4237  mpi_tag,
4238  tria->get_communicator(),
4239  &status);
4240  AssertThrowMPI(ierr);
4241 
4242  int len;
4243  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4244  AssertThrowMPI(ierr);
4245  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4246  ExcInternalError());
4247 
4248  const unsigned int n_cells =
4249  len / sizeof(typename CellId::binary_type);
4250  cell_data_to_send[idx].resize(n_cells);
4251 
4252  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4253  len,
4254  MPI_BYTE,
4255  status.MPI_SOURCE,
4256  status.MPI_TAG,
4257  tria->get_communicator(),
4258  &status);
4259  AssertThrowMPI(ierr);
4260 
4261  // store data for each cell
4262  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4263  {
4264  const auto cell =
4265  tria->create_cell_iterator(CellId(cell_data_to_send[idx][c]));
4266 
4267  MeshCellIteratorType mesh_it(tria,
4268  cell->level(),
4269  cell->index(),
4270  &mesh);
4271  const std_cxx17::optional<DataType> data = pack(mesh_it);
4272 
4273  if (data)
4274  {
4275  typename DestinationToBufferMap::iterator p =
4276  destination_to_data_buffer_map
4277  .insert(std::make_pair(
4278  idx,
4279  GridTools::CellDataTransferBuffer<dim, DataType>()))
4280  .first;
4281 
4282  p->second.cell_ids.emplace_back(cell->id());
4283  p->second.data.emplace_back(*data);
4284  }
4285  }
4286 
4287  // send reply
4288  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4289  destination_to_data_buffer_map[idx];
4290 
4291  sendbuffers[idx] =
4292  Utilities::pack(data, /*enable_compression*/ false);
4293  ierr = MPI_Isend(sendbuffers[idx].data(),
4294  sendbuffers[idx].size(),
4295  MPI_BYTE,
4296  status.MPI_SOURCE,
4297  mpi_tag_reply,
4298  tria->get_communicator(),
4299  &reply_requests[idx]);
4300  AssertThrowMPI(ierr);
4301  }
4302 
4303  // finally receive the replies
4304  std::vector<char> receive;
4305  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4306  {
4307  MPI_Status status;
4308  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4309  mpi_tag_reply,
4310  tria->get_communicator(),
4311  &status);
4312  AssertThrowMPI(ierr);
4313 
4314  int len;
4315  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4316  AssertThrowMPI(ierr);
4317 
4318  receive.resize(len);
4319 
4320  char *ptr = receive.data();
4321  ierr = MPI_Recv(ptr,
4322  len,
4323  MPI_BYTE,
4324  status.MPI_SOURCE,
4325  status.MPI_TAG,
4326  tria->get_communicator(),
4327  &status);
4328  AssertThrowMPI(ierr);
4329 
4330  auto cellinfo =
4331  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4332  receive, /*enable_compression*/ false);
4333 
4334  DataType *data = cellinfo.data.data();
4335  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4336  {
4338  tria_cell = tria->create_cell_iterator(cellinfo.cell_ids[c]);
4339 
4340  MeshCellIteratorType cell(tria,
4341  tria_cell->level(),
4342  tria_cell->index(),
4343  &mesh);
4344 
4345  unpack(cell, *data);
4346  }
4347  }
4348 
4349  // make sure that all communication is finished
4350  // when we leave this function.
4351  if (requests.size() > 0)
4352  {
4353  const int ierr =
4354  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4355  AssertThrowMPI(ierr);
4356  }
4357  if (reply_requests.size() > 0)
4358  {
4359  const int ierr = MPI_Waitall(reply_requests.size(),
4360  reply_requests.data(),
4361  MPI_STATUSES_IGNORE);
4362  AssertThrowMPI(ierr);
4363  }
4364 
4365 
4366 # endif // DEAL_II_WITH_MPI
4367  }
4368 
4369  } // namespace internal
4370 
4371  template <typename DataType, typename MeshType>
4372  inline void
4374  const MeshType & mesh,
4375  const std::function<std_cxx17::optional<DataType>(
4376  const typename MeshType::active_cell_iterator &)> &pack,
4377  const std::function<void(const typename MeshType::active_cell_iterator &,
4378  const DataType &)> & unpack,
4379  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4380  &cell_filter)
4381  {
4382 # ifndef DEAL_II_WITH_MPI
4383  (void)mesh;
4384  (void)pack;
4385  (void)unpack;
4386  (void)cell_filter;
4387  Assert(false, ExcNeedsMPI());
4388 # else
4389  internal::exchange_cell_data<DataType,
4390  MeshType,
4391  typename MeshType::active_cell_iterator>(
4392  mesh,
4393  pack,
4394  unpack,
4395  cell_filter,
4396  [&](const auto &process) {
4397  for (const auto &cell : mesh.active_cell_iterators())
4398  if (cell->is_ghost())
4399  process(cell, cell->subdomain_id());
4400  },
4401  [](const auto &tria) { return tria.ghost_owners(); });
4402 # endif
4403  }
4404 
4405 
4406 
4407  template <typename DataType, typename MeshType>
4408  inline void
4410  const MeshType & mesh,
4411  const std::function<std_cxx17::optional<DataType>(
4412  const typename MeshType::level_cell_iterator &)> &pack,
4413  const std::function<void(const typename MeshType::level_cell_iterator &,
4414  const DataType &)> & unpack,
4415  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4416  &cell_filter)
4417  {
4418 # ifndef DEAL_II_WITH_MPI
4419  (void)mesh;
4420  (void)pack;
4421  (void)unpack;
4422  (void)cell_filter;
4423  Assert(false, ExcNeedsMPI());
4424 # else
4425  internal::exchange_cell_data<DataType,
4426  MeshType,
4427  typename MeshType::level_cell_iterator>(
4428  mesh,
4429  pack,
4430  unpack,
4431  cell_filter,
4432  [&](const auto &process) {
4433  for (const auto &cell : mesh.cell_iterators())
4434  if (cell->level_subdomain_id() !=
4436  !cell->is_locally_owned_on_level())
4437  process(cell, cell->level_subdomain_id());
4438  },
4439  [](const auto &tria) { return tria.level_ghost_owners(); });
4440 # endif
4441  }
4442 } // namespace GridTools
4443 
4444 # endif
4445 
4447 
4448 /*---------------------------- grid_tools.h ---------------------------*/
4449 /* end of #ifndef dealii_grid_tools_h */
4450 #endif
4451 /*---------------------------- 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:4031
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:3816
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:5340
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:1168
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3791
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:3999
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:1917
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:1742
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:2194
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:4942
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:4093
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:1025
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2244
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:2293
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:2642
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:5180
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3077
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:5244
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:2680
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:2891
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:3256
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:1072
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:2709
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:3211
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:3283
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:1999
#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:4064
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:5403
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:3883
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:3749
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:1679
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:4444
#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:1200
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:1551
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:3199
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2646
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:2656
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:3911
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:1014
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:4408
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3229
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1097
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:2972
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:1489
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 & 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:1005
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:3213
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:5526