Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05: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 #include <deal.II/fe/mapping_q1.h>
36 
37 #include <deal.II/grid/manifold.h>
38 #include <deal.II/grid/tria.h>
41 
43 #include <deal.II/lac/la_vector.h>
47 
48 #include <deal.II/numerics/rtree.h>
49 
51 #include <boost/archive/binary_iarchive.hpp>
52 #include <boost/archive/binary_oarchive.hpp>
53 #include <boost/random/mersenne_twister.hpp>
54 #include <boost/serialization/array.hpp>
55 #include <boost/serialization/vector.hpp>
56 
57 #ifdef DEAL_II_WITH_ZLIB
58 # include <boost/iostreams/device/back_inserter.hpp>
59 # include <boost/iostreams/filter/gzip.hpp>
60 # include <boost/iostreams/filtering_stream.hpp>
61 # include <boost/iostreams/stream.hpp>
62 #endif
64 
65 #include <bitset>
66 #include <list>
67 #include <set>
68 
70 
71 // Forward declarations
72 #ifndef DOXYGEN
73 namespace parallel
74 {
75  namespace distributed
76  {
77  template <int, int>
78  class Triangulation;
79  }
80 } // namespace parallel
81 
82 namespace hp
83 {
84  template <int, int>
85  class MappingCollection;
86 }
87 
88 class SparsityPattern;
89 
90 namespace GridTools
91 {
92  template <int dim, int spacedim>
93  class Cache;
94 }
95 #endif
96 
97 namespace internal
98 {
99  template <int dim, int spacedim, class MeshType>
101  {
102  public:
103 #ifndef _MSC_VER
104  using type = typename MeshType::active_cell_iterator;
105 #else
107 #endif
108  };
109 
110 #ifdef _MSC_VER
111  template <int dim, int spacedim>
112  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
113  {
114  public:
115  using type =
117  };
118 #endif
119 } // namespace internal
120 
129 namespace GridTools
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 #ifndef _MSC_VER
178  .template get_default_linear_mapping<dim, spacedim>()
179 #else
180  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
181 #endif
182  ));
183 
194  template <int dim, int spacedim>
195  double
198  const Mapping<dim, spacedim> & mapping =
199  (ReferenceCells::get_hypercube<dim>()
200 #ifndef _MSC_VER
201  .template get_default_linear_mapping<dim, spacedim>()
202 #else
203  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
204 #endif
205  ));
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 
241  template <int dim>
242  DEAL_II_DEPRECATED double
244  const std::vector<Point<dim>> &all_vertices,
246 
266  template <int dim>
267  double
268  cell_measure(const std::vector<Point<dim>> & all_vertices,
270 
293  template <int dim, int spacedim>
294  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
296 
323  template <int dim>
327  const Quadrature<dim> & quadrature);
328 
336  template <int dim>
337  double
340  const Quadrature<dim> & quadrature);
341 
355  template <int dim, int spacedim>
358 
376  template <typename Iterator>
379  const Iterator & object,
381 
394  template <int dim, int spacedim>
395  std::
396  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
398 
404 
421  template <int dim, int spacedim>
422  void
424  std::vector<CellData<dim>> & cells,
425  SubCellData & subcelldata);
426 
445  template <int dim, int spacedim>
446  void
447  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
448  std::vector<CellData<dim>> & cells,
449  SubCellData & subcelldata,
450  std::vector<unsigned int> & considered_vertices,
451  const double tol = 1e-12);
452 
471  template <int dim, int spacedim>
472  void
474  const std::vector<Point<spacedim>> &all_vertices,
475  std::vector<CellData<dim>> & cells);
476 
486  template <int dim, int spacedim>
487  std::size_t
489  const std::vector<Point<spacedim>> &all_vertices,
490  std::vector<CellData<dim>> & cells);
491 
502  template <int dim>
503  void
504  consistently_order_cells(std::vector<CellData<dim>> &cells);
505 
511 
574  template <int dim, typename Transformation, int spacedim>
575  void
576  transform(const Transformation & transformation,
578 
585  template <int dim, int spacedim>
586  void
587  shift(const Tensor<1, spacedim> & shift_vector,
589 
590 
601  template <int dim>
602  void
604 
617  template <int dim>
618  void
619  rotate(const Tensor<1, 3, double> &axis,
620  const double angle,
622 
638  template <int dim>
639  DEAL_II_DEPRECATED void
640  rotate(const double angle,
641  const unsigned int axis,
643 
701  template <int dim>
702  void
703  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
705  const Function<dim, double> *coefficient = nullptr,
706  const bool solve_for_absolute_positions = false);
707 
713  template <int dim, int spacedim>
714  std::map<unsigned int, Point<spacedim>>
716 
724  template <int dim, int spacedim>
725  void
726  scale(const double scaling_factor,
728 
743  template <int dim, int spacedim>
744  void
746  const double factor,
748  const bool keep_boundary = true,
749  const unsigned int seed = boost::random::mt19937::default_seed);
750 
784  template <int dim, int spacedim>
785  void
787  const bool isotropic = false,
788  const unsigned int max_iterations = 100);
789 
814  template <int dim, int spacedim>
815  void
817  const double max_ratio = 1.6180339887,
818  const unsigned int max_iterations = 5);
819 
909  template <int dim, int spacedim>
910  void
912  const double limit_angle_fraction = .75);
913 
919 
974  template <int dim, int spacedim>
975 #ifndef DOXYGEN
976  std::tuple<
977  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
978  std::vector<std::vector<Point<dim>>>,
979  std::vector<std::vector<unsigned int>>>
980 #else
981  return_type
982 #endif
984  const Cache<dim, spacedim> & cache,
985  const std::vector<Point<spacedim>> &points,
987  &cell_hint =
989 
1023  template <int dim, int spacedim>
1024 #ifndef DOXYGEN
1025  std::tuple<
1026  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1027  std::vector<std::vector<Point<dim>>>,
1028  std::vector<std::vector<unsigned int>>,
1029  std::vector<unsigned int>>
1030 #else
1031  return_type
1032 #endif
1034  const Cache<dim, spacedim> & cache,
1035  const std::vector<Point<spacedim>> &points,
1037  &cell_hint =
1039 
1109  template <int dim, int spacedim>
1110 #ifndef DOXYGEN
1111  std::tuple<
1112  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1113  std::vector<std::vector<Point<dim>>>,
1114  std::vector<std::vector<unsigned int>>,
1115  std::vector<std::vector<Point<spacedim>>>,
1116  std::vector<std::vector<unsigned int>>>
1117 #else
1118  return_type
1119 #endif
1121  const GridTools::Cache<dim, spacedim> & cache,
1122  const std::vector<Point<spacedim>> & local_points,
1123  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1124  const double tolerance = 1e-10);
1125 
1126  namespace internal
1127  {
1142  template <int dim, int spacedim>
1144  {
1151  std::vector<std::tuple<std::pair<int, int>,
1152  unsigned int,
1153  unsigned int,
1154  Point<dim>,
1156  unsigned int>>
1158 
1162  std::vector<unsigned int> send_ranks;
1163 
1169  std::vector<unsigned int> send_ptrs;
1170 
1181  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1183 
1187  std::vector<unsigned int> recv_ranks;
1188 
1194  std::vector<unsigned int> recv_ptrs;
1195  };
1196 
1206  template <int dim, int spacedim>
1209  const GridTools::Cache<dim, spacedim> & cache,
1210  const std::vector<Point<spacedim>> & points,
1211  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1212  const std::vector<bool> & marked_vertices,
1213  const double tolerance,
1214  const bool perform_handshake,
1215  const bool enforce_unique_mapping = false);
1216 
1217  } // namespace internal
1218 
1255  template <int dim, int spacedim>
1256  std::map<unsigned int, Point<spacedim>>
1258  const Triangulation<dim, spacedim> &container,
1259  const Mapping<dim, spacedim> & mapping =
1260  (ReferenceCells::get_hypercube<dim>()
1261 #ifndef _MSC_VER
1262  .template get_default_linear_mapping<dim, spacedim>()
1263 #else
1264  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1265 #endif
1266  ));
1267 
1277  template <int spacedim>
1278  unsigned int
1279  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1280  const Point<spacedim> & p);
1281 
1305  template <int dim, template <int, int> class MeshType, int spacedim>
1306  unsigned int
1307  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1308  const Point<spacedim> & p,
1309  const std::vector<bool> & marked_vertices = {});
1310 
1334  template <int dim, template <int, int> class MeshType, int spacedim>
1335  unsigned int
1337  const MeshType<dim, spacedim> &mesh,
1338  const Point<spacedim> & p,
1339  const std::vector<bool> & marked_vertices = {});
1340 
1341 
1361  template <int dim, template <int, int> class MeshType, int spacedim>
1362 #ifndef _MSC_VER
1363  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1364 #else
1365  std::vector<
1366  typename ::internal::
1367  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1368 #endif
1369  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1370  const unsigned int vertex_index);
1371 
1434  template <int dim, template <int, int> class MeshType, int spacedim>
1435 #ifndef _MSC_VER
1436  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1437 #else
1438  std::pair<typename ::internal::
1439  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1440  Point<dim>>
1441 #endif
1443  const MeshType<dim, spacedim> &mesh,
1444  const Point<spacedim> & p,
1445  const std::vector<bool> &marked_vertices = {},
1446  const double tolerance = 1.e-10);
1447 
1455  template <int dim, template <int, int> class MeshType, int spacedim>
1456 #ifndef _MSC_VER
1457  typename MeshType<dim, spacedim>::active_cell_iterator
1458 #else
1459  typename ::internal::
1460  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1461 #endif
1462  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1463  const Point<spacedim> & p,
1464  const std::vector<bool> &marked_vertices = {},
1465  const double tolerance = 1.e-10);
1466 
1473  template <int dim, int spacedim>
1474  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1475  Point<dim>>
1477  const hp::MappingCollection<dim, spacedim> &mapping,
1478  const DoFHandler<dim, spacedim> & mesh,
1479  const Point<spacedim> & p,
1480  const double tolerance = 1.e-10);
1481 
1533  template <int dim, int spacedim>
1534  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1535  Point<dim>>
1537  const Cache<dim, spacedim> &cache,
1538  const Point<spacedim> & p,
1541  const std::vector<bool> &marked_vertices = {},
1542  const double tolerance = 1.e-10);
1543 
1557  template <int dim, template <int, int> class MeshType, int spacedim>
1558 #ifndef _MSC_VER
1559  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1560 #else
1561  std::pair<typename ::internal::
1562  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1563  Point<dim>>
1564 #endif
1566  const Mapping<dim, spacedim> & mapping,
1567  const MeshType<dim, spacedim> &mesh,
1568  const Point<spacedim> & p,
1569  const std::vector<
1570  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1572  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1573  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1574  typename MeshType<dim, spacedim>::active_cell_iterator(),
1575  const std::vector<bool> & marked_vertices = {},
1576  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1577  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1578  const double tolerance = 1.e-10,
1579  const RTree<
1580  std::pair<BoundingBox<spacedim>,
1582  *relevant_cell_bounding_boxes_rtree = nullptr);
1583 
1605  template <int dim, template <int, int> class MeshType, int spacedim>
1606 #ifndef _MSC_VER
1607  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1608  Point<dim>>>
1609 #else
1610  std::vector<std::pair<
1611  typename ::internal::
1612  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1613  Point<dim>>>
1614 #endif
1616  const Mapping<dim, spacedim> & mapping,
1617  const MeshType<dim, spacedim> &mesh,
1618  const Point<spacedim> & p,
1619  const double tolerance,
1620  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1621  Point<dim>> & first_cell);
1622 
1629  template <int dim, template <int, int> class MeshType, int spacedim>
1630 #ifndef _MSC_VER
1631  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1632  Point<dim>>>
1633 #else
1634  std::vector<std::pair<
1635  typename ::internal::
1636  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1637  Point<dim>>>
1638 #endif
1640  const Mapping<dim, spacedim> & mapping,
1641  const MeshType<dim, spacedim> &mesh,
1642  const Point<spacedim> & p,
1643  const double tolerance = 1e-10,
1644  const std::vector<bool> & marked_vertices = {});
1645 
1667  template <class MeshType>
1668  std::vector<typename MeshType::active_cell_iterator>
1669  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1670 
1695  template <class MeshType>
1696  void
1698  const typename MeshType::active_cell_iterator & cell,
1699  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1700 
1750  template <class MeshType>
1751  std::vector<typename MeshType::active_cell_iterator>
1753  const MeshType &mesh,
1754  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1755  &predicate);
1756 
1757 
1765  template <class MeshType>
1766  std::vector<typename MeshType::cell_iterator>
1768  const MeshType &mesh,
1769  const std::function<bool(const typename MeshType::cell_iterator &)>
1770  & predicate,
1771  const unsigned int level);
1772 
1773 
1786  template <class MeshType>
1787  std::vector<typename MeshType::active_cell_iterator>
1788  compute_ghost_cell_halo_layer(const MeshType &mesh);
1789 
1838  template <class MeshType>
1839  std::vector<typename MeshType::active_cell_iterator>
1841  const MeshType &mesh,
1842  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1843  & predicate,
1844  const double layer_thickness);
1845 
1868  template <class MeshType>
1869  std::vector<typename MeshType::active_cell_iterator>
1870  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1871  const double layer_thickness);
1872 
1888  template <class MeshType>
1889  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1891  const MeshType &mesh,
1892  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1893  &predicate);
1894 
1947  template <class MeshType>
1948  std::vector<BoundingBox<MeshType::space_dimension>>
1950  const MeshType &mesh,
1951  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1952  & predicate,
1953  const unsigned int refinement_level = 0,
1954  const bool allow_merge = false,
1955  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1956 
1984  template <int spacedim>
1985 #ifndef DOXYGEN
1986  std::tuple<std::vector<std::vector<unsigned int>>,
1987  std::map<unsigned int, unsigned int>,
1988  std::map<unsigned int, std::vector<unsigned int>>>
1989 #else
1990  return_type
1991 #endif
1993  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1994  const std::vector<Point<spacedim>> & points);
1995 
1996 
2031  template <int spacedim>
2032 #ifndef DOXYGEN
2033  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2034  std::map<unsigned int, unsigned int>,
2035  std::map<unsigned int, std::vector<unsigned int>>>
2036 #else
2037  return_type
2038 #endif
2040  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2041  const std::vector<Point<spacedim>> & points);
2042 
2043 
2052  template <int dim, int spacedim>
2053  std::vector<
2054  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2056 
2069  template <int dim, int spacedim>
2070  std::vector<std::vector<Tensor<1, spacedim>>>
2072  const Triangulation<dim, spacedim> &mesh,
2073  const std::vector<
2075  &vertex_to_cells);
2076 
2077 
2085  template <int dim, int spacedim>
2086  unsigned int
2089  const Point<spacedim> & position,
2090  const Mapping<dim, spacedim> & mapping =
2091  (ReferenceCells::get_hypercube<dim>()
2092 #ifndef _MSC_VER
2093  .template get_default_linear_mapping<dim, spacedim>()
2094 #else
2095  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2096 #endif
2097  ));
2098 
2110  template <int dim, int spacedim>
2111  std::map<unsigned int, types::global_vertex_index>
2114 
2126  template <int dim, int spacedim>
2127  std::pair<unsigned int, double>
2130 
2136 
2145  template <int dim, int spacedim>
2146  void
2149  DynamicSparsityPattern & connectivity);
2150 
2159  template <int dim, int spacedim>
2160  void
2163  DynamicSparsityPattern & connectivity);
2164 
2173  template <int dim, int spacedim>
2174  void
2177  const unsigned int level,
2178  DynamicSparsityPattern & connectivity);
2179 
2200  template <int dim, int spacedim>
2201  void
2202  partition_triangulation(const unsigned int n_partitions,
2204  const SparsityTools::Partitioner partitioner =
2206 
2217  template <int dim, int spacedim>
2218  void
2219  partition_triangulation(const unsigned int n_partitions,
2220  const std::vector<unsigned int> &cell_weights,
2222  const SparsityTools::Partitioner partitioner =
2224 
2270  template <int dim, int spacedim>
2271  void
2272  partition_triangulation(const unsigned int n_partitions,
2273  const SparsityPattern & cell_connection_graph,
2275  const SparsityTools::Partitioner partitioner =
2277 
2288  template <int dim, int spacedim>
2289  void
2290  partition_triangulation(const unsigned int n_partitions,
2291  const std::vector<unsigned int> &cell_weights,
2292  const SparsityPattern & cell_connection_graph,
2294  const SparsityTools::Partitioner partitioner =
2296 
2311  template <int dim, int spacedim>
2312  void
2313  partition_triangulation_zorder(const unsigned int n_partitions,
2315  const bool group_siblings = true);
2316 
2328  template <int dim, int spacedim>
2329  void
2331 
2339  template <int dim, int spacedim>
2340  std::vector<types::subdomain_id>
2342  const std::vector<CellId> & cell_ids);
2343 
2354  template <int dim, int spacedim>
2355  void
2357  std::vector<types::subdomain_id> & subdomain);
2358 
2373  template <int dim, int spacedim>
2374  unsigned int
2377  const types::subdomain_id subdomain);
2378 
2408  template <int dim, int spacedim>
2409  std::vector<bool>
2411 
2417 
2451  template <typename MeshType>
2452  std::list<std::pair<typename MeshType::cell_iterator,
2453  typename MeshType::cell_iterator>>
2454  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2455 
2465  template <int dim, int spacedim>
2466  bool
2468  const Triangulation<dim, spacedim> &mesh_2);
2469 
2479  template <typename MeshType>
2480  bool
2481  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2482 
2488 
2504  template <int dim, int spacedim>
2508  & distorted_cells,
2510 
2511 
2512 
2521 
2522 
2560  template <class MeshType>
2561  std::vector<typename MeshType::active_cell_iterator>
2562  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2563 
2564 
2586  template <class Container>
2587  std::vector<typename Container::cell_iterator>
2589  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2590 
2657  template <class Container>
2658  void
2660  const std::vector<typename Container::active_cell_iterator> &patch,
2662  &local_triangulation,
2663  std::map<
2664  typename Triangulation<Container::dimension,
2665  Container::space_dimension>::active_cell_iterator,
2666  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2667 
2699  template <int dim, int spacedim>
2700  std::map<
2702  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2704 
2705 
2712 
2718  template <typename CellIterator>
2720  {
2724  CellIterator cell[2];
2725 
2730  unsigned int face_idx[2];
2731 
2737  std::bitset<3> orientation;
2738 
2752  };
2753 
2754 
2818  template <typename FaceIterator>
2819  bool
2821  std::bitset<3> & orientation,
2822  const FaceIterator & face1,
2823  const FaceIterator & face2,
2824  const unsigned int direction,
2828 
2829 
2833  template <typename FaceIterator>
2834  bool
2836  const FaceIterator & face1,
2837  const FaceIterator & face2,
2838  const unsigned int direction,
2842 
2843 
2900  template <typename MeshType>
2901  void
2903  const MeshType & mesh,
2904  const types::boundary_id b_id1,
2905  const types::boundary_id b_id2,
2906  const unsigned int direction,
2908  & matched_pairs,
2909  const Tensor<1, MeshType::space_dimension> &offset =
2912 
2913 
2936  template <typename MeshType>
2937  void
2939  const MeshType & mesh,
2940  const types::boundary_id b_id,
2941  const unsigned int direction,
2943  & matched_pairs,
2944  const ::Tensor<1, MeshType::space_dimension> &offset =
2947 
2953 
2974  template <int dim, int spacedim>
2975  void
2977  const bool reset_boundary_ids = false);
2978 
3000  template <int dim, int spacedim>
3001  void
3003  const std::vector<types::boundary_id> &src_boundary_ids,
3004  const std::vector<types::manifold_id> &dst_manifold_ids,
3006  const std::vector<types::boundary_id> &reset_boundary_ids = {});
3007 
3037  template <int dim, int spacedim>
3038  void
3040  const bool compute_face_ids = false);
3041 
3066  template <int dim, int spacedim>
3067  void
3070  const std::function<types::manifold_id(
3071  const std::set<types::manifold_id> &)> &disambiguation_function =
3072  [](const std::set<types::manifold_id> &manifold_ids) {
3073  if (manifold_ids.size() == 1)
3074  return *manifold_ids.begin();
3075  else
3077  },
3078  bool overwrite_only_flat_manifold_ids = true);
3165  template <typename DataType, typename MeshType>
3166  void
3168  const MeshType & mesh,
3169  const std::function<std_cxx17::optional<DataType>(
3170  const typename MeshType::active_cell_iterator &)> &pack,
3171  const std::function<void(const typename MeshType::active_cell_iterator &,
3172  const DataType &)> & unpack,
3173  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3174  &cell_filter =
3176 
3187  template <typename DataType, typename MeshType>
3188  void
3190  const MeshType & mesh,
3191  const std::function<std_cxx17::optional<DataType>(
3192  const typename MeshType::level_cell_iterator &)> &pack,
3193  const std::function<void(const typename MeshType::level_cell_iterator &,
3194  const DataType &)> & unpack,
3195  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3197  true});
3198 
3199  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3200  * boxes @p local_bboxes.
3201  *
3202  * This function is meant to exchange bounding boxes describing the locally
3203  * owned cells in a distributed triangulation obtained with the function
3204  * GridTools::compute_mesh_predicate_bounding_box .
3205  *
3206  * The output vector's size is the number of processes of the MPI
3207  * communicator:
3208  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3209  */
3210  template <int spacedim>
3211  std::vector<std::vector<BoundingBox<spacedim>>>
3213  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3214  const MPI_Comm & mpi_communicator);
3215 
3248  template <int spacedim>
3251  const std::vector<BoundingBox<spacedim>> &local_description,
3252  const MPI_Comm & mpi_communicator);
3253 
3271  template <int dim, int spacedim>
3272  void
3275  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3276  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3277 
3284  template <int dim, int spacedim>
3285  std::map<unsigned int, std::set<::types::subdomain_id>>
3288 
3304  template <int dim, typename T>
3306  {
3310  std::vector<CellId> cell_ids;
3311 
3315  std::vector<T> data;
3316 
3325  template <class Archive>
3326  void
3327  save(Archive &ar, const unsigned int) const
3328  {
3329  // Implement the code inline as some compilers do otherwise complain
3330  // about the use of a deprecated interface.
3331  Assert(cell_ids.size() == data.size(),
3332  ExcDimensionMismatch(cell_ids.size(), data.size()));
3333  // archive the cellids in an efficient binary format
3334  const std::size_t n_cells = cell_ids.size();
3335  ar & n_cells;
3336  for (const auto &id : cell_ids)
3337  {
3338  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3339  ar & binary_cell_id;
3340  }
3341 
3342  ar &data;
3343  }
3344 
3351  template <class Archive>
3352  void
3353  load(Archive &ar, const unsigned int)
3354  {
3355  // Implement the code inline as some compilers do otherwise complain
3356  // about the use of a deprecated interface.
3357  std::size_t n_cells;
3358  ar & n_cells;
3359  cell_ids.clear();
3360  cell_ids.reserve(n_cells);
3361  for (unsigned int c = 0; c < n_cells; ++c)
3362  {
3363  CellId::binary_type value;
3364  ar & value;
3365  cell_ids.emplace_back(value);
3366  }
3367  ar &data;
3368  }
3369 
3370 
3371 #ifdef DOXYGEN
3377  template <class Archive>
3378  void
3379  serialize(Archive &archive, const unsigned int version);
3380 #else
3381  // This macro defines the serialize() method that is compatible with
3382  // the templated save() and load() method that have been implemented.
3383  BOOST_SERIALIZATION_SPLIT_MEMBER()
3384 #endif
3385  };
3386 
3387 
3388 
3406  template <int dim, typename VectorType>
3408  {
3409  public:
3413  using value_type = typename VectorType::value_type;
3414 
3418  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3419  const FiniteElement<dim, dim> &fe,
3420  const unsigned int n_subdivisions = 1,
3421  const double tolerance = 1e-10);
3422 
3427  void
3428  process(const DoFHandler<dim> & background_dof_handler,
3429  const VectorType & ls_vector,
3430  const double iso_level,
3431  std::vector<Point<dim>> & vertices,
3432  std::vector<CellData<dim - 1>> &cells) const;
3433 
3440  void
3442  const VectorType & ls_vector,
3443  const double iso_level,
3444  std::vector<Point<dim>> & vertices,
3445  std::vector<CellData<dim - 1>> &cells) const;
3446 
3447  private:
3452  static Quadrature<dim>
3453  create_quadrature_rule(const unsigned int n_subdivisions);
3454 
3458  void
3459  process_cell(std::vector<value_type> & ls_values,
3460  const std::vector<Point<dim>> & points,
3461  const double iso_level,
3462  std::vector<Point<dim>> & vertices,
3463  std::vector<CellData<dim - 1>> &cells) const;
3464 
3471  void
3472  process_sub_cell(const std::vector<value_type> & ls_values,
3473  const std::vector<Point<2>> & points,
3474  const std::vector<unsigned int> mask,
3475  const double iso_level,
3476  std::vector<Point<2>> & vertices,
3477  std::vector<CellData<1>> & cells) const;
3478 
3482  void
3483  process_sub_cell(const std::vector<value_type> & ls_values,
3484  const std::vector<Point<3>> & points,
3485  const std::vector<unsigned int> mask,
3486  const double iso_level,
3487  std::vector<Point<3>> & vertices,
3488  std::vector<CellData<2>> & cells) const;
3489 
3494  const unsigned int n_subdivisions;
3495 
3500  const double tolerance;
3501 
3507  };
3508 
3509 
3510 
3515 
3520  int,
3521  << "The number of partitions you gave is " << arg1
3522  << ", but must be greater than zero.");
3527  int,
3528  << "The subdomain id " << arg1
3529  << " has no cells associated with it.");
3534 
3539  double,
3540  << "The scaling factor must be positive, but it is " << arg1
3541  << '.');
3542 
3547  unsigned int,
3548  << "The given vertex with index " << arg1
3549  << " is not used in the given triangulation.");
3550 
3553 } /*namespace GridTools*/
3554 
3555 
3566  "The edges of the mesh are not consistently orientable.");
3567 
3568 
3569 
3570 /* ----------------- Template function --------------- */
3571 
3572 #ifndef DOXYGEN
3573 
3574 namespace GridTools
3575 {
3576  template <int dim>
3577  double
3578  cell_measure(
3579  const std::vector<Point<dim>> &all_vertices,
3580  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3581  {
3582  // We forward call to the ArrayView version:
3583  const ArrayView<const unsigned int> view(
3585  return cell_measure(all_vertices, view);
3586  }
3587 
3588 
3589 
3590  // This specialization is defined here so that the general template in the
3591  // source file doesn't need to have further 1D overloads for the internal
3592  // functions it calls.
3593  template <>
3597  {
3598  return {};
3599  }
3600 
3601 
3602 
3603  template <int dim, typename Predicate, int spacedim>
3604  void
3605  transform(const Predicate & predicate,
3607  {
3608  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3609 
3610  // loop over all active cells, and
3611  // transform those vertices that
3612  // have not yet been touched. note
3613  // that we get to all vertices in
3614  // the triangulation by only
3615  // visiting the active cells.
3617  cell = triangulation.begin_active(),
3618  endc = triangulation.end();
3619  for (; cell != endc; ++cell)
3620  for (const unsigned int v : cell->vertex_indices())
3621  if (treated_vertices[cell->vertex_index(v)] == false)
3622  {
3623  // transform this vertex
3624  cell->vertex(v) = predicate(cell->vertex(v));
3625  // and mark it as treated
3626  treated_vertices[cell->vertex_index(v)] = true;
3627  };
3628 
3629 
3630  // now fix any vertices on hanging nodes so that we don't create any holes
3631  if (dim == 2)
3632  {
3634  cell = triangulation.begin_active(),
3635  endc = triangulation.end();
3636  for (; cell != endc; ++cell)
3637  for (const unsigned int face : cell->face_indices())
3638  if (cell->face(face)->has_children() &&
3639  !cell->face(face)->at_boundary())
3640  {
3641  Assert(cell->reference_cell() ==
3642  ReferenceCells::get_hypercube<dim>(),
3643  ExcNotImplemented());
3644 
3645  // this line has children
3646  cell->face(face)->child(0)->vertex(1) =
3647  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3648  2;
3649  }
3650  }
3651  else if (dim == 3)
3652  {
3654  cell = triangulation.begin_active(),
3655  endc = triangulation.end();
3656  for (; cell != endc; ++cell)
3657  for (const unsigned int face : cell->face_indices())
3658  if (cell->face(face)->has_children() &&
3659  !cell->face(face)->at_boundary())
3660  {
3661  Assert(cell->reference_cell() ==
3662  ReferenceCells::get_hypercube<dim>(),
3663  ExcNotImplemented());
3664 
3665  // this face has hanging nodes
3666  cell->face(face)->child(0)->vertex(1) =
3667  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3668  2.0;
3669  cell->face(face)->child(0)->vertex(2) =
3670  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3671  2.0;
3672  cell->face(face)->child(1)->vertex(3) =
3673  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3674  2.0;
3675  cell->face(face)->child(2)->vertex(3) =
3676  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3677  2.0;
3678 
3679  // center of the face
3680  cell->face(face)->child(0)->vertex(3) =
3681  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3682  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3683  4.0;
3684  }
3685  }
3686 
3687  // Make sure FEValues notices that the mesh has changed
3688  triangulation.signals.mesh_movement();
3689  }
3690 
3691 
3692 
3693  template <class MeshType>
3694  std::vector<typename MeshType::active_cell_iterator>
3695  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3696  {
3697  std::vector<typename MeshType::active_cell_iterator> child_cells;
3698 
3699  if (cell->has_children())
3700  {
3701  for (unsigned int child = 0; child < cell->n_children(); ++child)
3702  if (cell->child(child)->has_children())
3703  {
3704  const std::vector<typename MeshType::active_cell_iterator>
3705  children = get_active_child_cells<MeshType>(cell->child(child));
3706  child_cells.insert(child_cells.end(),
3707  children.begin(),
3708  children.end());
3709  }
3710  else
3711  child_cells.push_back(cell->child(child));
3712  }
3713 
3714  return child_cells;
3715  }
3716 
3717 
3718 
3719  template <class MeshType>
3720  void
3722  const typename MeshType::active_cell_iterator & cell,
3723  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3724  {
3725  active_neighbors.clear();
3726  for (const unsigned int n : cell->face_indices())
3727  if (!cell->at_boundary(n))
3728  {
3729  if (MeshType::dimension == 1)
3730  {
3731  // check children of neighbor. note
3732  // that in 1d children of the neighbor
3733  // may be further refined. In 1d the
3734  // case is simple since we know what
3735  // children bound to the present cell
3736  typename MeshType::cell_iterator neighbor_child =
3737  cell->neighbor(n);
3738  if (!neighbor_child->is_active())
3739  {
3740  while (neighbor_child->has_children())
3741  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3742 
3743  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3744  ExcInternalError());
3745  }
3746  active_neighbors.push_back(neighbor_child);
3747  }
3748  else
3749  {
3750  if (cell->face(n)->has_children())
3751  // this neighbor has children. find
3752  // out which border to the present
3753  // cell
3754  for (unsigned int c = 0;
3755  c < cell->face(n)->n_active_descendants();
3756  ++c)
3757  active_neighbors.push_back(
3758  cell->neighbor_child_on_subface(n, c));
3759  else
3760  {
3761  // the neighbor must be active
3762  // himself
3763  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3764  active_neighbors.push_back(cell->neighbor(n));
3765  }
3766  }
3767  }
3768  }
3769 
3770 
3771 
3772  namespace internal
3773  {
3774  namespace ProjectToObject
3775  {
3788  struct CrossDerivative
3789  {
3790  const unsigned int direction_0;
3791  const unsigned int direction_1;
3792 
3793  CrossDerivative(const unsigned int d0, const unsigned int d1);
3794  };
3795 
3796  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3797  const unsigned int d1)
3798  : direction_0(d0)
3799  , direction_1(d1)
3800  {}
3801 
3802 
3803 
3808  template <typename F>
3809  inline auto
3810  centered_first_difference(const double center,
3811  const double step,
3812  const F &f) -> decltype(f(center) - f(center))
3813  {
3814  return (f(center + step) - f(center - step)) / (2.0 * step);
3815  }
3816 
3817 
3818 
3823  template <typename F>
3824  inline auto
3825  centered_second_difference(const double center,
3826  const double step,
3827  const F &f) -> decltype(f(center) - f(center))
3828  {
3829  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3830  (step * step);
3831  }
3832 
3833 
3834 
3844  template <int structdim, typename F>
3845  inline auto
3846  cross_stencil(
3847  const CrossDerivative cross_derivative,
3849  const double step,
3850  const F &f) -> decltype(f(center) - f(center))
3851  {
3853  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3854  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3855  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3856  1.0 / 3.0 * f(center - simplex_vector) +
3857  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3858  step;
3859  }
3860 
3861 
3862 
3869  template <int spacedim, int structdim, typename F>
3870  inline double
3871  gradient_entry(
3872  const unsigned int row_n,
3873  const unsigned int dependent_direction,
3874  const Point<spacedim> &p0,
3876  const double step,
3877  const F & f)
3878  {
3880  dependent_direction <
3882  ExcMessage("This function assumes that the last weight is a "
3883  "dependent variable (and hence we cannot take its "
3884  "derivative directly)."));
3885  Assert(row_n != dependent_direction,
3886  ExcMessage(
3887  "We cannot differentiate with respect to the variable "
3888  "that is assumed to be dependent."));
3889 
3890  const Point<spacedim> manifold_point = f(center);
3891  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3892  {row_n, dependent_direction}, center, step, f);
3893  double entry = 0.0;
3894  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3895  entry +=
3896  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3897  return entry;
3898  }
3899 
3905  template <typename Iterator, int spacedim, int structdim>
3907  project_to_d_linear_object(const Iterator & object,
3908  const Point<spacedim> &trial_point)
3909  {
3910  // let's look at this for simplicity for a quad (structdim==2) in a
3911  // space with spacedim>2 (notate trial_point by y): all points on the
3912  // surface are given by
3913  // x(\xi) = sum_i v_i phi_x(\xi)
3914  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3915  // reference coordinates of the quad. so what we are trying to do is
3916  // find a point x on the surface that is closest to the point y. there
3917  // are different ways to solve this problem, but in the end it's a
3918  // nonlinear problem and we have to find reference coordinates \xi so
3919  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3920  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3921  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3922  // have to use a Newton method to find the answer. This leads to the
3923  // following formulation of Newton steps:
3924  //
3925  // Given \xi_k, find \delta\xi_k so that
3926  // H_k \delta\xi_k = - F_k
3927  // where H_k is an approximation to the second derivatives of J at
3928  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3929  // number of times until the right hand side is small enough. As a
3930  // stopping criterion, we terminate if ||\delta\xi||<eps.
3931  //
3932  // As for the Hessian, the best choice would be
3933  // H_k = J''(\xi_k)
3934  // but we'll opt for the simpler Gauss-Newton form
3935  // H_k = A^T A
3936  // i.e.
3937  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3938  // \partial_n phi_i *
3939  // \partial_m phi_j
3940  // we start at xi=(0.5, 0.5).
3941  Point<structdim> xi;
3942  for (unsigned int d = 0; d < structdim; ++d)
3943  xi[d] = 0.5;
3944 
3945  Point<spacedim> x_k;
3946  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3947  x_k += object->vertex(i) *
3949 
3950  do
3951  {
3953  for (const unsigned int i :
3955  F_k +=
3956  (x_k - trial_point) * object->vertex(i) *
3958  i);
3959 
3961  for (const unsigned int i :
3963  for (const unsigned int j :
3965  {
3968  xi, i),
3970  xi, j));
3971  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3972  }
3973 
3974  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3975  xi += delta_xi;
3976 
3977  x_k = Point<spacedim>();
3978  for (const unsigned int i :
3980  x_k += object->vertex(i) *
3982 
3983  if (delta_xi.norm() < 1e-7)
3984  break;
3985  }
3986  while (true);
3987 
3988  return x_k;
3989  }
3990  } // namespace ProjectToObject
3991  } // namespace internal
3992 
3993 
3994 
3995  namespace internal
3996  {
3997  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3998  // inside the project_to_object function below.
3999  template <int structdim>
4000  inline bool
4001  weights_are_ok(
4003  {
4004  // clang has trouble figuring out structdim here, so define it
4005  // again:
4006  static const std::size_t n_vertices_per_cell =
4008  n_independent_components;
4009  std::array<double, n_vertices_per_cell> copied_weights;
4010  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4011  {
4012  copied_weights[i] = v[i];
4013  if (v[i] < 0.0 || v[i] > 1.0)
4014  return false;
4015  }
4016 
4017  // check the sum: try to avoid some roundoff errors by summing in order
4018  std::sort(copied_weights.begin(), copied_weights.end());
4019  const double sum =
4020  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4021  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4022  }
4023  } // namespace internal
4024 
4025  template <typename Iterator>
4028  const Iterator & object,
4030  {
4031  const int spacedim = Iterator::AccessorType::space_dimension;
4032  const int structdim = Iterator::AccessorType::structure_dimension;
4033 
4034  Point<spacedim> projected_point = trial_point;
4035 
4036  if (structdim >= spacedim)
4037  return projected_point;
4038  else if (structdim == 1 || structdim == 2)
4039  {
4040  using namespace internal::ProjectToObject;
4041  // Try to use the special flat algorithm for quads (this is better
4042  // than the general algorithm in 3D). This does not take into account
4043  // whether projected_point is outside the quad, but we optimize along
4044  // lines below anyway:
4045  const int dim = Iterator::AccessorType::dimension;
4046  const Manifold<dim, spacedim> &manifold = object->get_manifold();
4047  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4048  &manifold) != nullptr)
4049  {
4050  projected_point =
4051  project_to_d_linear_object<Iterator, spacedim, structdim>(
4052  object, trial_point);
4053  }
4054  else
4055  {
4056  // We want to find a point on the convex hull (defined by the
4057  // vertices of the object and the manifold description) that is
4058  // relatively close to the trial point. This has a few issues:
4059  //
4060  // 1. For a general convex hull we are not guaranteed that a unique
4061  // minimum exists.
4062  // 2. The independent variables in the optimization process are the
4063  // weights given to Manifold::get_new_point, which must sum to 1,
4064  // so we cannot use standard finite differences to approximate a
4065  // gradient.
4066  //
4067  // There is not much we can do about 1., but for 2. we can derive
4068  // finite difference stencils that work on a structdim-dimensional
4069  // simplex and rewrite the optimization problem to use those
4070  // instead. Consider the structdim 2 case and let
4071  //
4072  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4073  // c2, c3})
4074  //
4075  // where {c0, c1, c2, c3} are the weights for the four vertices on
4076  // the quadrilateral. We seek to minimize the Euclidean distance
4077  // between F(...) and trial_point. We can solve for c3 in terms of
4078  // the other weights and get, for one coordinate direction
4079  //
4080  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4081  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4082  //
4083  // where we substitute back in for c3 after taking the
4084  // derivative. We can compute a stencil for the cross derivative
4085  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4086  // (and gradient_entry computes the sum over the independent
4087  // variables). Below, we somewhat arbitrarily pick the last
4088  // component as the dependent one.
4089  //
4090  // Since we can now calculate derivatives of the objective
4091  // function we can use gradient descent to minimize it.
4092  //
4093  // Of course, this is much simpler in the structdim = 1 case (we
4094  // could rewrite the projection as a 1D optimization problem), but
4095  // to reduce the potential for bugs we use the same code in both
4096  // cases.
4097  const double step_size = object->diameter() / 64.0;
4098 
4099  constexpr unsigned int n_vertices_per_cell =
4101 
4102  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4103  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4104  ++vertex_n)
4105  vertices[vertex_n] = object->vertex(vertex_n);
4106 
4107  auto get_point_from_weights =
4108  [&](const Tensor<1, n_vertices_per_cell> &weights)
4109  -> Point<spacedim> {
4110  return object->get_manifold().get_new_point(
4111  make_array_view(vertices.begin(), vertices.end()),
4112  make_array_view(weights.begin_raw(), weights.end_raw()));
4113  };
4114 
4115  // pick the initial weights as (normalized) inverse distances from
4116  // the trial point:
4117  Tensor<1, n_vertices_per_cell> guess_weights;
4118  double guess_weights_sum = 0.0;
4119  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4120  ++vertex_n)
4121  {
4122  const double distance =
4123  vertices[vertex_n].distance(trial_point);
4124  if (distance == 0.0)
4125  {
4126  guess_weights = 0.0;
4127  guess_weights[vertex_n] = 1.0;
4128  guess_weights_sum = 1.0;
4129  break;
4130  }
4131  else
4132  {
4133  guess_weights[vertex_n] = 1.0 / distance;
4134  guess_weights_sum += guess_weights[vertex_n];
4135  }
4136  }
4137  guess_weights /= guess_weights_sum;
4138  Assert(internal::weights_are_ok<structdim>(guess_weights),
4139  ExcInternalError());
4140 
4141  // The optimization algorithm consists of two parts:
4142  //
4143  // 1. An outer loop where we apply the gradient descent algorithm.
4144  // 2. An inner loop where we do a line search to find the optimal
4145  // length of the step one should take in the gradient direction.
4146  //
4147  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4148  {
4149  const unsigned int dependent_direction =
4150  n_vertices_per_cell - 1;
4151  Tensor<1, n_vertices_per_cell> current_gradient;
4152  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4153  ++row_n)
4154  {
4155  if (row_n != dependent_direction)
4156  {
4157  current_gradient[row_n] =
4158  gradient_entry<spacedim, structdim>(
4159  row_n,
4160  dependent_direction,
4161  trial_point,
4162  guess_weights,
4163  step_size,
4164  get_point_from_weights);
4165 
4166  current_gradient[dependent_direction] -=
4167  current_gradient[row_n];
4168  }
4169  }
4170 
4171  // We need to travel in the -gradient direction, as noted
4172  // above, but we may not want to take a full step in that
4173  // direction; instead, guess that we will go -0.5*gradient and
4174  // do quasi-Newton iteration to pick the best multiplier. The
4175  // goal is to find a scalar alpha such that
4176  //
4177  // F(x - alpha g)
4178  //
4179  // is minimized, where g is the gradient and F is the
4180  // objective function. To find the optimal value we find roots
4181  // of the derivative of the objective function with respect to
4182  // alpha by Newton iteration, where we approximate the first
4183  // and second derivatives of F(x - alpha g) with centered
4184  // finite differences.
4185  double gradient_weight = -0.5;
4186  auto gradient_weight_objective_function =
4187  [&](const double gradient_weight_guess) -> double {
4188  return (trial_point -
4189  get_point_from_weights(guess_weights +
4190  gradient_weight_guess *
4191  current_gradient))
4192  .norm_square();
4193  };
4194 
4195  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4196  {
4197  const double update_numerator = centered_first_difference(
4198  gradient_weight,
4199  step_size,
4200  gradient_weight_objective_function);
4201  const double update_denominator =
4202  centered_second_difference(
4203  gradient_weight,
4204  step_size,
4205  gradient_weight_objective_function);
4206 
4207  // avoid division by zero. Note that we limit the gradient
4208  // weight below
4209  if (std::abs(update_denominator) == 0.0)
4210  break;
4211  gradient_weight =
4212  gradient_weight - update_numerator / update_denominator;
4213 
4214  // Put a fairly lenient bound on the largest possible
4215  // gradient (things tend to be locally flat, so the gradient
4216  // itself is usually small)
4217  if (std::abs(gradient_weight) > 10)
4218  {
4219  gradient_weight = -10.0;
4220  break;
4221  }
4222  }
4223 
4224  // It only makes sense to take convex combinations with weights
4225  // between zero and one. If the update takes us outside of this
4226  // region then rescale the update to stay within the region and
4227  // try again
4228  Tensor<1, n_vertices_per_cell> tentative_weights =
4229  guess_weights + gradient_weight * current_gradient;
4230 
4231  double new_gradient_weight = gradient_weight;
4232  for (unsigned int iteration_count = 0; iteration_count < 40;
4233  ++iteration_count)
4234  {
4235  if (internal::weights_are_ok<structdim>(tentative_weights))
4236  break;
4237 
4238  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4239  {
4240  if (tentative_weights[i] < 0.0)
4241  {
4242  tentative_weights -=
4243  (tentative_weights[i] / current_gradient[i]) *
4244  current_gradient;
4245  }
4246  if (tentative_weights[i] < 0.0 ||
4247  1.0 < tentative_weights[i])
4248  {
4249  new_gradient_weight /= 2.0;
4250  tentative_weights =
4251  guess_weights +
4252  new_gradient_weight * current_gradient;
4253  }
4254  }
4255  }
4256 
4257  // the update might still send us outside the valid region, so
4258  // check again and quit if the update is still not valid
4259  if (!internal::weights_are_ok<structdim>(tentative_weights))
4260  break;
4261 
4262  // if we cannot get closer by traveling in the gradient
4263  // direction then quit
4264  if (get_point_from_weights(tentative_weights)
4265  .distance(trial_point) <
4266  get_point_from_weights(guess_weights).distance(trial_point))
4267  guess_weights = tentative_weights;
4268  else
4269  break;
4270  Assert(internal::weights_are_ok<structdim>(guess_weights),
4271  ExcInternalError());
4272  }
4273  Assert(internal::weights_are_ok<structdim>(guess_weights),
4274  ExcInternalError());
4275  projected_point = get_point_from_weights(guess_weights);
4276  }
4277 
4278  // if structdim == 2 and the optimal point is not on the interior then
4279  // we may be able to get a more accurate result by projecting onto the
4280  // lines.
4281  if (structdim == 2)
4282  {
4283  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4284  line_projections;
4285  for (unsigned int line_n = 0;
4286  line_n < GeometryInfo<structdim>::lines_per_cell;
4287  ++line_n)
4288  {
4289  line_projections[line_n] =
4290  project_to_object(object->line(line_n), trial_point);
4291  }
4292  std::sort(line_projections.begin(),
4293  line_projections.end(),
4294  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4295  return a.distance(trial_point) <
4296  b.distance(trial_point);
4297  });
4298  if (line_projections[0].distance(trial_point) <
4299  projected_point.distance(trial_point))
4300  projected_point = line_projections[0];
4301  }
4302  }
4303  else
4304  {
4305  Assert(false, ExcNotImplemented());
4306  return projected_point;
4307  }
4308 
4309  return projected_point;
4310  }
4311 
4312 
4313 
4314  namespace internal
4315  {
4316  template <typename DataType,
4317  typename MeshType,
4318  typename MeshCellIteratorType>
4319  inline void
4320  exchange_cell_data(
4321  const MeshType &mesh,
4322  const std::function<
4323  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4324  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4325  & unpack,
4326  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4327  const std::function<void(
4328  const std::function<void(const MeshCellIteratorType &,
4329  const types::subdomain_id)> &)> &process_cells,
4330  const std::function<std::set<types::subdomain_id>(
4331  const parallel::TriangulationBase<MeshType::dimension,
4332  MeshType::space_dimension> &)>
4333  &compute_ghost_owners)
4334  {
4335 # ifndef DEAL_II_WITH_MPI
4336  (void)mesh;
4337  (void)pack;
4338  (void)unpack;
4339  (void)cell_filter;
4340  (void)process_cells;
4341  (void)compute_ghost_owners;
4342  Assert(false, ExcNeedsMPI());
4343 # else
4344  constexpr int dim = MeshType::dimension;
4345  constexpr int spacedim = MeshType::space_dimension;
4346  auto tria =
4347  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4348  &mesh.get_triangulation());
4349  Assert(
4350  tria != nullptr,
4351  ExcMessage(
4352  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4353 
4354  if (const auto tria = dynamic_cast<
4356  &mesh.get_triangulation()))
4357  {
4358  Assert(
4359  tria->with_artificial_cells(),
4360  ExcMessage(
4361  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4362  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4363  "operate on a single layer of ghost cells. However, you have "
4364  "given a Triangulation object of type "
4365  "parallel::shared::Triangulation without artificial cells "
4366  "resulting in arbitrary numbers of ghost layers."));
4367  }
4368 
4369  // build list of cells to request for each neighbor
4370  std::set<::types::subdomain_id> ghost_owners =
4371  compute_ghost_owners(*tria);
4372  std::map<::types::subdomain_id,
4373  std::vector<typename CellId::binary_type>>
4374  neighbor_cell_list;
4375 
4376  for (const auto ghost_owner : ghost_owners)
4377  neighbor_cell_list[ghost_owner] = {};
4378 
4379  process_cells([&](const auto &cell, const auto key) -> void {
4380  if (cell_filter(cell))
4381  {
4382  constexpr int spacedim = MeshType::space_dimension;
4383  neighbor_cell_list[key].emplace_back(
4384  cell->id().template to_binary<spacedim>());
4385  }
4386  });
4387 
4388  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4389  ExcInternalError());
4390 
4391 
4392  // Before sending & receiving, make sure we protect this section with
4393  // a mutex:
4394  static Utilities::MPI::CollectiveMutex mutex;
4396  mutex, tria->get_communicator());
4397 
4398  const int mpi_tag =
4400  const int mpi_tag_reply =
4402 
4403  // send our requests
4404  std::vector<MPI_Request> requests(ghost_owners.size());
4405  {
4406  unsigned int idx = 0;
4407  for (const auto &it : neighbor_cell_list)
4408  {
4409  // send the data about the relevant cells
4410  const int ierr = MPI_Isend(it.second.data(),
4411  it.second.size() * sizeof(it.second[0]),
4412  MPI_BYTE,
4413  it.first,
4414  mpi_tag,
4416  &requests[idx]);
4417  AssertThrowMPI(ierr);
4418  ++idx;
4419  }
4420  }
4421 
4422  // receive requests and reply with the results
4423  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4424  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4425 
4426  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4427  {
4428  MPI_Status status;
4429  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4430  mpi_tag,
4432  &status);
4433  AssertThrowMPI(ierr);
4434 
4435  int len;
4436  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4437  AssertThrowMPI(ierr);
4438  Assert(len % sizeof(typename CellId::binary_type) == 0,
4439  ExcInternalError());
4440 
4441  const unsigned int n_cells =
4442  len / sizeof(typename CellId::binary_type);
4443  std::vector<typename CellId::binary_type> cells_with_requests(
4444  n_cells);
4445  std::vector<DataType> data_to_send;
4446  data_to_send.reserve(n_cells);
4447  std::vector<bool> cell_carries_data(n_cells, false);
4448 
4449  ierr = MPI_Recv(cells_with_requests.data(),
4450  len,
4451  MPI_BYTE,
4452  status.MPI_SOURCE,
4453  status.MPI_TAG,
4455  &status);
4456  AssertThrowMPI(ierr);
4457 
4458  // store data for each cell
4459  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4460  {
4461  const auto cell =
4462  tria->create_cell_iterator(CellId(cells_with_requests[c]));
4463 
4464  MeshCellIteratorType mesh_it(tria,
4465  cell->level(),
4466  cell->index(),
4467  &mesh);
4468 
4469  const std_cxx17::optional<DataType> data = pack(mesh_it);
4470  if (data)
4471  {
4472  data_to_send.emplace_back(*data);
4473  cell_carries_data[c] = true;
4474  }
4475  }
4476 
4477  // collect data for sending the reply in a buffer
4478 
4479  // (a) make room for storing the local offsets in case we receive
4480  // other data
4481  sendbuffers[idx].resize(sizeof(std::size_t));
4482 
4483  // (b) append the actual data and store how much memory it
4484  // corresponds to, which we then insert into the leading position of
4485  // the sendbuffer
4486  std::size_t size_of_send =
4487  Utilities::pack(data_to_send,
4488  sendbuffers[idx],
4489  /*enable_compression*/ false);
4490  std::memcpy(sendbuffers[idx].data(),
4491  &size_of_send,
4492  sizeof(std::size_t));
4493 
4494  // (c) append information of certain cells that got left out in case
4495  // we need it
4496  if (data_to_send.size() < n_cells)
4497  Utilities::pack(cell_carries_data,
4498  sendbuffers[idx],
4499  /*enable_compression*/ false);
4500 
4501  // send data
4502  ierr = MPI_Isend(sendbuffers[idx].data(),
4503  sendbuffers[idx].size(),
4504  MPI_BYTE,
4505  status.MPI_SOURCE,
4506  mpi_tag_reply,
4508  &reply_requests[idx]);
4509  AssertThrowMPI(ierr);
4510  }
4511 
4512  // finally receive the replies
4513  std::vector<char> receive;
4514  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4515  {
4516  MPI_Status status;
4517  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4518  mpi_tag_reply,
4520  &status);
4521  AssertThrowMPI(ierr);
4522 
4523  int len;
4524  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4525  AssertThrowMPI(ierr);
4526 
4527  receive.resize(len);
4528 
4529  ierr = MPI_Recv(receive.data(),
4530  len,
4531  MPI_BYTE,
4532  status.MPI_SOURCE,
4533  status.MPI_TAG,
4535  &status);
4536  AssertThrowMPI(ierr);
4537 
4538  // (a) first determine the length of the data section in the
4539  // received buffer
4540  auto data_iterator = receive.begin();
4541  std::size_t size_of_received_data =
4542  Utilities::unpack<std::size_t>(data_iterator,
4543  data_iterator + sizeof(std::size_t));
4544  data_iterator += sizeof(std::size_t);
4545 
4546  // (b) unpack the data section in the indicated region
4547  auto received_data = Utilities::unpack<std::vector<DataType>>(
4548  data_iterator,
4549  data_iterator + size_of_received_data,
4550  /*enable_compression*/ false);
4551  data_iterator += size_of_received_data;
4552 
4553  // (c) check if the received data contained fewer entries than the
4554  // number of cells we identified in the beginning, in which case we
4555  // need to extract the boolean vector with the relevant information
4556  const std::vector<typename CellId::binary_type> &this_cell_list =
4557  neighbor_cell_list[status.MPI_SOURCE];
4558  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4559  std::vector<bool> cells_with_data;
4560  if (received_data.size() < this_cell_list.size())
4561  {
4562  cells_with_data = Utilities::unpack<std::vector<bool>>(
4563  data_iterator, receive.end(), /*enable_compression*/ false);
4564  AssertDimension(cells_with_data.size(), this_cell_list.size());
4565  }
4566 
4567  // (d) go through the received data and call the user-provided
4568  // unpack function
4569  auto received_data_iterator = received_data.begin();
4570  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4571  if (cells_with_data.empty() || cells_with_data[c])
4572  {
4574  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4575 
4576  MeshCellIteratorType cell(tria,
4577  tria_cell->level(),
4578  tria_cell->index(),
4579  &mesh);
4580 
4581  unpack(cell, *received_data_iterator);
4582  ++received_data_iterator;
4583  }
4584  }
4585 
4586  // make sure that all communication is finished
4587  // when we leave this function.
4588  if (requests.size() > 0)
4589  {
4590  const int ierr =
4591  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4592  AssertThrowMPI(ierr);
4593  }
4594  if (reply_requests.size() > 0)
4595  {
4596  const int ierr = MPI_Waitall(reply_requests.size(),
4597  reply_requests.data(),
4598  MPI_STATUSES_IGNORE);
4599  AssertThrowMPI(ierr);
4600  }
4601 
4602 
4603 # endif // DEAL_II_WITH_MPI
4604  }
4605 
4606  } // namespace internal
4607 
4608  template <typename DataType, typename MeshType>
4609  inline void
4611  const MeshType & mesh,
4612  const std::function<std_cxx17::optional<DataType>(
4613  const typename MeshType::active_cell_iterator &)> &pack,
4614  const std::function<void(const typename MeshType::active_cell_iterator &,
4615  const DataType &)> & unpack,
4616  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4617  &cell_filter)
4618  {
4619 # ifndef DEAL_II_WITH_MPI
4620  (void)mesh;
4621  (void)pack;
4622  (void)unpack;
4623  (void)cell_filter;
4624  Assert(false, ExcNeedsMPI());
4625 # else
4626  internal::exchange_cell_data<DataType,
4627  MeshType,
4628  typename MeshType::active_cell_iterator>(
4629  mesh,
4630  pack,
4631  unpack,
4632  cell_filter,
4633  [&](const auto &process) {
4634  for (const auto &cell : mesh.active_cell_iterators())
4635  if (cell->is_ghost())
4636  process(cell, cell->subdomain_id());
4637  },
4638  [](const auto &tria) { return tria.ghost_owners(); });
4639 # endif
4640  }
4641 
4642 
4643 
4644  template <typename DataType, typename MeshType>
4645  inline void
4647  const MeshType & mesh,
4648  const std::function<std_cxx17::optional<DataType>(
4649  const typename MeshType::level_cell_iterator &)> &pack,
4650  const std::function<void(const typename MeshType::level_cell_iterator &,
4651  const DataType &)> & unpack,
4652  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4653  &cell_filter)
4654  {
4655 # ifndef DEAL_II_WITH_MPI
4656  (void)mesh;
4657  (void)pack;
4658  (void)unpack;
4659  (void)cell_filter;
4660  Assert(false, ExcNeedsMPI());
4661 # else
4662  internal::exchange_cell_data<DataType,
4663  MeshType,
4664  typename MeshType::level_cell_iterator>(
4665  mesh,
4666  pack,
4667  unpack,
4668  cell_filter,
4669  [&](const auto &process) {
4670  for (const auto &cell : mesh.cell_iterators())
4671  if (cell->is_ghost_on_level())
4672  process(cell, cell->level_subdomain_id());
4673  },
4674  [](const auto &tria) { return tria.level_ghost_owners(); });
4675 # endif
4676  }
4677 } // namespace GridTools
4678 
4679 #endif // DOXYGEN
4680 
4682 
4683 #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:699
Definition: cell_id.h:71
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
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 >> &cells) const
Definition: grid_tools.cc:6694
typename VectorType::value_type value_type
Definition: grid_tools.h:3413
const unsigned int n_subdivisions
Definition: grid_tools.h:3494
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:6631
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:6648
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:6678
void process_sub_cell(const std::vector< value_type > &ls_values, const std::vector< Point< 2 >> &points, const std::vector< unsigned int > mask, const double iso_level, std::vector< Point< 2 >> &vertices, std::vector< CellData< 1 >> &cells) const
Definition: grid_tools.cc:6859
Definition: point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
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:104
#define DEAL_II_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
Point< 3 > center
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4607
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1296
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:464
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:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
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:6450
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:4977
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:5069
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:5002
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6573
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< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:6291
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})
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:6387
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:5097
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:5982
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3784
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:289
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:3137
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2087
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4124
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2230
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:4935
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:6246
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:939
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:3055
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:741
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:6227
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:4415
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:5279
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5250
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:4442
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2053
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:5624
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5217
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:4230
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:2741
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::size_t invert_cells_with_negative_measure(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:872
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3284
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4399
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1962
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:5595
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 rotate(const double angle, Triangulation< dim > &triangulation)
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:3883
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:313
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3382
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:139
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:636
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:379
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 get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3818
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
std::vector< types::subdomain_id > get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const std::vector< CellId > &cell_ids)
Definition: grid_tools.cc:4329
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:2262
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)
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:3847
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:2613
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:3431
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:395
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:4469
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5185
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:535
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:5792
@ 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:1483
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1648
Definition: hp.h:118
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13717
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::manifold_id flat_manifold_id
Definition: types.h:269
unsigned int manifold_id
Definition: types.h:141
unsigned int global_dof_index
Definition: types.h:76
unsigned int subdomain_id
Definition: types.h:43
unsigned int boundary_id
Definition: types.h:129
::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:146
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:3353
void serialize(Archive &archive, const unsigned int version)
void save(Archive &ar, const unsigned int) const
Definition: grid_tools.h:3327
std::vector< CellId > cell_ids
Definition: grid_tools.h:3310
std::bitset< 3 > orientation
Definition: grid_tools.h:2737
FullMatrix< double > matrix
Definition: grid_tools.h:2751
unsigned int face_idx[2]
Definition: grid_tools.h:2730
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1182
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1157
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