Reference documentation for deal.II version Git 4e68a80cad 2021-10-22 15:50:12 -0600
\(\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 - 2021 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 #include <deal.II/fe/mapping_q1.h>
36 
37 #include <deal.II/grid/manifold.h>
38 #include <deal.II/grid/tria.h>
41 
42 #include <deal.II/hp/dof_handler.h>
43 
45 #include <deal.II/lac/la_vector.h>
49 
50 #include <deal.II/numerics/rtree.h>
51 
53 #include <boost/archive/binary_iarchive.hpp>
54 #include <boost/archive/binary_oarchive.hpp>
55 #include <boost/geometry/index/rtree.hpp>
56 #include <boost/random/mersenne_twister.hpp>
57 #include <boost/serialization/array.hpp>
58 #include <boost/serialization/vector.hpp>
59 
60 #ifdef DEAL_II_WITH_ZLIB
61 # include <boost/iostreams/device/back_inserter.hpp>
62 # include <boost/iostreams/filter/gzip.hpp>
63 # include <boost/iostreams/filtering_stream.hpp>
64 # include <boost/iostreams/stream.hpp>
65 #endif
67 
68 #include <bitset>
69 #include <list>
70 #include <set>
71 
73 
74 // Forward declarations
75 #ifndef DOXYGEN
76 namespace parallel
77 {
78  namespace distributed
79  {
80  template <int, int>
81  class Triangulation;
82  }
83 } // namespace parallel
84 
85 namespace hp
86 {
87  template <int, int>
88  class MappingCollection;
89 }
90 
91 class SparsityPattern;
92 #endif
93 
94 namespace internal
95 {
96  template <int dim, int spacedim, class MeshType>
98  {
99  public:
100 #ifndef _MSC_VER
101  using type = typename MeshType::active_cell_iterator;
102 #else
104 #endif
105  };
106 
107 #ifdef _MSC_VER
108  template <int dim, int spacedim>
109  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
110  {
111  public:
112  using type =
114  };
115 #endif
116 } // namespace internal
117 
126 namespace GridTools
127 {
128  template <int dim, int spacedim>
129  class Cache;
130 
135 
142  template <int dim, int spacedim>
143  double
145 
172  template <int dim, int spacedim>
173  double
175  const Mapping<dim, spacedim> & mapping =
176  (ReferenceCells::get_hypercube<dim>()
177  .template get_default_linear_mapping<dim, spacedim>()));
178 
189  template <int dim, int spacedim>
190  double
193  const Mapping<dim, spacedim> & mapping =
194  (ReferenceCells::get_hypercube<dim>()
195  .template get_default_linear_mapping<dim, spacedim>()));
196 
207  template <int dim, int spacedim>
208  double
210  const Triangulation<dim, spacedim> &triangulation,
211  const Mapping<dim, spacedim> & mapping =
212  (ReferenceCells::get_hypercube<dim>()
213  .template get_default_linear_mapping<dim, spacedim>()));
214 
226  template <int dim>
227  DEAL_II_DEPRECATED double
228  cell_measure(
229  const std::vector<Point<dim>> &all_vertices,
231 
248  template <int dim>
249  double
250  cell_measure(const std::vector<Point<dim>> & all_vertices,
252 
275  template <int dim, int spacedim>
276  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
278 
305  template <int dim>
306  Vector<double>
308  const Triangulation<dim> &triangulation,
309  const Quadrature<dim> & quadrature);
310 
318  template <int dim>
319  double
321  const Triangulation<dim> &triangulation,
322  const Quadrature<dim> & quadrature);
323 
337  template <int dim, int spacedim>
340 
358  template <typename Iterator>
361  const Iterator & object,
363 
376  template <int dim, int spacedim>
377  std::
378  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
380 
386 
403  template <int dim, int spacedim>
404  void
406  std::vector<CellData<dim>> & cells,
407  SubCellData & subcelldata);
408 
426  template <int dim, int spacedim>
427  void
428  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
429  std::vector<CellData<dim>> & cells,
430  SubCellData & subcelldata,
431  std::vector<unsigned int> & considered_vertices,
432  const double tol = 1e-12);
433 
452  template <int dim, int spacedim>
453  void
455  const std::vector<Point<spacedim>> &all_vertices,
456  std::vector<CellData<dim>> & cells);
457 
467  template <int dim>
468  void
469  consistently_order_cells(std::vector<CellData<dim>> &cells);
470 
476 
530  template <int dim, typename Transformation, int spacedim>
531  void
532  transform(const Transformation & transformation,
533  Triangulation<dim, spacedim> &triangulation);
534 
540  template <int dim, int spacedim>
541  void
542  shift(const Tensor<1, spacedim> & shift_vector,
543  Triangulation<dim, spacedim> &triangulation);
544 
545 
555  template <int dim>
556  void
557  rotate(const double angle, Triangulation<dim> &triangulation);
558 
570  template <int dim>
571  void
572  rotate(const Tensor<1, 3, double> &axis,
573  const double angle,
574  Triangulation<dim, 3> & triangulation);
575 
590  template <int dim>
591  DEAL_II_DEPRECATED_EARLY void
592  rotate(const double angle,
593  const unsigned int axis,
594  Triangulation<dim, 3> &triangulation);
595 
653  template <int dim>
654  void
655  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
656  Triangulation<dim> & tria,
657  const Function<dim, double> *coefficient = nullptr,
658  const bool solve_for_absolute_positions = false);
659 
665  template <int dim, int spacedim>
666  std::map<unsigned int, Point<spacedim>>
668 
676  template <int dim, int spacedim>
677  void
678  scale(const double scaling_factor,
679  Triangulation<dim, spacedim> &triangulation);
680 
695  template <int dim, int spacedim>
696  void
698  const double factor,
699  Triangulation<dim, spacedim> &triangulation,
700  const bool keep_boundary = true,
701  const unsigned int seed = boost::random::mt19937::default_seed);
702 
736  template <int dim, int spacedim>
737  void
739  const bool isotropic = false,
740  const unsigned int max_iterations = 100);
741 
766  template <int dim, int spacedim>
767  void
769  const double max_ratio = 1.6180339887,
770  const unsigned int max_iterations = 5);
771 
861  template <int dim, int spacedim>
862  void
864  const double limit_angle_fraction = .75);
865 
871 
926  template <int dim, int spacedim>
927 #ifndef DOXYGEN
928  std::tuple<
929  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
930  std::vector<std::vector<Point<dim>>>,
931  std::vector<std::vector<unsigned int>>>
932 #else
933  return_type
934 #endif
936  const Cache<dim, spacedim> & cache,
937  const std::vector<Point<spacedim>> &points,
939  &cell_hint =
941 
975  template <int dim, int spacedim>
976 #ifndef DOXYGEN
977  std::tuple<
978  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
979  std::vector<std::vector<Point<dim>>>,
980  std::vector<std::vector<unsigned int>>,
981  std::vector<unsigned int>>
982 #else
983  return_type
984 #endif
986  const Cache<dim, spacedim> & cache,
987  const std::vector<Point<spacedim>> &points,
989  &cell_hint =
991 
1061  template <int dim, int spacedim>
1062 #ifndef DOXYGEN
1063  std::tuple<
1064  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1065  std::vector<std::vector<Point<dim>>>,
1066  std::vector<std::vector<unsigned int>>,
1067  std::vector<std::vector<Point<spacedim>>>,
1068  std::vector<std::vector<unsigned int>>>
1069 #else
1070  return_type
1071 #endif
1073  const GridTools::Cache<dim, spacedim> & cache,
1074  const std::vector<Point<spacedim>> & local_points,
1075  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1076  const double tolerance = 1e-10);
1077 
1078  namespace internal
1079  {
1094  template <int dim, int spacedim>
1096  {
1103  std::vector<std::tuple<std::pair<int, int>,
1104  unsigned int,
1105  unsigned int,
1106  Point<dim>,
1108  unsigned int>>
1110 
1114  std::vector<unsigned int> send_ranks;
1115 
1121  std::vector<unsigned int> send_ptrs;
1122 
1133  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1135 
1139  std::vector<unsigned int> recv_ranks;
1140 
1146  std::vector<unsigned int> recv_ptrs;
1147  };
1148 
1158  template <int dim, int spacedim>
1161  const GridTools::Cache<dim, spacedim> & cache,
1162  const std::vector<Point<spacedim>> & points,
1163  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1164  const double tolerance,
1165  const bool perform_handshake,
1166  const bool enforce_unique_mapping = false);
1167 
1168  } // namespace internal
1169 
1206  template <int dim, int spacedim>
1207  std::map<unsigned int, Point<spacedim>>
1209  const Triangulation<dim, spacedim> &container,
1210  const Mapping<dim, spacedim> & mapping =
1211  (ReferenceCells::get_hypercube<dim>()
1212  .template get_default_linear_mapping<dim, spacedim>()));
1213 
1223  template <int spacedim>
1224  unsigned int
1225  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1226  const Point<spacedim> & p);
1227 
1251  template <int dim, template <int, int> class MeshType, int spacedim>
1252  unsigned int
1253  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1254  const Point<spacedim> & p,
1255  const std::vector<bool> & marked_vertices = {});
1256 
1280  template <int dim, template <int, int> class MeshType, int spacedim>
1281  unsigned int
1283  const MeshType<dim, spacedim> &mesh,
1284  const Point<spacedim> & p,
1285  const std::vector<bool> & marked_vertices = {});
1286 
1287 
1307  template <int dim, template <int, int> class MeshType, int spacedim>
1308 #ifndef _MSC_VER
1309  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1310 #else
1311  std::vector<
1312  typename ::internal::
1313  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1314 #endif
1315  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1316  const unsigned int vertex_index);
1317 
1380  template <int dim, template <int, int> class MeshType, int spacedim>
1381 #ifndef _MSC_VER
1382  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1383 #else
1384  std::pair<typename ::internal::
1385  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1386  Point<dim>>
1387 #endif
1389  const MeshType<dim, spacedim> &mesh,
1390  const Point<spacedim> & p,
1391  const std::vector<bool> &marked_vertices = {},
1392  const double tolerance = 1.e-10);
1393 
1401  template <int dim, template <int, int> class MeshType, int spacedim>
1402 #ifndef _MSC_VER
1403  typename MeshType<dim, spacedim>::active_cell_iterator
1404 #else
1405  typename ::internal::
1406  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1407 #endif
1408  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1409  const Point<spacedim> & p,
1410  const std::vector<bool> &marked_vertices = {},
1411  const double tolerance = 1.e-10);
1412 
1419  template <int dim, int spacedim>
1420  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1421  Point<dim>>
1423  const hp::MappingCollection<dim, spacedim> &mapping,
1424  const DoFHandler<dim, spacedim> & mesh,
1425  const Point<spacedim> & p,
1426  const double tolerance = 1.e-10);
1427 
1479  template <int dim, int spacedim>
1480  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1481  Point<dim>>
1483  const Cache<dim, spacedim> &cache,
1484  const Point<spacedim> & p,
1487  const std::vector<bool> &marked_vertices = {},
1488  const double tolerance = 1.e-10);
1489 
1503  template <int dim, template <int, int> class MeshType, int spacedim>
1504 #ifndef _MSC_VER
1505  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1506 #else
1507  std::pair<typename ::internal::
1508  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1509  Point<dim>>
1510 #endif
1512  const Mapping<dim, spacedim> & mapping,
1513  const MeshType<dim, spacedim> &mesh,
1514  const Point<spacedim> & p,
1515  const std::vector<
1516  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1518  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1519  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1520  typename MeshType<dim, spacedim>::active_cell_iterator(),
1521  const std::vector<bool> & marked_vertices = {},
1522  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1523  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1524  const double tolerance = 1.e-10,
1525  const RTree<
1526  std::pair<BoundingBox<spacedim>,
1528  *relevant_cell_bounding_boxes_rtree = nullptr);
1529 
1551  template <int dim, template <int, int> class MeshType, int spacedim>
1552 #ifndef _MSC_VER
1553  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1554  Point<dim>>>
1555 #else
1556  std::vector<std::pair<
1557  typename ::internal::
1558  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1559  Point<dim>>>
1560 #endif
1562  const Mapping<dim, spacedim> & mapping,
1563  const MeshType<dim, spacedim> &mesh,
1564  const Point<spacedim> & p,
1565  const double tolerance,
1566  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1567  Point<dim>> & first_cell);
1568 
1575  template <int dim, template <int, int> class MeshType, int spacedim>
1576 #ifndef _MSC_VER
1577  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1578  Point<dim>>>
1579 #else
1580  std::vector<std::pair<
1581  typename ::internal::
1582  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1583  Point<dim>>>
1584 #endif
1586  const Mapping<dim, spacedim> & mapping,
1587  const MeshType<dim, spacedim> &mesh,
1588  const Point<spacedim> & p,
1589  const double tolerance = 1e-10,
1590  const std::vector<bool> & marked_vertices = {});
1591 
1613  template <class MeshType>
1614  std::vector<typename MeshType::active_cell_iterator>
1615  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1616 
1641  template <class MeshType>
1642  void
1644  const typename MeshType::active_cell_iterator & cell,
1645  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1646 
1696  template <class MeshType>
1697  std::vector<typename MeshType::active_cell_iterator>
1699  const MeshType &mesh,
1700  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1701  &predicate);
1702 
1703 
1711  template <class MeshType>
1712  std::vector<typename MeshType::cell_iterator>
1714  const MeshType &mesh,
1715  const std::function<bool(const typename MeshType::cell_iterator &)>
1716  & predicate,
1717  const unsigned int level);
1718 
1719 
1732  template <class MeshType>
1733  std::vector<typename MeshType::active_cell_iterator>
1734  compute_ghost_cell_halo_layer(const MeshType &mesh);
1735 
1784  template <class MeshType>
1785  std::vector<typename MeshType::active_cell_iterator>
1787  const MeshType &mesh,
1788  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1789  & predicate,
1790  const double layer_thickness);
1791 
1814  template <class MeshType>
1815  std::vector<typename MeshType::active_cell_iterator>
1816  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1817  const double layer_thickness);
1818 
1834  template <class MeshType>
1835  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1837  const MeshType &mesh,
1838  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1839  &predicate);
1840 
1893  template <class MeshType>
1894  std::vector<BoundingBox<MeshType::space_dimension>>
1896  const MeshType &mesh,
1897  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1898  & predicate,
1899  const unsigned int refinement_level = 0,
1900  const bool allow_merge = false,
1901  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1902 
1930  template <int spacedim>
1931 #ifndef DOXYGEN
1932  std::tuple<std::vector<std::vector<unsigned int>>,
1933  std::map<unsigned int, unsigned int>,
1934  std::map<unsigned int, std::vector<unsigned int>>>
1935 #else
1936  return_type
1937 #endif
1939  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1940  const std::vector<Point<spacedim>> & points);
1941 
1942 
1977  template <int spacedim>
1978 #ifndef DOXYGEN
1979  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1980  std::map<unsigned int, unsigned int>,
1981  std::map<unsigned int, std::vector<unsigned int>>>
1982 #else
1983  return_type
1984 #endif
1986  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1987  const std::vector<Point<spacedim>> & points);
1988 
1989 
1998  template <int dim, int spacedim>
1999  std::vector<
2000  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2001  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
2002 
2015  template <int dim, int spacedim>
2016  std::vector<std::vector<Tensor<1, spacedim>>>
2018  const Triangulation<dim, spacedim> &mesh,
2019  const std::vector<
2021  &vertex_to_cells);
2022 
2023 
2031  template <int dim, int spacedim>
2032  unsigned int
2035  const Point<spacedim> & position,
2036  const Mapping<dim, spacedim> & mapping =
2037  (ReferenceCells::get_hypercube<dim>()
2038  .template get_default_linear_mapping<dim, spacedim>()));
2039 
2051  template <int dim, int spacedim>
2052  std::map<unsigned int, types::global_vertex_index>
2055 
2067  template <int dim, int spacedim>
2068  std::pair<unsigned int, double>
2071 
2077 
2086  template <int dim, int spacedim>
2087  void
2089  const Triangulation<dim, spacedim> &triangulation,
2090  DynamicSparsityPattern & connectivity);
2091 
2100  template <int dim, int spacedim>
2101  void
2103  const Triangulation<dim, spacedim> &triangulation,
2104  DynamicSparsityPattern & connectivity);
2105 
2114  template <int dim, int spacedim>
2115  void
2117  const Triangulation<dim, spacedim> &triangulation,
2118  const unsigned int level,
2119  DynamicSparsityPattern & connectivity);
2120 
2141  template <int dim, int spacedim>
2142  void
2143  partition_triangulation(const unsigned int n_partitions,
2144  Triangulation<dim, spacedim> & triangulation,
2145  const SparsityTools::Partitioner partitioner =
2147 
2158  template <int dim, int spacedim>
2159  void
2160  partition_triangulation(const unsigned int n_partitions,
2161  const std::vector<unsigned int> &cell_weights,
2162  Triangulation<dim, spacedim> & triangulation,
2163  const SparsityTools::Partitioner partitioner =
2165 
2211  template <int dim, int spacedim>
2212  void
2213  partition_triangulation(const unsigned int n_partitions,
2214  const SparsityPattern & cell_connection_graph,
2215  Triangulation<dim, spacedim> &triangulation,
2216  const SparsityTools::Partitioner partitioner =
2218 
2229  template <int dim, int spacedim>
2230  void
2231  partition_triangulation(const unsigned int n_partitions,
2232  const std::vector<unsigned int> &cell_weights,
2233  const SparsityPattern & cell_connection_graph,
2234  Triangulation<dim, spacedim> &triangulation,
2235  const SparsityTools::Partitioner partitioner =
2237 
2252  template <int dim, int spacedim>
2253  void
2254  partition_triangulation_zorder(const unsigned int n_partitions,
2255  Triangulation<dim, spacedim> &triangulation,
2256  const bool group_siblings = true);
2257 
2269  template <int dim, int spacedim>
2270  void
2272 
2280  template <int dim, int spacedim>
2281  std::vector<types::subdomain_id>
2283  const std::vector<CellId> & cell_ids);
2284 
2295  template <int dim, int spacedim>
2296  void
2298  std::vector<types::subdomain_id> & subdomain);
2299 
2314  template <int dim, int spacedim>
2315  unsigned int
2317  const Triangulation<dim, spacedim> &triangulation,
2318  const types::subdomain_id subdomain);
2319 
2349  template <int dim, int spacedim>
2350  std::vector<bool>
2352 
2358 
2392  template <typename MeshType>
2393  std::list<std::pair<typename MeshType::cell_iterator,
2394  typename MeshType::cell_iterator>>
2395  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2396 
2406  template <int dim, int spacedim>
2407  bool
2409  const Triangulation<dim, spacedim> &mesh_2);
2410 
2420  template <typename MeshType>
2421  bool
2422  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2423 
2429 
2445  template <int dim, int spacedim>
2449  & distorted_cells,
2450  Triangulation<dim, spacedim> &triangulation);
2451 
2452 
2453 
2462 
2463 
2501  template <class MeshType>
2502  std::vector<typename MeshType::active_cell_iterator>
2503  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2504 
2505 
2527  template <class Container>
2528  std::vector<typename Container::cell_iterator>
2530  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2531 
2598  template <class Container>
2599  void
2601  const std::vector<typename Container::active_cell_iterator> &patch,
2603  &local_triangulation,
2604  std::map<
2605  typename Triangulation<Container::dimension,
2606  Container::space_dimension>::active_cell_iterator,
2607  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2608 
2640  template <int dim, int spacedim>
2641  std::map<
2643  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2645 
2646 
2653 
2659  template <typename CellIterator>
2661  {
2665  CellIterator cell[2];
2666 
2671  unsigned int face_idx[2];
2672 
2678  std::bitset<3> orientation;
2679 
2693  };
2694 
2695 
2759  template <typename FaceIterator>
2760  bool
2762  std::bitset<3> & orientation,
2763  const FaceIterator & face1,
2764  const FaceIterator & face2,
2765  const int direction,
2769 
2770 
2774  template <typename FaceIterator>
2775  bool
2777  const FaceIterator & face1,
2778  const FaceIterator & face2,
2779  const int direction,
2783 
2784 
2841  template <typename MeshType>
2842  void
2844  const MeshType & mesh,
2845  const types::boundary_id b_id1,
2846  const types::boundary_id b_id2,
2847  const int direction,
2849  & matched_pairs,
2850  const Tensor<1, MeshType::space_dimension> &offset =
2853 
2854 
2877  template <typename MeshType>
2878  void
2880  const MeshType & mesh,
2881  const types::boundary_id b_id,
2882  const int direction,
2884  & matched_pairs,
2885  const ::Tensor<1, MeshType::space_dimension> &offset =
2888 
2894 
2915  template <int dim, int spacedim>
2916  void
2918  const bool reset_boundary_ids = false);
2919 
2941  template <int dim, int spacedim>
2942  void
2944  const std::vector<types::boundary_id> &src_boundary_ids,
2945  const std::vector<types::manifold_id> &dst_manifold_ids,
2947  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2948 
2978  template <int dim, int spacedim>
2979  void
2981  const bool compute_face_ids = false);
2982 
3007  template <int dim, int spacedim>
3008  void
3011  const std::function<types::manifold_id(
3012  const std::set<types::manifold_id> &)> &disambiguation_function =
3013  [](const std::set<types::manifold_id> &manifold_ids) {
3014  if (manifold_ids.size() == 1)
3015  return *manifold_ids.begin();
3016  else
3018  },
3019  bool overwrite_only_flat_manifold_ids = true);
3106  template <typename DataType, typename MeshType>
3107  void
3109  const MeshType & mesh,
3110  const std::function<std_cxx17::optional<DataType>(
3111  const typename MeshType::active_cell_iterator &)> &pack,
3112  const std::function<void(const typename MeshType::active_cell_iterator &,
3113  const DataType &)> & unpack,
3114  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3115  &cell_filter =
3117 
3128  template <typename DataType, typename MeshType>
3129  void
3131  const MeshType & mesh,
3132  const std::function<std_cxx17::optional<DataType>(
3133  const typename MeshType::level_cell_iterator &)> &pack,
3134  const std::function<void(const typename MeshType::level_cell_iterator &,
3135  const DataType &)> & unpack,
3136  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3138  true});
3139 
3140  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3141  * boxes @p local_bboxes.
3142  *
3143  * This function is meant to exchange bounding boxes describing the locally
3144  * owned cells in a distributed triangulation obtained with the function
3145  * GridTools::compute_mesh_predicate_bounding_box .
3146  *
3147  * The output vector's size is the number of processes of the MPI
3148  * communicator:
3149  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3150  */
3151  template <int spacedim>
3152  std::vector<std::vector<BoundingBox<spacedim>>>
3154  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3155  const MPI_Comm & mpi_communicator);
3156 
3189  template <int spacedim>
3192  const std::vector<BoundingBox<spacedim>> &local_description,
3193  const MPI_Comm & mpi_communicator);
3194 
3212  template <int dim, int spacedim>
3213  void
3215  const Triangulation<dim, spacedim> & tria,
3216  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3217  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3218 
3225  template <int dim, int spacedim>
3226  std::map<unsigned int, std::set<::types::subdomain_id>>
3228  const Triangulation<dim, spacedim> &tria);
3229 
3242  template <int dim, typename T>
3244  {
3248  std::vector<CellId> cell_ids;
3249 
3253  std::vector<T> data;
3254 
3263  template <class Archive>
3264  void
3265  save(Archive &ar, const unsigned int version) const;
3266 
3273  template <class Archive>
3274  void
3275  load(Archive &ar, const unsigned int version);
3276 
3277 #ifdef DOXYGEN
3278 
3283  template <class Archive>
3284  void
3285  serialize(Archive &archive, const unsigned int version);
3286 #else
3287  // This macro defines the serialize() method that is compatible with
3288  // the templated save() and load() method that have been implemented.
3289  BOOST_SERIALIZATION_SPLIT_MEMBER()
3290 #endif
3291  };
3292 
3293 
3294 
3312  template <int dim, typename VectorType>
3314  {
3315  public:
3319  using value_type = typename VectorType::value_type;
3320 
3324  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3325  const FiniteElement<dim, dim> &fe,
3326  const unsigned int n_subdivisions = 1,
3327  const double tolerance = 1e-10);
3328 
3333  void
3334  process(const DoFHandler<dim> & background_dof_handler,
3335  const VectorType & ls_vector,
3336  const double iso_level,
3337  std::vector<Point<dim>> & vertices,
3338  std::vector<CellData<dim - 1>> &cells) const;
3339 
3346  void
3347  process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
3348  const VectorType & ls_vector,
3349  const double iso_level,
3350  std::vector<Point<dim>> & vertices,
3351  std::vector<CellData<dim - 1>> &cells) const;
3352 
3353  private:
3358  static Quadrature<dim>
3359  create_quadrature_rule(const unsigned int n_subdivisions);
3360 
3364  void
3365  process_cell(std::vector<value_type> & ls_values,
3366  const std::vector<Point<dim>> & points,
3367  const double iso_level,
3368  std::vector<Point<dim>> & vertices,
3369  std::vector<CellData<dim - 1>> &cells) const;
3370 
3377  void
3378  process_sub_cell(const std::vector<value_type> & ls_values,
3379  const std::vector<Point<2>> & points,
3380  const std::vector<unsigned int> mask,
3381  const double iso_level,
3382  std::vector<Point<2>> & vertices,
3383  std::vector<CellData<1>> & cells) const;
3384 
3388  void
3389  process_sub_cell(const std::vector<value_type> & ls_values,
3390  const std::vector<Point<3>> & points,
3391  const std::vector<unsigned int> mask,
3392  const double iso_level,
3393  std::vector<Point<3>> & vertices,
3394  std::vector<CellData<2>> & cells) const;
3395 
3400  const unsigned int n_subdivisions;
3401 
3406  const double tolerance;
3407 
3413  };
3414 
3415 
3416 
3421 
3426  int,
3427  << "The number of partitions you gave is " << arg1
3428  << ", but must be greater than zero.");
3433  int,
3434  << "The subdomain id " << arg1
3435  << " has no cells associated with it.");
3440 
3445  double,
3446  << "The scaling factor must be positive, but it is " << arg1
3447  << ".");
3448 
3453  unsigned int,
3454  << "The given vertex with index " << arg1
3455  << " is not used in the given triangulation.");
3456 
3459 } /*namespace GridTools*/
3460 
3461 
3472  "The edges of the mesh are not consistently orientable.");
3473 
3474 
3475 
3476 /* ----------------- Template function --------------- */
3477 
3478 #ifndef DOXYGEN
3479 
3480 namespace GridTools
3481 {
3482  template <int dim>
3483  double
3484  cell_measure(
3485  const std::vector<Point<dim>> &all_vertices,
3486  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3487  {
3488  // We forward call to the ArrayView version:
3489  const ArrayView<const unsigned int> view(
3490  indices, GeometryInfo<dim>::vertices_per_cell);
3491  return cell_measure(all_vertices, view);
3492  }
3493 
3494 
3495 
3496  // This specialization is defined here so that the general template in the
3497  // source file doesn't need to have further 1D overloads for the internal
3498  // functions it calls.
3499  template <>
3503  {
3504  return {};
3505  }
3506 
3507 
3508 
3509  template <int dim, typename Predicate, int spacedim>
3510  void
3511  transform(const Predicate & predicate,
3512  Triangulation<dim, spacedim> &triangulation)
3513  {
3514  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3515 
3516  // loop over all active cells, and
3517  // transform those vertices that
3518  // have not yet been touched. note
3519  // that we get to all vertices in
3520  // the triangulation by only
3521  // visiting the active cells.
3523  cell = triangulation.begin_active(),
3524  endc = triangulation.end();
3525  for (; cell != endc; ++cell)
3526  for (const unsigned int v : cell->vertex_indices())
3527  if (treated_vertices[cell->vertex_index(v)] == false)
3528  {
3529  // transform this vertex
3530  cell->vertex(v) = predicate(cell->vertex(v));
3531  // and mark it as treated
3532  treated_vertices[cell->vertex_index(v)] = true;
3533  };
3534 
3535 
3536  // now fix any vertices on hanging nodes so that we don't create any holes
3537  if (dim == 2)
3538  {
3540  cell = triangulation.begin_active(),
3541  endc = triangulation.end();
3542  for (; cell != endc; ++cell)
3543  for (const unsigned int face : cell->face_indices())
3544  if (cell->face(face)->has_children() &&
3545  !cell->face(face)->at_boundary())
3546  {
3547  Assert(cell->reference_cell() ==
3548  ReferenceCells::get_hypercube<dim>(),
3549  ExcNotImplemented());
3550 
3551  // this line has children
3552  cell->face(face)->child(0)->vertex(1) =
3553  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3554  2;
3555  }
3556  }
3557  else if (dim == 3)
3558  {
3560  cell = triangulation.begin_active(),
3561  endc = triangulation.end();
3562  for (; cell != endc; ++cell)
3563  for (const unsigned int face : cell->face_indices())
3564  if (cell->face(face)->has_children() &&
3565  !cell->face(face)->at_boundary())
3566  {
3567  Assert(cell->reference_cell() ==
3568  ReferenceCells::get_hypercube<dim>(),
3569  ExcNotImplemented());
3570 
3571  // this face has hanging nodes
3572  cell->face(face)->child(0)->vertex(1) =
3573  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3574  2.0;
3575  cell->face(face)->child(0)->vertex(2) =
3576  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3577  2.0;
3578  cell->face(face)->child(1)->vertex(3) =
3579  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3580  2.0;
3581  cell->face(face)->child(2)->vertex(3) =
3582  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3583  2.0;
3584 
3585  // center of the face
3586  cell->face(face)->child(0)->vertex(3) =
3587  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3588  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3589  4.0;
3590  }
3591  }
3592 
3593  // Make sure FEValues notices that the mesh has changed
3594  triangulation.signals.mesh_movement();
3595  }
3596 
3597 
3598 
3599  template <class MeshType>
3600  std::vector<typename MeshType::active_cell_iterator>
3601  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3602  {
3603  std::vector<typename MeshType::active_cell_iterator> child_cells;
3604 
3605  if (cell->has_children())
3606  {
3607  for (unsigned int child = 0; child < cell->n_children(); ++child)
3608  if (cell->child(child)->has_children())
3609  {
3610  const std::vector<typename MeshType::active_cell_iterator>
3611  children = get_active_child_cells<MeshType>(cell->child(child));
3612  child_cells.insert(child_cells.end(),
3613  children.begin(),
3614  children.end());
3615  }
3616  else
3617  child_cells.push_back(cell->child(child));
3618  }
3619 
3620  return child_cells;
3621  }
3622 
3623 
3624 
3625  template <class MeshType>
3626  void
3628  const typename MeshType::active_cell_iterator & cell,
3629  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3630  {
3631  active_neighbors.clear();
3632  for (const unsigned int n : cell->face_indices())
3633  if (!cell->at_boundary(n))
3634  {
3635  if (MeshType::dimension == 1)
3636  {
3637  // check children of neighbor. note
3638  // that in 1d children of the neighbor
3639  // may be further refined. In 1d the
3640  // case is simple since we know what
3641  // children bound to the present cell
3642  typename MeshType::cell_iterator neighbor_child =
3643  cell->neighbor(n);
3644  if (!neighbor_child->is_active())
3645  {
3646  while (neighbor_child->has_children())
3647  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3648 
3649  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3650  ExcInternalError());
3651  }
3652  active_neighbors.push_back(neighbor_child);
3653  }
3654  else
3655  {
3656  if (cell->face(n)->has_children())
3657  // this neighbor has children. find
3658  // out which border to the present
3659  // cell
3660  for (unsigned int c = 0;
3661  c < cell->face(n)->n_active_descendants();
3662  ++c)
3663  active_neighbors.push_back(
3664  cell->neighbor_child_on_subface(n, c));
3665  else
3666  {
3667  // the neighbor must be active
3668  // himself
3669  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3670  active_neighbors.push_back(cell->neighbor(n));
3671  }
3672  }
3673  }
3674  }
3675 
3676 
3677 
3678  namespace internal
3679  {
3680  namespace ProjectToObject
3681  {
3694  struct CrossDerivative
3695  {
3696  const unsigned int direction_0;
3697  const unsigned int direction_1;
3698 
3699  CrossDerivative(const unsigned int d0, const unsigned int d1);
3700  };
3701 
3702  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3703  const unsigned int d1)
3704  : direction_0(d0)
3705  , direction_1(d1)
3706  {}
3707 
3708 
3709 
3714  template <typename F>
3715  inline auto
3716  centered_first_difference(const double center,
3717  const double step,
3718  const F &f) -> decltype(f(center) - f(center))
3719  {
3720  return (f(center + step) - f(center - step)) / (2.0 * step);
3721  }
3722 
3723 
3724 
3729  template <typename F>
3730  inline auto
3731  centered_second_difference(const double center,
3732  const double step,
3733  const F &f) -> decltype(f(center) - f(center))
3734  {
3735  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3736  (step * step);
3737  }
3738 
3739 
3740 
3750  template <int structdim, typename F>
3751  inline auto
3752  cross_stencil(
3753  const CrossDerivative cross_derivative,
3755  const double step,
3756  const F &f) -> decltype(f(center) - f(center))
3757  {
3759  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3760  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3761  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3762  1.0 / 3.0 * f(center - simplex_vector) +
3763  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3764  step;
3765  }
3766 
3767 
3768 
3775  template <int spacedim, int structdim, typename F>
3776  inline double
3777  gradient_entry(
3778  const unsigned int row_n,
3779  const unsigned int dependent_direction,
3780  const Point<spacedim> &p0,
3782  const double step,
3783  const F & f)
3784  {
3786  dependent_direction <
3788  ExcMessage("This function assumes that the last weight is a "
3789  "dependent variable (and hence we cannot take its "
3790  "derivative directly)."));
3791  Assert(row_n != dependent_direction,
3792  ExcMessage(
3793  "We cannot differentiate with respect to the variable "
3794  "that is assumed to be dependent."));
3795 
3796  const Point<spacedim> manifold_point = f(center);
3797  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3798  {row_n, dependent_direction}, center, step, f);
3799  double entry = 0.0;
3800  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3801  entry +=
3802  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3803  return entry;
3804  }
3805 
3811  template <typename Iterator, int spacedim, int structdim>
3813  project_to_d_linear_object(const Iterator & object,
3814  const Point<spacedim> &trial_point)
3815  {
3816  // let's look at this for simplicity for a quad (structdim==2) in a
3817  // space with spacedim>2 (notate trial_point by y): all points on the
3818  // surface are given by
3819  // x(\xi) = sum_i v_i phi_x(\xi)
3820  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3821  // reference coordinates of the quad. so what we are trying to do is
3822  // find a point x on the surface that is closest to the point y. there
3823  // are different ways to solve this problem, but in the end it's a
3824  // nonlinear problem and we have to find reference coordinates \xi so
3825  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3826  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3827  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3828  // have to use a Newton method to find the answer. This leads to the
3829  // following formulation of Newton steps:
3830  //
3831  // Given \xi_k, find \delta\xi_k so that
3832  // H_k \delta\xi_k = - F_k
3833  // where H_k is an approximation to the second derivatives of J at
3834  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3835  // number of times until the right hand side is small enough. As a
3836  // stopping criterion, we terminate if ||\delta\xi||<eps.
3837  //
3838  // As for the Hessian, the best choice would be
3839  // H_k = J''(\xi_k)
3840  // but we'll opt for the simpler Gauss-Newton form
3841  // H_k = A^T A
3842  // i.e.
3843  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3844  // \partial_n phi_i *
3845  // \partial_m phi_j
3846  // we start at xi=(0.5, 0.5).
3847  Point<structdim> xi;
3848  for (unsigned int d = 0; d < structdim; ++d)
3849  xi[d] = 0.5;
3850 
3851  Point<spacedim> x_k;
3852  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3853  x_k += object->vertex(i) *
3855 
3856  do
3857  {
3859  for (const unsigned int i :
3861  F_k +=
3862  (x_k - trial_point) * object->vertex(i) *
3864  i);
3865 
3867  for (const unsigned int i :
3869  for (const unsigned int j :
3871  {
3874  xi, i),
3876  xi, j));
3877  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3878  }
3879 
3880  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3881  xi += delta_xi;
3882 
3883  x_k = Point<spacedim>();
3884  for (const unsigned int i :
3886  x_k += object->vertex(i) *
3888 
3889  if (delta_xi.norm() < 1e-7)
3890  break;
3891  }
3892  while (true);
3893 
3894  return x_k;
3895  }
3896  } // namespace ProjectToObject
3897  } // namespace internal
3898 
3899 
3900 
3901  namespace internal
3902  {
3903  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3904  // inside the project_to_object function below.
3905  template <int structdim>
3906  inline bool
3907  weights_are_ok(
3909  {
3910  // clang has trouble figuring out structdim here, so define it
3911  // again:
3912  static const std::size_t n_vertices_per_cell =
3914  n_independent_components;
3915  std::array<double, n_vertices_per_cell> copied_weights;
3916  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3917  {
3918  copied_weights[i] = v[i];
3919  if (v[i] < 0.0 || v[i] > 1.0)
3920  return false;
3921  }
3922 
3923  // check the sum: try to avoid some roundoff errors by summing in order
3924  std::sort(copied_weights.begin(), copied_weights.end());
3925  const double sum =
3926  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3927  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3928  }
3929  } // namespace internal
3930 
3931  template <typename Iterator>
3934  const Iterator & object,
3936  {
3937  const int spacedim = Iterator::AccessorType::space_dimension;
3938  const int structdim = Iterator::AccessorType::structure_dimension;
3939 
3940  Point<spacedim> projected_point = trial_point;
3941 
3942  if (structdim >= spacedim)
3943  return projected_point;
3944  else if (structdim == 1 || structdim == 2)
3945  {
3946  using namespace internal::ProjectToObject;
3947  // Try to use the special flat algorithm for quads (this is better
3948  // than the general algorithm in 3D). This does not take into account
3949  // whether projected_point is outside the quad, but we optimize along
3950  // lines below anyway:
3951  const int dim = Iterator::AccessorType::dimension;
3952  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3953  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3954  &manifold) != nullptr)
3955  {
3956  projected_point =
3957  project_to_d_linear_object<Iterator, spacedim, structdim>(
3958  object, trial_point);
3959  }
3960  else
3961  {
3962  // We want to find a point on the convex hull (defined by the
3963  // vertices of the object and the manifold description) that is
3964  // relatively close to the trial point. This has a few issues:
3965  //
3966  // 1. For a general convex hull we are not guaranteed that a unique
3967  // minimum exists.
3968  // 2. The independent variables in the optimization process are the
3969  // weights given to Manifold::get_new_point, which must sum to 1,
3970  // so we cannot use standard finite differences to approximate a
3971  // gradient.
3972  //
3973  // There is not much we can do about 1., but for 2. we can derive
3974  // finite difference stencils that work on a structdim-dimensional
3975  // simplex and rewrite the optimization problem to use those
3976  // instead. Consider the structdim 2 case and let
3977  //
3978  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3979  // c2, c3})
3980  //
3981  // where {c0, c1, c2, c3} are the weights for the four vertices on
3982  // the quadrilateral. We seek to minimize the Euclidean distance
3983  // between F(...) and trial_point. We can solve for c3 in terms of
3984  // the other weights and get, for one coordinate direction
3985  //
3986  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3987  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3988  //
3989  // where we substitute back in for c3 after taking the
3990  // derivative. We can compute a stencil for the cross derivative
3991  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3992  // (and gradient_entry computes the sum over the independent
3993  // variables). Below, we somewhat arbitrarily pick the last
3994  // component as the dependent one.
3995  //
3996  // Since we can now calculate derivatives of the objective
3997  // function we can use gradient descent to minimize it.
3998  //
3999  // Of course, this is much simpler in the structdim = 1 case (we
4000  // could rewrite the projection as a 1D optimization problem), but
4001  // to reduce the potential for bugs we use the same code in both
4002  // cases.
4003  const double step_size = object->diameter() / 64.0;
4004 
4005  constexpr unsigned int n_vertices_per_cell =
4007 
4008  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4009  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4010  ++vertex_n)
4011  vertices[vertex_n] = object->vertex(vertex_n);
4012 
4013  auto get_point_from_weights =
4014  [&](const Tensor<1, n_vertices_per_cell> &weights)
4015  -> Point<spacedim> {
4016  return object->get_manifold().get_new_point(
4017  make_array_view(vertices.begin(), vertices.end()),
4018  make_array_view(weights.begin_raw(), weights.end_raw()));
4019  };
4020 
4021  // pick the initial weights as (normalized) inverse distances from
4022  // the trial point:
4023  Tensor<1, n_vertices_per_cell> guess_weights;
4024  double guess_weights_sum = 0.0;
4025  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4026  ++vertex_n)
4027  {
4028  const double distance =
4029  vertices[vertex_n].distance(trial_point);
4030  if (distance == 0.0)
4031  {
4032  guess_weights = 0.0;
4033  guess_weights[vertex_n] = 1.0;
4034  guess_weights_sum = 1.0;
4035  break;
4036  }
4037  else
4038  {
4039  guess_weights[vertex_n] = 1.0 / distance;
4040  guess_weights_sum += guess_weights[vertex_n];
4041  }
4042  }
4043  guess_weights /= guess_weights_sum;
4044  Assert(internal::weights_are_ok<structdim>(guess_weights),
4045  ExcInternalError());
4046 
4047  // The optimization algorithm consists of two parts:
4048  //
4049  // 1. An outer loop where we apply the gradient descent algorithm.
4050  // 2. An inner loop where we do a line search to find the optimal
4051  // length of the step one should take in the gradient direction.
4052  //
4053  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4054  {
4055  const unsigned int dependent_direction =
4056  n_vertices_per_cell - 1;
4057  Tensor<1, n_vertices_per_cell> current_gradient;
4058  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4059  ++row_n)
4060  {
4061  if (row_n != dependent_direction)
4062  {
4063  current_gradient[row_n] =
4064  gradient_entry<spacedim, structdim>(
4065  row_n,
4066  dependent_direction,
4067  trial_point,
4068  guess_weights,
4069  step_size,
4070  get_point_from_weights);
4071 
4072  current_gradient[dependent_direction] -=
4073  current_gradient[row_n];
4074  }
4075  }
4076 
4077  // We need to travel in the -gradient direction, as noted
4078  // above, but we may not want to take a full step in that
4079  // direction; instead, guess that we will go -0.5*gradient and
4080  // do quasi-Newton iteration to pick the best multiplier. The
4081  // goal is to find a scalar alpha such that
4082  //
4083  // F(x - alpha g)
4084  //
4085  // is minimized, where g is the gradient and F is the
4086  // objective function. To find the optimal value we find roots
4087  // of the derivative of the objective function with respect to
4088  // alpha by Newton iteration, where we approximate the first
4089  // and second derivatives of F(x - alpha g) with centered
4090  // finite differences.
4091  double gradient_weight = -0.5;
4092  auto gradient_weight_objective_function =
4093  [&](const double gradient_weight_guess) -> double {
4094  return (trial_point -
4095  get_point_from_weights(guess_weights +
4096  gradient_weight_guess *
4097  current_gradient))
4098  .norm_square();
4099  };
4100 
4101  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4102  {
4103  const double update_numerator = centered_first_difference(
4104  gradient_weight,
4105  step_size,
4106  gradient_weight_objective_function);
4107  const double update_denominator =
4108  centered_second_difference(
4109  gradient_weight,
4110  step_size,
4111  gradient_weight_objective_function);
4112 
4113  // avoid division by zero. Note that we limit the gradient
4114  // weight below
4115  if (std::abs(update_denominator) == 0.0)
4116  break;
4117  gradient_weight =
4118  gradient_weight - update_numerator / update_denominator;
4119 
4120  // Put a fairly lenient bound on the largest possible
4121  // gradient (things tend to be locally flat, so the gradient
4122  // itself is usually small)
4123  if (std::abs(gradient_weight) > 10)
4124  {
4125  gradient_weight = -10.0;
4126  break;
4127  }
4128  }
4129 
4130  // It only makes sense to take convex combinations with weights
4131  // between zero and one. If the update takes us outside of this
4132  // region then rescale the update to stay within the region and
4133  // try again
4134  Tensor<1, n_vertices_per_cell> tentative_weights =
4135  guess_weights + gradient_weight * current_gradient;
4136 
4137  double new_gradient_weight = gradient_weight;
4138  for (unsigned int iteration_count = 0; iteration_count < 40;
4139  ++iteration_count)
4140  {
4141  if (internal::weights_are_ok<structdim>(tentative_weights))
4142  break;
4143 
4144  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4145  {
4146  if (tentative_weights[i] < 0.0)
4147  {
4148  tentative_weights -=
4149  (tentative_weights[i] / current_gradient[i]) *
4150  current_gradient;
4151  }
4152  if (tentative_weights[i] < 0.0 ||
4153  1.0 < tentative_weights[i])
4154  {
4155  new_gradient_weight /= 2.0;
4156  tentative_weights =
4157  guess_weights +
4158  new_gradient_weight * current_gradient;
4159  }
4160  }
4161  }
4162 
4163  // the update might still send us outside the valid region, so
4164  // check again and quit if the update is still not valid
4165  if (!internal::weights_are_ok<structdim>(tentative_weights))
4166  break;
4167 
4168  // if we cannot get closer by traveling in the gradient
4169  // direction then quit
4170  if (get_point_from_weights(tentative_weights)
4171  .distance(trial_point) <
4172  get_point_from_weights(guess_weights).distance(trial_point))
4173  guess_weights = tentative_weights;
4174  else
4175  break;
4176  Assert(internal::weights_are_ok<structdim>(guess_weights),
4177  ExcInternalError());
4178  }
4179  Assert(internal::weights_are_ok<structdim>(guess_weights),
4180  ExcInternalError());
4181  projected_point = get_point_from_weights(guess_weights);
4182  }
4183 
4184  // if structdim == 2 and the optimal point is not on the interior then
4185  // we may be able to get a more accurate result by projecting onto the
4186  // lines.
4187  if (structdim == 2)
4188  {
4189  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4190  line_projections;
4191  for (unsigned int line_n = 0;
4192  line_n < GeometryInfo<structdim>::lines_per_cell;
4193  ++line_n)
4194  {
4195  line_projections[line_n] =
4196  project_to_object(object->line(line_n), trial_point);
4197  }
4198  std::sort(line_projections.begin(),
4199  line_projections.end(),
4200  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4201  return a.distance(trial_point) <
4202  b.distance(trial_point);
4203  });
4204  if (line_projections[0].distance(trial_point) <
4205  projected_point.distance(trial_point))
4206  projected_point = line_projections[0];
4207  }
4208  }
4209  else
4210  {
4211  Assert(false, ExcNotImplemented());
4212  return projected_point;
4213  }
4214 
4215  return projected_point;
4216  }
4217 
4218 
4219 
4220  template <int dim, typename T>
4221  template <class Archive>
4222  void
4224  const unsigned int /*version*/) const
4225  {
4226  Assert(cell_ids.size() == data.size(),
4227  ExcDimensionMismatch(cell_ids.size(), data.size()));
4228  // archive the cellids in an efficient binary format
4229  const std::size_t n_cells = cell_ids.size();
4230  ar & n_cells;
4231  for (const auto &id : cell_ids)
4232  {
4233  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
4234  ar & binary_cell_id;
4235  }
4236 
4237  ar &data;
4238  }
4239 
4240 
4241 
4242  template <int dim, typename T>
4243  template <class Archive>
4244  void
4246  const unsigned int /*version*/)
4247  {
4248  std::size_t n_cells;
4249  ar & n_cells;
4250  cell_ids.clear();
4251  cell_ids.reserve(n_cells);
4252  for (unsigned int c = 0; c < n_cells; ++c)
4253  {
4254  CellId::binary_type value;
4255  ar & value;
4256  cell_ids.emplace_back(value);
4257  }
4258  ar &data;
4259  }
4260 
4261 
4262  namespace internal
4263  {
4264  template <typename DataType,
4265  typename MeshType,
4266  typename MeshCellIteratorType>
4267  inline void
4268  exchange_cell_data(
4269  const MeshType &mesh,
4270  const std::function<
4271  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4272  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4273  & unpack,
4274  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4275  const std::function<void(
4276  const std::function<void(const MeshCellIteratorType &,
4277  const types::subdomain_id)> &)> &process_cells,
4278  const std::function<std::set<types::subdomain_id>(
4279  const parallel::TriangulationBase<MeshType::dimension,
4280  MeshType::space_dimension> &)>
4281  &compute_ghost_owners)
4282  {
4283 # ifndef DEAL_II_WITH_MPI
4284  (void)mesh;
4285  (void)pack;
4286  (void)unpack;
4287  (void)cell_filter;
4288  (void)process_cells;
4289  (void)compute_ghost_owners;
4290  Assert(false, ExcNeedsMPI());
4291 # else
4292  constexpr int dim = MeshType::dimension;
4293  constexpr int spacedim = MeshType::space_dimension;
4294  auto tria =
4295  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4296  &mesh.get_triangulation());
4297  Assert(
4298  tria != nullptr,
4299  ExcMessage(
4300  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4301 
4302  if (const auto tria = dynamic_cast<
4304  &mesh.get_triangulation()))
4305  {
4306  Assert(
4307  tria->with_artificial_cells(),
4308  ExcMessage(
4309  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4310  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4311  "operate on a single layer ghost cells. However, you have "
4312  "given a Triangulation object of type "
4313  "parallel::shared::Triangulation without artificial cells "
4314  "resulting in arbitrary numbers of ghost layers."));
4315  }
4316 
4317  // build list of cells to request for each neighbor
4318  std::set<::types::subdomain_id> ghost_owners =
4319  compute_ghost_owners(*tria);
4320  std::map<::types::subdomain_id,
4321  std::vector<typename CellId::binary_type>>
4322  neighbor_cell_list;
4323 
4324  for (const auto ghost_owner : ghost_owners)
4325  neighbor_cell_list[ghost_owner] = {};
4326 
4327  process_cells([&](const auto &cell, const auto key) {
4328  if (cell_filter(cell))
4329  neighbor_cell_list[key].emplace_back(
4330  cell->id().template to_binary<spacedim>());
4331  });
4332 
4333  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4334  ExcInternalError());
4335 
4336 
4337  // Before sending & receiving, make sure we protect this section with
4338  // a mutex:
4339  static Utilities::MPI::CollectiveMutex mutex;
4341  mutex, tria->get_communicator());
4342 
4343  const int mpi_tag =
4345  const int mpi_tag_reply =
4347 
4348  // send our requests:
4349  std::vector<MPI_Request> requests(ghost_owners.size());
4350  {
4351  unsigned int idx = 0;
4352  for (const auto &it : neighbor_cell_list)
4353  {
4354  // send the data about the relevant cells
4355  const int ierr = MPI_Isend(it.second.data(),
4356  it.second.size() * sizeof(it.second[0]),
4357  MPI_BYTE,
4358  it.first,
4359  mpi_tag,
4360  tria->get_communicator(),
4361  &requests[idx]);
4362  AssertThrowMPI(ierr);
4363  ++idx;
4364  }
4365  }
4366 
4367  using DestinationToBufferMap =
4368  std::map<::types::subdomain_id,
4370  DestinationToBufferMap destination_to_data_buffer_map;
4371 
4372  // receive requests and reply with the ghost indices
4373  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4374  ghost_owners.size());
4375  std::vector<std::vector<types::global_dof_index>>
4376  send_dof_numbers_and_indices(ghost_owners.size());
4377  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4378  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4379 
4380  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4381  {
4382  MPI_Status status;
4383  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4384  mpi_tag,
4385  tria->get_communicator(),
4386  &status);
4387  AssertThrowMPI(ierr);
4388 
4389  int len;
4390  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4391  AssertThrowMPI(ierr);
4392  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4393  ExcInternalError());
4394 
4395  const unsigned int n_cells =
4396  len / sizeof(typename CellId::binary_type);
4397  cell_data_to_send[idx].resize(n_cells);
4398 
4399  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4400  len,
4401  MPI_BYTE,
4402  status.MPI_SOURCE,
4403  status.MPI_TAG,
4404  tria->get_communicator(),
4405  &status);
4406  AssertThrowMPI(ierr);
4407 
4408  // store data for each cell
4409  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4410  {
4411  const auto cell =
4412  tria->create_cell_iterator(CellId(cell_data_to_send[idx][c]));
4413 
4414  MeshCellIteratorType mesh_it(tria,
4415  cell->level(),
4416  cell->index(),
4417  &mesh);
4418  const std_cxx17::optional<DataType> data = pack(mesh_it);
4419 
4420  if (data)
4421  {
4422  typename DestinationToBufferMap::iterator p =
4423  destination_to_data_buffer_map
4424  .insert(std::make_pair(
4425  idx,
4426  GridTools::CellDataTransferBuffer<dim, DataType>()))
4427  .first;
4428 
4429  p->second.cell_ids.emplace_back(cell->id());
4430  p->second.data.emplace_back(*data);
4431  }
4432  }
4433 
4434  // send reply
4435  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4436  destination_to_data_buffer_map[idx];
4437 
4438  sendbuffers[idx] =
4439  Utilities::pack(data, /*enable_compression*/ false);
4440  ierr = MPI_Isend(sendbuffers[idx].data(),
4441  sendbuffers[idx].size(),
4442  MPI_BYTE,
4443  status.MPI_SOURCE,
4444  mpi_tag_reply,
4445  tria->get_communicator(),
4446  &reply_requests[idx]);
4447  AssertThrowMPI(ierr);
4448  }
4449 
4450  // finally receive the replies
4451  std::vector<char> receive;
4452  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4453  {
4454  MPI_Status status;
4455  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4456  mpi_tag_reply,
4457  tria->get_communicator(),
4458  &status);
4459  AssertThrowMPI(ierr);
4460 
4461  int len;
4462  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4463  AssertThrowMPI(ierr);
4464 
4465  receive.resize(len);
4466 
4467  char *ptr = receive.data();
4468  ierr = MPI_Recv(ptr,
4469  len,
4470  MPI_BYTE,
4471  status.MPI_SOURCE,
4472  status.MPI_TAG,
4473  tria->get_communicator(),
4474  &status);
4475  AssertThrowMPI(ierr);
4476 
4477  auto cellinfo =
4478  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4479  receive, /*enable_compression*/ false);
4480 
4481  DataType *data = cellinfo.data.data();
4482  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4483  {
4485  tria_cell = tria->create_cell_iterator(cellinfo.cell_ids[c]);
4486 
4487  MeshCellIteratorType cell(tria,
4488  tria_cell->level(),
4489  tria_cell->index(),
4490  &mesh);
4491 
4492  unpack(cell, *data);
4493  }
4494  }
4495 
4496  // make sure that all communication is finished
4497  // when we leave this function.
4498  if (requests.size() > 0)
4499  {
4500  const int ierr =
4501  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4502  AssertThrowMPI(ierr);
4503  }
4504  if (reply_requests.size() > 0)
4505  {
4506  const int ierr = MPI_Waitall(reply_requests.size(),
4507  reply_requests.data(),
4508  MPI_STATUSES_IGNORE);
4509  AssertThrowMPI(ierr);
4510  }
4511 
4512 
4513 # endif // DEAL_II_WITH_MPI
4514  }
4515 
4516  } // namespace internal
4517 
4518  template <typename DataType, typename MeshType>
4519  inline void
4521  const MeshType & mesh,
4522  const std::function<std_cxx17::optional<DataType>(
4523  const typename MeshType::active_cell_iterator &)> &pack,
4524  const std::function<void(const typename MeshType::active_cell_iterator &,
4525  const DataType &)> & unpack,
4526  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4527  &cell_filter)
4528  {
4529 # ifndef DEAL_II_WITH_MPI
4530  (void)mesh;
4531  (void)pack;
4532  (void)unpack;
4533  (void)cell_filter;
4534  Assert(false, ExcNeedsMPI());
4535 # else
4536  internal::exchange_cell_data<DataType,
4537  MeshType,
4538  typename MeshType::active_cell_iterator>(
4539  mesh,
4540  pack,
4541  unpack,
4542  cell_filter,
4543  [&](const auto &process) {
4544  for (const auto &cell : mesh.active_cell_iterators())
4545  if (cell->is_ghost())
4546  process(cell, cell->subdomain_id());
4547  },
4548  [](const auto &tria) { return tria.ghost_owners(); });
4549 # endif
4550  }
4551 
4552 
4553 
4554  template <typename DataType, typename MeshType>
4555  inline void
4557  const MeshType & mesh,
4558  const std::function<std_cxx17::optional<DataType>(
4559  const typename MeshType::level_cell_iterator &)> &pack,
4560  const std::function<void(const typename MeshType::level_cell_iterator &,
4561  const DataType &)> & unpack,
4562  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4563  &cell_filter)
4564  {
4565 # ifndef DEAL_II_WITH_MPI
4566  (void)mesh;
4567  (void)pack;
4568  (void)unpack;
4569  (void)cell_filter;
4570  Assert(false, ExcNeedsMPI());
4571 # else
4572  internal::exchange_cell_data<DataType,
4573  MeshType,
4574  typename MeshType::level_cell_iterator>(
4575  mesh,
4576  pack,
4577  unpack,
4578  cell_filter,
4579  [&](const auto &process) {
4580  for (const auto &cell : mesh.cell_iterators())
4581  if (cell->level_subdomain_id() !=
4583  !cell->is_locally_owned_on_level())
4584  process(cell, cell->level_subdomain_id());
4585  },
4586  [](const auto &tria) { return tria.level_ghost_owners(); });
4587 # endif
4588  }
4589 } // namespace GridTools
4590 
4591 #endif // DOXYGEN
4592 
4594 
4595 #endif
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5116
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:4901
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:6279
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:213
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
unsigned int n_vertices() const
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static ::ExceptionBase & ExcNeedsMPI()
const unsigned int n_subdivisions
Definition: grid_tools.h:3400
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2190
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:4876
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1922
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:82
virtual MPI_Comm get_communicator() const
Definition: tria.cc:11233
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:312
Contents is actually a matrix.
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5084
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:2966
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void load(Archive &ar, const unsigned int version)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
return_type guess_point_owner(const RTree< std::pair< BoundingBox< spacedim >, unsigned int >> &covering_rtree, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3243
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:5178
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2047
const ::Triangulation< dim, spacedim > & tria
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3293
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:13067
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:3342
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
std::pair< Point< MeshType::space_dimension >, Point< MeshType::space_dimension > > compute_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
std::bitset< 3 > orientation
Definition: grid_tools.h:2678
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:6119
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4126
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:150
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:414
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:6183
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:626
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3729
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
void partition_triangulation(const unsigned int n_partitions, const std::vector< unsigned int > &cell_weights, const SparsityPattern &cell_connection_graph, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:3940
cell_iterator end() const
Definition: tria.cc:13158
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:534
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:4341
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
void process_sub_cell(const std::array< unsigned int, n_configurations > &cut_line_table, const ndarray< unsigned int, n_configurations, n_cols > &new_line_table, const ndarray< unsigned int, n_lines, 2 > &line_to_vertex_table, const std::vector< value_type > &ls_values, const std::vector< Point< dim >> &points, const std::vector< unsigned int > &mask, const double iso_level, const double tolerance, std::vector< Point< dim >> &vertices, std::vector< CellData< dim - 1 >> &cells)
Definition: grid_tools.cc:6663
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1109
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:101
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< std::set< typename MeshType< dim, spacedim >::active_cell_iterator >> &vertex_to_cell_map, const std::vector< std::vector< Tensor< 1, spacedim >>> &vertex_to_cell_centers, const typename MeshType< dim, spacedim >::active_cell_iterator &cell_hint=typename MeshType< dim, spacedim >::active_cell_iterator(), const std::vector< bool > &marked_vertices={}, const RTree< std::pair< Point< spacedim >, unsigned int >> &used_vertices_rtree=RTree< std::pair< Point< spacedim >, unsigned int >>{}, const double tolerance=1.e-10, const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator >> *relevant_cell_bounding_boxes_rtree=nullptr)
Definition: grid_tools.cc:2764
std::map< types::global_dof_index, std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandler< dim, spacedim > &dof_handler)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
T sum(const T &t, const MPI_Comm &mpi_communicator)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3758
typename VectorType::value_type value_type
Definition: grid_tools.h:3319
bool orthogonal_equality(const FaceIterator &face1, const FaceIterator &face2, const int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
#define Assert(cond, exc)
Definition: exceptions.h:1461
Signals signals
Definition: tria.h:2289
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3248
Abstract base class for mapping classes.
Definition: mapping.h:303
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:4368
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1256
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:698
#define DeclException0(Exception0)
Definition: exceptions.h:464
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:3048
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 >> &cells) const
Definition: grid_tools.cc:6570
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
void save(Archive &ar, const unsigned int version) const
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5149
Point< 3 > vertices[4]
void build_triangulation_from_patch(const std::vector< typename Container::active_cell_iterator > &patch, Triangulation< Container::dimension, Container::space_dimension > &local_triangulation, std::map< typename Triangulation< Container::dimension, Container::space_dimension >::active_cell_iterator, typename Container::active_cell_iterator > &patch_to_global_tria_map)
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:6342
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:288
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:4968
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:4834
cell_iterator create_cell_iterator(const CellId &cell_id) const
Definition: tria.cc:13133
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1219
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: hp.h:117
std::vector< std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_active_cells_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const double tolerance=1e-10, const std::vector< bool > &marked_vertices={})
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13697
Definition: cell_id.h:70
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator >> &vertex_to_cells)
Definition: grid_tools.cc:2701
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
unsigned int global_dof_index
Definition: types.h:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
__global__ void set(Number *val, const Number s, const size_type N)
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:5522
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1778
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
double cell_measure(const std::vector< Point< dim >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:452
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
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:2222
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:138
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:2573
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:378
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:4284
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3695
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition: grid_tools.cc:5880
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1322
Point< 3 > center
FullMatrix< double > matrix
Definition: grid_tools.h:2692
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:4996
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:2032
static ::ExceptionBase & ExcNotImplemented()
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:5493
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4314
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1134
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4021
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
numbers::NumberTraits< Number >::real_type norm() const
unsigned int find_closest_vertex(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={})
Definition: grid_tools.cc:2511
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
#define DEAL_II_DEPRECATED
Definition: config.h:160
void invert_all_negative_measure_cells(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:863
static ::ExceptionBase & ExcMeshNotOrientable()
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id, const int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator >> &matched_pairs, const ::Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2013
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_layer_within_distance(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const double layer_thickness)
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4298
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
static ::ExceptionBase & ExcInternalError()
void delete_duplicated_vertices(std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
Definition: grid_tools.cc:731
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6465