Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55:02+00:00
\(\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 - 2022 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 
24 #include <deal.II/base/point.h>
26 
28 
30 
32 
33 #include <deal.II/fe/fe_values.h>
34 #include <deal.II/fe/mapping.h>
35 
36 #include <deal.II/grid/manifold.h>
37 #include <deal.II/grid/tria.h>
40 
42 #include <deal.II/lac/la_vector.h>
46 
47 #include <deal.II/numerics/rtree.h>
48 
49 #include <boost/archive/binary_iarchive.hpp>
50 #include <boost/archive/binary_oarchive.hpp>
51 #include <boost/random/mersenne_twister.hpp>
52 #include <boost/serialization/array.hpp>
53 #include <boost/serialization/vector.hpp>
54 
55 #ifdef DEAL_II_WITH_ZLIB
56 # include <boost/iostreams/device/back_inserter.hpp>
57 # include <boost/iostreams/filter/gzip.hpp>
58 # include <boost/iostreams/filtering_stream.hpp>
59 # include <boost/iostreams/stream.hpp>
60 #endif
61 
62 #include <bitset>
63 #include <list>
64 #include <set>
65 
66 #ifdef DEAL_II_HAVE_CXX20
67 # include <concepts>
68 #endif
69 
70 
72 
73 // Forward declarations
74 #ifndef DOXYGEN
75 namespace parallel
76 {
77  namespace distributed
78  {
79  template <int dim, int spacedim>
80  class Triangulation;
81  }
82 } // namespace parallel
83 
84 namespace hp
85 {
86  template <int, int>
87  class MappingCollection;
88 }
89 
90 class SparsityPattern;
91 
92 namespace GridTools
93 {
94  template <int dim, int spacedim>
95  class Cache;
96 }
97 #endif
98 
99 namespace internal
100 {
101  template <int dim, int spacedim, class MeshType>
102  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
104  {
105  public:
106 #ifndef _MSC_VER
107  using type = typename MeshType::active_cell_iterator;
108 #else
110 #endif
111  };
112 
113 #ifdef _MSC_VER
114  template <int dim, int spacedim>
115  class ActiveCellIterator<dim, spacedim, DoFHandler<dim, spacedim>>
116  {
117  public:
118  using type =
120  };
121 #endif
122 } // namespace internal
123 
132 namespace GridTools
133 {
145  template <int dim, int spacedim>
146  double
148 
172  template <int dim, int spacedim>
173  double
175 
203  template <int dim, int spacedim>
204  double
206  const Mapping<dim, spacedim> & mapping);
207 
218  template <int dim, int spacedim>
219  double
222  const Mapping<dim, spacedim> & mapping =
223  (ReferenceCells::get_hypercube<dim>()
224 #ifndef _MSC_VER
225  .template get_default_linear_mapping<dim, spacedim>()
226 #else
227  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
228 #endif
229  ));
230 
241  template <int dim, int spacedim>
242  double
245  const Mapping<dim, spacedim> & mapping =
246  (ReferenceCells::get_hypercube<dim>()
247 #ifndef _MSC_VER
248  .template get_default_linear_mapping<dim, spacedim>()
249 #else
250  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
251 #endif
252  ));
253 
265  template <int dim>
266  DEAL_II_DEPRECATED double
268  const std::vector<Point<dim>> &all_vertices,
270 
290  template <int dim>
291  double
292  cell_measure(const std::vector<Point<dim>> & all_vertices,
294 
317  template <int dim, int spacedim>
318  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
320 
347  template <int dim>
351  const Quadrature<dim> & quadrature);
352 
360  template <int dim>
361  double
364  const Quadrature<dim> & quadrature);
365 
379  template <int dim, int spacedim>
382 
400  template <typename Iterator>
403  const Iterator & object,
405 
418  template <int dim, int spacedim>
419  std::
420  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
422 
445  template <int dim, int spacedim>
446  void
448  std::vector<CellData<dim>> & cells,
449  SubCellData & subcelldata);
450 
469  template <int dim, int spacedim>
470  void
471  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
472  std::vector<CellData<dim>> & cells,
473  SubCellData & subcelldata,
474  std::vector<unsigned int> & considered_vertices,
475  const double tol = 1e-12);
476 
484  template <int dim>
485  void
487  const double tol = 1e-12);
488 
507  template <int dim, int spacedim>
508  void
510  const std::vector<Point<spacedim>> &all_vertices,
511  std::vector<CellData<dim>> & cells);
512 
522  template <int dim, int spacedim>
523  std::size_t
525  const std::vector<Point<spacedim>> &all_vertices,
526  std::vector<CellData<dim>> & cells);
527 
538  template <int dim>
539  void
540  consistently_order_cells(std::vector<CellData<dim>> &cells);
541 
610  template <int dim, typename Transformation, int spacedim>
612  (std::invocable<Transformation, Point<spacedim>> &&
613  std::assignable_from<
614  Point<spacedim> &,
615  std::invoke_result_t<Transformation, Point<spacedim>>>))
616  void transform(const Transformation & transformation,
617  Triangulation<dim, spacedim> &triangulation);
618 
625  template <int dim, int spacedim>
626  void
627  shift(const Tensor<1, spacedim> & shift_vector,
628  Triangulation<dim, spacedim> &triangulation);
629 
630 
641  template <int dim, int spacedim>
642  void
643  rotate(const double angle, Triangulation<dim, spacedim> &triangulation);
644 
657  template <int dim>
658  void
659  rotate(const Tensor<1, 3, double> &axis,
660  const double angle,
661  Triangulation<dim, 3> & triangulation);
662 
678  template <int dim>
679  DEAL_II_DEPRECATED void
680  rotate(const double angle,
681  const unsigned int axis,
682  Triangulation<dim, 3> &triangulation);
683 
741  template <int dim>
742  void
743  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
744  Triangulation<dim> & tria,
745  const Function<dim, double> *coefficient = nullptr,
746  const bool solve_for_absolute_positions = false);
747 
753  template <int dim, int spacedim>
754  std::map<unsigned int, Point<spacedim>>
755  get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria);
756 
764  template <int dim, int spacedim>
765  void
766  scale(const double scaling_factor,
767  Triangulation<dim, spacedim> &triangulation);
768 
783  template <int dim, int spacedim>
784  void
786  const double factor,
787  Triangulation<dim, spacedim> &triangulation,
788  const bool keep_boundary = true,
789  const unsigned int seed = boost::random::mt19937::default_seed);
790 
824  template <int dim, int spacedim>
825  void
826  remove_hanging_nodes(Triangulation<dim, spacedim> &tria,
827  const bool isotropic = false,
828  const unsigned int max_iterations = 100);
829 
854  template <int dim, int spacedim>
855  void
856  remove_anisotropy(Triangulation<dim, spacedim> &tria,
857  const double max_ratio = 1.6180339887,
858  const unsigned int max_iterations = 5);
859 
949  template <int dim, int spacedim>
950  void
952  const double limit_angle_fraction = .75);
953 
1014  template <int dim, int spacedim>
1015 #ifndef DOXYGEN
1016  std::tuple<
1017  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1018  std::vector<std::vector<Point<dim>>>,
1019  std::vector<std::vector<unsigned int>>>
1020 #else
1021  return_type
1022 #endif
1024  const Cache<dim, spacedim> & cache,
1025  const std::vector<Point<spacedim>> &points,
1027  &cell_hint =
1029 
1063  template <int dim, int spacedim>
1064 #ifndef DOXYGEN
1065  std::tuple<
1066  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1067  std::vector<std::vector<Point<dim>>>,
1068  std::vector<std::vector<unsigned int>>,
1069  std::vector<unsigned int>>
1070 #else
1071  return_type
1072 #endif
1074  const Cache<dim, spacedim> & cache,
1075  const std::vector<Point<spacedim>> &points,
1077  &cell_hint =
1079 
1149  template <int dim, int spacedim>
1150 #ifndef DOXYGEN
1151  std::tuple<
1152  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1153  std::vector<std::vector<Point<dim>>>,
1154  std::vector<std::vector<unsigned int>>,
1155  std::vector<std::vector<Point<spacedim>>>,
1156  std::vector<std::vector<unsigned int>>>
1157 #else
1158  return_type
1159 #endif
1161  const GridTools::Cache<dim, spacedim> & cache,
1162  const std::vector<Point<spacedim>> & local_points,
1163  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1164  const double tolerance = 1e-10);
1165 
1166  namespace internal
1167  {
1182  template <int dim, int spacedim>
1184  {
1191  std::vector<std::tuple<std::pair<int, int>,
1192  unsigned int,
1193  unsigned int,
1194  Point<dim>,
1196  unsigned int>>
1198 
1202  std::vector<unsigned int> send_ranks;
1203 
1209  std::vector<unsigned int> send_ptrs;
1210 
1221  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1223 
1227  std::vector<unsigned int> recv_ranks;
1228 
1234  std::vector<unsigned int> recv_ptrs;
1235  };
1236 
1246  template <int dim, int spacedim>
1249  const GridTools::Cache<dim, spacedim> & cache,
1250  const std::vector<Point<spacedim>> & points,
1251  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1252  const std::vector<bool> & marked_vertices,
1253  const double tolerance,
1254  const bool perform_handshake,
1255  const bool enforce_unique_mapping = false);
1256 
1257  } // namespace internal
1258 
1295  template <int dim, int spacedim>
1296  std::map<unsigned int, Point<spacedim>>
1298  const Triangulation<dim, spacedim> &container,
1299  const Mapping<dim, spacedim> & mapping =
1300  (ReferenceCells::get_hypercube<dim>()
1301 #ifndef _MSC_VER
1302  .template get_default_linear_mapping<dim, spacedim>()
1303 #else
1304  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1305 #endif
1306  ));
1307 
1317  template <int spacedim>
1318  unsigned int
1319  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1320  const Point<spacedim> & p);
1321 
1345  template <int dim, template <int, int> class MeshType, int spacedim>
1347  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1348  unsigned int find_closest_vertex(
1349  const MeshType<dim, spacedim> &mesh,
1350  const Point<spacedim> & p,
1351  const std::vector<bool> & marked_vertices = {});
1352 
1376  template <int dim, template <int, int> class MeshType, int spacedim>
1378  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1379  unsigned int find_closest_vertex(
1380  const Mapping<dim, spacedim> & mapping,
1381  const MeshType<dim, spacedim> &mesh,
1382  const Point<spacedim> & p,
1383  const std::vector<bool> & marked_vertices = {});
1384 
1385 
1405  template <int dim, template <int, int> class MeshType, int spacedim>
1407  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1408 #ifndef _MSC_VER
1409  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1410 #else
1411  std::vector<
1412  typename ::internal::
1413  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1414 #endif
1415  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1416  const unsigned int vertex_index);
1417 
1480  template <int dim, template <int, int> class MeshType, int spacedim>
1482  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1483 #ifndef _MSC_VER
1484  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1485 #else
1486  std::pair<typename ::internal::
1487  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1488  Point<dim>>
1489 #endif
1491  const MeshType<dim, spacedim> &mesh,
1492  const Point<spacedim> & p,
1493  const std::vector<bool> &marked_vertices = {},
1494  const double tolerance = 1.e-10);
1495 
1503  template <int dim, template <int, int> class MeshType, int spacedim>
1505  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1506 #ifndef _MSC_VER
1507  typename MeshType<dim, spacedim>::active_cell_iterator
1508 #else
1509  typename ::internal::
1510  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1511 #endif
1512  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1513  const Point<spacedim> & p,
1514  const std::vector<bool> &marked_vertices = {},
1515  const double tolerance = 1.e-10);
1516 
1523  template <int dim, int spacedim>
1524  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1525  Point<dim>>
1527  const hp::MappingCollection<dim, spacedim> &mapping,
1528  const DoFHandler<dim, spacedim> & mesh,
1529  const Point<spacedim> & p,
1530  const double tolerance = 1.e-10);
1531 
1583  template <int dim, int spacedim>
1584  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1585  Point<dim>>
1587  const Cache<dim, spacedim> &cache,
1588  const Point<spacedim> & p,
1591  const std::vector<bool> &marked_vertices = {},
1592  const double tolerance = 1.e-10);
1593 
1607  template <int dim, template <int, int> class MeshType, int spacedim>
1609  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1610 #ifndef _MSC_VER
1611  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1612 #else
1613  std::pair<typename ::internal::
1614  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1615  Point<dim>>
1616 #endif
1618  const Mapping<dim, spacedim> & mapping,
1619  const MeshType<dim, spacedim> &mesh,
1620  const Point<spacedim> & p,
1621  const std::vector<
1622  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1624  const std::vector<std::vector<Tensor<1, spacedim>>>
1625  &vertex_to_cell_centers,
1626  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1627  typename MeshType<dim, spacedim>::active_cell_iterator(),
1628  const std::vector<bool> &marked_vertices = {},
1629  const RTree<std::pair<Point<spacedim>, unsigned int>> &
1630  used_vertices_rtree = RTree<std::pair<Point<spacedim>, unsigned int>>{},
1631  const double tolerance = 1.e-10,
1632  const RTree<
1633  std::pair<BoundingBox<spacedim>,
1635  *relevant_cell_bounding_boxes_rtree = nullptr);
1636 
1658  template <int dim, template <int, int> class MeshType, int spacedim>
1660  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1661 #ifndef _MSC_VER
1662  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1663  Point<dim>>>
1664 #else
1665  std::vector<std::pair<
1666  typename ::internal::
1667  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1668  Point<dim>>>
1669 #endif
1671  const Mapping<dim, spacedim> & mapping,
1672  const MeshType<dim, spacedim> &mesh,
1673  const Point<spacedim> & p,
1674  const double tolerance,
1675  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1676  Point<dim>> & first_cell);
1677 
1684  template <int dim, template <int, int> class MeshType, int spacedim>
1686  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1687 #ifndef _MSC_VER
1688  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1689  Point<dim>>>
1690 #else
1691  std::vector<std::pair<
1692  typename ::internal::
1693  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1694  Point<dim>>>
1695 #endif
1697  const Mapping<dim, spacedim> & mapping,
1698  const MeshType<dim, spacedim> &mesh,
1699  const Point<spacedim> & p,
1700  const double tolerance = 1e-10,
1701  const std::vector<bool> & marked_vertices = {});
1702 
1724  template <class MeshType>
1725  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1726  std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
1727  const typename MeshType::cell_iterator &cell);
1728 
1753  template <class MeshType>
1756  const typename MeshType::active_cell_iterator & cell,
1757  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1758 
1808  template <class MeshType>
1810  std::
1811  vector<typename MeshType::active_cell_iterator> compute_active_cell_halo_layer(
1812  const MeshType &mesh,
1813  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1814  &predicate);
1815 
1816 
1824  template <class MeshType>
1826  std::
1827  vector<typename MeshType::cell_iterator> compute_cell_halo_layer_on_level(
1828  const MeshType &mesh,
1829  const std::function<bool(const typename MeshType::cell_iterator &)>
1830  & predicate,
1831  const unsigned int level);
1832 
1833 
1846  template <class MeshType>
1848  std::vector<
1849  typename MeshType::
1850  active_cell_iterator> compute_ghost_cell_halo_layer(const MeshType &mesh);
1851 
1900  template <class MeshType>
1902  std::
1903  vector<typename MeshType::active_cell_iterator> compute_active_cell_layer_within_distance(
1904  const MeshType &mesh,
1905  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1906  & predicate,
1907  const double layer_thickness);
1908 
1931  template <class MeshType>
1933  std::
1934  vector<typename MeshType::active_cell_iterator> compute_ghost_cell_layer_within_distance(
1935  const MeshType &mesh,
1936  const double layer_thickness);
1937 
1953  template <class MeshType>
1955  std::pair<
1956  Point<MeshType::space_dimension>,
1957  Point<MeshType::
1958  space_dimension>> compute_bounding_box(const MeshType &mesh,
1959  const std::function<bool(
1960  const typename MeshType::
1961  active_cell_iterator &)>
1962  &predicate);
1963 
2034  template <class MeshType>
2036  std::
2037  vector<BoundingBox<MeshType::space_dimension>> compute_mesh_predicate_bounding_box(
2038  const MeshType &mesh,
2039  const std::function<bool(const typename MeshType::active_cell_iterator &)>
2040  & predicate,
2041  const unsigned int refinement_level = 0,
2042  const bool allow_merge = false,
2043  const unsigned int max_boxes = numbers::invalid_unsigned_int);
2044 
2072  template <int spacedim>
2073 #ifndef DOXYGEN
2074  std::tuple<std::vector<std::vector<unsigned int>>,
2075  std::map<unsigned int, unsigned int>,
2076  std::map<unsigned int, std::vector<unsigned int>>>
2077 #else
2078  return_type
2079 #endif
2081  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
2082  const std::vector<Point<spacedim>> & points);
2083 
2084 
2119  template <int spacedim>
2120 #ifndef DOXYGEN
2121  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2122  std::map<unsigned int, unsigned int>,
2123  std::map<unsigned int, std::vector<unsigned int>>>
2124 #else
2125  return_type
2126 #endif
2128  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2129  const std::vector<Point<spacedim>> & points);
2130 
2131 
2140  template <int dim, int spacedim>
2141  std::vector<
2142  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2144 
2157  template <int dim, int spacedim>
2158  std::vector<std::vector<Tensor<1, spacedim>>>
2160  const Triangulation<dim, spacedim> &mesh,
2161  const std::vector<
2163  &vertex_to_cells);
2164 
2165 
2173  template <int dim, int spacedim>
2174  unsigned int
2177  const Point<spacedim> & position,
2178  const Mapping<dim, spacedim> & mapping =
2179  (ReferenceCells::get_hypercube<dim>()
2180 #ifndef _MSC_VER
2181  .template get_default_linear_mapping<dim, spacedim>()
2182 #else
2183  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2184 #endif
2185  ));
2186 
2198  template <int dim, int spacedim>
2199  std::map<unsigned int, types::global_vertex_index>
2202 
2214  template <int dim, int spacedim>
2215  std::pair<unsigned int, double>
2218 
2233  template <int dim, int spacedim>
2234  void
2237  DynamicSparsityPattern & connectivity);
2238 
2247  template <int dim, int spacedim>
2248  void
2251  DynamicSparsityPattern & connectivity);
2252 
2261  template <int dim, int spacedim>
2262  void
2265  const unsigned int level,
2266  DynamicSparsityPattern & connectivity);
2267 
2288  template <int dim, int spacedim>
2289  void
2290  partition_triangulation(const unsigned int n_partitions,
2292  const SparsityTools::Partitioner partitioner =
2294 
2305  template <int dim, int spacedim>
2306  void
2307  partition_triangulation(const unsigned int n_partitions,
2308  const std::vector<unsigned int> &cell_weights,
2310  const SparsityTools::Partitioner partitioner =
2312 
2358  template <int dim, int spacedim>
2359  void
2360  partition_triangulation(const unsigned int n_partitions,
2361  const SparsityPattern & cell_connection_graph,
2363  const SparsityTools::Partitioner partitioner =
2365 
2376  template <int dim, int spacedim>
2377  void
2378  partition_triangulation(const unsigned int n_partitions,
2379  const std::vector<unsigned int> &cell_weights,
2380  const SparsityPattern & cell_connection_graph,
2382  const SparsityTools::Partitioner partitioner =
2384 
2399  template <int dim, int spacedim>
2400  void
2401  partition_triangulation_zorder(const unsigned int n_partitions,
2403  const bool group_siblings = true);
2404 
2416  template <int dim, int spacedim>
2417  void
2419 
2427  template <int dim, int spacedim>
2428  std::vector<types::subdomain_id>
2430  const std::vector<CellId> & cell_ids);
2431 
2442  template <int dim, int spacedim>
2443  void
2445  std::vector<types::subdomain_id> & subdomain);
2446 
2461  template <int dim, int spacedim>
2462  unsigned int
2465  const types::subdomain_id subdomain);
2466 
2496  template <int dim, int spacedim>
2497  std::vector<bool>
2499 
2539  template <typename MeshType>
2540  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2541  std::list<std::pair<
2542  typename MeshType::cell_iterator,
2543  typename MeshType::cell_iterator>> get_finest_common_cells(const MeshType
2544  &mesh_1,
2545  const MeshType
2546  &mesh_2);
2547 
2557  template <int dim, int spacedim>
2558  bool
2560  const Triangulation<dim, spacedim> &mesh_2);
2561 
2571  template <typename MeshType>
2572  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2573  bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2574 
2596  template <int dim, int spacedim>
2600  & distorted_cells,
2602 
2603 
2604 
2652  template <class MeshType>
2653  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2654  std::vector<typename MeshType::active_cell_iterator> get_patch_around_cell(
2655  const typename MeshType::active_cell_iterator &cell);
2656 
2657 
2679  template <class Container>
2680  std::vector<typename Container::cell_iterator>
2682  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2683 
2750  template <class Container>
2751  void
2753  const std::vector<typename Container::active_cell_iterator> &patch,
2755  &local_triangulation,
2756  std::map<
2757  typename Triangulation<Container::dimension,
2758  Container::space_dimension>::active_cell_iterator,
2759  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2760 
2792  template <int dim, int spacedim>
2793  std::map<
2795  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2797 
2798 
2811  template <typename CellIterator>
2813  {
2817  CellIterator cell[2];
2818 
2823  unsigned int face_idx[2];
2824 
2830  std::bitset<3> orientation;
2831 
2845 
2849  std::size_t
2851  };
2852 
2853 
2917  template <typename FaceIterator>
2918  bool
2920  std::bitset<3> & orientation,
2921  const FaceIterator & face1,
2922  const FaceIterator & face2,
2923  const unsigned int direction,
2927 
2928 
2932  template <typename FaceIterator>
2933  bool
2935  const FaceIterator & face1,
2936  const FaceIterator & face2,
2937  const unsigned int direction,
2941 
2942 
2999  template <typename MeshType>
3000  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3002  const MeshType & mesh,
3003  const types::boundary_id b_id1,
3004  const types::boundary_id b_id2,
3005  const unsigned int direction,
3007  & matched_pairs,
3008  const Tensor<1, MeshType::space_dimension> &offset =
3011 
3012 
3035  template <typename MeshType>
3036  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3038  const MeshType & mesh,
3039  const types::boundary_id b_id,
3040  const unsigned int direction,
3041  std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
3042  & matched_pairs,
3043  const ::Tensor<1, MeshType::space_dimension> &offset =
3044  ::Tensor<1, MeshType::space_dimension>(),
3045  const FullMatrix<double> &matrix = FullMatrix<double>());
3046 
3073  template <int dim, int spacedim>
3074  void
3076  const bool reset_boundary_ids = false);
3077 
3099  template <int dim, int spacedim>
3100  void
3102  const std::vector<types::boundary_id> &src_boundary_ids,
3103  const std::vector<types::manifold_id> &dst_manifold_ids,
3104  Triangulation<dim, spacedim> & tria,
3105  const std::vector<types::boundary_id> &reset_boundary_ids = {});
3106 
3136  template <int dim, int spacedim>
3137  void
3139  const bool compute_face_ids = false);
3140 
3165  template <int dim, int spacedim>
3166  void
3169  const std::function<types::manifold_id(
3170  const std::set<types::manifold_id> &)> &disambiguation_function =
3171  [](const std::set<types::manifold_id> &manifold_ids) {
3172  if (manifold_ids.size() == 1)
3173  return *manifold_ids.begin();
3174  else
3176  },
3177  bool overwrite_only_flat_manifold_ids = true);
3264  template <typename DataType, typename MeshType>
3265  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3267  const MeshType & mesh,
3268  const std::function<std_cxx17::optional<DataType>(
3269  const typename MeshType::active_cell_iterator &)> &pack,
3270  const std::function<void(const typename MeshType::active_cell_iterator &,
3271  const DataType &)> & unpack,
3272  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3273  &cell_filter =
3274  always_return<typename MeshType::active_cell_iterator, bool>{true});
3275 
3286  template <typename DataType, typename MeshType>
3287  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3289  const MeshType & mesh,
3290  const std::function<std_cxx17::optional<DataType>(
3291  const typename MeshType::level_cell_iterator &)> &pack,
3292  const std::function<void(const typename MeshType::level_cell_iterator &,
3293  const DataType &)> & unpack,
3294  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3295  cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
3296  true});
3297 
3298  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3299  * boxes @p local_bboxes.
3300  *
3301  * This function is meant to exchange bounding boxes describing the locally
3302  * owned cells in a distributed triangulation obtained with the function
3303  * GridTools::compute_mesh_predicate_bounding_box .
3304  *
3305  * The output vector's size is the number of processes of the MPI
3306  * communicator:
3307  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3308  */
3309  template <int spacedim>
3310  std::vector<std::vector<BoundingBox<spacedim>>>
3312  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3313  const MPI_Comm & mpi_communicator);
3314 
3347  template <int spacedim>
3350  const std::vector<BoundingBox<spacedim>> &local_description,
3351  const MPI_Comm & mpi_communicator);
3352 
3370  template <int dim, int spacedim>
3371  void
3374  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3375  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3376 
3396  template <int dim, int spacedim>
3397  std::map<unsigned int, std::set<::types::subdomain_id>>
3400 
3416  template <int dim, typename T>
3418  {
3422  std::vector<CellId> cell_ids;
3423 
3427  std::vector<T> data;
3428 
3437  template <class Archive>
3438  void
3439  save(Archive &ar, const unsigned int) const
3440  {
3441  // Implement the code inline as some compilers do otherwise complain
3442  // about the use of a deprecated interface.
3443  Assert(cell_ids.size() == data.size(),
3444  ExcDimensionMismatch(cell_ids.size(), data.size()));
3445  // archive the cellids in an efficient binary format
3446  const std::size_t n_cells = cell_ids.size();
3447  ar & n_cells;
3448  for (const auto &id : cell_ids)
3449  {
3450  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3451  ar & binary_cell_id;
3452  }
3453 
3454  ar &data;
3455  }
3456 
3463  template <class Archive>
3464  void
3465  load(Archive &ar, const unsigned int)
3466  {
3467  // Implement the code inline as some compilers do otherwise complain
3468  // about the use of a deprecated interface.
3469  std::size_t n_cells;
3470  ar & n_cells;
3471  cell_ids.clear();
3472  cell_ids.reserve(n_cells);
3473  for (unsigned int c = 0; c < n_cells; ++c)
3474  {
3475  CellId::binary_type value;
3476  ar & value;
3477  cell_ids.emplace_back(value);
3478  }
3479  ar &data;
3480  }
3481 
3482 
3483 #ifdef DOXYGEN
3489  template <class Archive>
3490  void
3491  serialize(Archive &archive, const unsigned int version);
3492 #else
3493  // This macro defines the serialize() method that is compatible with
3494  // the templated save() and load() method that have been implemented.
3495  BOOST_SERIALIZATION_SPLIT_MEMBER()
3496 #endif
3497  };
3498 
3499 
3500 
3521  template <int dim, typename VectorType>
3523  {
3524  public:
3528  using value_type = typename VectorType::value_type;
3529 
3533  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3534  const FiniteElement<dim, dim> &fe,
3535  const unsigned int n_subdivisions = 1,
3536  const double tolerance = 1e-10);
3537 
3548  void
3549  process(const DoFHandler<dim> & background_dof_handler,
3550  const VectorType & ls_vector,
3551  const double iso_level,
3552  std::vector<Point<dim>> &vertices,
3553  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3554 
3559  void
3560  process(const DoFHandler<dim> & background_dof_handler,
3561  const VectorType & ls_vector,
3562  const double iso_level,
3563  std::vector<Point<dim>> &vertices) const;
3564 
3574  void
3576  const VectorType & ls_vector,
3577  const double iso_level,
3578  std::vector<Point<dim>> & vertices,
3579  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3585  void
3587  const VectorType & ls_vector,
3588  const double iso_level,
3589  std::vector<Point<dim>> &vertices) const;
3590 
3591  private:
3596  static Quadrature<dim>
3597  create_quadrature_rule(const unsigned int n_subdivisions);
3598 
3602  void
3603  process_cell(std::vector<value_type> & ls_values,
3604  const std::vector<Point<dim>> & points,
3605  const double iso_level,
3606  std::vector<Point<dim>> & vertices,
3607  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
3608  const bool write_back_cell_data = true) const;
3609 
3613  void
3614  process_sub_cell(const std::vector<value_type> &,
3615  const std::vector<Point<1>> &,
3616  const std::vector<unsigned int> &,
3617  const double,
3618  std::vector<Point<1>> &,
3619  std::vector<CellData<1>> &,
3620  const bool) const
3621  {
3622  AssertThrow(false, ExcNotImplemented());
3623  }
3624 
3631  void
3632  process_sub_cell(const std::vector<value_type> & ls_values,
3633  const std::vector<Point<2>> & points,
3634  const std::vector<unsigned int> &mask,
3635  const double iso_level,
3636  std::vector<Point<2>> & vertices,
3637  std::vector<CellData<1>> & cells,
3638  const bool write_back_cell_data) const;
3639 
3643  void
3644  process_sub_cell(const std::vector<value_type> & ls_values,
3645  const std::vector<Point<3>> & points,
3646  const std::vector<unsigned int> &mask,
3647  const double iso_level,
3648  std::vector<Point<3>> & vertices,
3649  std::vector<CellData<2>> & cells,
3650  const bool write_back_cell_data) const;
3651 
3656  const unsigned int n_subdivisions;
3657 
3662  const double tolerance;
3663 
3669  };
3670 
3671 
3672 
3682  int,
3683  << "The number of partitions you gave is " << arg1
3684  << ", but must be greater than zero.");
3689  int,
3690  << "The subdomain id " << arg1
3691  << " has no cells associated with it.");
3696 
3701  double,
3702  << "The scaling factor must be positive, but it is " << arg1
3703  << '.');
3704 
3709  unsigned int,
3710  << "The given vertex with index " << arg1
3711  << " is not used in the given triangulation.");
3712 
3715 } /*namespace GridTools*/
3716 
3717 
3728  "The edges of the mesh are not consistently orientable.");
3729 
3730 
3731 
3732 /* ----------------- Template function --------------- */
3733 
3734 #ifndef DOXYGEN
3735 
3736 namespace GridTools
3737 {
3738  template <int dim>
3739  double
3740  cell_measure(
3741  const std::vector<Point<dim>> &all_vertices,
3742  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3743  {
3744  // We forward call to the ArrayView version:
3745  const ArrayView<const unsigned int> view(
3747  return cell_measure(all_vertices, view);
3748  }
3749 
3750 
3751 
3752  // This specialization is defined here so that the general template in the
3753  // source file doesn't need to have further 1d overloads for the internal
3754  // functions it calls.
3755  template <>
3759  {
3760  return {};
3761  }
3762 
3763 
3764 
3765  template <int dim, typename Transformation, int spacedim>
3767  (std::invocable<Transformation, Point<spacedim>> &&
3768  std::assignable_from<
3769  Point<spacedim> &,
3770  std::invoke_result_t<Transformation, Point<spacedim>>>))
3771  void transform(const Transformation & transformation,
3772  Triangulation<dim, spacedim> &triangulation)
3773  {
3774  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3775 
3776  // loop over all active cells, and
3777  // transform those vertices that
3778  // have not yet been touched. note
3779  // that we get to all vertices in
3780  // the triangulation by only
3781  // visiting the active cells.
3783  cell = triangulation.begin_active(),
3784  endc = triangulation.end();
3785  for (; cell != endc; ++cell)
3786  for (const unsigned int v : cell->vertex_indices())
3787  if (treated_vertices[cell->vertex_index(v)] == false)
3788  {
3789  // transform this vertex
3790  cell->vertex(v) = transformation(cell->vertex(v));
3791  // and mark it as treated
3792  treated_vertices[cell->vertex_index(v)] = true;
3793  };
3794 
3795 
3796  // now fix any vertices on hanging nodes so that we don't create any holes
3797  if (dim == 2)
3798  {
3800  cell = triangulation.begin_active(),
3801  endc = triangulation.end();
3802  for (; cell != endc; ++cell)
3803  for (const unsigned int face : cell->face_indices())
3804  if (cell->face(face)->has_children() &&
3805  !cell->face(face)->at_boundary())
3806  {
3807  Assert(cell->reference_cell() ==
3808  ReferenceCells::get_hypercube<dim>(),
3809  ExcNotImplemented());
3810 
3811  // this line has children
3812  cell->face(face)->child(0)->vertex(1) =
3813  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3814  2;
3815  }
3816  }
3817  else if (dim == 3)
3818  {
3820  cell = triangulation.begin_active(),
3821  endc = triangulation.end();
3822  for (; cell != endc; ++cell)
3823  for (const unsigned int face : cell->face_indices())
3824  if (cell->face(face)->has_children() &&
3825  !cell->face(face)->at_boundary())
3826  {
3827  Assert(cell->reference_cell() ==
3828  ReferenceCells::get_hypercube<dim>(),
3829  ExcNotImplemented());
3830 
3831  // this face has hanging nodes
3832  cell->face(face)->child(0)->vertex(1) =
3833  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3834  2.0;
3835  cell->face(face)->child(0)->vertex(2) =
3836  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3837  2.0;
3838  cell->face(face)->child(1)->vertex(3) =
3839  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3840  2.0;
3841  cell->face(face)->child(2)->vertex(3) =
3842  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3843  2.0;
3844 
3845  // center of the face
3846  cell->face(face)->child(0)->vertex(3) =
3847  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3848  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3849  4.0;
3850  }
3851  }
3852 
3853  // Make sure FEValues notices that the mesh has changed
3854  triangulation.signals.mesh_movement();
3855  }
3856 
3857 
3858 
3859  template <class MeshType>
3860  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3861  std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
3862  const typename MeshType::cell_iterator &cell)
3863  {
3864  std::vector<typename MeshType::active_cell_iterator> child_cells;
3865 
3866  if (cell->has_children())
3867  {
3868  for (unsigned int child = 0; child < cell->n_children(); ++child)
3869  if (cell->child(child)->has_children())
3870  {
3871  const std::vector<typename MeshType::active_cell_iterator>
3872  children = get_active_child_cells<MeshType>(cell->child(child));
3873  child_cells.insert(child_cells.end(),
3874  children.begin(),
3875  children.end());
3876  }
3877  else
3878  child_cells.push_back(cell->child(child));
3879  }
3880 
3881  return child_cells;
3882  }
3883 
3884 
3885 
3886  template <class MeshType>
3887  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3888  void get_active_neighbors(
3889  const typename MeshType::active_cell_iterator & cell,
3890  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3891  {
3892  active_neighbors.clear();
3893  for (const unsigned int n : cell->face_indices())
3894  if (!cell->at_boundary(n))
3895  {
3896  if (MeshType::dimension == 1)
3897  {
3898  // check children of neighbor. note
3899  // that in 1d children of the neighbor
3900  // may be further refined. In 1d the
3901  // case is simple since we know what
3902  // children bound to the present cell
3903  typename MeshType::cell_iterator neighbor_child =
3904  cell->neighbor(n);
3905  if (!neighbor_child->is_active())
3906  {
3907  while (neighbor_child->has_children())
3908  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3909 
3910  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3911  ExcInternalError());
3912  }
3913  active_neighbors.push_back(neighbor_child);
3914  }
3915  else
3916  {
3917  if (cell->face(n)->has_children())
3918  // this neighbor has children. find
3919  // out which border to the present
3920  // cell
3921  for (unsigned int c = 0;
3922  c < cell->face(n)->n_active_descendants();
3923  ++c)
3924  active_neighbors.push_back(
3925  cell->neighbor_child_on_subface(n, c));
3926  else
3927  {
3928  // the neighbor must be active
3929  // himself
3930  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3931  active_neighbors.push_back(cell->neighbor(n));
3932  }
3933  }
3934  }
3935  }
3936 
3937 
3938 
3939  template <typename CellIterator>
3940  std::size_t
3942  {
3943  return sizeof(*this) + matrix.memory_consumption();
3944  }
3945 
3946 
3947 
3948  namespace internal
3949  {
3950  namespace ProjectToObject
3951  {
3964  struct CrossDerivative
3965  {
3966  const unsigned int direction_0;
3967  const unsigned int direction_1;
3968 
3969  CrossDerivative(const unsigned int d0, const unsigned int d1);
3970  };
3971 
3972  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3973  const unsigned int d1)
3974  : direction_0(d0)
3975  , direction_1(d1)
3976  {}
3977 
3978 
3979 
3984  template <typename F>
3985  inline auto
3986  centered_first_difference(const double center,
3987  const double step,
3988  const F &f) -> decltype(f(center) - f(center))
3989  {
3990  return (f(center + step) - f(center - step)) / (2.0 * step);
3991  }
3992 
3993 
3994 
3999  template <typename F>
4000  inline auto
4001  centered_second_difference(const double center,
4002  const double step,
4003  const F &f) -> decltype(f(center) - f(center))
4004  {
4005  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
4006  (step * step);
4007  }
4008 
4009 
4010 
4020  template <int structdim, typename F>
4021  inline auto
4022  cross_stencil(
4023  const CrossDerivative cross_derivative,
4025  const double step,
4026  const F &f) -> decltype(f(center) - f(center))
4027  {
4029  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
4030  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
4031  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
4032  1.0 / 3.0 * f(center - simplex_vector) +
4033  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
4034  step;
4035  }
4036 
4037 
4038 
4045  template <int spacedim, int structdim, typename F>
4046  inline double
4047  gradient_entry(
4048  const unsigned int row_n,
4049  const unsigned int dependent_direction,
4050  const Point<spacedim> &p0,
4052  const double step,
4053  const F & f)
4054  {
4056  dependent_direction <
4058  ExcMessage("This function assumes that the last weight is a "
4059  "dependent variable (and hence we cannot take its "
4060  "derivative directly)."));
4061  Assert(row_n != dependent_direction,
4062  ExcMessage(
4063  "We cannot differentiate with respect to the variable "
4064  "that is assumed to be dependent."));
4065 
4066  const Point<spacedim> manifold_point = f(center);
4067  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
4068  {row_n, dependent_direction}, center, step, f);
4069  double entry = 0.0;
4070  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
4071  entry +=
4072  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
4073  return entry;
4074  }
4075 
4081  template <typename Iterator, int spacedim, int structdim>
4083  project_to_d_linear_object(const Iterator & object,
4084  const Point<spacedim> &trial_point)
4085  {
4086  // let's look at this for simplicity for a quad (structdim==2) in a
4087  // space with spacedim>2 (notate trial_point by y): all points on the
4088  // surface are given by
4089  // x(\xi) = sum_i v_i phi_x(\xi)
4090  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
4091  // reference coordinates of the quad. so what we are trying to do is
4092  // find a point x on the surface that is closest to the point y. there
4093  // are different ways to solve this problem, but in the end it's a
4094  // nonlinear problem and we have to find reference coordinates \xi so
4095  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
4096  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
4097  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
4098  // have to use a Newton method to find the answer. This leads to the
4099  // following formulation of Newton steps:
4100  //
4101  // Given \xi_k, find \delta\xi_k so that
4102  // H_k \delta\xi_k = - F_k
4103  // where H_k is an approximation to the second derivatives of J at
4104  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
4105  // number of times until the right hand side is small enough. As a
4106  // stopping criterion, we terminate if ||\delta\xi||<eps.
4107  //
4108  // As for the Hessian, the best choice would be
4109  // H_k = J''(\xi_k)
4110  // but we'll opt for the simpler Gauss-Newton form
4111  // H_k = A^T A
4112  // i.e.
4113  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
4114  // \partial_n phi_i *
4115  // \partial_m phi_j
4116  // we start at xi=(0.5, 0.5).
4117  Point<structdim> xi;
4118  for (unsigned int d = 0; d < structdim; ++d)
4119  xi[d] = 0.5;
4120 
4121  Point<spacedim> x_k;
4122  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4123  x_k += object->vertex(i) *
4125 
4126  do
4127  {
4129  for (const unsigned int i :
4131  F_k +=
4132  (x_k - trial_point) * object->vertex(i) *
4134  i);
4135 
4137  for (const unsigned int i :
4139  for (const unsigned int j :
4141  {
4144  xi, i),
4146  xi, j));
4147  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
4148  }
4149 
4150  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
4151  xi += delta_xi;
4152 
4153  x_k = Point<spacedim>();
4154  for (const unsigned int i :
4156  x_k += object->vertex(i) *
4158 
4159  if (delta_xi.norm() < 1e-7)
4160  break;
4161  }
4162  while (true);
4163 
4164  return x_k;
4165  }
4166  } // namespace ProjectToObject
4167  } // namespace internal
4168 
4169 
4170 
4171  namespace internal
4172  {
4173  // We hit an internal compiler error in ICC 15 if we define this as a lambda
4174  // inside the project_to_object function below.
4175  template <int structdim>
4176  inline bool
4177  weights_are_ok(
4179  {
4180  // clang has trouble figuring out structdim here, so define it
4181  // again:
4182  static const std::size_t n_vertices_per_cell =
4184  n_independent_components;
4185  std::array<double, n_vertices_per_cell> copied_weights;
4186  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4187  {
4188  copied_weights[i] = v[i];
4189  if (v[i] < 0.0 || v[i] > 1.0)
4190  return false;
4191  }
4192 
4193  // check the sum: try to avoid some roundoff errors by summing in order
4194  std::sort(copied_weights.begin(), copied_weights.end());
4195  const double sum =
4196  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4197  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4198  }
4199  } // namespace internal
4200 
4201  template <typename Iterator>
4204  const Iterator & object,
4206  {
4207  const int spacedim = Iterator::AccessorType::space_dimension;
4208  const int structdim = Iterator::AccessorType::structure_dimension;
4209 
4210  Point<spacedim> projected_point = trial_point;
4211 
4212  if (structdim >= spacedim)
4213  return projected_point;
4214  else if (structdim == 1 || structdim == 2)
4215  {
4216  using namespace internal::ProjectToObject;
4217  // Try to use the special flat algorithm for quads (this is better
4218  // than the general algorithm in 3d). This does not take into account
4219  // whether projected_point is outside the quad, but we optimize along
4220  // lines below anyway:
4221  const int dim = Iterator::AccessorType::dimension;
4222  const Manifold<dim, spacedim> &manifold = object->get_manifold();
4223  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4224  &manifold) != nullptr)
4225  {
4226  projected_point =
4227  project_to_d_linear_object<Iterator, spacedim, structdim>(
4228  object, trial_point);
4229  }
4230  else
4231  {
4232  // We want to find a point on the convex hull (defined by the
4233  // vertices of the object and the manifold description) that is
4234  // relatively close to the trial point. This has a few issues:
4235  //
4236  // 1. For a general convex hull we are not guaranteed that a unique
4237  // minimum exists.
4238  // 2. The independent variables in the optimization process are the
4239  // weights given to Manifold::get_new_point, which must sum to 1,
4240  // so we cannot use standard finite differences to approximate a
4241  // gradient.
4242  //
4243  // There is not much we can do about 1., but for 2. we can derive
4244  // finite difference stencils that work on a structdim-dimensional
4245  // simplex and rewrite the optimization problem to use those
4246  // instead. Consider the structdim 2 case and let
4247  //
4248  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4249  // c2, c3})
4250  //
4251  // where {c0, c1, c2, c3} are the weights for the four vertices on
4252  // the quadrilateral. We seek to minimize the Euclidean distance
4253  // between F(...) and trial_point. We can solve for c3 in terms of
4254  // the other weights and get, for one coordinate direction
4255  //
4256  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4257  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4258  //
4259  // where we substitute back in for c3 after taking the
4260  // derivative. We can compute a stencil for the cross derivative
4261  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4262  // (and gradient_entry computes the sum over the independent
4263  // variables). Below, we somewhat arbitrarily pick the last
4264  // component as the dependent one.
4265  //
4266  // Since we can now calculate derivatives of the objective
4267  // function we can use gradient descent to minimize it.
4268  //
4269  // Of course, this is much simpler in the structdim = 1 case (we
4270  // could rewrite the projection as a 1d optimization problem), but
4271  // to reduce the potential for bugs we use the same code in both
4272  // cases.
4273  const double step_size = object->diameter() / 64.0;
4274 
4275  constexpr unsigned int n_vertices_per_cell =
4277 
4278  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4279  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4280  ++vertex_n)
4281  vertices[vertex_n] = object->vertex(vertex_n);
4282 
4283  auto get_point_from_weights =
4284  [&](const Tensor<1, n_vertices_per_cell> &weights)
4285  -> Point<spacedim> {
4286  return object->get_manifold().get_new_point(
4287  make_array_view(vertices.begin(), vertices.end()),
4288  make_array_view(weights.begin_raw(), weights.end_raw()));
4289  };
4290 
4291  // pick the initial weights as (normalized) inverse distances from
4292  // the trial point:
4293  Tensor<1, n_vertices_per_cell> guess_weights;
4294  double guess_weights_sum = 0.0;
4295  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4296  ++vertex_n)
4297  {
4298  const double distance =
4299  vertices[vertex_n].distance(trial_point);
4300  if (distance == 0.0)
4301  {
4302  guess_weights = 0.0;
4303  guess_weights[vertex_n] = 1.0;
4304  guess_weights_sum = 1.0;
4305  break;
4306  }
4307  else
4308  {
4309  guess_weights[vertex_n] = 1.0 / distance;
4310  guess_weights_sum += guess_weights[vertex_n];
4311  }
4312  }
4313  guess_weights /= guess_weights_sum;
4314  Assert(internal::weights_are_ok<structdim>(guess_weights),
4315  ExcInternalError());
4316 
4317  // The optimization algorithm consists of two parts:
4318  //
4319  // 1. An outer loop where we apply the gradient descent algorithm.
4320  // 2. An inner loop where we do a line search to find the optimal
4321  // length of the step one should take in the gradient direction.
4322  //
4323  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4324  {
4325  const unsigned int dependent_direction =
4326  n_vertices_per_cell - 1;
4327  Tensor<1, n_vertices_per_cell> current_gradient;
4328  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4329  ++row_n)
4330  {
4331  if (row_n != dependent_direction)
4332  {
4333  current_gradient[row_n] =
4334  gradient_entry<spacedim, structdim>(
4335  row_n,
4336  dependent_direction,
4337  trial_point,
4338  guess_weights,
4339  step_size,
4340  get_point_from_weights);
4341 
4342  current_gradient[dependent_direction] -=
4343  current_gradient[row_n];
4344  }
4345  }
4346 
4347  // We need to travel in the -gradient direction, as noted
4348  // above, but we may not want to take a full step in that
4349  // direction; instead, guess that we will go -0.5*gradient and
4350  // do quasi-Newton iteration to pick the best multiplier. The
4351  // goal is to find a scalar alpha such that
4352  //
4353  // F(x - alpha g)
4354  //
4355  // is minimized, where g is the gradient and F is the
4356  // objective function. To find the optimal value we find roots
4357  // of the derivative of the objective function with respect to
4358  // alpha by Newton iteration, where we approximate the first
4359  // and second derivatives of F(x - alpha g) with centered
4360  // finite differences.
4361  double gradient_weight = -0.5;
4362  auto gradient_weight_objective_function =
4363  [&](const double gradient_weight_guess) -> double {
4364  return (trial_point -
4365  get_point_from_weights(guess_weights +
4366  gradient_weight_guess *
4367  current_gradient))
4368  .norm_square();
4369  };
4370 
4371  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4372  {
4373  const double update_numerator = centered_first_difference(
4374  gradient_weight,
4375  step_size,
4376  gradient_weight_objective_function);
4377  const double update_denominator =
4378  centered_second_difference(
4379  gradient_weight,
4380  step_size,
4381  gradient_weight_objective_function);
4382 
4383  // avoid division by zero. Note that we limit the gradient
4384  // weight below
4385  if (std::abs(update_denominator) == 0.0)
4386  break;
4387  gradient_weight =
4388  gradient_weight - update_numerator / update_denominator;
4389 
4390  // Put a fairly lenient bound on the largest possible
4391  // gradient (things tend to be locally flat, so the gradient
4392  // itself is usually small)
4393  if (std::abs(gradient_weight) > 10)
4394  {
4395  gradient_weight = -10.0;
4396  break;
4397  }
4398  }
4399 
4400  // It only makes sense to take convex combinations with weights
4401  // between zero and one. If the update takes us outside of this
4402  // region then rescale the update to stay within the region and
4403  // try again
4404  Tensor<1, n_vertices_per_cell> tentative_weights =
4405  guess_weights + gradient_weight * current_gradient;
4406 
4407  double new_gradient_weight = gradient_weight;
4408  for (unsigned int iteration_count = 0; iteration_count < 40;
4409  ++iteration_count)
4410  {
4411  if (internal::weights_are_ok<structdim>(tentative_weights))
4412  break;
4413 
4414  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4415  {
4416  if (tentative_weights[i] < 0.0)
4417  {
4418  tentative_weights -=
4419  (tentative_weights[i] / current_gradient[i]) *
4420  current_gradient;
4421  }
4422  if (tentative_weights[i] < 0.0 ||
4423  1.0 < tentative_weights[i])
4424  {
4425  new_gradient_weight /= 2.0;
4426  tentative_weights =
4427  guess_weights +
4428  new_gradient_weight * current_gradient;
4429  }
4430  }
4431  }
4432 
4433  // the update might still send us outside the valid region, so
4434  // check again and quit if the update is still not valid
4435  if (!internal::weights_are_ok<structdim>(tentative_weights))
4436  break;
4437 
4438  // if we cannot get closer by traveling in the gradient
4439  // direction then quit
4440  if (get_point_from_weights(tentative_weights)
4441  .distance(trial_point) <
4442  get_point_from_weights(guess_weights).distance(trial_point))
4443  guess_weights = tentative_weights;
4444  else
4445  break;
4446  Assert(internal::weights_are_ok<structdim>(guess_weights),
4447  ExcInternalError());
4448  }
4449  Assert(internal::weights_are_ok<structdim>(guess_weights),
4450  ExcInternalError());
4451  projected_point = get_point_from_weights(guess_weights);
4452  }
4453 
4454  // if structdim == 2 and the optimal point is not on the interior then
4455  // we may be able to get a more accurate result by projecting onto the
4456  // lines.
4457  if (structdim == 2)
4458  {
4459  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4460  line_projections;
4461  for (unsigned int line_n = 0;
4462  line_n < GeometryInfo<structdim>::lines_per_cell;
4463  ++line_n)
4464  {
4465  line_projections[line_n] =
4466  project_to_object(object->line(line_n), trial_point);
4467  }
4468  std::sort(line_projections.begin(),
4469  line_projections.end(),
4470  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4471  return a.distance(trial_point) <
4472  b.distance(trial_point);
4473  });
4474  if (line_projections[0].distance(trial_point) <
4475  projected_point.distance(trial_point))
4476  projected_point = line_projections[0];
4477  }
4478  }
4479  else
4480  {
4481  Assert(false, ExcNotImplemented());
4482  return projected_point;
4483  }
4484 
4485  return projected_point;
4486  }
4487 
4488 
4489 
4490  namespace internal
4491  {
4492  template <typename DataType,
4493  typename MeshType,
4494  typename MeshCellIteratorType>
4495  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
4496  inline void exchange_cell_data(
4497  const MeshType &mesh,
4498  const std::function<
4499  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4500  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4501  & unpack,
4502  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4503  const std::function<void(
4504  const std::function<void(const MeshCellIteratorType &,
4505  const types::subdomain_id)> &)> &process_cells,
4506  const std::function<std::set<types::subdomain_id>(
4507  const parallel::TriangulationBase<MeshType::dimension,
4508  MeshType::space_dimension> &)>
4509  &compute_ghost_owners)
4510  {
4511 # ifndef DEAL_II_WITH_MPI
4512  (void)mesh;
4513  (void)pack;
4514  (void)unpack;
4515  (void)cell_filter;
4516  (void)process_cells;
4517  (void)compute_ghost_owners;
4518  Assert(false, ExcNeedsMPI());
4519 # else
4520  constexpr int dim = MeshType::dimension;
4521  constexpr int spacedim = MeshType::space_dimension;
4522  auto tria =
4523  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4524  &mesh.get_triangulation());
4525  Assert(
4526  tria != nullptr,
4527  ExcMessage(
4528  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4529 
4530  if (const auto tria = dynamic_cast<
4532  &mesh.get_triangulation()))
4533  {
4534  Assert(
4535  tria->with_artificial_cells(),
4536  ExcMessage(
4537  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4538  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4539  "operate on a single layer of ghost cells. However, you have "
4540  "given a Triangulation object of type "
4541  "parallel::shared::Triangulation without artificial cells "
4542  "resulting in arbitrary numbers of ghost layers."));
4543  }
4544 
4545  // build list of cells to request for each neighbor
4546  std::set<::types::subdomain_id> ghost_owners =
4547  compute_ghost_owners(*tria);
4548  std::map<::types::subdomain_id,
4549  std::vector<typename CellId::binary_type>>
4550  neighbor_cell_list;
4551 
4552  for (const auto ghost_owner : ghost_owners)
4553  neighbor_cell_list[ghost_owner] = {};
4554 
4555  process_cells([&](const auto &cell, const auto key) -> void {
4556  if (cell_filter(cell))
4557  {
4558  constexpr int spacedim = MeshType::space_dimension;
4559  neighbor_cell_list[key].emplace_back(
4560  cell->id().template to_binary<spacedim>());
4561  }
4562  });
4563 
4564  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4565  ExcInternalError());
4566 
4567 
4568  // Before sending & receiving, make sure we protect this section with
4569  // a mutex:
4570  static Utilities::MPI::CollectiveMutex mutex;
4572  mutex, tria->get_communicator());
4573 
4574  const int mpi_tag =
4576  const int mpi_tag_reply =
4578 
4579  // send our requests
4580  std::vector<MPI_Request> requests(ghost_owners.size());
4581  {
4582  unsigned int idx = 0;
4583  for (const auto &it : neighbor_cell_list)
4584  {
4585  // send the data about the relevant cells
4586  const int ierr = MPI_Isend(it.second.data(),
4587  it.second.size() * sizeof(it.second[0]),
4588  MPI_BYTE,
4589  it.first,
4590  mpi_tag,
4592  &requests[idx]);
4593  AssertThrowMPI(ierr);
4594  ++idx;
4595  }
4596  }
4597 
4598  // receive requests and reply with the results
4599  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4600  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4601 
4602  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4603  {
4604  MPI_Status status;
4605  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4606  mpi_tag,
4608  &status);
4609  AssertThrowMPI(ierr);
4610 
4611  int len;
4612  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4613  AssertThrowMPI(ierr);
4614  Assert(len % sizeof(typename CellId::binary_type) == 0,
4615  ExcInternalError());
4616 
4617  const unsigned int n_cells =
4618  len / sizeof(typename CellId::binary_type);
4619  std::vector<typename CellId::binary_type> cells_with_requests(
4620  n_cells);
4621  std::vector<DataType> data_to_send;
4622  data_to_send.reserve(n_cells);
4623  std::vector<bool> cell_carries_data(n_cells, false);
4624 
4625  ierr = MPI_Recv(cells_with_requests.data(),
4626  len,
4627  MPI_BYTE,
4628  status.MPI_SOURCE,
4629  status.MPI_TAG,
4631  &status);
4632  AssertThrowMPI(ierr);
4633 
4634  // store data for each cell
4635  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4636  {
4637  const auto cell =
4638  tria->create_cell_iterator(CellId(cells_with_requests[c]));
4639 
4640  MeshCellIteratorType mesh_it(tria,
4641  cell->level(),
4642  cell->index(),
4643  &mesh);
4644 
4645  const std_cxx17::optional<DataType> data = pack(mesh_it);
4646  if (data)
4647  {
4648  data_to_send.emplace_back(*data);
4649  cell_carries_data[c] = true;
4650  }
4651  }
4652 
4653  // collect data for sending the reply in a buffer
4654 
4655  // (a) make room for storing the local offsets in case we receive
4656  // other data
4657  sendbuffers[idx].resize(sizeof(std::size_t));
4658 
4659  // (b) append the actual data and store how much memory it
4660  // corresponds to, which we then insert into the leading position of
4661  // the sendbuffer
4662  std::size_t size_of_send =
4663  Utilities::pack(data_to_send,
4664  sendbuffers[idx],
4665  /*enable_compression*/ false);
4666  std::memcpy(sendbuffers[idx].data(),
4667  &size_of_send,
4668  sizeof(std::size_t));
4669 
4670  // (c) append information of certain cells that got left out in case
4671  // we need it
4672  if (data_to_send.size() < n_cells)
4673  Utilities::pack(cell_carries_data,
4674  sendbuffers[idx],
4675  /*enable_compression*/ false);
4676 
4677  // send data
4678  ierr = MPI_Isend(sendbuffers[idx].data(),
4679  sendbuffers[idx].size(),
4680  MPI_BYTE,
4681  status.MPI_SOURCE,
4682  mpi_tag_reply,
4684  &reply_requests[idx]);
4685  AssertThrowMPI(ierr);
4686  }
4687 
4688  // finally receive the replies
4689  std::vector<char> receive;
4690  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4691  {
4692  MPI_Status status;
4693  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4694  mpi_tag_reply,
4696  &status);
4697  AssertThrowMPI(ierr);
4698 
4699  int len;
4700  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4701  AssertThrowMPI(ierr);
4702 
4703  receive.resize(len);
4704 
4705  ierr = MPI_Recv(receive.data(),
4706  len,
4707  MPI_BYTE,
4708  status.MPI_SOURCE,
4709  status.MPI_TAG,
4711  &status);
4712  AssertThrowMPI(ierr);
4713 
4714  // (a) first determine the length of the data section in the
4715  // received buffer
4716  auto data_iterator = receive.begin();
4717  std::size_t size_of_received_data =
4718  Utilities::unpack<std::size_t>(data_iterator,
4719  data_iterator + sizeof(std::size_t));
4720  data_iterator += sizeof(std::size_t);
4721 
4722  // (b) unpack the data section in the indicated region
4723  auto received_data = Utilities::unpack<std::vector<DataType>>(
4724  data_iterator,
4725  data_iterator + size_of_received_data,
4726  /*enable_compression*/ false);
4727  data_iterator += size_of_received_data;
4728 
4729  // (c) check if the received data contained fewer entries than the
4730  // number of cells we identified in the beginning, in which case we
4731  // need to extract the boolean vector with the relevant information
4732  const std::vector<typename CellId::binary_type> &this_cell_list =
4733  neighbor_cell_list[status.MPI_SOURCE];
4734  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4735  std::vector<bool> cells_with_data;
4736  if (received_data.size() < this_cell_list.size())
4737  {
4738  cells_with_data = Utilities::unpack<std::vector<bool>>(
4739  data_iterator, receive.end(), /*enable_compression*/ false);
4740  AssertDimension(cells_with_data.size(), this_cell_list.size());
4741  }
4742 
4743  // (d) go through the received data and call the user-provided
4744  // unpack function
4745  auto received_data_iterator = received_data.begin();
4746  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4747  if (cells_with_data.empty() || cells_with_data[c])
4748  {
4750  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4751 
4752  MeshCellIteratorType cell(tria,
4753  tria_cell->level(),
4754  tria_cell->index(),
4755  &mesh);
4756 
4757  unpack(cell, *received_data_iterator);
4758  ++received_data_iterator;
4759  }
4760  }
4761 
4762  // make sure that all communication is finished
4763  // when we leave this function.
4764  if (requests.size() > 0)
4765  {
4766  const int ierr =
4767  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4768  AssertThrowMPI(ierr);
4769  }
4770  if (reply_requests.size() > 0)
4771  {
4772  const int ierr = MPI_Waitall(reply_requests.size(),
4773  reply_requests.data(),
4774  MPI_STATUSES_IGNORE);
4775  AssertThrowMPI(ierr);
4776  }
4777 
4778 
4779 # endif // DEAL_II_WITH_MPI
4780  }
4781 
4782  } // namespace internal
4783 
4784  template <typename DataType, typename MeshType>
4785  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
4786  inline void exchange_cell_data_to_ghosts(
4787  const MeshType & mesh,
4788  const std::function<std_cxx17::optional<DataType>(
4789  const typename MeshType::active_cell_iterator &)> &pack,
4790  const std::function<void(const typename MeshType::active_cell_iterator &,
4791  const DataType &)> & unpack,
4792  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4793  &cell_filter)
4794  {
4795 # ifndef DEAL_II_WITH_MPI
4796  (void)mesh;
4797  (void)pack;
4798  (void)unpack;
4799  (void)cell_filter;
4800  Assert(false, ExcNeedsMPI());
4801 # else
4802  internal::exchange_cell_data<DataType,
4803  MeshType,
4804  typename MeshType::active_cell_iterator>(
4805  mesh,
4806  pack,
4807  unpack,
4808  cell_filter,
4809  [&](const auto &process) {
4810  for (const auto &cell : mesh.active_cell_iterators())
4811  if (cell->is_ghost())
4812  process(cell, cell->subdomain_id());
4813  },
4814  [](const auto &tria) { return tria.ghost_owners(); });
4815 # endif
4816  }
4817 
4818 
4819 
4820  template <typename DataType, typename MeshType>
4821  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
4823  const MeshType & mesh,
4824  const std::function<std_cxx17::optional<DataType>(
4825  const typename MeshType::level_cell_iterator &)> &pack,
4826  const std::function<void(const typename MeshType::level_cell_iterator &,
4827  const DataType &)> & unpack,
4828  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4829  &cell_filter)
4830  {
4831 # ifndef DEAL_II_WITH_MPI
4832  (void)mesh;
4833  (void)pack;
4834  (void)unpack;
4835  (void)cell_filter;
4836  Assert(false, ExcNeedsMPI());
4837 # else
4838  internal::exchange_cell_data<DataType,
4839  MeshType,
4840  typename MeshType::level_cell_iterator>(
4841  mesh,
4842  pack,
4843  unpack,
4844  cell_filter,
4845  [&](const auto &process) {
4846  for (const auto &cell : mesh.cell_iterators())
4847  if (cell->is_ghost_on_level())
4848  process(cell, cell->level_subdomain_id());
4849  },
4850  [](const auto &tria) { return tria.level_ghost_owners(); });
4851 # endif
4852  }
4853 } // namespace GridTools
4854 
4855 #endif // DOXYGEN
4856 
4858 
4859 #endif
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:704
Definition: cell_id.h:72
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:81
typename VectorType::value_type value_type
Definition: grid_tools.h:3528
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:6926
void process_cell(const typename DoFHandler< dim >::active_cell_iterator &cell, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:6962
void process_sub_cell(const std::vector< value_type > &, const std::vector< Point< 1 >> &, const std::vector< unsigned int > &, const double, std::vector< Point< 1 >> &, std::vector< CellData< 1 >> &, const bool) const
Definition: grid_tools.h:3614
const unsigned int n_subdivisions
Definition: grid_tools.h:3656
MarchingCubeAlgorithm(const Mapping< dim, dim > &mapping, const FiniteElement< dim, dim > &fe, const unsigned int n_subdivisions=1, const double tolerance=1e-10)
Definition: grid_tools.cc:6874
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:6891
Definition: point.h:110
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
virtual MPI_Comm get_communicator() const
cell_iterator create_cell_iterator(const CellId &cell_id) const
cell_iterator begin(const unsigned int level=0) const
Definition: vector.h:109
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:107
#define DEAL_II_DEPRECATED
Definition: config.h:174
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:162
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 3 > center
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4618
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1355
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcMeshNotOrientable()
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1886
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:488
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:510
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:439
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:5108
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:5200
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:5133
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:5228
void random(DoFHandler< dim, spacedim > &dof_handler)
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 std::vector< bool > &marked_vertices, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition: grid_tools.cc:6113
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3915
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:304
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:3277
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2197
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4255
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:6593
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2347
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:5066
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2141
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:6389
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
void invert_all_negative_measure_cells(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:998
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:3189
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< bool > &marked_vertices={}, const double tolerance=1.e-10)
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:756
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:6370
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4546
double cell_measure(const std::vector< Point< dim >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< dim >::vertices_per_cell])
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:5410
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5381
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:4573
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2131
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 collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const unsigned 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 >())
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)
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:5755
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5348
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4361
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:2862
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)
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6716
std::size_t invert_cells_with_negative_measure(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:918
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3424
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4530
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:2021
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:5726
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)
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})
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)
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:4014
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:328
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3522
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:651
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:394
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:6434
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
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})
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3949
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
double volume(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:139
std::vector< types::subdomain_id > get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const std::vector< CellId > &cell_ids)
Definition: grid_tools.cc:4460
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:2379
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)
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:6530
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:83
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3978
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:2734
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:3571
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:410
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:4600
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5316
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, const std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim >> &first_cell)
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:550
bool orthogonal_equality(std::bitset< 3 > &orientation, const FaceIterator &face1, const FaceIterator &face2, const unsigned int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &local_points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const double tolerance=1e-10)
Definition: grid_tools.cc:5923
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
@ exchange_cell_data_request
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
@ exchange_cell_data_reply
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
T sum(const T &t, const MPI_Comm &mpi_communicator)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1334
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1493
constexpr bool is_triangulation_or_dof_handler
concept is_triangulation_or_dof_handler
Definition: hp.h:118
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13819
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const types::manifold_id flat_manifold_id
Definition: types.h:286
Definition: c++.h:141
Definition: types.h:33
unsigned int manifold_id
Definition: types.h:153
unsigned int global_dof_index
Definition: types.h:82
unsigned int subdomain_id
Definition: types.h:44
unsigned int boundary_id
Definition: types.h:141
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:144
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
void load(Archive &ar, const unsigned int)
Definition: grid_tools.h:3465
void serialize(Archive &archive, const unsigned int version)
void save(Archive &ar, const unsigned int) const
Definition: grid_tools.h:3439
std::vector< CellId > cell_ids
Definition: grid_tools.h:3422
std::bitset< 3 > orientation
Definition: grid_tools.h:2830
std::size_t memory_consumption() const
FullMatrix< double > matrix
Definition: grid_tools.h:2844
unsigned int face_idx[2]
Definition: grid_tools.h:2823
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1222
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1197
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
const ::Triangulation< dim, spacedim > & tria