Reference documentation for deal.II version GIT 20f059c89a 2022-12-04 14:55:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_grid_tools_h
17 #define dealii_grid_tools_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 #include <deal.II/base/point.h>
26 
28 
30 
32 
33 #include <deal.II/fe/fe_values.h>
34 #include <deal.II/fe/mapping.h>
35 
36 #include <deal.II/grid/manifold.h>
37 #include <deal.II/grid/tria.h>
40 
42 #include <deal.II/lac/la_vector.h>
46 
47 #include <deal.II/numerics/rtree.h>
48 
50 #include <boost/archive/binary_iarchive.hpp>
51 #include <boost/archive/binary_oarchive.hpp>
52 #include <boost/random/mersenne_twister.hpp>
53 #include <boost/serialization/array.hpp>
54 #include <boost/serialization/vector.hpp>
55 
56 #ifdef DEAL_II_WITH_ZLIB
57 # include <boost/iostreams/device/back_inserter.hpp>
58 # include <boost/iostreams/filter/gzip.hpp>
59 # include <boost/iostreams/filtering_stream.hpp>
60 # include <boost/iostreams/stream.hpp>
61 #endif
63 
64 #include <bitset>
65 #include <list>
66 #include <set>
67 
69 
70 // Forward declarations
71 #ifndef DOXYGEN
72 namespace parallel
73 {
74  namespace distributed
75  {
76  template <int, int>
77  class Triangulation;
78  }
79 } // namespace parallel
80 
81 namespace hp
82 {
83  template <int, int>
84  class MappingCollection;
85 }
86 
87 class SparsityPattern;
88 
89 namespace GridTools
90 {
91  template <int dim, int spacedim>
92  class Cache;
93 }
94 #endif
95 
96 namespace internal
97 {
98  template <int dim, int spacedim, class MeshType>
100  {
101  public:
102 #ifndef _MSC_VER
103  using type = typename MeshType::active_cell_iterator;
104 #else
106 #endif
107  };
108 
109 #ifdef _MSC_VER
110  template <int dim, int spacedim>
111  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
112  {
113  public:
114  using type =
116  };
117 #endif
118 } // namespace internal
119 
128 namespace GridTools
129 {
141  template <int dim, int spacedim>
142  double
144 
171  template <int dim, int spacedim>
172  double
174  const Mapping<dim, spacedim> & mapping =
175  (ReferenceCells::get_hypercube<dim>()
176 #ifndef _MSC_VER
177  .template get_default_linear_mapping<dim, spacedim>()
178 #else
179  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
180 #endif
181  ));
182 
193  template <int dim, int spacedim>
194  double
197  const Mapping<dim, spacedim> & mapping =
198  (ReferenceCells::get_hypercube<dim>()
199 #ifndef _MSC_VER
200  .template get_default_linear_mapping<dim, spacedim>()
201 #else
202  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
203 #endif
204  ));
205 
216  template <int dim, int spacedim>
217  double
220  const Mapping<dim, spacedim> & mapping =
221  (ReferenceCells::get_hypercube<dim>()
222 #ifndef _MSC_VER
223  .template get_default_linear_mapping<dim, spacedim>()
224 #else
225  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
226 #endif
227  ));
228 
240  template <int dim>
241  DEAL_II_DEPRECATED double
243  const std::vector<Point<dim>> &all_vertices,
245 
265  template <int dim>
266  double
267  cell_measure(const std::vector<Point<dim>> & all_vertices,
269 
292  template <int dim, int spacedim>
293  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
295 
322  template <int dim>
326  const Quadrature<dim> & quadrature);
327 
335  template <int dim>
336  double
339  const Quadrature<dim> & quadrature);
340 
354  template <int dim, int spacedim>
357 
375  template <typename Iterator>
378  const Iterator & object,
380 
393  template <int dim, int spacedim>
394  std::
395  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
397 
420  template <int dim, int spacedim>
421  void
423  std::vector<CellData<dim>> & cells,
424  SubCellData & subcelldata);
425 
444  template <int dim, int spacedim>
445  void
446  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
447  std::vector<CellData<dim>> & cells,
448  SubCellData & subcelldata,
449  std::vector<unsigned int> & considered_vertices,
450  const double tol = 1e-12);
451 
459  template <int dim>
460  void
462  const double tol = 1e-12);
463 
482  template <int dim, int spacedim>
483  void
485  const std::vector<Point<spacedim>> &all_vertices,
486  std::vector<CellData<dim>> & cells);
487 
497  template <int dim, int spacedim>
498  std::size_t
500  const std::vector<Point<spacedim>> &all_vertices,
501  std::vector<CellData<dim>> & cells);
502 
513  template <int dim>
514  void
515  consistently_order_cells(std::vector<CellData<dim>> &cells);
516 
585  template <int dim, typename Transformation, int spacedim>
586  void
587  transform(const Transformation & transformation,
589 
596  template <int dim, int spacedim>
597  void
598  shift(const Tensor<1, spacedim> & shift_vector,
600 
601 
612  template <int dim>
613  void
615 
628  template <int dim>
629  void
630  rotate(const Tensor<1, 3, double> &axis,
631  const double angle,
633 
649  template <int dim>
650  DEAL_II_DEPRECATED void
651  rotate(const double angle,
652  const unsigned int axis,
654 
712  template <int dim>
713  void
714  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
716  const Function<dim, double> *coefficient = nullptr,
717  const bool solve_for_absolute_positions = false);
718 
724  template <int dim, int spacedim>
725  std::map<unsigned int, Point<spacedim>>
727 
735  template <int dim, int spacedim>
736  void
737  scale(const double scaling_factor,
739 
754  template <int dim, int spacedim>
755  void
757  const double factor,
759  const bool keep_boundary = true,
760  const unsigned int seed = boost::random::mt19937::default_seed);
761 
795  template <int dim, int spacedim>
796  void
798  const bool isotropic = false,
799  const unsigned int max_iterations = 100);
800 
825  template <int dim, int spacedim>
826  void
828  const double max_ratio = 1.6180339887,
829  const unsigned int max_iterations = 5);
830 
920  template <int dim, int spacedim>
921  void
923  const double limit_angle_fraction = .75);
924 
985  template <int dim, int spacedim>
986 #ifndef DOXYGEN
987  std::tuple<
988  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
989  std::vector<std::vector<Point<dim>>>,
990  std::vector<std::vector<unsigned int>>>
991 #else
992  return_type
993 #endif
995  const Cache<dim, spacedim> & cache,
996  const std::vector<Point<spacedim>> &points,
998  &cell_hint =
1000 
1034  template <int dim, int spacedim>
1035 #ifndef DOXYGEN
1036  std::tuple<
1037  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1038  std::vector<std::vector<Point<dim>>>,
1039  std::vector<std::vector<unsigned int>>,
1040  std::vector<unsigned int>>
1041 #else
1042  return_type
1043 #endif
1045  const Cache<dim, spacedim> & cache,
1046  const std::vector<Point<spacedim>> &points,
1048  &cell_hint =
1050 
1120  template <int dim, int spacedim>
1121 #ifndef DOXYGEN
1122  std::tuple<
1123  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1124  std::vector<std::vector<Point<dim>>>,
1125  std::vector<std::vector<unsigned int>>,
1126  std::vector<std::vector<Point<spacedim>>>,
1127  std::vector<std::vector<unsigned int>>>
1128 #else
1129  return_type
1130 #endif
1132  const GridTools::Cache<dim, spacedim> & cache,
1133  const std::vector<Point<spacedim>> & local_points,
1134  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1135  const double tolerance = 1e-10);
1136 
1137  namespace internal
1138  {
1153  template <int dim, int spacedim>
1155  {
1162  std::vector<std::tuple<std::pair<int, int>,
1163  unsigned int,
1164  unsigned int,
1165  Point<dim>,
1167  unsigned int>>
1169 
1173  std::vector<unsigned int> send_ranks;
1174 
1180  std::vector<unsigned int> send_ptrs;
1181 
1192  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1194 
1198  std::vector<unsigned int> recv_ranks;
1199 
1205  std::vector<unsigned int> recv_ptrs;
1206  };
1207 
1217  template <int dim, int spacedim>
1220  const GridTools::Cache<dim, spacedim> & cache,
1221  const std::vector<Point<spacedim>> & points,
1222  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1223  const std::vector<bool> & marked_vertices,
1224  const double tolerance,
1225  const bool perform_handshake,
1226  const bool enforce_unique_mapping = false);
1227 
1228  } // namespace internal
1229 
1266  template <int dim, int spacedim>
1267  std::map<unsigned int, Point<spacedim>>
1269  const Triangulation<dim, spacedim> &container,
1270  const Mapping<dim, spacedim> & mapping =
1271  (ReferenceCells::get_hypercube<dim>()
1272 #ifndef _MSC_VER
1273  .template get_default_linear_mapping<dim, spacedim>()
1274 #else
1275  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1276 #endif
1277  ));
1278 
1288  template <int spacedim>
1289  unsigned int
1290  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1291  const Point<spacedim> & p);
1292 
1316  template <int dim, template <int, int> class MeshType, int spacedim>
1317  unsigned int
1318  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1319  const Point<spacedim> & p,
1320  const std::vector<bool> & marked_vertices = {});
1321 
1345  template <int dim, template <int, int> class MeshType, int spacedim>
1346  unsigned int
1348  const MeshType<dim, spacedim> &mesh,
1349  const Point<spacedim> & p,
1350  const std::vector<bool> & marked_vertices = {});
1351 
1352 
1372  template <int dim, template <int, int> class MeshType, int spacedim>
1373 #ifndef _MSC_VER
1374  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1375 #else
1376  std::vector<
1377  typename ::internal::
1378  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1379 #endif
1380  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1381  const unsigned int vertex_index);
1382 
1445  template <int dim, template <int, int> class MeshType, int spacedim>
1446 #ifndef _MSC_VER
1447  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1448 #else
1449  std::pair<typename ::internal::
1450  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1451  Point<dim>>
1452 #endif
1454  const MeshType<dim, spacedim> &mesh,
1455  const Point<spacedim> & p,
1456  const std::vector<bool> &marked_vertices = {},
1457  const double tolerance = 1.e-10);
1458 
1466  template <int dim, template <int, int> class MeshType, int spacedim>
1467 #ifndef _MSC_VER
1468  typename MeshType<dim, spacedim>::active_cell_iterator
1469 #else
1470  typename ::internal::
1471  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1472 #endif
1473  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1474  const Point<spacedim> & p,
1475  const std::vector<bool> &marked_vertices = {},
1476  const double tolerance = 1.e-10);
1477 
1484  template <int dim, int spacedim>
1485  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1486  Point<dim>>
1488  const hp::MappingCollection<dim, spacedim> &mapping,
1489  const DoFHandler<dim, spacedim> & mesh,
1490  const Point<spacedim> & p,
1491  const double tolerance = 1.e-10);
1492 
1544  template <int dim, int spacedim>
1545  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1546  Point<dim>>
1548  const Cache<dim, spacedim> &cache,
1549  const Point<spacedim> & p,
1552  const std::vector<bool> &marked_vertices = {},
1553  const double tolerance = 1.e-10);
1554 
1568  template <int dim, template <int, int> class MeshType, int spacedim>
1569 #ifndef _MSC_VER
1570  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1571 #else
1572  std::pair<typename ::internal::
1573  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1574  Point<dim>>
1575 #endif
1577  const Mapping<dim, spacedim> & mapping,
1578  const MeshType<dim, spacedim> &mesh,
1579  const Point<spacedim> & p,
1580  const std::vector<
1581  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1583  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1584  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1585  typename MeshType<dim, spacedim>::active_cell_iterator(),
1586  const std::vector<bool> & marked_vertices = {},
1587  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1588  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1589  const double tolerance = 1.e-10,
1590  const RTree<
1591  std::pair<BoundingBox<spacedim>,
1593  *relevant_cell_bounding_boxes_rtree = nullptr);
1594 
1616  template <int dim, template <int, int> class MeshType, int spacedim>
1617 #ifndef _MSC_VER
1618  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1619  Point<dim>>>
1620 #else
1621  std::vector<std::pair<
1622  typename ::internal::
1623  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1624  Point<dim>>>
1625 #endif
1627  const Mapping<dim, spacedim> & mapping,
1628  const MeshType<dim, spacedim> &mesh,
1629  const Point<spacedim> & p,
1630  const double tolerance,
1631  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1632  Point<dim>> & first_cell);
1633 
1640  template <int dim, template <int, int> class MeshType, int spacedim>
1641 #ifndef _MSC_VER
1642  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1643  Point<dim>>>
1644 #else
1645  std::vector<std::pair<
1646  typename ::internal::
1647  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1648  Point<dim>>>
1649 #endif
1651  const Mapping<dim, spacedim> & mapping,
1652  const MeshType<dim, spacedim> &mesh,
1653  const Point<spacedim> & p,
1654  const double tolerance = 1e-10,
1655  const std::vector<bool> & marked_vertices = {});
1656 
1678  template <class MeshType>
1679  std::vector<typename MeshType::active_cell_iterator>
1680  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1681 
1706  template <class MeshType>
1707  void
1709  const typename MeshType::active_cell_iterator & cell,
1710  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1711 
1761  template <class MeshType>
1762  std::vector<typename MeshType::active_cell_iterator>
1764  const MeshType &mesh,
1765  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1766  &predicate);
1767 
1768 
1776  template <class MeshType>
1777  std::vector<typename MeshType::cell_iterator>
1779  const MeshType &mesh,
1780  const std::function<bool(const typename MeshType::cell_iterator &)>
1781  & predicate,
1782  const unsigned int level);
1783 
1784 
1797  template <class MeshType>
1798  std::vector<typename MeshType::active_cell_iterator>
1799  compute_ghost_cell_halo_layer(const MeshType &mesh);
1800 
1849  template <class MeshType>
1850  std::vector<typename MeshType::active_cell_iterator>
1852  const MeshType &mesh,
1853  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1854  & predicate,
1855  const double layer_thickness);
1856 
1879  template <class MeshType>
1880  std::vector<typename MeshType::active_cell_iterator>
1881  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1882  const double layer_thickness);
1883 
1899  template <class MeshType>
1900  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1902  const MeshType &mesh,
1903  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1904  &predicate);
1905 
1976  template <class MeshType>
1977  std::vector<BoundingBox<MeshType::space_dimension>>
1979  const MeshType &mesh,
1980  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1981  & predicate,
1982  const unsigned int refinement_level = 0,
1983  const bool allow_merge = false,
1984  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1985 
2013  template <int spacedim>
2014 #ifndef DOXYGEN
2015  std::tuple<std::vector<std::vector<unsigned int>>,
2016  std::map<unsigned int, unsigned int>,
2017  std::map<unsigned int, std::vector<unsigned int>>>
2018 #else
2019  return_type
2020 #endif
2022  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
2023  const std::vector<Point<spacedim>> & points);
2024 
2025 
2060  template <int spacedim>
2061 #ifndef DOXYGEN
2062  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2063  std::map<unsigned int, unsigned int>,
2064  std::map<unsigned int, std::vector<unsigned int>>>
2065 #else
2066  return_type
2067 #endif
2069  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2070  const std::vector<Point<spacedim>> & points);
2071 
2072 
2081  template <int dim, int spacedim>
2082  std::vector<
2083  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2085 
2098  template <int dim, int spacedim>
2099  std::vector<std::vector<Tensor<1, spacedim>>>
2101  const Triangulation<dim, spacedim> &mesh,
2102  const std::vector<
2104  &vertex_to_cells);
2105 
2106 
2114  template <int dim, int spacedim>
2115  unsigned int
2118  const Point<spacedim> & position,
2119  const Mapping<dim, spacedim> & mapping =
2120  (ReferenceCells::get_hypercube<dim>()
2121 #ifndef _MSC_VER
2122  .template get_default_linear_mapping<dim, spacedim>()
2123 #else
2124  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2125 #endif
2126  ));
2127 
2139  template <int dim, int spacedim>
2140  std::map<unsigned int, types::global_vertex_index>
2143 
2155  template <int dim, int spacedim>
2156  std::pair<unsigned int, double>
2159 
2174  template <int dim, int spacedim>
2175  void
2178  DynamicSparsityPattern & connectivity);
2179 
2188  template <int dim, int spacedim>
2189  void
2192  DynamicSparsityPattern & connectivity);
2193 
2202  template <int dim, int spacedim>
2203  void
2206  const unsigned int level,
2207  DynamicSparsityPattern & connectivity);
2208 
2229  template <int dim, int spacedim>
2230  void
2231  partition_triangulation(const unsigned int n_partitions,
2233  const SparsityTools::Partitioner partitioner =
2235 
2246  template <int dim, int spacedim>
2247  void
2248  partition_triangulation(const unsigned int n_partitions,
2249  const std::vector<unsigned int> &cell_weights,
2251  const SparsityTools::Partitioner partitioner =
2253 
2299  template <int dim, int spacedim>
2300  void
2301  partition_triangulation(const unsigned int n_partitions,
2302  const SparsityPattern & cell_connection_graph,
2304  const SparsityTools::Partitioner partitioner =
2306 
2317  template <int dim, int spacedim>
2318  void
2319  partition_triangulation(const unsigned int n_partitions,
2320  const std::vector<unsigned int> &cell_weights,
2321  const SparsityPattern & cell_connection_graph,
2323  const SparsityTools::Partitioner partitioner =
2325 
2340  template <int dim, int spacedim>
2341  void
2342  partition_triangulation_zorder(const unsigned int n_partitions,
2344  const bool group_siblings = true);
2345 
2357  template <int dim, int spacedim>
2358  void
2360 
2368  template <int dim, int spacedim>
2369  std::vector<types::subdomain_id>
2371  const std::vector<CellId> & cell_ids);
2372 
2383  template <int dim, int spacedim>
2384  void
2386  std::vector<types::subdomain_id> & subdomain);
2387 
2402  template <int dim, int spacedim>
2403  unsigned int
2406  const types::subdomain_id subdomain);
2407 
2437  template <int dim, int spacedim>
2438  std::vector<bool>
2440 
2480  template <typename MeshType>
2481  std::list<std::pair<typename MeshType::cell_iterator,
2482  typename MeshType::cell_iterator>>
2483  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2484 
2494  template <int dim, int spacedim>
2495  bool
2497  const Triangulation<dim, spacedim> &mesh_2);
2498 
2508  template <typename MeshType>
2509  bool
2510  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2511 
2533  template <int dim, int spacedim>
2537  & distorted_cells,
2539 
2540 
2541 
2589  template <class MeshType>
2590  std::vector<typename MeshType::active_cell_iterator>
2591  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2592 
2593 
2615  template <class Container>
2616  std::vector<typename Container::cell_iterator>
2618  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2619 
2686  template <class Container>
2687  void
2689  const std::vector<typename Container::active_cell_iterator> &patch,
2691  &local_triangulation,
2692  std::map<
2693  typename Triangulation<Container::dimension,
2694  Container::space_dimension>::active_cell_iterator,
2695  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2696 
2728  template <int dim, int spacedim>
2729  std::map<
2731  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2733 
2734 
2747  template <typename CellIterator>
2749  {
2753  CellIterator cell[2];
2754 
2759  unsigned int face_idx[2];
2760 
2766  std::bitset<3> orientation;
2767 
2781  };
2782 
2783 
2847  template <typename FaceIterator>
2848  bool
2850  std::bitset<3> & orientation,
2851  const FaceIterator & face1,
2852  const FaceIterator & face2,
2853  const unsigned int direction,
2857 
2858 
2862  template <typename FaceIterator>
2863  bool
2865  const FaceIterator & face1,
2866  const FaceIterator & face2,
2867  const unsigned int direction,
2871 
2872 
2929  template <typename MeshType>
2930  void
2932  const MeshType & mesh,
2933  const types::boundary_id b_id1,
2934  const types::boundary_id b_id2,
2935  const unsigned int direction,
2937  & matched_pairs,
2938  const Tensor<1, MeshType::space_dimension> &offset =
2941 
2942 
2965  template <typename MeshType>
2966  void
2968  const MeshType & mesh,
2969  const types::boundary_id b_id,
2970  const unsigned int direction,
2972  & matched_pairs,
2973  const ::Tensor<1, MeshType::space_dimension> &offset =
2976 
3003  template <int dim, int spacedim>
3004  void
3006  const bool reset_boundary_ids = false);
3007 
3029  template <int dim, int spacedim>
3030  void
3032  const std::vector<types::boundary_id> &src_boundary_ids,
3033  const std::vector<types::manifold_id> &dst_manifold_ids,
3035  const std::vector<types::boundary_id> &reset_boundary_ids = {});
3036 
3066  template <int dim, int spacedim>
3067  void
3069  const bool compute_face_ids = false);
3070 
3095  template <int dim, int spacedim>
3096  void
3099  const std::function<types::manifold_id(
3100  const std::set<types::manifold_id> &)> &disambiguation_function =
3101  [](const std::set<types::manifold_id> &manifold_ids) {
3102  if (manifold_ids.size() == 1)
3103  return *manifold_ids.begin();
3104  else
3106  },
3107  bool overwrite_only_flat_manifold_ids = true);
3194  template <typename DataType, typename MeshType>
3195  void
3197  const MeshType & mesh,
3198  const std::function<std_cxx17::optional<DataType>(
3199  const typename MeshType::active_cell_iterator &)> &pack,
3200  const std::function<void(const typename MeshType::active_cell_iterator &,
3201  const DataType &)> & unpack,
3202  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3203  &cell_filter =
3205 
3216  template <typename DataType, typename MeshType>
3217  void
3219  const MeshType & mesh,
3220  const std::function<std_cxx17::optional<DataType>(
3221  const typename MeshType::level_cell_iterator &)> &pack,
3222  const std::function<void(const typename MeshType::level_cell_iterator &,
3223  const DataType &)> & unpack,
3224  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3226  true});
3227 
3228  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3229  * boxes @p local_bboxes.
3230  *
3231  * This function is meant to exchange bounding boxes describing the locally
3232  * owned cells in a distributed triangulation obtained with the function
3233  * GridTools::compute_mesh_predicate_bounding_box .
3234  *
3235  * The output vector's size is the number of processes of the MPI
3236  * communicator:
3237  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3238  */
3239  template <int spacedim>
3240  std::vector<std::vector<BoundingBox<spacedim>>>
3242  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3243  const MPI_Comm & mpi_communicator);
3244 
3277  template <int spacedim>
3280  const std::vector<BoundingBox<spacedim>> &local_description,
3281  const MPI_Comm & mpi_communicator);
3282 
3300  template <int dim, int spacedim>
3301  void
3304  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3305  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3306 
3313  template <int dim, int spacedim>
3314  std::map<unsigned int, std::set<::types::subdomain_id>>
3317 
3333  template <int dim, typename T>
3335  {
3339  std::vector<CellId> cell_ids;
3340 
3344  std::vector<T> data;
3345 
3354  template <class Archive>
3355  void
3356  save(Archive &ar, const unsigned int) const
3357  {
3358  // Implement the code inline as some compilers do otherwise complain
3359  // about the use of a deprecated interface.
3360  Assert(cell_ids.size() == data.size(),
3361  ExcDimensionMismatch(cell_ids.size(), data.size()));
3362  // archive the cellids in an efficient binary format
3363  const std::size_t n_cells = cell_ids.size();
3364  ar & n_cells;
3365  for (const auto &id : cell_ids)
3366  {
3367  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3368  ar & binary_cell_id;
3369  }
3370 
3371  ar &data;
3372  }
3373 
3380  template <class Archive>
3381  void
3382  load(Archive &ar, const unsigned int)
3383  {
3384  // Implement the code inline as some compilers do otherwise complain
3385  // about the use of a deprecated interface.
3386  std::size_t n_cells;
3387  ar & n_cells;
3388  cell_ids.clear();
3389  cell_ids.reserve(n_cells);
3390  for (unsigned int c = 0; c < n_cells; ++c)
3391  {
3392  CellId::binary_type value;
3393  ar & value;
3394  cell_ids.emplace_back(value);
3395  }
3396  ar &data;
3397  }
3398 
3399 
3400 #ifdef DOXYGEN
3406  template <class Archive>
3407  void
3408  serialize(Archive &archive, const unsigned int version);
3409 #else
3410  // This macro defines the serialize() method that is compatible with
3411  // the templated save() and load() method that have been implemented.
3412  BOOST_SERIALIZATION_SPLIT_MEMBER()
3413 #endif
3414  };
3415 
3416 
3417 
3435  template <int dim, typename VectorType>
3437  {
3438  public:
3442  using value_type = typename VectorType::value_type;
3443 
3447  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3448  const FiniteElement<dim, dim> &fe,
3449  const unsigned int n_subdivisions = 1,
3450  const double tolerance = 1e-10);
3451 
3462  void
3463  process(const DoFHandler<dim> & background_dof_handler,
3464  const VectorType & ls_vector,
3465  const double iso_level,
3466  std::vector<Point<dim>> &vertices,
3467  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3468 
3473  void
3474  process(const DoFHandler<dim> & background_dof_handler,
3475  const VectorType & ls_vector,
3476  const double iso_level,
3477  std::vector<Point<dim>> &vertices) const;
3478 
3488  void
3490  const VectorType & ls_vector,
3491  const double iso_level,
3492  std::vector<Point<dim>> & vertices,
3493  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3499  void
3501  const VectorType & ls_vector,
3502  const double iso_level,
3503  std::vector<Point<dim>> &vertices) const;
3504 
3505  private:
3510  static Quadrature<dim>
3511  create_quadrature_rule(const unsigned int n_subdivisions);
3512 
3516  void
3517  process_cell(std::vector<value_type> & ls_values,
3518  const std::vector<Point<dim>> & points,
3519  const double iso_level,
3520  std::vector<Point<dim>> & vertices,
3521  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
3522  const bool write_back_cell_data = true) const;
3523 
3527  void
3528  process_sub_cell(const std::vector<value_type> &,
3529  const std::vector<Point<1>> &,
3530  const std::vector<unsigned int> &,
3531  const double,
3532  std::vector<Point<1>> &,
3533  std::vector<CellData<1>> &,
3534  const bool) const
3535  {
3536  AssertThrow(false, ExcNotImplemented());
3537  }
3538 
3545  void
3546  process_sub_cell(const std::vector<value_type> & ls_values,
3547  const std::vector<Point<2>> & points,
3548  const std::vector<unsigned int> &mask,
3549  const double iso_level,
3550  std::vector<Point<2>> & vertices,
3551  std::vector<CellData<1>> & cells,
3552  const bool write_back_cell_data) const;
3553 
3557  void
3558  process_sub_cell(const std::vector<value_type> & ls_values,
3559  const std::vector<Point<3>> & points,
3560  const std::vector<unsigned int> &mask,
3561  const double iso_level,
3562  std::vector<Point<3>> & vertices,
3563  std::vector<CellData<2>> & cells,
3564  const bool write_back_cell_data) const;
3565 
3570  const unsigned int n_subdivisions;
3571 
3576  const double tolerance;
3577 
3583  };
3584 
3585 
3586 
3596  int,
3597  << "The number of partitions you gave is " << arg1
3598  << ", but must be greater than zero.");
3603  int,
3604  << "The subdomain id " << arg1
3605  << " has no cells associated with it.");
3610 
3615  double,
3616  << "The scaling factor must be positive, but it is " << arg1
3617  << '.');
3618 
3623  unsigned int,
3624  << "The given vertex with index " << arg1
3625  << " is not used in the given triangulation.");
3626 
3629 } /*namespace GridTools*/
3630 
3631 
3642  "The edges of the mesh are not consistently orientable.");
3643 
3644 
3645 
3646 /* ----------------- Template function --------------- */
3647 
3648 #ifndef DOXYGEN
3649 
3650 namespace GridTools
3651 {
3652  template <int dim>
3653  double
3654  cell_measure(
3655  const std::vector<Point<dim>> &all_vertices,
3656  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3657  {
3658  // We forward call to the ArrayView version:
3659  const ArrayView<const unsigned int> view(
3661  return cell_measure(all_vertices, view);
3662  }
3663 
3664 
3665 
3666  // This specialization is defined here so that the general template in the
3667  // source file doesn't need to have further 1D overloads for the internal
3668  // functions it calls.
3669  template <>
3673  {
3674  return {};
3675  }
3676 
3677 
3678 
3679  template <int dim, typename Predicate, int spacedim>
3680  void
3681  transform(const Predicate & predicate,
3683  {
3684  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3685 
3686  // loop over all active cells, and
3687  // transform those vertices that
3688  // have not yet been touched. note
3689  // that we get to all vertices in
3690  // the triangulation by only
3691  // visiting the active cells.
3693  cell = triangulation.begin_active(),
3694  endc = triangulation.end();
3695  for (; cell != endc; ++cell)
3696  for (const unsigned int v : cell->vertex_indices())
3697  if (treated_vertices[cell->vertex_index(v)] == false)
3698  {
3699  // transform this vertex
3700  cell->vertex(v) = predicate(cell->vertex(v));
3701  // and mark it as treated
3702  treated_vertices[cell->vertex_index(v)] = true;
3703  };
3704 
3705 
3706  // now fix any vertices on hanging nodes so that we don't create any holes
3707  if (dim == 2)
3708  {
3710  cell = triangulation.begin_active(),
3711  endc = triangulation.end();
3712  for (; cell != endc; ++cell)
3713  for (const unsigned int face : cell->face_indices())
3714  if (cell->face(face)->has_children() &&
3715  !cell->face(face)->at_boundary())
3716  {
3717  Assert(cell->reference_cell() ==
3718  ReferenceCells::get_hypercube<dim>(),
3719  ExcNotImplemented());
3720 
3721  // this line has children
3722  cell->face(face)->child(0)->vertex(1) =
3723  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3724  2;
3725  }
3726  }
3727  else if (dim == 3)
3728  {
3730  cell = triangulation.begin_active(),
3731  endc = triangulation.end();
3732  for (; cell != endc; ++cell)
3733  for (const unsigned int face : cell->face_indices())
3734  if (cell->face(face)->has_children() &&
3735  !cell->face(face)->at_boundary())
3736  {
3737  Assert(cell->reference_cell() ==
3738  ReferenceCells::get_hypercube<dim>(),
3739  ExcNotImplemented());
3740 
3741  // this face has hanging nodes
3742  cell->face(face)->child(0)->vertex(1) =
3743  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3744  2.0;
3745  cell->face(face)->child(0)->vertex(2) =
3746  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3747  2.0;
3748  cell->face(face)->child(1)->vertex(3) =
3749  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3750  2.0;
3751  cell->face(face)->child(2)->vertex(3) =
3752  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3753  2.0;
3754 
3755  // center of the face
3756  cell->face(face)->child(0)->vertex(3) =
3757  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3758  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3759  4.0;
3760  }
3761  }
3762 
3763  // Make sure FEValues notices that the mesh has changed
3764  triangulation.signals.mesh_movement();
3765  }
3766 
3767 
3768 
3769  template <class MeshType>
3770  std::vector<typename MeshType::active_cell_iterator>
3771  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3772  {
3773  std::vector<typename MeshType::active_cell_iterator> child_cells;
3774 
3775  if (cell->has_children())
3776  {
3777  for (unsigned int child = 0; child < cell->n_children(); ++child)
3778  if (cell->child(child)->has_children())
3779  {
3780  const std::vector<typename MeshType::active_cell_iterator>
3781  children = get_active_child_cells<MeshType>(cell->child(child));
3782  child_cells.insert(child_cells.end(),
3783  children.begin(),
3784  children.end());
3785  }
3786  else
3787  child_cells.push_back(cell->child(child));
3788  }
3789 
3790  return child_cells;
3791  }
3792 
3793 
3794 
3795  template <class MeshType>
3796  void
3798  const typename MeshType::active_cell_iterator & cell,
3799  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3800  {
3801  active_neighbors.clear();
3802  for (const unsigned int n : cell->face_indices())
3803  if (!cell->at_boundary(n))
3804  {
3805  if (MeshType::dimension == 1)
3806  {
3807  // check children of neighbor. note
3808  // that in 1d children of the neighbor
3809  // may be further refined. In 1d the
3810  // case is simple since we know what
3811  // children bound to the present cell
3812  typename MeshType::cell_iterator neighbor_child =
3813  cell->neighbor(n);
3814  if (!neighbor_child->is_active())
3815  {
3816  while (neighbor_child->has_children())
3817  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3818 
3819  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3820  ExcInternalError());
3821  }
3822  active_neighbors.push_back(neighbor_child);
3823  }
3824  else
3825  {
3826  if (cell->face(n)->has_children())
3827  // this neighbor has children. find
3828  // out which border to the present
3829  // cell
3830  for (unsigned int c = 0;
3831  c < cell->face(n)->n_active_descendants();
3832  ++c)
3833  active_neighbors.push_back(
3834  cell->neighbor_child_on_subface(n, c));
3835  else
3836  {
3837  // the neighbor must be active
3838  // himself
3839  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3840  active_neighbors.push_back(cell->neighbor(n));
3841  }
3842  }
3843  }
3844  }
3845 
3846 
3847 
3848  namespace internal
3849  {
3850  namespace ProjectToObject
3851  {
3864  struct CrossDerivative
3865  {
3866  const unsigned int direction_0;
3867  const unsigned int direction_1;
3868 
3869  CrossDerivative(const unsigned int d0, const unsigned int d1);
3870  };
3871 
3872  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3873  const unsigned int d1)
3874  : direction_0(d0)
3875  , direction_1(d1)
3876  {}
3877 
3878 
3879 
3884  template <typename F>
3885  inline auto
3886  centered_first_difference(const double center,
3887  const double step,
3888  const F &f) -> decltype(f(center) - f(center))
3889  {
3890  return (f(center + step) - f(center - step)) / (2.0 * step);
3891  }
3892 
3893 
3894 
3899  template <typename F>
3900  inline auto
3901  centered_second_difference(const double center,
3902  const double step,
3903  const F &f) -> decltype(f(center) - f(center))
3904  {
3905  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3906  (step * step);
3907  }
3908 
3909 
3910 
3920  template <int structdim, typename F>
3921  inline auto
3922  cross_stencil(
3923  const CrossDerivative cross_derivative,
3925  const double step,
3926  const F &f) -> decltype(f(center) - f(center))
3927  {
3929  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3930  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3931  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3932  1.0 / 3.0 * f(center - simplex_vector) +
3933  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3934  step;
3935  }
3936 
3937 
3938 
3945  template <int spacedim, int structdim, typename F>
3946  inline double
3947  gradient_entry(
3948  const unsigned int row_n,
3949  const unsigned int dependent_direction,
3950  const Point<spacedim> &p0,
3952  const double step,
3953  const F & f)
3954  {
3956  dependent_direction <
3958  ExcMessage("This function assumes that the last weight is a "
3959  "dependent variable (and hence we cannot take its "
3960  "derivative directly)."));
3961  Assert(row_n != dependent_direction,
3962  ExcMessage(
3963  "We cannot differentiate with respect to the variable "
3964  "that is assumed to be dependent."));
3965 
3966  const Point<spacedim> manifold_point = f(center);
3967  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3968  {row_n, dependent_direction}, center, step, f);
3969  double entry = 0.0;
3970  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3971  entry +=
3972  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3973  return entry;
3974  }
3975 
3981  template <typename Iterator, int spacedim, int structdim>
3983  project_to_d_linear_object(const Iterator & object,
3984  const Point<spacedim> &trial_point)
3985  {
3986  // let's look at this for simplicity for a quad (structdim==2) in a
3987  // space with spacedim>2 (notate trial_point by y): all points on the
3988  // surface are given by
3989  // x(\xi) = sum_i v_i phi_x(\xi)
3990  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3991  // reference coordinates of the quad. so what we are trying to do is
3992  // find a point x on the surface that is closest to the point y. there
3993  // are different ways to solve this problem, but in the end it's a
3994  // nonlinear problem and we have to find reference coordinates \xi so
3995  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3996  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3997  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3998  // have to use a Newton method to find the answer. This leads to the
3999  // following formulation of Newton steps:
4000  //
4001  // Given \xi_k, find \delta\xi_k so that
4002  // H_k \delta\xi_k = - F_k
4003  // where H_k is an approximation to the second derivatives of J at
4004  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
4005  // number of times until the right hand side is small enough. As a
4006  // stopping criterion, we terminate if ||\delta\xi||<eps.
4007  //
4008  // As for the Hessian, the best choice would be
4009  // H_k = J''(\xi_k)
4010  // but we'll opt for the simpler Gauss-Newton form
4011  // H_k = A^T A
4012  // i.e.
4013  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
4014  // \partial_n phi_i *
4015  // \partial_m phi_j
4016  // we start at xi=(0.5, 0.5).
4017  Point<structdim> xi;
4018  for (unsigned int d = 0; d < structdim; ++d)
4019  xi[d] = 0.5;
4020 
4021  Point<spacedim> x_k;
4022  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4023  x_k += object->vertex(i) *
4025 
4026  do
4027  {
4029  for (const unsigned int i :
4031  F_k +=
4032  (x_k - trial_point) * object->vertex(i) *
4034  i);
4035 
4037  for (const unsigned int i :
4039  for (const unsigned int j :
4041  {
4044  xi, i),
4046  xi, j));
4047  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
4048  }
4049 
4050  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
4051  xi += delta_xi;
4052 
4053  x_k = Point<spacedim>();
4054  for (const unsigned int i :
4056  x_k += object->vertex(i) *
4058 
4059  if (delta_xi.norm() < 1e-7)
4060  break;
4061  }
4062  while (true);
4063 
4064  return x_k;
4065  }
4066  } // namespace ProjectToObject
4067  } // namespace internal
4068 
4069 
4070 
4071  namespace internal
4072  {
4073  // We hit an internal compiler error in ICC 15 if we define this as a lambda
4074  // inside the project_to_object function below.
4075  template <int structdim>
4076  inline bool
4077  weights_are_ok(
4079  {
4080  // clang has trouble figuring out structdim here, so define it
4081  // again:
4082  static const std::size_t n_vertices_per_cell =
4084  n_independent_components;
4085  std::array<double, n_vertices_per_cell> copied_weights;
4086  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4087  {
4088  copied_weights[i] = v[i];
4089  if (v[i] < 0.0 || v[i] > 1.0)
4090  return false;
4091  }
4092 
4093  // check the sum: try to avoid some roundoff errors by summing in order
4094  std::sort(copied_weights.begin(), copied_weights.end());
4095  const double sum =
4096  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4097  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4098  }
4099  } // namespace internal
4100 
4101  template <typename Iterator>
4104  const Iterator & object,
4106  {
4107  const int spacedim = Iterator::AccessorType::space_dimension;
4108  const int structdim = Iterator::AccessorType::structure_dimension;
4109 
4110  Point<spacedim> projected_point = trial_point;
4111 
4112  if (structdim >= spacedim)
4113  return projected_point;
4114  else if (structdim == 1 || structdim == 2)
4115  {
4116  using namespace internal::ProjectToObject;
4117  // Try to use the special flat algorithm for quads (this is better
4118  // than the general algorithm in 3D). This does not take into account
4119  // whether projected_point is outside the quad, but we optimize along
4120  // lines below anyway:
4121  const int dim = Iterator::AccessorType::dimension;
4122  const Manifold<dim, spacedim> &manifold = object->get_manifold();
4123  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4124  &manifold) != nullptr)
4125  {
4126  projected_point =
4127  project_to_d_linear_object<Iterator, spacedim, structdim>(
4128  object, trial_point);
4129  }
4130  else
4131  {
4132  // We want to find a point on the convex hull (defined by the
4133  // vertices of the object and the manifold description) that is
4134  // relatively close to the trial point. This has a few issues:
4135  //
4136  // 1. For a general convex hull we are not guaranteed that a unique
4137  // minimum exists.
4138  // 2. The independent variables in the optimization process are the
4139  // weights given to Manifold::get_new_point, which must sum to 1,
4140  // so we cannot use standard finite differences to approximate a
4141  // gradient.
4142  //
4143  // There is not much we can do about 1., but for 2. we can derive
4144  // finite difference stencils that work on a structdim-dimensional
4145  // simplex and rewrite the optimization problem to use those
4146  // instead. Consider the structdim 2 case and let
4147  //
4148  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4149  // c2, c3})
4150  //
4151  // where {c0, c1, c2, c3} are the weights for the four vertices on
4152  // the quadrilateral. We seek to minimize the Euclidean distance
4153  // between F(...) and trial_point. We can solve for c3 in terms of
4154  // the other weights and get, for one coordinate direction
4155  //
4156  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4157  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4158  //
4159  // where we substitute back in for c3 after taking the
4160  // derivative. We can compute a stencil for the cross derivative
4161  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4162  // (and gradient_entry computes the sum over the independent
4163  // variables). Below, we somewhat arbitrarily pick the last
4164  // component as the dependent one.
4165  //
4166  // Since we can now calculate derivatives of the objective
4167  // function we can use gradient descent to minimize it.
4168  //
4169  // Of course, this is much simpler in the structdim = 1 case (we
4170  // could rewrite the projection as a 1D optimization problem), but
4171  // to reduce the potential for bugs we use the same code in both
4172  // cases.
4173  const double step_size = object->diameter() / 64.0;
4174 
4175  constexpr unsigned int n_vertices_per_cell =
4177 
4178  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4179  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4180  ++vertex_n)
4181  vertices[vertex_n] = object->vertex(vertex_n);
4182 
4183  auto get_point_from_weights =
4184  [&](const Tensor<1, n_vertices_per_cell> &weights)
4185  -> Point<spacedim> {
4186  return object->get_manifold().get_new_point(
4187  make_array_view(vertices.begin(), vertices.end()),
4188  make_array_view(weights.begin_raw(), weights.end_raw()));
4189  };
4190 
4191  // pick the initial weights as (normalized) inverse distances from
4192  // the trial point:
4193  Tensor<1, n_vertices_per_cell> guess_weights;
4194  double guess_weights_sum = 0.0;
4195  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4196  ++vertex_n)
4197  {
4198  const double distance =
4199  vertices[vertex_n].distance(trial_point);
4200  if (distance == 0.0)
4201  {
4202  guess_weights = 0.0;
4203  guess_weights[vertex_n] = 1.0;
4204  guess_weights_sum = 1.0;
4205  break;
4206  }
4207  else
4208  {
4209  guess_weights[vertex_n] = 1.0 / distance;
4210  guess_weights_sum += guess_weights[vertex_n];
4211  }
4212  }
4213  guess_weights /= guess_weights_sum;
4214  Assert(internal::weights_are_ok<structdim>(guess_weights),
4215  ExcInternalError());
4216 
4217  // The optimization algorithm consists of two parts:
4218  //
4219  // 1. An outer loop where we apply the gradient descent algorithm.
4220  // 2. An inner loop where we do a line search to find the optimal
4221  // length of the step one should take in the gradient direction.
4222  //
4223  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4224  {
4225  const unsigned int dependent_direction =
4226  n_vertices_per_cell - 1;
4227  Tensor<1, n_vertices_per_cell> current_gradient;
4228  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4229  ++row_n)
4230  {
4231  if (row_n != dependent_direction)
4232  {
4233  current_gradient[row_n] =
4234  gradient_entry<spacedim, structdim>(
4235  row_n,
4236  dependent_direction,
4237  trial_point,
4238  guess_weights,
4239  step_size,
4240  get_point_from_weights);
4241 
4242  current_gradient[dependent_direction] -=
4243  current_gradient[row_n];
4244  }
4245  }
4246 
4247  // We need to travel in the -gradient direction, as noted
4248  // above, but we may not want to take a full step in that
4249  // direction; instead, guess that we will go -0.5*gradient and
4250  // do quasi-Newton iteration to pick the best multiplier. The
4251  // goal is to find a scalar alpha such that
4252  //
4253  // F(x - alpha g)
4254  //
4255  // is minimized, where g is the gradient and F is the
4256  // objective function. To find the optimal value we find roots
4257  // of the derivative of the objective function with respect to
4258  // alpha by Newton iteration, where we approximate the first
4259  // and second derivatives of F(x - alpha g) with centered
4260  // finite differences.
4261  double gradient_weight = -0.5;
4262  auto gradient_weight_objective_function =
4263  [&](const double gradient_weight_guess) -> double {
4264  return (trial_point -
4265  get_point_from_weights(guess_weights +
4266  gradient_weight_guess *
4267  current_gradient))
4268  .norm_square();
4269  };
4270 
4271  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4272  {
4273  const double update_numerator = centered_first_difference(
4274  gradient_weight,
4275  step_size,
4276  gradient_weight_objective_function);
4277  const double update_denominator =
4278  centered_second_difference(
4279  gradient_weight,
4280  step_size,
4281  gradient_weight_objective_function);
4282 
4283  // avoid division by zero. Note that we limit the gradient
4284  // weight below
4285  if (std::abs(update_denominator) == 0.0)
4286  break;
4287  gradient_weight =
4288  gradient_weight - update_numerator / update_denominator;
4289 
4290  // Put a fairly lenient bound on the largest possible
4291  // gradient (things tend to be locally flat, so the gradient
4292  // itself is usually small)
4293  if (std::abs(gradient_weight) > 10)
4294  {
4295  gradient_weight = -10.0;
4296  break;
4297  }
4298  }
4299 
4300  // It only makes sense to take convex combinations with weights
4301  // between zero and one. If the update takes us outside of this
4302  // region then rescale the update to stay within the region and
4303  // try again
4304  Tensor<1, n_vertices_per_cell> tentative_weights =
4305  guess_weights + gradient_weight * current_gradient;
4306 
4307  double new_gradient_weight = gradient_weight;
4308  for (unsigned int iteration_count = 0; iteration_count < 40;
4309  ++iteration_count)
4310  {
4311  if (internal::weights_are_ok<structdim>(tentative_weights))
4312  break;
4313 
4314  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4315  {
4316  if (tentative_weights[i] < 0.0)
4317  {
4318  tentative_weights -=
4319  (tentative_weights[i] / current_gradient[i]) *
4320  current_gradient;
4321  }
4322  if (tentative_weights[i] < 0.0 ||
4323  1.0 < tentative_weights[i])
4324  {
4325  new_gradient_weight /= 2.0;
4326  tentative_weights =
4327  guess_weights +
4328  new_gradient_weight * current_gradient;
4329  }
4330  }
4331  }
4332 
4333  // the update might still send us outside the valid region, so
4334  // check again and quit if the update is still not valid
4335  if (!internal::weights_are_ok<structdim>(tentative_weights))
4336  break;
4337 
4338  // if we cannot get closer by traveling in the gradient
4339  // direction then quit
4340  if (get_point_from_weights(tentative_weights)
4341  .distance(trial_point) <
4342  get_point_from_weights(guess_weights).distance(trial_point))
4343  guess_weights = tentative_weights;
4344  else
4345  break;
4346  Assert(internal::weights_are_ok<structdim>(guess_weights),
4347  ExcInternalError());
4348  }
4349  Assert(internal::weights_are_ok<structdim>(guess_weights),
4350  ExcInternalError());
4351  projected_point = get_point_from_weights(guess_weights);
4352  }
4353 
4354  // if structdim == 2 and the optimal point is not on the interior then
4355  // we may be able to get a more accurate result by projecting onto the
4356  // lines.
4357  if (structdim == 2)
4358  {
4359  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4360  line_projections;
4361  for (unsigned int line_n = 0;
4362  line_n < GeometryInfo<structdim>::lines_per_cell;
4363  ++line_n)
4364  {
4365  line_projections[line_n] =
4366  project_to_object(object->line(line_n), trial_point);
4367  }
4368  std::sort(line_projections.begin(),
4369  line_projections.end(),
4370  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4371  return a.distance(trial_point) <
4372  b.distance(trial_point);
4373  });
4374  if (line_projections[0].distance(trial_point) <
4375  projected_point.distance(trial_point))
4376  projected_point = line_projections[0];
4377  }
4378  }
4379  else
4380  {
4381  Assert(false, ExcNotImplemented());
4382  return projected_point;
4383  }
4384 
4385  return projected_point;
4386  }
4387 
4388 
4389 
4390  namespace internal
4391  {
4392  template <typename DataType,
4393  typename MeshType,
4394  typename MeshCellIteratorType>
4395  inline void
4396  exchange_cell_data(
4397  const MeshType &mesh,
4398  const std::function<
4399  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4400  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4401  & unpack,
4402  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4403  const std::function<void(
4404  const std::function<void(const MeshCellIteratorType &,
4405  const types::subdomain_id)> &)> &process_cells,
4406  const std::function<std::set<types::subdomain_id>(
4407  const parallel::TriangulationBase<MeshType::dimension,
4408  MeshType::space_dimension> &)>
4409  &compute_ghost_owners)
4410  {
4411 # ifndef DEAL_II_WITH_MPI
4412  (void)mesh;
4413  (void)pack;
4414  (void)unpack;
4415  (void)cell_filter;
4416  (void)process_cells;
4417  (void)compute_ghost_owners;
4418  Assert(false, ExcNeedsMPI());
4419 # else
4420  constexpr int dim = MeshType::dimension;
4421  constexpr int spacedim = MeshType::space_dimension;
4422  auto tria =
4423  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4424  &mesh.get_triangulation());
4425  Assert(
4426  tria != nullptr,
4427  ExcMessage(
4428  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4429 
4430  if (const auto tria = dynamic_cast<
4432  &mesh.get_triangulation()))
4433  {
4434  Assert(
4435  tria->with_artificial_cells(),
4436  ExcMessage(
4437  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4438  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4439  "operate on a single layer of ghost cells. However, you have "
4440  "given a Triangulation object of type "
4441  "parallel::shared::Triangulation without artificial cells "
4442  "resulting in arbitrary numbers of ghost layers."));
4443  }
4444 
4445  // build list of cells to request for each neighbor
4446  std::set<::types::subdomain_id> ghost_owners =
4447  compute_ghost_owners(*tria);
4448  std::map<::types::subdomain_id,
4449  std::vector<typename CellId::binary_type>>
4450  neighbor_cell_list;
4451 
4452  for (const auto ghost_owner : ghost_owners)
4453  neighbor_cell_list[ghost_owner] = {};
4454 
4455  process_cells([&](const auto &cell, const auto key) -> void {
4456  if (cell_filter(cell))
4457  {
4458  constexpr int spacedim = MeshType::space_dimension;
4459  neighbor_cell_list[key].emplace_back(
4460  cell->id().template to_binary<spacedim>());
4461  }
4462  });
4463 
4464  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4465  ExcInternalError());
4466 
4467 
4468  // Before sending & receiving, make sure we protect this section with
4469  // a mutex:
4470  static Utilities::MPI::CollectiveMutex mutex;
4472  mutex, tria->get_communicator());
4473 
4474  const int mpi_tag =
4476  const int mpi_tag_reply =
4478 
4479  // send our requests
4480  std::vector<MPI_Request> requests(ghost_owners.size());
4481  {
4482  unsigned int idx = 0;
4483  for (const auto &it : neighbor_cell_list)
4484  {
4485  // send the data about the relevant cells
4486  const int ierr = MPI_Isend(it.second.data(),
4487  it.second.size() * sizeof(it.second[0]),
4488  MPI_BYTE,
4489  it.first,
4490  mpi_tag,
4492  &requests[idx]);
4493  AssertThrowMPI(ierr);
4494  ++idx;
4495  }
4496  }
4497 
4498  // receive requests and reply with the results
4499  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4500  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4501 
4502  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4503  {
4504  MPI_Status status;
4505  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4506  mpi_tag,
4508  &status);
4509  AssertThrowMPI(ierr);
4510 
4511  int len;
4512  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4513  AssertThrowMPI(ierr);
4514  Assert(len % sizeof(typename CellId::binary_type) == 0,
4515  ExcInternalError());
4516 
4517  const unsigned int n_cells =
4518  len / sizeof(typename CellId::binary_type);
4519  std::vector<typename CellId::binary_type> cells_with_requests(
4520  n_cells);
4521  std::vector<DataType> data_to_send;
4522  data_to_send.reserve(n_cells);
4523  std::vector<bool> cell_carries_data(n_cells, false);
4524 
4525  ierr = MPI_Recv(cells_with_requests.data(),
4526  len,
4527  MPI_BYTE,
4528  status.MPI_SOURCE,
4529  status.MPI_TAG,
4531  &status);
4532  AssertThrowMPI(ierr);
4533 
4534  // store data for each cell
4535  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4536  {
4537  const auto cell =
4538  tria->create_cell_iterator(CellId(cells_with_requests[c]));
4539 
4540  MeshCellIteratorType mesh_it(tria,
4541  cell->level(),
4542  cell->index(),
4543  &mesh);
4544 
4545  const std_cxx17::optional<DataType> data = pack(mesh_it);
4546  if (data)
4547  {
4548  data_to_send.emplace_back(*data);
4549  cell_carries_data[c] = true;
4550  }
4551  }
4552 
4553  // collect data for sending the reply in a buffer
4554 
4555  // (a) make room for storing the local offsets in case we receive
4556  // other data
4557  sendbuffers[idx].resize(sizeof(std::size_t));
4558 
4559  // (b) append the actual data and store how much memory it
4560  // corresponds to, which we then insert into the leading position of
4561  // the sendbuffer
4562  std::size_t size_of_send =
4563  Utilities::pack(data_to_send,
4564  sendbuffers[idx],
4565  /*enable_compression*/ false);
4566  std::memcpy(sendbuffers[idx].data(),
4567  &size_of_send,
4568  sizeof(std::size_t));
4569 
4570  // (c) append information of certain cells that got left out in case
4571  // we need it
4572  if (data_to_send.size() < n_cells)
4573  Utilities::pack(cell_carries_data,
4574  sendbuffers[idx],
4575  /*enable_compression*/ false);
4576 
4577  // send data
4578  ierr = MPI_Isend(sendbuffers[idx].data(),
4579  sendbuffers[idx].size(),
4580  MPI_BYTE,
4581  status.MPI_SOURCE,
4582  mpi_tag_reply,
4584  &reply_requests[idx]);
4585  AssertThrowMPI(ierr);
4586  }
4587 
4588  // finally receive the replies
4589  std::vector<char> receive;
4590  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4591  {
4592  MPI_Status status;
4593  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4594  mpi_tag_reply,
4596  &status);
4597  AssertThrowMPI(ierr);
4598 
4599  int len;
4600  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4601  AssertThrowMPI(ierr);
4602 
4603  receive.resize(len);
4604 
4605  ierr = MPI_Recv(receive.data(),
4606  len,
4607  MPI_BYTE,
4608  status.MPI_SOURCE,
4609  status.MPI_TAG,
4611  &status);
4612  AssertThrowMPI(ierr);
4613 
4614  // (a) first determine the length of the data section in the
4615  // received buffer
4616  auto data_iterator = receive.begin();
4617  std::size_t size_of_received_data =
4618  Utilities::unpack<std::size_t>(data_iterator,
4619  data_iterator + sizeof(std::size_t));
4620  data_iterator += sizeof(std::size_t);
4621 
4622  // (b) unpack the data section in the indicated region
4623  auto received_data = Utilities::unpack<std::vector<DataType>>(
4624  data_iterator,
4625  data_iterator + size_of_received_data,
4626  /*enable_compression*/ false);
4627  data_iterator += size_of_received_data;
4628 
4629  // (c) check if the received data contained fewer entries than the
4630  // number of cells we identified in the beginning, in which case we
4631  // need to extract the boolean vector with the relevant information
4632  const std::vector<typename CellId::binary_type> &this_cell_list =
4633  neighbor_cell_list[status.MPI_SOURCE];
4634  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4635  std::vector<bool> cells_with_data;
4636  if (received_data.size() < this_cell_list.size())
4637  {
4638  cells_with_data = Utilities::unpack<std::vector<bool>>(
4639  data_iterator, receive.end(), /*enable_compression*/ false);
4640  AssertDimension(cells_with_data.size(), this_cell_list.size());
4641  }
4642 
4643  // (d) go through the received data and call the user-provided
4644  // unpack function
4645  auto received_data_iterator = received_data.begin();
4646  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4647  if (cells_with_data.empty() || cells_with_data[c])
4648  {
4650  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4651 
4652  MeshCellIteratorType cell(tria,
4653  tria_cell->level(),
4654  tria_cell->index(),
4655  &mesh);
4656 
4657  unpack(cell, *received_data_iterator);
4658  ++received_data_iterator;
4659  }
4660  }
4661 
4662  // make sure that all communication is finished
4663  // when we leave this function.
4664  if (requests.size() > 0)
4665  {
4666  const int ierr =
4667  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4668  AssertThrowMPI(ierr);
4669  }
4670  if (reply_requests.size() > 0)
4671  {
4672  const int ierr = MPI_Waitall(reply_requests.size(),
4673  reply_requests.data(),
4674  MPI_STATUSES_IGNORE);
4675  AssertThrowMPI(ierr);
4676  }
4677 
4678 
4679 # endif // DEAL_II_WITH_MPI
4680  }
4681 
4682  } // namespace internal
4683 
4684  template <typename DataType, typename MeshType>
4685  inline void
4687  const MeshType & mesh,
4688  const std::function<std_cxx17::optional<DataType>(
4689  const typename MeshType::active_cell_iterator &)> &pack,
4690  const std::function<void(const typename MeshType::active_cell_iterator &,
4691  const DataType &)> & unpack,
4692  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4693  &cell_filter)
4694  {
4695 # ifndef DEAL_II_WITH_MPI
4696  (void)mesh;
4697  (void)pack;
4698  (void)unpack;
4699  (void)cell_filter;
4700  Assert(false, ExcNeedsMPI());
4701 # else
4702  internal::exchange_cell_data<DataType,
4703  MeshType,
4704  typename MeshType::active_cell_iterator>(
4705  mesh,
4706  pack,
4707  unpack,
4708  cell_filter,
4709  [&](const auto &process) {
4710  for (const auto &cell : mesh.active_cell_iterators())
4711  if (cell->is_ghost())
4712  process(cell, cell->subdomain_id());
4713  },
4714  [](const auto &tria) { return tria.ghost_owners(); });
4715 # endif
4716  }
4717 
4718 
4719 
4720  template <typename DataType, typename MeshType>
4721  inline void
4723  const MeshType & mesh,
4724  const std::function<std_cxx17::optional<DataType>(
4725  const typename MeshType::level_cell_iterator &)> &pack,
4726  const std::function<void(const typename MeshType::level_cell_iterator &,
4727  const DataType &)> & unpack,
4728  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4729  &cell_filter)
4730  {
4731 # ifndef DEAL_II_WITH_MPI
4732  (void)mesh;
4733  (void)pack;
4734  (void)unpack;
4735  (void)cell_filter;
4736  Assert(false, ExcNeedsMPI());
4737 # else
4738  internal::exchange_cell_data<DataType,
4739  MeshType,
4740  typename MeshType::level_cell_iterator>(
4741  mesh,
4742  pack,
4743  unpack,
4744  cell_filter,
4745  [&](const auto &process) {
4746  for (const auto &cell : mesh.cell_iterators())
4747  if (cell->is_ghost_on_level())
4748  process(cell, cell->level_subdomain_id());
4749  },
4750  [](const auto &tria) { return tria.level_ghost_owners(); });
4751 # endif
4752  }
4753 } // namespace GridTools
4754 
4755 #endif // DOXYGEN
4756 
4758 
4759 #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
typename VectorType::value_type value_type
Definition: grid_tools.h:3442
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:6835
void process_cell(const typename DoFHandler< dim >::active_cell_iterator &cell, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:6871
void process_sub_cell(const std::vector< value_type > &, const std::vector< Point< 1 >> &, const std::vector< unsigned int > &, const double, std::vector< Point< 1 >> &, std::vector< CellData< 1 >> &, const bool) const
Definition: grid_tools.h:3528
const unsigned int n_subdivisions
Definition: grid_tools.h:3570
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:6783
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:6800
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:103
#define DEAL_II_DEPRECATED
Definition: config.h:164
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:512
Point< 3 > center
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4608
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1341
__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:1501
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1818
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#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)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:437
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:5020
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:5112
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:5045
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:5140
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:6025
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3827
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:290
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:3189
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2132
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4167
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:6502
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2282
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:4978
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:6298
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:984
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:3107
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:742
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:6279
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:4458
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:5322
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5293
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:4485
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2098
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:5667
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5260
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:4273
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:2793
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6625
std::size_t invert_cells_with_negative_measure(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:904
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3336
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4442
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:2007
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:5638
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)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:3926
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:314
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3434
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:140
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:637
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:380
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:6343
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3861
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:4372
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:2314
void build_triangulation_from_patch(const std::vector< typename Container::active_cell_iterator > &patch, Triangulation< Container::dimension, Container::space_dimension > &local_triangulation, std::map< typename Triangulation< Container::dimension, Container::space_dimension >::active_cell_iterator, typename Container::active_cell_iterator > &patch_to_global_tria_map)
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:6439
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:84
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3890
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:2665
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:3483
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:396
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:4512
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5228
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:536
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:5835
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
@ exchange_cell_data_request
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
@ exchange_cell_data_reply
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
T sum(const T &t, const MPI_Comm &mpi_communicator)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1334
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1499
Definition: hp.h:118
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13750
static const unsigned int invalid_unsigned_int
Definition: types.h:206
const types::manifold_id flat_manifold_id
Definition: types.h:279
unsigned int manifold_id
Definition: types.h:146
unsigned int global_dof_index
Definition: types.h:81
unsigned int subdomain_id
Definition: types.h:43
unsigned int boundary_id
Definition: types.h:134
::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:3382
void serialize(Archive &archive, const unsigned int version)
void save(Archive &ar, const unsigned int) const
Definition: grid_tools.h:3356
std::vector< CellId > cell_ids
Definition: grid_tools.h:3339
std::bitset< 3 > orientation
Definition: grid_tools.h:2766
FullMatrix< double > matrix
Definition: grid_tools.h:2780
unsigned int face_idx[2]
Definition: grid_tools.h:2759
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1193
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1168
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