Reference documentation for deal.II version GIT 92344eb39a 2023-03-23 02:40: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>
103  {
104  public:
105 #ifndef _MSC_VER
106  using type = typename MeshType::active_cell_iterator;
107 #else
109 #endif
110  };
111 
112 #ifdef _MSC_VER
113  template <int dim, int spacedim>
114  class ActiveCellIterator<dim, spacedim, DoFHandler<dim, spacedim>>
115  {
116  public:
117  using type =
119  };
120 #endif
121 } // namespace internal
122 
131 namespace GridTools
132 {
144  template <int dim, int spacedim>
145  double
147 
171  template <int dim, int spacedim>
172  double
174 
202  template <int dim, int spacedim>
203  double
205  const Mapping<dim, spacedim> & mapping);
206 
217  template <int dim, int spacedim>
218  double
221  const Mapping<dim, spacedim> & mapping =
222  (ReferenceCells::get_hypercube<dim>()
223 #ifndef _MSC_VER
224  .template get_default_linear_mapping<dim, spacedim>()
225 #else
226  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
227 #endif
228  ));
229 
240  template <int dim, int spacedim>
241  double
244  const Mapping<dim, spacedim> & mapping =
245  (ReferenceCells::get_hypercube<dim>()
246 #ifndef _MSC_VER
247  .template get_default_linear_mapping<dim, spacedim>()
248 #else
249  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
250 #endif
251  ));
252 
264  template <int dim>
265  DEAL_II_DEPRECATED double
267  const std::vector<Point<dim>> &all_vertices,
269 
289  template <int dim>
290  double
291  cell_measure(const std::vector<Point<dim>> & all_vertices,
293 
316  template <int dim, int spacedim>
317  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
319 
346  template <int dim>
350  const Quadrature<dim> & quadrature);
351 
359  template <int dim>
360  double
363  const Quadrature<dim> & quadrature);
364 
378  template <int dim, int spacedim>
381 
399  template <typename Iterator>
402  const Iterator & object,
404 
417  template <int dim, int spacedim>
418  std::
419  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
421 
444  template <int dim, int spacedim>
445  void
447  std::vector<CellData<dim>> & cells,
448  SubCellData & subcelldata);
449 
468  template <int dim, int spacedim>
469  void
470  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
471  std::vector<CellData<dim>> & cells,
472  SubCellData & subcelldata,
473  std::vector<unsigned int> & considered_vertices,
474  const double tol = 1e-12);
475 
483  template <int dim>
484  void
486  const double tol = 1e-12);
487 
506  template <int dim, int spacedim>
507  void
509  const std::vector<Point<spacedim>> &all_vertices,
510  std::vector<CellData<dim>> & cells);
511 
521  template <int dim, int spacedim>
522  std::size_t
524  const std::vector<Point<spacedim>> &all_vertices,
525  std::vector<CellData<dim>> & cells);
526 
537  template <int dim>
538  void
539  consistently_order_cells(std::vector<CellData<dim>> &cells);
540 
609  template <int dim, typename Transformation, int spacedim>
611  (std::invocable<Transformation, Point<spacedim>> &&
612  std::assignable_from<
613  Point<spacedim> &,
614  std::invoke_result_t<Transformation, Point<spacedim>>>))
615  void transform(const Transformation & transformation,
616  Triangulation<dim, spacedim> &triangulation);
617 
624  template <int dim, int spacedim>
625  void
626  shift(const Tensor<1, spacedim> & shift_vector,
627  Triangulation<dim, spacedim> &triangulation);
628 
629 
640  template <int dim, int spacedim>
641  void
642  rotate(const double angle, Triangulation<dim, spacedim> &triangulation);
643 
656  template <int dim>
657  void
658  rotate(const Tensor<1, 3, double> &axis,
659  const double angle,
660  Triangulation<dim, 3> & triangulation);
661 
677  template <int dim>
678  DEAL_II_DEPRECATED void
679  rotate(const double angle,
680  const unsigned int axis,
681  Triangulation<dim, 3> &triangulation);
682 
740  template <int dim>
741  void
742  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
743  Triangulation<dim> & tria,
744  const Function<dim, double> *coefficient = nullptr,
745  const bool solve_for_absolute_positions = false);
746 
752  template <int dim, int spacedim>
753  std::map<unsigned int, Point<spacedim>>
754  get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria);
755 
763  template <int dim, int spacedim>
764  void
765  scale(const double scaling_factor,
766  Triangulation<dim, spacedim> &triangulation);
767 
782  template <int dim, int spacedim>
783  void
785  const double factor,
786  Triangulation<dim, spacedim> &triangulation,
787  const bool keep_boundary = true,
788  const unsigned int seed = boost::random::mt19937::default_seed);
789 
823  template <int dim, int spacedim>
824  void
825  remove_hanging_nodes(Triangulation<dim, spacedim> &tria,
826  const bool isotropic = false,
827  const unsigned int max_iterations = 100);
828 
853  template <int dim, int spacedim>
854  void
855  remove_anisotropy(Triangulation<dim, spacedim> &tria,
856  const double max_ratio = 1.6180339887,
857  const unsigned int max_iterations = 5);
858 
948  template <int dim, int spacedim>
949  void
951  const double limit_angle_fraction = .75);
952 
1013  template <int dim, int spacedim>
1014 #ifndef DOXYGEN
1015  std::tuple<
1016  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1017  std::vector<std::vector<Point<dim>>>,
1018  std::vector<std::vector<unsigned int>>>
1019 #else
1020  return_type
1021 #endif
1023  const Cache<dim, spacedim> & cache,
1024  const std::vector<Point<spacedim>> &points,
1026  &cell_hint =
1028 
1062  template <int dim, int spacedim>
1063 #ifndef DOXYGEN
1064  std::tuple<
1065  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1066  std::vector<std::vector<Point<dim>>>,
1067  std::vector<std::vector<unsigned int>>,
1068  std::vector<unsigned int>>
1069 #else
1070  return_type
1071 #endif
1073  const Cache<dim, spacedim> & cache,
1074  const std::vector<Point<spacedim>> &points,
1076  &cell_hint =
1078 
1148  template <int dim, int spacedim>
1149 #ifndef DOXYGEN
1150  std::tuple<
1151  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1152  std::vector<std::vector<Point<dim>>>,
1153  std::vector<std::vector<unsigned int>>,
1154  std::vector<std::vector<Point<spacedim>>>,
1155  std::vector<std::vector<unsigned int>>>
1156 #else
1157  return_type
1158 #endif
1160  const GridTools::Cache<dim, spacedim> & cache,
1161  const std::vector<Point<spacedim>> & local_points,
1162  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1163  const double tolerance = 1e-10);
1164 
1165  namespace internal
1166  {
1181  template <int dim, int spacedim>
1183  {
1190  std::vector<std::tuple<std::pair<int, int>,
1191  unsigned int,
1192  unsigned int,
1193  Point<dim>,
1195  unsigned int>>
1197 
1201  std::vector<unsigned int> send_ranks;
1202 
1208  std::vector<unsigned int> send_ptrs;
1209 
1220  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1222 
1226  std::vector<unsigned int> recv_ranks;
1227 
1233  std::vector<unsigned int> recv_ptrs;
1234  };
1235 
1245  template <int dim, int spacedim>
1248  const GridTools::Cache<dim, spacedim> & cache,
1249  const std::vector<Point<spacedim>> & points,
1250  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1251  const std::vector<bool> & marked_vertices,
1252  const double tolerance,
1253  const bool perform_handshake,
1254  const bool enforce_unique_mapping = false);
1255 
1256  } // namespace internal
1257 
1294  template <int dim, int spacedim>
1295  std::map<unsigned int, Point<spacedim>>
1297  const Triangulation<dim, spacedim> &container,
1298  const Mapping<dim, spacedim> & mapping =
1299  (ReferenceCells::get_hypercube<dim>()
1300 #ifndef _MSC_VER
1301  .template get_default_linear_mapping<dim, spacedim>()
1302 #else
1303  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1304 #endif
1305  ));
1306 
1316  template <int spacedim>
1317  unsigned int
1318  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1319  const Point<spacedim> & p);
1320 
1344  template <int dim, template <int, int> class MeshType, int spacedim>
1345  unsigned int
1346  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1347  const Point<spacedim> & p,
1348  const std::vector<bool> & marked_vertices = {});
1349 
1373  template <int dim, template <int, int> class MeshType, int spacedim>
1374  unsigned int
1376  const MeshType<dim, spacedim> &mesh,
1377  const Point<spacedim> & p,
1378  const std::vector<bool> & marked_vertices = {});
1379 
1380 
1400  template <int dim, template <int, int> class MeshType, int spacedim>
1401 #ifndef _MSC_VER
1402  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1403 #else
1404  std::vector<
1405  typename ::internal::
1406  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1407 #endif
1408  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1409  const unsigned int vertex_index);
1410 
1473  template <int dim, template <int, int> class MeshType, int spacedim>
1474 #ifndef _MSC_VER
1475  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1476 #else
1477  std::pair<typename ::internal::
1478  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1479  Point<dim>>
1480 #endif
1482  const MeshType<dim, spacedim> &mesh,
1483  const Point<spacedim> & p,
1484  const std::vector<bool> &marked_vertices = {},
1485  const double tolerance = 1.e-10);
1486 
1494  template <int dim, template <int, int> class MeshType, int spacedim>
1495 #ifndef _MSC_VER
1496  typename MeshType<dim, spacedim>::active_cell_iterator
1497 #else
1498  typename ::internal::
1499  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1500 #endif
1501  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1502  const Point<spacedim> & p,
1503  const std::vector<bool> &marked_vertices = {},
1504  const double tolerance = 1.e-10);
1505 
1512  template <int dim, int spacedim>
1513  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1514  Point<dim>>
1516  const hp::MappingCollection<dim, spacedim> &mapping,
1517  const DoFHandler<dim, spacedim> & mesh,
1518  const Point<spacedim> & p,
1519  const double tolerance = 1.e-10);
1520 
1572  template <int dim, int spacedim>
1573  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1574  Point<dim>>
1576  const Cache<dim, spacedim> &cache,
1577  const Point<spacedim> & p,
1580  const std::vector<bool> &marked_vertices = {},
1581  const double tolerance = 1.e-10);
1582 
1596  template <int dim, template <int, int> class MeshType, int spacedim>
1597 #ifndef _MSC_VER
1598  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1599 #else
1600  std::pair<typename ::internal::
1601  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1602  Point<dim>>
1603 #endif
1605  const Mapping<dim, spacedim> & mapping,
1606  const MeshType<dim, spacedim> &mesh,
1607  const Point<spacedim> & p,
1608  const std::vector<
1609  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1611  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1612  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1613  typename MeshType<dim, spacedim>::active_cell_iterator(),
1614  const std::vector<bool> & marked_vertices = {},
1615  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1616  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1617  const double tolerance = 1.e-10,
1618  const RTree<
1619  std::pair<BoundingBox<spacedim>,
1621  *relevant_cell_bounding_boxes_rtree = nullptr);
1622 
1644  template <int dim, template <int, int> class MeshType, int spacedim>
1645 #ifndef _MSC_VER
1646  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1647  Point<dim>>>
1648 #else
1649  std::vector<std::pair<
1650  typename ::internal::
1651  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1652  Point<dim>>>
1653 #endif
1655  const Mapping<dim, spacedim> & mapping,
1656  const MeshType<dim, spacedim> &mesh,
1657  const Point<spacedim> & p,
1658  const double tolerance,
1659  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1660  Point<dim>> & first_cell);
1661 
1668  template <int dim, template <int, int> class MeshType, int spacedim>
1669 #ifndef _MSC_VER
1670  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1671  Point<dim>>>
1672 #else
1673  std::vector<std::pair<
1674  typename ::internal::
1675  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1676  Point<dim>>>
1677 #endif
1679  const Mapping<dim, spacedim> & mapping,
1680  const MeshType<dim, spacedim> &mesh,
1681  const Point<spacedim> & p,
1682  const double tolerance = 1e-10,
1683  const std::vector<bool> & marked_vertices = {});
1684 
1706  template <class MeshType>
1707  std::vector<typename MeshType::active_cell_iterator>
1708  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1709 
1734  template <class MeshType>
1735  void
1737  const typename MeshType::active_cell_iterator & cell,
1738  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1739 
1789  template <class MeshType>
1790  std::vector<typename MeshType::active_cell_iterator>
1792  const MeshType &mesh,
1793  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1794  &predicate);
1795 
1796 
1804  template <class MeshType>
1805  std::vector<typename MeshType::cell_iterator>
1807  const MeshType &mesh,
1808  const std::function<bool(const typename MeshType::cell_iterator &)>
1809  & predicate,
1810  const unsigned int level);
1811 
1812 
1825  template <class MeshType>
1826  std::vector<typename MeshType::active_cell_iterator>
1827  compute_ghost_cell_halo_layer(const MeshType &mesh);
1828 
1877  template <class MeshType>
1878  std::vector<typename MeshType::active_cell_iterator>
1880  const MeshType &mesh,
1881  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1882  & predicate,
1883  const double layer_thickness);
1884 
1907  template <class MeshType>
1908  std::vector<typename MeshType::active_cell_iterator>
1909  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1910  const double layer_thickness);
1911 
1927  template <class MeshType>
1928  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1930  const MeshType &mesh,
1931  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1932  &predicate);
1933 
2004  template <class MeshType>
2005  std::vector<BoundingBox<MeshType::space_dimension>>
2007  const MeshType &mesh,
2008  const std::function<bool(const typename MeshType::active_cell_iterator &)>
2009  & predicate,
2010  const unsigned int refinement_level = 0,
2011  const bool allow_merge = false,
2012  const unsigned int max_boxes = numbers::invalid_unsigned_int);
2013 
2041  template <int spacedim>
2042 #ifndef DOXYGEN
2043  std::tuple<std::vector<std::vector<unsigned int>>,
2044  std::map<unsigned int, unsigned int>,
2045  std::map<unsigned int, std::vector<unsigned int>>>
2046 #else
2047  return_type
2048 #endif
2050  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
2051  const std::vector<Point<spacedim>> & points);
2052 
2053 
2088  template <int spacedim>
2089 #ifndef DOXYGEN
2090  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2091  std::map<unsigned int, unsigned int>,
2092  std::map<unsigned int, std::vector<unsigned int>>>
2093 #else
2094  return_type
2095 #endif
2097  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2098  const std::vector<Point<spacedim>> & points);
2099 
2100 
2109  template <int dim, int spacedim>
2110  std::vector<
2111  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2113 
2126  template <int dim, int spacedim>
2127  std::vector<std::vector<Tensor<1, spacedim>>>
2129  const Triangulation<dim, spacedim> &mesh,
2130  const std::vector<
2132  &vertex_to_cells);
2133 
2134 
2142  template <int dim, int spacedim>
2143  unsigned int
2146  const Point<spacedim> & position,
2147  const Mapping<dim, spacedim> & mapping =
2148  (ReferenceCells::get_hypercube<dim>()
2149 #ifndef _MSC_VER
2150  .template get_default_linear_mapping<dim, spacedim>()
2151 #else
2152  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2153 #endif
2154  ));
2155 
2167  template <int dim, int spacedim>
2168  std::map<unsigned int, types::global_vertex_index>
2171 
2183  template <int dim, int spacedim>
2184  std::pair<unsigned int, double>
2187 
2202  template <int dim, int spacedim>
2203  void
2206  DynamicSparsityPattern & connectivity);
2207 
2216  template <int dim, int spacedim>
2217  void
2220  DynamicSparsityPattern & connectivity);
2221 
2230  template <int dim, int spacedim>
2231  void
2234  const unsigned int level,
2235  DynamicSparsityPattern & connectivity);
2236 
2257  template <int dim, int spacedim>
2258  void
2259  partition_triangulation(const unsigned int n_partitions,
2261  const SparsityTools::Partitioner partitioner =
2263 
2274  template <int dim, int spacedim>
2275  void
2276  partition_triangulation(const unsigned int n_partitions,
2277  const std::vector<unsigned int> &cell_weights,
2279  const SparsityTools::Partitioner partitioner =
2281 
2327  template <int dim, int spacedim>
2328  void
2329  partition_triangulation(const unsigned int n_partitions,
2330  const SparsityPattern & cell_connection_graph,
2332  const SparsityTools::Partitioner partitioner =
2334 
2345  template <int dim, int spacedim>
2346  void
2347  partition_triangulation(const unsigned int n_partitions,
2348  const std::vector<unsigned int> &cell_weights,
2349  const SparsityPattern & cell_connection_graph,
2351  const SparsityTools::Partitioner partitioner =
2353 
2368  template <int dim, int spacedim>
2369  void
2370  partition_triangulation_zorder(const unsigned int n_partitions,
2372  const bool group_siblings = true);
2373 
2385  template <int dim, int spacedim>
2386  void
2388 
2396  template <int dim, int spacedim>
2397  std::vector<types::subdomain_id>
2399  const std::vector<CellId> & cell_ids);
2400 
2411  template <int dim, int spacedim>
2412  void
2414  std::vector<types::subdomain_id> & subdomain);
2415 
2430  template <int dim, int spacedim>
2431  unsigned int
2434  const types::subdomain_id subdomain);
2435 
2465  template <int dim, int spacedim>
2466  std::vector<bool>
2468 
2508  template <typename MeshType>
2509  std::list<std::pair<typename MeshType::cell_iterator,
2510  typename MeshType::cell_iterator>>
2511  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2512 
2522  template <int dim, int spacedim>
2523  bool
2525  const Triangulation<dim, spacedim> &mesh_2);
2526 
2536  template <typename MeshType>
2537  bool
2538  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2539 
2561  template <int dim, int spacedim>
2565  & distorted_cells,
2567 
2568 
2569 
2617  template <class MeshType>
2618  std::vector<typename MeshType::active_cell_iterator>
2619  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2620 
2621 
2643  template <class Container>
2644  std::vector<typename Container::cell_iterator>
2646  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2647 
2714  template <class Container>
2715  void
2717  const std::vector<typename Container::active_cell_iterator> &patch,
2719  &local_triangulation,
2720  std::map<
2721  typename Triangulation<Container::dimension,
2722  Container::space_dimension>::active_cell_iterator,
2723  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2724 
2756  template <int dim, int spacedim>
2757  std::map<
2759  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2761 
2762 
2775  template <typename CellIterator>
2777  {
2781  CellIterator cell[2];
2782 
2787  unsigned int face_idx[2];
2788 
2794  std::bitset<3> orientation;
2795 
2809 
2813  std::size_t
2815  };
2816 
2817 
2881  template <typename FaceIterator>
2882  bool
2884  std::bitset<3> & orientation,
2885  const FaceIterator & face1,
2886  const FaceIterator & face2,
2887  const unsigned int direction,
2891 
2892 
2896  template <typename FaceIterator>
2897  bool
2899  const FaceIterator & face1,
2900  const FaceIterator & face2,
2901  const unsigned int direction,
2905 
2906 
2963  template <typename MeshType>
2964  void
2966  const MeshType & mesh,
2967  const types::boundary_id b_id1,
2968  const types::boundary_id b_id2,
2969  const unsigned int direction,
2971  & matched_pairs,
2972  const Tensor<1, MeshType::space_dimension> &offset =
2975 
2976 
2999  template <typename MeshType>
3000  void
3002  const MeshType & mesh,
3003  const types::boundary_id b_id,
3004  const unsigned int direction,
3006  & matched_pairs,
3007  const ::Tensor<1, MeshType::space_dimension> &offset =
3010 
3037  template <int dim, int spacedim>
3038  void
3040  const bool reset_boundary_ids = false);
3041 
3063  template <int dim, int spacedim>
3064  void
3066  const std::vector<types::boundary_id> &src_boundary_ids,
3067  const std::vector<types::manifold_id> &dst_manifold_ids,
3069  const std::vector<types::boundary_id> &reset_boundary_ids = {});
3070 
3100  template <int dim, int spacedim>
3101  void
3103  const bool compute_face_ids = false);
3104 
3129  template <int dim, int spacedim>
3130  void
3133  const std::function<types::manifold_id(
3134  const std::set<types::manifold_id> &)> &disambiguation_function =
3135  [](const std::set<types::manifold_id> &manifold_ids) {
3136  if (manifold_ids.size() == 1)
3137  return *manifold_ids.begin();
3138  else
3140  },
3141  bool overwrite_only_flat_manifold_ids = true);
3228  template <typename DataType, typename MeshType>
3229  void
3231  const MeshType & mesh,
3232  const std::function<std_cxx17::optional<DataType>(
3233  const typename MeshType::active_cell_iterator &)> &pack,
3234  const std::function<void(const typename MeshType::active_cell_iterator &,
3235  const DataType &)> & unpack,
3236  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3237  &cell_filter =
3239 
3250  template <typename DataType, typename MeshType>
3251  void
3253  const MeshType & mesh,
3254  const std::function<std_cxx17::optional<DataType>(
3255  const typename MeshType::level_cell_iterator &)> &pack,
3256  const std::function<void(const typename MeshType::level_cell_iterator &,
3257  const DataType &)> & unpack,
3258  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3260  true});
3261 
3262  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3263  * boxes @p local_bboxes.
3264  *
3265  * This function is meant to exchange bounding boxes describing the locally
3266  * owned cells in a distributed triangulation obtained with the function
3267  * GridTools::compute_mesh_predicate_bounding_box .
3268  *
3269  * The output vector's size is the number of processes of the MPI
3270  * communicator:
3271  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3272  */
3273  template <int spacedim>
3274  std::vector<std::vector<BoundingBox<spacedim>>>
3276  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3277  const MPI_Comm & mpi_communicator);
3278 
3311  template <int spacedim>
3314  const std::vector<BoundingBox<spacedim>> &local_description,
3315  const MPI_Comm & mpi_communicator);
3316 
3334  template <int dim, int spacedim>
3335  void
3338  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3339  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3340 
3347  template <int dim, int spacedim>
3348  std::map<unsigned int, std::set<::types::subdomain_id>>
3351 
3367  template <int dim, typename T>
3369  {
3373  std::vector<CellId> cell_ids;
3374 
3378  std::vector<T> data;
3379 
3388  template <class Archive>
3389  void
3390  save(Archive &ar, const unsigned int) const
3391  {
3392  // Implement the code inline as some compilers do otherwise complain
3393  // about the use of a deprecated interface.
3394  Assert(cell_ids.size() == data.size(),
3395  ExcDimensionMismatch(cell_ids.size(), data.size()));
3396  // archive the cellids in an efficient binary format
3397  const std::size_t n_cells = cell_ids.size();
3398  ar & n_cells;
3399  for (const auto &id : cell_ids)
3400  {
3401  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3402  ar & binary_cell_id;
3403  }
3404 
3405  ar &data;
3406  }
3407 
3414  template <class Archive>
3415  void
3416  load(Archive &ar, const unsigned int)
3417  {
3418  // Implement the code inline as some compilers do otherwise complain
3419  // about the use of a deprecated interface.
3420  std::size_t n_cells;
3421  ar & n_cells;
3422  cell_ids.clear();
3423  cell_ids.reserve(n_cells);
3424  for (unsigned int c = 0; c < n_cells; ++c)
3425  {
3426  CellId::binary_type value;
3427  ar & value;
3428  cell_ids.emplace_back(value);
3429  }
3430  ar &data;
3431  }
3432 
3433 
3434 #ifdef DOXYGEN
3440  template <class Archive>
3441  void
3442  serialize(Archive &archive, const unsigned int version);
3443 #else
3444  // This macro defines the serialize() method that is compatible with
3445  // the templated save() and load() method that have been implemented.
3446  BOOST_SERIALIZATION_SPLIT_MEMBER()
3447 #endif
3448  };
3449 
3450 
3451 
3472  template <int dim, typename VectorType>
3474  {
3475  public:
3479  using value_type = typename VectorType::value_type;
3480 
3484  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3485  const FiniteElement<dim, dim> &fe,
3486  const unsigned int n_subdivisions = 1,
3487  const double tolerance = 1e-10);
3488 
3499  void
3500  process(const DoFHandler<dim> & background_dof_handler,
3501  const VectorType & ls_vector,
3502  const double iso_level,
3503  std::vector<Point<dim>> &vertices,
3504  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3505 
3510  void
3511  process(const DoFHandler<dim> & background_dof_handler,
3512  const VectorType & ls_vector,
3513  const double iso_level,
3514  std::vector<Point<dim>> &vertices) const;
3515 
3525  void
3527  const VectorType & ls_vector,
3528  const double iso_level,
3529  std::vector<Point<dim>> & vertices,
3530  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3536  void
3538  const VectorType & ls_vector,
3539  const double iso_level,
3540  std::vector<Point<dim>> &vertices) const;
3541 
3542  private:
3547  static Quadrature<dim>
3548  create_quadrature_rule(const unsigned int n_subdivisions);
3549 
3553  void
3554  process_cell(std::vector<value_type> & ls_values,
3555  const std::vector<Point<dim>> & points,
3556  const double iso_level,
3557  std::vector<Point<dim>> & vertices,
3558  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
3559  const bool write_back_cell_data = true) const;
3560 
3564  void
3565  process_sub_cell(const std::vector<value_type> &,
3566  const std::vector<Point<1>> &,
3567  const std::vector<unsigned int> &,
3568  const double,
3569  std::vector<Point<1>> &,
3570  std::vector<CellData<1>> &,
3571  const bool) const
3572  {
3573  AssertThrow(false, ExcNotImplemented());
3574  }
3575 
3582  void
3583  process_sub_cell(const std::vector<value_type> & ls_values,
3584  const std::vector<Point<2>> & points,
3585  const std::vector<unsigned int> &mask,
3586  const double iso_level,
3587  std::vector<Point<2>> & vertices,
3588  std::vector<CellData<1>> & cells,
3589  const bool write_back_cell_data) const;
3590 
3594  void
3595  process_sub_cell(const std::vector<value_type> & ls_values,
3596  const std::vector<Point<3>> & points,
3597  const std::vector<unsigned int> &mask,
3598  const double iso_level,
3599  std::vector<Point<3>> & vertices,
3600  std::vector<CellData<2>> & cells,
3601  const bool write_back_cell_data) const;
3602 
3607  const unsigned int n_subdivisions;
3608 
3613  const double tolerance;
3614 
3620  };
3621 
3622 
3623 
3633  int,
3634  << "The number of partitions you gave is " << arg1
3635  << ", but must be greater than zero.");
3640  int,
3641  << "The subdomain id " << arg1
3642  << " has no cells associated with it.");
3647 
3652  double,
3653  << "The scaling factor must be positive, but it is " << arg1
3654  << '.');
3655 
3660  unsigned int,
3661  << "The given vertex with index " << arg1
3662  << " is not used in the given triangulation.");
3663 
3666 } /*namespace GridTools*/
3667 
3668 
3679  "The edges of the mesh are not consistently orientable.");
3680 
3681 
3682 
3683 /* ----------------- Template function --------------- */
3684 
3685 #ifndef DOXYGEN
3686 
3687 namespace GridTools
3688 {
3689  template <int dim>
3690  double
3691  cell_measure(
3692  const std::vector<Point<dim>> &all_vertices,
3693  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3694  {
3695  // We forward call to the ArrayView version:
3696  const ArrayView<const unsigned int> view(
3698  return cell_measure(all_vertices, view);
3699  }
3700 
3701 
3702 
3703  // This specialization is defined here so that the general template in the
3704  // source file doesn't need to have further 1d overloads for the internal
3705  // functions it calls.
3706  template <>
3710  {
3711  return {};
3712  }
3713 
3714 
3715 
3716  template <int dim, typename Transformation, int spacedim>
3718  (std::invocable<Transformation, Point<spacedim>> &&
3719  std::assignable_from<
3720  Point<spacedim> &,
3721  std::invoke_result_t<Transformation, Point<spacedim>>>))
3722  void transform(const Transformation & transformation,
3723  Triangulation<dim, spacedim> &triangulation)
3724  {
3725  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3726 
3727  // loop over all active cells, and
3728  // transform those vertices that
3729  // have not yet been touched. note
3730  // that we get to all vertices in
3731  // the triangulation by only
3732  // visiting the active cells.
3734  cell = triangulation.begin_active(),
3735  endc = triangulation.end();
3736  for (; cell != endc; ++cell)
3737  for (const unsigned int v : cell->vertex_indices())
3738  if (treated_vertices[cell->vertex_index(v)] == false)
3739  {
3740  // transform this vertex
3741  cell->vertex(v) = transformation(cell->vertex(v));
3742  // and mark it as treated
3743  treated_vertices[cell->vertex_index(v)] = true;
3744  };
3745 
3746 
3747  // now fix any vertices on hanging nodes so that we don't create any holes
3748  if (dim == 2)
3749  {
3751  cell = triangulation.begin_active(),
3752  endc = triangulation.end();
3753  for (; cell != endc; ++cell)
3754  for (const unsigned int face : cell->face_indices())
3755  if (cell->face(face)->has_children() &&
3756  !cell->face(face)->at_boundary())
3757  {
3758  Assert(cell->reference_cell() ==
3759  ReferenceCells::get_hypercube<dim>(),
3760  ExcNotImplemented());
3761 
3762  // this line has children
3763  cell->face(face)->child(0)->vertex(1) =
3764  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3765  2;
3766  }
3767  }
3768  else if (dim == 3)
3769  {
3771  cell = triangulation.begin_active(),
3772  endc = triangulation.end();
3773  for (; cell != endc; ++cell)
3774  for (const unsigned int face : cell->face_indices())
3775  if (cell->face(face)->has_children() &&
3776  !cell->face(face)->at_boundary())
3777  {
3778  Assert(cell->reference_cell() ==
3779  ReferenceCells::get_hypercube<dim>(),
3780  ExcNotImplemented());
3781 
3782  // this face has hanging nodes
3783  cell->face(face)->child(0)->vertex(1) =
3784  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3785  2.0;
3786  cell->face(face)->child(0)->vertex(2) =
3787  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3788  2.0;
3789  cell->face(face)->child(1)->vertex(3) =
3790  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3791  2.0;
3792  cell->face(face)->child(2)->vertex(3) =
3793  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3794  2.0;
3795 
3796  // center of the face
3797  cell->face(face)->child(0)->vertex(3) =
3798  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3799  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3800  4.0;
3801  }
3802  }
3803 
3804  // Make sure FEValues notices that the mesh has changed
3805  triangulation.signals.mesh_movement();
3806  }
3807 
3808 
3809 
3810  template <class MeshType>
3811  std::vector<typename MeshType::active_cell_iterator>
3812  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3813  {
3814  std::vector<typename MeshType::active_cell_iterator> child_cells;
3815 
3816  if (cell->has_children())
3817  {
3818  for (unsigned int child = 0; child < cell->n_children(); ++child)
3819  if (cell->child(child)->has_children())
3820  {
3821  const std::vector<typename MeshType::active_cell_iterator>
3822  children = get_active_child_cells<MeshType>(cell->child(child));
3823  child_cells.insert(child_cells.end(),
3824  children.begin(),
3825  children.end());
3826  }
3827  else
3828  child_cells.push_back(cell->child(child));
3829  }
3830 
3831  return child_cells;
3832  }
3833 
3834 
3835 
3836  template <class MeshType>
3837  void
3839  const typename MeshType::active_cell_iterator & cell,
3840  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3841  {
3842  active_neighbors.clear();
3843  for (const unsigned int n : cell->face_indices())
3844  if (!cell->at_boundary(n))
3845  {
3846  if (MeshType::dimension == 1)
3847  {
3848  // check children of neighbor. note
3849  // that in 1d children of the neighbor
3850  // may be further refined. In 1d the
3851  // case is simple since we know what
3852  // children bound to the present cell
3853  typename MeshType::cell_iterator neighbor_child =
3854  cell->neighbor(n);
3855  if (!neighbor_child->is_active())
3856  {
3857  while (neighbor_child->has_children())
3858  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3859 
3860  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3861  ExcInternalError());
3862  }
3863  active_neighbors.push_back(neighbor_child);
3864  }
3865  else
3866  {
3867  if (cell->face(n)->has_children())
3868  // this neighbor has children. find
3869  // out which border to the present
3870  // cell
3871  for (unsigned int c = 0;
3872  c < cell->face(n)->n_active_descendants();
3873  ++c)
3874  active_neighbors.push_back(
3875  cell->neighbor_child_on_subface(n, c));
3876  else
3877  {
3878  // the neighbor must be active
3879  // himself
3880  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3881  active_neighbors.push_back(cell->neighbor(n));
3882  }
3883  }
3884  }
3885  }
3886 
3887 
3888 
3889  template <typename CellIterator>
3890  std::size_t
3892  {
3893  return sizeof(*this) + matrix.memory_consumption();
3894  }
3895 
3896 
3897 
3898  namespace internal
3899  {
3900  namespace ProjectToObject
3901  {
3914  struct CrossDerivative
3915  {
3916  const unsigned int direction_0;
3917  const unsigned int direction_1;
3918 
3919  CrossDerivative(const unsigned int d0, const unsigned int d1);
3920  };
3921 
3922  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3923  const unsigned int d1)
3924  : direction_0(d0)
3925  , direction_1(d1)
3926  {}
3927 
3928 
3929 
3934  template <typename F>
3935  inline auto
3936  centered_first_difference(const double center,
3937  const double step,
3938  const F &f) -> decltype(f(center) - f(center))
3939  {
3940  return (f(center + step) - f(center - step)) / (2.0 * step);
3941  }
3942 
3943 
3944 
3949  template <typename F>
3950  inline auto
3951  centered_second_difference(const double center,
3952  const double step,
3953  const F &f) -> decltype(f(center) - f(center))
3954  {
3955  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3956  (step * step);
3957  }
3958 
3959 
3960 
3970  template <int structdim, typename F>
3971  inline auto
3972  cross_stencil(
3973  const CrossDerivative cross_derivative,
3975  const double step,
3976  const F &f) -> decltype(f(center) - f(center))
3977  {
3979  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3980  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3981  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3982  1.0 / 3.0 * f(center - simplex_vector) +
3983  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3984  step;
3985  }
3986 
3987 
3988 
3995  template <int spacedim, int structdim, typename F>
3996  inline double
3997  gradient_entry(
3998  const unsigned int row_n,
3999  const unsigned int dependent_direction,
4000  const Point<spacedim> &p0,
4002  const double step,
4003  const F & f)
4004  {
4006  dependent_direction <
4008  ExcMessage("This function assumes that the last weight is a "
4009  "dependent variable (and hence we cannot take its "
4010  "derivative directly)."));
4011  Assert(row_n != dependent_direction,
4012  ExcMessage(
4013  "We cannot differentiate with respect to the variable "
4014  "that is assumed to be dependent."));
4015 
4016  const Point<spacedim> manifold_point = f(center);
4017  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
4018  {row_n, dependent_direction}, center, step, f);
4019  double entry = 0.0;
4020  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
4021  entry +=
4022  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
4023  return entry;
4024  }
4025 
4031  template <typename Iterator, int spacedim, int structdim>
4033  project_to_d_linear_object(const Iterator & object,
4034  const Point<spacedim> &trial_point)
4035  {
4036  // let's look at this for simplicity for a quad (structdim==2) in a
4037  // space with spacedim>2 (notate trial_point by y): all points on the
4038  // surface are given by
4039  // x(\xi) = sum_i v_i phi_x(\xi)
4040  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
4041  // reference coordinates of the quad. so what we are trying to do is
4042  // find a point x on the surface that is closest to the point y. there
4043  // are different ways to solve this problem, but in the end it's a
4044  // nonlinear problem and we have to find reference coordinates \xi so
4045  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
4046  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
4047  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
4048  // have to use a Newton method to find the answer. This leads to the
4049  // following formulation of Newton steps:
4050  //
4051  // Given \xi_k, find \delta\xi_k so that
4052  // H_k \delta\xi_k = - F_k
4053  // where H_k is an approximation to the second derivatives of J at
4054  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
4055  // number of times until the right hand side is small enough. As a
4056  // stopping criterion, we terminate if ||\delta\xi||<eps.
4057  //
4058  // As for the Hessian, the best choice would be
4059  // H_k = J''(\xi_k)
4060  // but we'll opt for the simpler Gauss-Newton form
4061  // H_k = A^T A
4062  // i.e.
4063  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
4064  // \partial_n phi_i *
4065  // \partial_m phi_j
4066  // we start at xi=(0.5, 0.5).
4067  Point<structdim> xi;
4068  for (unsigned int d = 0; d < structdim; ++d)
4069  xi[d] = 0.5;
4070 
4071  Point<spacedim> x_k;
4072  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4073  x_k += object->vertex(i) *
4075 
4076  do
4077  {
4079  for (const unsigned int i :
4081  F_k +=
4082  (x_k - trial_point) * object->vertex(i) *
4084  i);
4085 
4087  for (const unsigned int i :
4089  for (const unsigned int j :
4091  {
4094  xi, i),
4096  xi, j));
4097  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
4098  }
4099 
4100  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
4101  xi += delta_xi;
4102 
4103  x_k = Point<spacedim>();
4104  for (const unsigned int i :
4106  x_k += object->vertex(i) *
4108 
4109  if (delta_xi.norm() < 1e-7)
4110  break;
4111  }
4112  while (true);
4113 
4114  return x_k;
4115  }
4116  } // namespace ProjectToObject
4117  } // namespace internal
4118 
4119 
4120 
4121  namespace internal
4122  {
4123  // We hit an internal compiler error in ICC 15 if we define this as a lambda
4124  // inside the project_to_object function below.
4125  template <int structdim>
4126  inline bool
4127  weights_are_ok(
4129  {
4130  // clang has trouble figuring out structdim here, so define it
4131  // again:
4132  static const std::size_t n_vertices_per_cell =
4134  n_independent_components;
4135  std::array<double, n_vertices_per_cell> copied_weights;
4136  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4137  {
4138  copied_weights[i] = v[i];
4139  if (v[i] < 0.0 || v[i] > 1.0)
4140  return false;
4141  }
4142 
4143  // check the sum: try to avoid some roundoff errors by summing in order
4144  std::sort(copied_weights.begin(), copied_weights.end());
4145  const double sum =
4146  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4147  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4148  }
4149  } // namespace internal
4150 
4151  template <typename Iterator>
4154  const Iterator & object,
4156  {
4157  const int spacedim = Iterator::AccessorType::space_dimension;
4158  const int structdim = Iterator::AccessorType::structure_dimension;
4159 
4160  Point<spacedim> projected_point = trial_point;
4161 
4162  if (structdim >= spacedim)
4163  return projected_point;
4164  else if (structdim == 1 || structdim == 2)
4165  {
4166  using namespace internal::ProjectToObject;
4167  // Try to use the special flat algorithm for quads (this is better
4168  // than the general algorithm in 3d). This does not take into account
4169  // whether projected_point is outside the quad, but we optimize along
4170  // lines below anyway:
4171  const int dim = Iterator::AccessorType::dimension;
4172  const Manifold<dim, spacedim> &manifold = object->get_manifold();
4173  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4174  &manifold) != nullptr)
4175  {
4176  projected_point =
4177  project_to_d_linear_object<Iterator, spacedim, structdim>(
4178  object, trial_point);
4179  }
4180  else
4181  {
4182  // We want to find a point on the convex hull (defined by the
4183  // vertices of the object and the manifold description) that is
4184  // relatively close to the trial point. This has a few issues:
4185  //
4186  // 1. For a general convex hull we are not guaranteed that a unique
4187  // minimum exists.
4188  // 2. The independent variables in the optimization process are the
4189  // weights given to Manifold::get_new_point, which must sum to 1,
4190  // so we cannot use standard finite differences to approximate a
4191  // gradient.
4192  //
4193  // There is not much we can do about 1., but for 2. we can derive
4194  // finite difference stencils that work on a structdim-dimensional
4195  // simplex and rewrite the optimization problem to use those
4196  // instead. Consider the structdim 2 case and let
4197  //
4198  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4199  // c2, c3})
4200  //
4201  // where {c0, c1, c2, c3} are the weights for the four vertices on
4202  // the quadrilateral. We seek to minimize the Euclidean distance
4203  // between F(...) and trial_point. We can solve for c3 in terms of
4204  // the other weights and get, for one coordinate direction
4205  //
4206  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4207  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4208  //
4209  // where we substitute back in for c3 after taking the
4210  // derivative. We can compute a stencil for the cross derivative
4211  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4212  // (and gradient_entry computes the sum over the independent
4213  // variables). Below, we somewhat arbitrarily pick the last
4214  // component as the dependent one.
4215  //
4216  // Since we can now calculate derivatives of the objective
4217  // function we can use gradient descent to minimize it.
4218  //
4219  // Of course, this is much simpler in the structdim = 1 case (we
4220  // could rewrite the projection as a 1d optimization problem), but
4221  // to reduce the potential for bugs we use the same code in both
4222  // cases.
4223  const double step_size = object->diameter() / 64.0;
4224 
4225  constexpr unsigned int n_vertices_per_cell =
4227 
4228  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4229  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4230  ++vertex_n)
4231  vertices[vertex_n] = object->vertex(vertex_n);
4232 
4233  auto get_point_from_weights =
4234  [&](const Tensor<1, n_vertices_per_cell> &weights)
4235  -> Point<spacedim> {
4236  return object->get_manifold().get_new_point(
4237  make_array_view(vertices.begin(), vertices.end()),
4238  make_array_view(weights.begin_raw(), weights.end_raw()));
4239  };
4240 
4241  // pick the initial weights as (normalized) inverse distances from
4242  // the trial point:
4243  Tensor<1, n_vertices_per_cell> guess_weights;
4244  double guess_weights_sum = 0.0;
4245  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4246  ++vertex_n)
4247  {
4248  const double distance =
4249  vertices[vertex_n].distance(trial_point);
4250  if (distance == 0.0)
4251  {
4252  guess_weights = 0.0;
4253  guess_weights[vertex_n] = 1.0;
4254  guess_weights_sum = 1.0;
4255  break;
4256  }
4257  else
4258  {
4259  guess_weights[vertex_n] = 1.0 / distance;
4260  guess_weights_sum += guess_weights[vertex_n];
4261  }
4262  }
4263  guess_weights /= guess_weights_sum;
4264  Assert(internal::weights_are_ok<structdim>(guess_weights),
4265  ExcInternalError());
4266 
4267  // The optimization algorithm consists of two parts:
4268  //
4269  // 1. An outer loop where we apply the gradient descent algorithm.
4270  // 2. An inner loop where we do a line search to find the optimal
4271  // length of the step one should take in the gradient direction.
4272  //
4273  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4274  {
4275  const unsigned int dependent_direction =
4276  n_vertices_per_cell - 1;
4277  Tensor<1, n_vertices_per_cell> current_gradient;
4278  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4279  ++row_n)
4280  {
4281  if (row_n != dependent_direction)
4282  {
4283  current_gradient[row_n] =
4284  gradient_entry<spacedim, structdim>(
4285  row_n,
4286  dependent_direction,
4287  trial_point,
4288  guess_weights,
4289  step_size,
4290  get_point_from_weights);
4291 
4292  current_gradient[dependent_direction] -=
4293  current_gradient[row_n];
4294  }
4295  }
4296 
4297  // We need to travel in the -gradient direction, as noted
4298  // above, but we may not want to take a full step in that
4299  // direction; instead, guess that we will go -0.5*gradient and
4300  // do quasi-Newton iteration to pick the best multiplier. The
4301  // goal is to find a scalar alpha such that
4302  //
4303  // F(x - alpha g)
4304  //
4305  // is minimized, where g is the gradient and F is the
4306  // objective function. To find the optimal value we find roots
4307  // of the derivative of the objective function with respect to
4308  // alpha by Newton iteration, where we approximate the first
4309  // and second derivatives of F(x - alpha g) with centered
4310  // finite differences.
4311  double gradient_weight = -0.5;
4312  auto gradient_weight_objective_function =
4313  [&](const double gradient_weight_guess) -> double {
4314  return (trial_point -
4315  get_point_from_weights(guess_weights +
4316  gradient_weight_guess *
4317  current_gradient))
4318  .norm_square();
4319  };
4320 
4321  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4322  {
4323  const double update_numerator = centered_first_difference(
4324  gradient_weight,
4325  step_size,
4326  gradient_weight_objective_function);
4327  const double update_denominator =
4328  centered_second_difference(
4329  gradient_weight,
4330  step_size,
4331  gradient_weight_objective_function);
4332 
4333  // avoid division by zero. Note that we limit the gradient
4334  // weight below
4335  if (std::abs(update_denominator) == 0.0)
4336  break;
4337  gradient_weight =
4338  gradient_weight - update_numerator / update_denominator;
4339 
4340  // Put a fairly lenient bound on the largest possible
4341  // gradient (things tend to be locally flat, so the gradient
4342  // itself is usually small)
4343  if (std::abs(gradient_weight) > 10)
4344  {
4345  gradient_weight = -10.0;
4346  break;
4347  }
4348  }
4349 
4350  // It only makes sense to take convex combinations with weights
4351  // between zero and one. If the update takes us outside of this
4352  // region then rescale the update to stay within the region and
4353  // try again
4354  Tensor<1, n_vertices_per_cell> tentative_weights =
4355  guess_weights + gradient_weight * current_gradient;
4356 
4357  double new_gradient_weight = gradient_weight;
4358  for (unsigned int iteration_count = 0; iteration_count < 40;
4359  ++iteration_count)
4360  {
4361  if (internal::weights_are_ok<structdim>(tentative_weights))
4362  break;
4363 
4364  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4365  {
4366  if (tentative_weights[i] < 0.0)
4367  {
4368  tentative_weights -=
4369  (tentative_weights[i] / current_gradient[i]) *
4370  current_gradient;
4371  }
4372  if (tentative_weights[i] < 0.0 ||
4373  1.0 < tentative_weights[i])
4374  {
4375  new_gradient_weight /= 2.0;
4376  tentative_weights =
4377  guess_weights +
4378  new_gradient_weight * current_gradient;
4379  }
4380  }
4381  }
4382 
4383  // the update might still send us outside the valid region, so
4384  // check again and quit if the update is still not valid
4385  if (!internal::weights_are_ok<structdim>(tentative_weights))
4386  break;
4387 
4388  // if we cannot get closer by traveling in the gradient
4389  // direction then quit
4390  if (get_point_from_weights(tentative_weights)
4391  .distance(trial_point) <
4392  get_point_from_weights(guess_weights).distance(trial_point))
4393  guess_weights = tentative_weights;
4394  else
4395  break;
4396  Assert(internal::weights_are_ok<structdim>(guess_weights),
4397  ExcInternalError());
4398  }
4399  Assert(internal::weights_are_ok<structdim>(guess_weights),
4400  ExcInternalError());
4401  projected_point = get_point_from_weights(guess_weights);
4402  }
4403 
4404  // if structdim == 2 and the optimal point is not on the interior then
4405  // we may be able to get a more accurate result by projecting onto the
4406  // lines.
4407  if (structdim == 2)
4408  {
4409  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4410  line_projections;
4411  for (unsigned int line_n = 0;
4412  line_n < GeometryInfo<structdim>::lines_per_cell;
4413  ++line_n)
4414  {
4415  line_projections[line_n] =
4416  project_to_object(object->line(line_n), trial_point);
4417  }
4418  std::sort(line_projections.begin(),
4419  line_projections.end(),
4420  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4421  return a.distance(trial_point) <
4422  b.distance(trial_point);
4423  });
4424  if (line_projections[0].distance(trial_point) <
4425  projected_point.distance(trial_point))
4426  projected_point = line_projections[0];
4427  }
4428  }
4429  else
4430  {
4431  Assert(false, ExcNotImplemented());
4432  return projected_point;
4433  }
4434 
4435  return projected_point;
4436  }
4437 
4438 
4439 
4440  namespace internal
4441  {
4442  template <typename DataType,
4443  typename MeshType,
4444  typename MeshCellIteratorType>
4445  inline void
4446  exchange_cell_data(
4447  const MeshType &mesh,
4448  const std::function<
4449  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4450  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4451  & unpack,
4452  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4453  const std::function<void(
4454  const std::function<void(const MeshCellIteratorType &,
4455  const types::subdomain_id)> &)> &process_cells,
4456  const std::function<std::set<types::subdomain_id>(
4457  const parallel::TriangulationBase<MeshType::dimension,
4458  MeshType::space_dimension> &)>
4459  &compute_ghost_owners)
4460  {
4461 # ifndef DEAL_II_WITH_MPI
4462  (void)mesh;
4463  (void)pack;
4464  (void)unpack;
4465  (void)cell_filter;
4466  (void)process_cells;
4467  (void)compute_ghost_owners;
4468  Assert(false, ExcNeedsMPI());
4469 # else
4470  constexpr int dim = MeshType::dimension;
4471  constexpr int spacedim = MeshType::space_dimension;
4472  auto tria =
4473  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4474  &mesh.get_triangulation());
4475  Assert(
4476  tria != nullptr,
4477  ExcMessage(
4478  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4479 
4480  if (const auto tria = dynamic_cast<
4482  &mesh.get_triangulation()))
4483  {
4484  Assert(
4485  tria->with_artificial_cells(),
4486  ExcMessage(
4487  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4488  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4489  "operate on a single layer of ghost cells. However, you have "
4490  "given a Triangulation object of type "
4491  "parallel::shared::Triangulation without artificial cells "
4492  "resulting in arbitrary numbers of ghost layers."));
4493  }
4494 
4495  // build list of cells to request for each neighbor
4496  std::set<::types::subdomain_id> ghost_owners =
4497  compute_ghost_owners(*tria);
4498  std::map<::types::subdomain_id,
4499  std::vector<typename CellId::binary_type>>
4500  neighbor_cell_list;
4501 
4502  for (const auto ghost_owner : ghost_owners)
4503  neighbor_cell_list[ghost_owner] = {};
4504 
4505  process_cells([&](const auto &cell, const auto key) -> void {
4506  if (cell_filter(cell))
4507  {
4508  constexpr int spacedim = MeshType::space_dimension;
4509  neighbor_cell_list[key].emplace_back(
4510  cell->id().template to_binary<spacedim>());
4511  }
4512  });
4513 
4514  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4515  ExcInternalError());
4516 
4517 
4518  // Before sending & receiving, make sure we protect this section with
4519  // a mutex:
4520  static Utilities::MPI::CollectiveMutex mutex;
4522  mutex, tria->get_communicator());
4523 
4524  const int mpi_tag =
4526  const int mpi_tag_reply =
4528 
4529  // send our requests
4530  std::vector<MPI_Request> requests(ghost_owners.size());
4531  {
4532  unsigned int idx = 0;
4533  for (const auto &it : neighbor_cell_list)
4534  {
4535  // send the data about the relevant cells
4536  const int ierr = MPI_Isend(it.second.data(),
4537  it.second.size() * sizeof(it.second[0]),
4538  MPI_BYTE,
4539  it.first,
4540  mpi_tag,
4542  &requests[idx]);
4543  AssertThrowMPI(ierr);
4544  ++idx;
4545  }
4546  }
4547 
4548  // receive requests and reply with the results
4549  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4550  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4551 
4552  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4553  {
4554  MPI_Status status;
4555  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4556  mpi_tag,
4558  &status);
4559  AssertThrowMPI(ierr);
4560 
4561  int len;
4562  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4563  AssertThrowMPI(ierr);
4564  Assert(len % sizeof(typename CellId::binary_type) == 0,
4565  ExcInternalError());
4566 
4567  const unsigned int n_cells =
4568  len / sizeof(typename CellId::binary_type);
4569  std::vector<typename CellId::binary_type> cells_with_requests(
4570  n_cells);
4571  std::vector<DataType> data_to_send;
4572  data_to_send.reserve(n_cells);
4573  std::vector<bool> cell_carries_data(n_cells, false);
4574 
4575  ierr = MPI_Recv(cells_with_requests.data(),
4576  len,
4577  MPI_BYTE,
4578  status.MPI_SOURCE,
4579  status.MPI_TAG,
4581  &status);
4582  AssertThrowMPI(ierr);
4583 
4584  // store data for each cell
4585  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4586  {
4587  const auto cell =
4588  tria->create_cell_iterator(CellId(cells_with_requests[c]));
4589 
4590  MeshCellIteratorType mesh_it(tria,
4591  cell->level(),
4592  cell->index(),
4593  &mesh);
4594 
4595  const std_cxx17::optional<DataType> data = pack(mesh_it);
4596  if (data)
4597  {
4598  data_to_send.emplace_back(*data);
4599  cell_carries_data[c] = true;
4600  }
4601  }
4602 
4603  // collect data for sending the reply in a buffer
4604 
4605  // (a) make room for storing the local offsets in case we receive
4606  // other data
4607  sendbuffers[idx].resize(sizeof(std::size_t));
4608 
4609  // (b) append the actual data and store how much memory it
4610  // corresponds to, which we then insert into the leading position of
4611  // the sendbuffer
4612  std::size_t size_of_send =
4613  Utilities::pack(data_to_send,
4614  sendbuffers[idx],
4615  /*enable_compression*/ false);
4616  std::memcpy(sendbuffers[idx].data(),
4617  &size_of_send,
4618  sizeof(std::size_t));
4619 
4620  // (c) append information of certain cells that got left out in case
4621  // we need it
4622  if (data_to_send.size() < n_cells)
4623  Utilities::pack(cell_carries_data,
4624  sendbuffers[idx],
4625  /*enable_compression*/ false);
4626 
4627  // send data
4628  ierr = MPI_Isend(sendbuffers[idx].data(),
4629  sendbuffers[idx].size(),
4630  MPI_BYTE,
4631  status.MPI_SOURCE,
4632  mpi_tag_reply,
4634  &reply_requests[idx]);
4635  AssertThrowMPI(ierr);
4636  }
4637 
4638  // finally receive the replies
4639  std::vector<char> receive;
4640  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4641  {
4642  MPI_Status status;
4643  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4644  mpi_tag_reply,
4646  &status);
4647  AssertThrowMPI(ierr);
4648 
4649  int len;
4650  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4651  AssertThrowMPI(ierr);
4652 
4653  receive.resize(len);
4654 
4655  ierr = MPI_Recv(receive.data(),
4656  len,
4657  MPI_BYTE,
4658  status.MPI_SOURCE,
4659  status.MPI_TAG,
4661  &status);
4662  AssertThrowMPI(ierr);
4663 
4664  // (a) first determine the length of the data section in the
4665  // received buffer
4666  auto data_iterator = receive.begin();
4667  std::size_t size_of_received_data =
4668  Utilities::unpack<std::size_t>(data_iterator,
4669  data_iterator + sizeof(std::size_t));
4670  data_iterator += sizeof(std::size_t);
4671 
4672  // (b) unpack the data section in the indicated region
4673  auto received_data = Utilities::unpack<std::vector<DataType>>(
4674  data_iterator,
4675  data_iterator + size_of_received_data,
4676  /*enable_compression*/ false);
4677  data_iterator += size_of_received_data;
4678 
4679  // (c) check if the received data contained fewer entries than the
4680  // number of cells we identified in the beginning, in which case we
4681  // need to extract the boolean vector with the relevant information
4682  const std::vector<typename CellId::binary_type> &this_cell_list =
4683  neighbor_cell_list[status.MPI_SOURCE];
4684  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4685  std::vector<bool> cells_with_data;
4686  if (received_data.size() < this_cell_list.size())
4687  {
4688  cells_with_data = Utilities::unpack<std::vector<bool>>(
4689  data_iterator, receive.end(), /*enable_compression*/ false);
4690  AssertDimension(cells_with_data.size(), this_cell_list.size());
4691  }
4692 
4693  // (d) go through the received data and call the user-provided
4694  // unpack function
4695  auto received_data_iterator = received_data.begin();
4696  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4697  if (cells_with_data.empty() || cells_with_data[c])
4698  {
4700  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4701 
4702  MeshCellIteratorType cell(tria,
4703  tria_cell->level(),
4704  tria_cell->index(),
4705  &mesh);
4706 
4707  unpack(cell, *received_data_iterator);
4708  ++received_data_iterator;
4709  }
4710  }
4711 
4712  // make sure that all communication is finished
4713  // when we leave this function.
4714  if (requests.size() > 0)
4715  {
4716  const int ierr =
4717  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4718  AssertThrowMPI(ierr);
4719  }
4720  if (reply_requests.size() > 0)
4721  {
4722  const int ierr = MPI_Waitall(reply_requests.size(),
4723  reply_requests.data(),
4724  MPI_STATUSES_IGNORE);
4725  AssertThrowMPI(ierr);
4726  }
4727 
4728 
4729 # endif // DEAL_II_WITH_MPI
4730  }
4731 
4732  } // namespace internal
4733 
4734  template <typename DataType, typename MeshType>
4735  inline void
4737  const MeshType & mesh,
4738  const std::function<std_cxx17::optional<DataType>(
4739  const typename MeshType::active_cell_iterator &)> &pack,
4740  const std::function<void(const typename MeshType::active_cell_iterator &,
4741  const DataType &)> & unpack,
4742  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4743  &cell_filter)
4744  {
4745 # ifndef DEAL_II_WITH_MPI
4746  (void)mesh;
4747  (void)pack;
4748  (void)unpack;
4749  (void)cell_filter;
4750  Assert(false, ExcNeedsMPI());
4751 # else
4752  internal::exchange_cell_data<DataType,
4753  MeshType,
4754  typename MeshType::active_cell_iterator>(
4755  mesh,
4756  pack,
4757  unpack,
4758  cell_filter,
4759  [&](const auto &process) {
4760  for (const auto &cell : mesh.active_cell_iterators())
4761  if (cell->is_ghost())
4762  process(cell, cell->subdomain_id());
4763  },
4764  [](const auto &tria) { return tria.ghost_owners(); });
4765 # endif
4766  }
4767 
4768 
4769 
4770  template <typename DataType, typename MeshType>
4771  inline void
4773  const MeshType & mesh,
4774  const std::function<std_cxx17::optional<DataType>(
4775  const typename MeshType::level_cell_iterator &)> &pack,
4776  const std::function<void(const typename MeshType::level_cell_iterator &,
4777  const DataType &)> & unpack,
4778  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4779  &cell_filter)
4780  {
4781 # ifndef DEAL_II_WITH_MPI
4782  (void)mesh;
4783  (void)pack;
4784  (void)unpack;
4785  (void)cell_filter;
4786  Assert(false, ExcNeedsMPI());
4787 # else
4788  internal::exchange_cell_data<DataType,
4789  MeshType,
4790  typename MeshType::level_cell_iterator>(
4791  mesh,
4792  pack,
4793  unpack,
4794  cell_filter,
4795  [&](const auto &process) {
4796  for (const auto &cell : mesh.cell_iterators())
4797  if (cell->is_ghost_on_level())
4798  process(cell, cell->level_subdomain_id());
4799  },
4800  [](const auto &tria) { return tria.level_ghost_owners(); });
4801 # endif
4802  }
4803 } // namespace GridTools
4804 
4805 #endif // DOXYGEN
4806 
4808 
4809 #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:3479
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:6912
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:6948
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:3565
const unsigned int n_subdivisions
Definition: grid_tools.h:3607
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:6860
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:6877
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:106
#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:5094
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:5186
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:5119
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:5214
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:6099
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3901
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:3263
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:4241
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:6579
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:5052
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:6375
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:3181
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:6356
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:4532
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:5396
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5367
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:4559
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:5741
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5334
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:4347
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:2858
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:6702
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:3410
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4516
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:5712
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:4000
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:3508
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:6420
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:3935
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:4446
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:6516
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:3964
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:2730
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:3557
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:4586
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5302
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:5909
@ 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
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
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:3416
void serialize(Archive &archive, const unsigned int version)
void save(Archive &ar, const unsigned int) const
Definition: grid_tools.h:3390
std::vector< CellId > cell_ids
Definition: grid_tools.h:3373
std::bitset< 3 > orientation
Definition: grid_tools.h:2794
std::size_t memory_consumption() const
FullMatrix< double > matrix
Definition: grid_tools.h:2808
unsigned int face_idx[2]
Definition: grid_tools.h:2787
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1221
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1196
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