Reference documentation for deal.II version GIT b9c422a7ef 2022-07-03 21:00: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 
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 
470  template <int dim, int spacedim>
471  void
473  const std::vector<Point<spacedim>> &all_vertices,
474  std::vector<CellData<dim>> & cells);
475 
485  template <int dim, int spacedim>
486  std::size_t
488  const std::vector<Point<spacedim>> &all_vertices,
489  std::vector<CellData<dim>> & cells);
490 
501  template <int dim>
502  void
503  consistently_order_cells(std::vector<CellData<dim>> &cells);
504 
510 
573  template <int dim, typename Transformation, int spacedim>
574  void
575  transform(const Transformation & transformation,
577 
584  template <int dim, int spacedim>
585  void
586  shift(const Tensor<1, spacedim> & shift_vector,
588 
589 
600  template <int dim>
601  void
603 
616  template <int dim>
617  void
618  rotate(const Tensor<1, 3, double> &axis,
619  const double angle,
621 
637  template <int dim>
638  DEAL_II_DEPRECATED void
639  rotate(const double angle,
640  const unsigned int axis,
642 
700  template <int dim>
701  void
702  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
704  const Function<dim, double> *coefficient = nullptr,
705  const bool solve_for_absolute_positions = false);
706 
712  template <int dim, int spacedim>
713  std::map<unsigned int, Point<spacedim>>
715 
723  template <int dim, int spacedim>
724  void
725  scale(const double scaling_factor,
727 
742  template <int dim, int spacedim>
743  void
745  const double factor,
747  const bool keep_boundary = true,
748  const unsigned int seed = boost::random::mt19937::default_seed);
749 
783  template <int dim, int spacedim>
784  void
786  const bool isotropic = false,
787  const unsigned int max_iterations = 100);
788 
813  template <int dim, int spacedim>
814  void
816  const double max_ratio = 1.6180339887,
817  const unsigned int max_iterations = 5);
818 
908  template <int dim, int spacedim>
909  void
911  const double limit_angle_fraction = .75);
912 
918 
973  template <int dim, int spacedim>
974 #ifndef DOXYGEN
975  std::tuple<
976  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
977  std::vector<std::vector<Point<dim>>>,
978  std::vector<std::vector<unsigned int>>>
979 #else
980  return_type
981 #endif
983  const Cache<dim, spacedim> & cache,
984  const std::vector<Point<spacedim>> &points,
986  &cell_hint =
988 
1022  template <int dim, int spacedim>
1023 #ifndef DOXYGEN
1024  std::tuple<
1025  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1026  std::vector<std::vector<Point<dim>>>,
1027  std::vector<std::vector<unsigned int>>,
1028  std::vector<unsigned int>>
1029 #else
1030  return_type
1031 #endif
1033  const Cache<dim, spacedim> & cache,
1034  const std::vector<Point<spacedim>> &points,
1036  &cell_hint =
1038 
1108  template <int dim, int spacedim>
1109 #ifndef DOXYGEN
1110  std::tuple<
1111  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1112  std::vector<std::vector<Point<dim>>>,
1113  std::vector<std::vector<unsigned int>>,
1114  std::vector<std::vector<Point<spacedim>>>,
1115  std::vector<std::vector<unsigned int>>>
1116 #else
1117  return_type
1118 #endif
1120  const GridTools::Cache<dim, spacedim> & cache,
1121  const std::vector<Point<spacedim>> & local_points,
1122  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1123  const double tolerance = 1e-10);
1124 
1125  namespace internal
1126  {
1141  template <int dim, int spacedim>
1143  {
1150  std::vector<std::tuple<std::pair<int, int>,
1151  unsigned int,
1152  unsigned int,
1153  Point<dim>,
1155  unsigned int>>
1157 
1161  std::vector<unsigned int> send_ranks;
1162 
1168  std::vector<unsigned int> send_ptrs;
1169 
1180  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1182 
1186  std::vector<unsigned int> recv_ranks;
1187 
1193  std::vector<unsigned int> recv_ptrs;
1194  };
1195 
1205  template <int dim, int spacedim>
1208  const GridTools::Cache<dim, spacedim> & cache,
1209  const std::vector<Point<spacedim>> & points,
1210  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1211  const std::vector<bool> & marked_vertices,
1212  const double tolerance,
1213  const bool perform_handshake,
1214  const bool enforce_unique_mapping = false);
1215 
1216  } // namespace internal
1217 
1254  template <int dim, int spacedim>
1255  std::map<unsigned int, Point<spacedim>>
1257  const Triangulation<dim, spacedim> &container,
1258  const Mapping<dim, spacedim> & mapping =
1259  (ReferenceCells::get_hypercube<dim>()
1260 #ifndef _MSC_VER
1261  .template get_default_linear_mapping<dim, spacedim>()
1262 #else
1263  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1264 #endif
1265  ));
1266 
1276  template <int spacedim>
1277  unsigned int
1278  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1279  const Point<spacedim> & p);
1280 
1304  template <int dim, template <int, int> class MeshType, int spacedim>
1305  unsigned int
1306  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1307  const Point<spacedim> & p,
1308  const std::vector<bool> & marked_vertices = {});
1309 
1333  template <int dim, template <int, int> class MeshType, int spacedim>
1334  unsigned int
1336  const MeshType<dim, spacedim> &mesh,
1337  const Point<spacedim> & p,
1338  const std::vector<bool> & marked_vertices = {});
1339 
1340 
1360  template <int dim, template <int, int> class MeshType, int spacedim>
1361 #ifndef _MSC_VER
1362  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1363 #else
1364  std::vector<
1365  typename ::internal::
1366  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1367 #endif
1368  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1369  const unsigned int vertex_index);
1370 
1433  template <int dim, template <int, int> class MeshType, int spacedim>
1434 #ifndef _MSC_VER
1435  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1436 #else
1437  std::pair<typename ::internal::
1438  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1439  Point<dim>>
1440 #endif
1442  const MeshType<dim, spacedim> &mesh,
1443  const Point<spacedim> & p,
1444  const std::vector<bool> &marked_vertices = {},
1445  const double tolerance = 1.e-10);
1446 
1454  template <int dim, template <int, int> class MeshType, int spacedim>
1455 #ifndef _MSC_VER
1456  typename MeshType<dim, spacedim>::active_cell_iterator
1457 #else
1458  typename ::internal::
1459  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1460 #endif
1461  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1462  const Point<spacedim> & p,
1463  const std::vector<bool> &marked_vertices = {},
1464  const double tolerance = 1.e-10);
1465 
1472  template <int dim, int spacedim>
1473  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1474  Point<dim>>
1476  const hp::MappingCollection<dim, spacedim> &mapping,
1477  const DoFHandler<dim, spacedim> & mesh,
1478  const Point<spacedim> & p,
1479  const double tolerance = 1.e-10);
1480 
1532  template <int dim, int spacedim>
1533  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1534  Point<dim>>
1536  const Cache<dim, spacedim> &cache,
1537  const Point<spacedim> & p,
1540  const std::vector<bool> &marked_vertices = {},
1541  const double tolerance = 1.e-10);
1542 
1556  template <int dim, template <int, int> class MeshType, int spacedim>
1557 #ifndef _MSC_VER
1558  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1559 #else
1560  std::pair<typename ::internal::
1561  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1562  Point<dim>>
1563 #endif
1565  const Mapping<dim, spacedim> & mapping,
1566  const MeshType<dim, spacedim> &mesh,
1567  const Point<spacedim> & p,
1568  const std::vector<
1569  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1571  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1572  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1573  typename MeshType<dim, spacedim>::active_cell_iterator(),
1574  const std::vector<bool> & marked_vertices = {},
1575  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1576  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1577  const double tolerance = 1.e-10,
1578  const RTree<
1579  std::pair<BoundingBox<spacedim>,
1581  *relevant_cell_bounding_boxes_rtree = nullptr);
1582 
1604  template <int dim, template <int, int> class MeshType, int spacedim>
1605 #ifndef _MSC_VER
1606  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1607  Point<dim>>>
1608 #else
1609  std::vector<std::pair<
1610  typename ::internal::
1611  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1612  Point<dim>>>
1613 #endif
1615  const Mapping<dim, spacedim> & mapping,
1616  const MeshType<dim, spacedim> &mesh,
1617  const Point<spacedim> & p,
1618  const double tolerance,
1619  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1620  Point<dim>> & first_cell);
1621 
1628  template <int dim, template <int, int> class MeshType, int spacedim>
1629 #ifndef _MSC_VER
1630  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1631  Point<dim>>>
1632 #else
1633  std::vector<std::pair<
1634  typename ::internal::
1635  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1636  Point<dim>>>
1637 #endif
1639  const Mapping<dim, spacedim> & mapping,
1640  const MeshType<dim, spacedim> &mesh,
1641  const Point<spacedim> & p,
1642  const double tolerance = 1e-10,
1643  const std::vector<bool> & marked_vertices = {});
1644 
1666  template <class MeshType>
1667  std::vector<typename MeshType::active_cell_iterator>
1668  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1669 
1694  template <class MeshType>
1695  void
1697  const typename MeshType::active_cell_iterator & cell,
1698  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1699 
1749  template <class MeshType>
1750  std::vector<typename MeshType::active_cell_iterator>
1752  const MeshType &mesh,
1753  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1754  &predicate);
1755 
1756 
1764  template <class MeshType>
1765  std::vector<typename MeshType::cell_iterator>
1767  const MeshType &mesh,
1768  const std::function<bool(const typename MeshType::cell_iterator &)>
1769  & predicate,
1770  const unsigned int level);
1771 
1772 
1785  template <class MeshType>
1786  std::vector<typename MeshType::active_cell_iterator>
1787  compute_ghost_cell_halo_layer(const MeshType &mesh);
1788 
1837  template <class MeshType>
1838  std::vector<typename MeshType::active_cell_iterator>
1840  const MeshType &mesh,
1841  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1842  & predicate,
1843  const double layer_thickness);
1844 
1867  template <class MeshType>
1868  std::vector<typename MeshType::active_cell_iterator>
1869  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1870  const double layer_thickness);
1871 
1887  template <class MeshType>
1888  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1890  const MeshType &mesh,
1891  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1892  &predicate);
1893 
1946  template <class MeshType>
1947  std::vector<BoundingBox<MeshType::space_dimension>>
1949  const MeshType &mesh,
1950  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1951  & predicate,
1952  const unsigned int refinement_level = 0,
1953  const bool allow_merge = false,
1954  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1955 
1983  template <int spacedim>
1984 #ifndef DOXYGEN
1985  std::tuple<std::vector<std::vector<unsigned int>>,
1986  std::map<unsigned int, unsigned int>,
1987  std::map<unsigned int, std::vector<unsigned int>>>
1988 #else
1989  return_type
1990 #endif
1992  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1993  const std::vector<Point<spacedim>> & points);
1994 
1995 
2030  template <int spacedim>
2031 #ifndef DOXYGEN
2032  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2033  std::map<unsigned int, unsigned int>,
2034  std::map<unsigned int, std::vector<unsigned int>>>
2035 #else
2036  return_type
2037 #endif
2039  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2040  const std::vector<Point<spacedim>> & points);
2041 
2042 
2051  template <int dim, int spacedim>
2052  std::vector<
2053  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2055 
2068  template <int dim, int spacedim>
2069  std::vector<std::vector<Tensor<1, spacedim>>>
2071  const Triangulation<dim, spacedim> &mesh,
2072  const std::vector<
2074  &vertex_to_cells);
2075 
2076 
2084  template <int dim, int spacedim>
2085  unsigned int
2088  const Point<spacedim> & position,
2089  const Mapping<dim, spacedim> & mapping =
2090  (ReferenceCells::get_hypercube<dim>()
2091 #ifndef _MSC_VER
2092  .template get_default_linear_mapping<dim, spacedim>()
2093 #else
2094  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2095 #endif
2096  ));
2097 
2109  template <int dim, int spacedim>
2110  std::map<unsigned int, types::global_vertex_index>
2113 
2125  template <int dim, int spacedim>
2126  std::pair<unsigned int, double>
2129 
2135 
2144  template <int dim, int spacedim>
2145  void
2148  DynamicSparsityPattern & connectivity);
2149 
2158  template <int dim, int spacedim>
2159  void
2162  DynamicSparsityPattern & connectivity);
2163 
2172  template <int dim, int spacedim>
2173  void
2176  const unsigned int level,
2177  DynamicSparsityPattern & connectivity);
2178 
2199  template <int dim, int spacedim>
2200  void
2201  partition_triangulation(const unsigned int n_partitions,
2203  const SparsityTools::Partitioner partitioner =
2205 
2216  template <int dim, int spacedim>
2217  void
2218  partition_triangulation(const unsigned int n_partitions,
2219  const std::vector<unsigned int> &cell_weights,
2221  const SparsityTools::Partitioner partitioner =
2223 
2269  template <int dim, int spacedim>
2270  void
2271  partition_triangulation(const unsigned int n_partitions,
2272  const SparsityPattern & cell_connection_graph,
2274  const SparsityTools::Partitioner partitioner =
2276 
2287  template <int dim, int spacedim>
2288  void
2289  partition_triangulation(const unsigned int n_partitions,
2290  const std::vector<unsigned int> &cell_weights,
2291  const SparsityPattern & cell_connection_graph,
2293  const SparsityTools::Partitioner partitioner =
2295 
2310  template <int dim, int spacedim>
2311  void
2312  partition_triangulation_zorder(const unsigned int n_partitions,
2314  const bool group_siblings = true);
2315 
2327  template <int dim, int spacedim>
2328  void
2330 
2338  template <int dim, int spacedim>
2339  std::vector<types::subdomain_id>
2341  const std::vector<CellId> & cell_ids);
2342 
2353  template <int dim, int spacedim>
2354  void
2356  std::vector<types::subdomain_id> & subdomain);
2357 
2372  template <int dim, int spacedim>
2373  unsigned int
2376  const types::subdomain_id subdomain);
2377 
2407  template <int dim, int spacedim>
2408  std::vector<bool>
2410 
2416 
2450  template <typename MeshType>
2451  std::list<std::pair<typename MeshType::cell_iterator,
2452  typename MeshType::cell_iterator>>
2453  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2454 
2464  template <int dim, int spacedim>
2465  bool
2467  const Triangulation<dim, spacedim> &mesh_2);
2468 
2478  template <typename MeshType>
2479  bool
2480  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2481 
2487 
2503  template <int dim, int spacedim>
2507  & distorted_cells,
2509 
2510 
2511 
2520 
2521 
2559  template <class MeshType>
2560  std::vector<typename MeshType::active_cell_iterator>
2561  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2562 
2563 
2585  template <class Container>
2586  std::vector<typename Container::cell_iterator>
2588  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2589 
2656  template <class Container>
2657  void
2659  const std::vector<typename Container::active_cell_iterator> &patch,
2661  &local_triangulation,
2662  std::map<
2663  typename Triangulation<Container::dimension,
2664  Container::space_dimension>::active_cell_iterator,
2665  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2666 
2698  template <int dim, int spacedim>
2699  std::map<
2701  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2703 
2704 
2711 
2717  template <typename CellIterator>
2719  {
2723  CellIterator cell[2];
2724 
2729  unsigned int face_idx[2];
2730 
2736  std::bitset<3> orientation;
2737 
2751  };
2752 
2753 
2817  template <typename FaceIterator>
2818  bool
2820  std::bitset<3> & orientation,
2821  const FaceIterator & face1,
2822  const FaceIterator & face2,
2823  const unsigned int direction,
2827 
2828 
2832  template <typename FaceIterator>
2833  bool
2835  const FaceIterator & face1,
2836  const FaceIterator & face2,
2837  const unsigned int direction,
2841 
2842 
2899  template <typename MeshType>
2900  void
2902  const MeshType & mesh,
2903  const types::boundary_id b_id1,
2904  const types::boundary_id b_id2,
2905  const unsigned int direction,
2907  & matched_pairs,
2908  const Tensor<1, MeshType::space_dimension> &offset =
2911 
2912 
2935  template <typename MeshType>
2936  void
2938  const MeshType & mesh,
2939  const types::boundary_id b_id,
2940  const unsigned int direction,
2942  & matched_pairs,
2943  const ::Tensor<1, MeshType::space_dimension> &offset =
2946 
2952 
2973  template <int dim, int spacedim>
2974  void
2976  const bool reset_boundary_ids = false);
2977 
2999  template <int dim, int spacedim>
3000  void
3002  const std::vector<types::boundary_id> &src_boundary_ids,
3003  const std::vector<types::manifold_id> &dst_manifold_ids,
3005  const std::vector<types::boundary_id> &reset_boundary_ids = {});
3006 
3036  template <int dim, int spacedim>
3037  void
3039  const bool compute_face_ids = false);
3040 
3065  template <int dim, int spacedim>
3066  void
3069  const std::function<types::manifold_id(
3070  const std::set<types::manifold_id> &)> &disambiguation_function =
3071  [](const std::set<types::manifold_id> &manifold_ids) {
3072  if (manifold_ids.size() == 1)
3073  return *manifold_ids.begin();
3074  else
3076  },
3077  bool overwrite_only_flat_manifold_ids = true);
3164  template <typename DataType, typename MeshType>
3165  void
3167  const MeshType & mesh,
3168  const std::function<std_cxx17::optional<DataType>(
3169  const typename MeshType::active_cell_iterator &)> &pack,
3170  const std::function<void(const typename MeshType::active_cell_iterator &,
3171  const DataType &)> & unpack,
3172  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3173  &cell_filter =
3175 
3186  template <typename DataType, typename MeshType>
3187  void
3189  const MeshType & mesh,
3190  const std::function<std_cxx17::optional<DataType>(
3191  const typename MeshType::level_cell_iterator &)> &pack,
3192  const std::function<void(const typename MeshType::level_cell_iterator &,
3193  const DataType &)> & unpack,
3194  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3196  true});
3197 
3198  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3199  * boxes @p local_bboxes.
3200  *
3201  * This function is meant to exchange bounding boxes describing the locally
3202  * owned cells in a distributed triangulation obtained with the function
3203  * GridTools::compute_mesh_predicate_bounding_box .
3204  *
3205  * The output vector's size is the number of processes of the MPI
3206  * communicator:
3207  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3208  */
3209  template <int spacedim>
3210  std::vector<std::vector<BoundingBox<spacedim>>>
3212  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3213  const MPI_Comm & mpi_communicator);
3214 
3247  template <int spacedim>
3250  const std::vector<BoundingBox<spacedim>> &local_description,
3251  const MPI_Comm & mpi_communicator);
3252 
3270  template <int dim, int spacedim>
3271  void
3274  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3275  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3276 
3283  template <int dim, int spacedim>
3284  std::map<unsigned int, std::set<::types::subdomain_id>>
3287 
3303  template <int dim, typename T>
3305  {
3309  std::vector<CellId> cell_ids;
3310 
3314  std::vector<T> data;
3315 
3324  template <class Archive>
3325  void
3326  save(Archive &ar, const unsigned int) const
3327  {
3328  // Implement the code inline as some compilers do otherwise complain
3329  // about the use of a deprecated interface.
3330  Assert(cell_ids.size() == data.size(),
3331  ExcDimensionMismatch(cell_ids.size(), data.size()));
3332  // archive the cellids in an efficient binary format
3333  const std::size_t n_cells = cell_ids.size();
3334  ar & n_cells;
3335  for (const auto &id : cell_ids)
3336  {
3337  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3338  ar & binary_cell_id;
3339  }
3340 
3341  ar &data;
3342  }
3343 
3350  template <class Archive>
3351  void
3352  load(Archive &ar, const unsigned int)
3353  {
3354  // Implement the code inline as some compilers do otherwise complain
3355  // about the use of a deprecated interface.
3356  std::size_t n_cells;
3357  ar & n_cells;
3358  cell_ids.clear();
3359  cell_ids.reserve(n_cells);
3360  for (unsigned int c = 0; c < n_cells; ++c)
3361  {
3362  CellId::binary_type value;
3363  ar & value;
3364  cell_ids.emplace_back(value);
3365  }
3366  ar &data;
3367  }
3368 
3369 
3370 #ifdef DOXYGEN
3376  template <class Archive>
3377  void
3378  serialize(Archive &archive, const unsigned int version);
3379 #else
3380  // This macro defines the serialize() method that is compatible with
3381  // the templated save() and load() method that have been implemented.
3382  BOOST_SERIALIZATION_SPLIT_MEMBER()
3383 #endif
3384  };
3385 
3386 
3387 
3405  template <int dim, typename VectorType>
3407  {
3408  public:
3412  using value_type = typename VectorType::value_type;
3413 
3417  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3418  const FiniteElement<dim, dim> &fe,
3419  const unsigned int n_subdivisions = 1,
3420  const double tolerance = 1e-10);
3421 
3426  void
3427  process(const DoFHandler<dim> & background_dof_handler,
3428  const VectorType & ls_vector,
3429  const double iso_level,
3430  std::vector<Point<dim>> & vertices,
3431  std::vector<CellData<dim - 1>> &cells) const;
3432 
3439  void
3441  const VectorType & ls_vector,
3442  const double iso_level,
3443  std::vector<Point<dim>> & vertices,
3444  std::vector<CellData<dim - 1>> &cells) const;
3445 
3446  private:
3451  static Quadrature<dim>
3452  create_quadrature_rule(const unsigned int n_subdivisions);
3453 
3457  void
3458  process_cell(std::vector<value_type> & ls_values,
3459  const std::vector<Point<dim>> & points,
3460  const double iso_level,
3461  std::vector<Point<dim>> & vertices,
3462  std::vector<CellData<dim - 1>> &cells) const;
3463 
3470  void
3471  process_sub_cell(const std::vector<value_type> & ls_values,
3472  const std::vector<Point<2>> & points,
3473  const std::vector<unsigned int> mask,
3474  const double iso_level,
3475  std::vector<Point<2>> & vertices,
3476  std::vector<CellData<1>> & cells) const;
3477 
3481  void
3482  process_sub_cell(const std::vector<value_type> & ls_values,
3483  const std::vector<Point<3>> & points,
3484  const std::vector<unsigned int> mask,
3485  const double iso_level,
3486  std::vector<Point<3>> & vertices,
3487  std::vector<CellData<2>> & cells) const;
3488 
3493  const unsigned int n_subdivisions;
3494 
3499  const double tolerance;
3500 
3506  };
3507 
3508 
3509 
3514 
3519  int,
3520  << "The number of partitions you gave is " << arg1
3521  << ", but must be greater than zero.");
3526  int,
3527  << "The subdomain id " << arg1
3528  << " has no cells associated with it.");
3533 
3538  double,
3539  << "The scaling factor must be positive, but it is " << arg1
3540  << '.');
3541 
3546  unsigned int,
3547  << "The given vertex with index " << arg1
3548  << " is not used in the given triangulation.");
3549 
3552 } /*namespace GridTools*/
3553 
3554 
3565  "The edges of the mesh are not consistently orientable.");
3566 
3567 
3568 
3569 /* ----------------- Template function --------------- */
3570 
3571 #ifndef DOXYGEN
3572 
3573 namespace GridTools
3574 {
3575  template <int dim>
3576  double
3577  cell_measure(
3578  const std::vector<Point<dim>> &all_vertices,
3579  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3580  {
3581  // We forward call to the ArrayView version:
3582  const ArrayView<const unsigned int> view(
3584  return cell_measure(all_vertices, view);
3585  }
3586 
3587 
3588 
3589  // This specialization is defined here so that the general template in the
3590  // source file doesn't need to have further 1D overloads for the internal
3591  // functions it calls.
3592  template <>
3596  {
3597  return {};
3598  }
3599 
3600 
3601 
3602  template <int dim, typename Predicate, int spacedim>
3603  void
3604  transform(const Predicate & predicate,
3606  {
3607  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3608 
3609  // loop over all active cells, and
3610  // transform those vertices that
3611  // have not yet been touched. note
3612  // that we get to all vertices in
3613  // the triangulation by only
3614  // visiting the active cells.
3616  cell = triangulation.begin_active(),
3617  endc = triangulation.end();
3618  for (; cell != endc; ++cell)
3619  for (const unsigned int v : cell->vertex_indices())
3620  if (treated_vertices[cell->vertex_index(v)] == false)
3621  {
3622  // transform this vertex
3623  cell->vertex(v) = predicate(cell->vertex(v));
3624  // and mark it as treated
3625  treated_vertices[cell->vertex_index(v)] = true;
3626  };
3627 
3628 
3629  // now fix any vertices on hanging nodes so that we don't create any holes
3630  if (dim == 2)
3631  {
3633  cell = triangulation.begin_active(),
3634  endc = triangulation.end();
3635  for (; cell != endc; ++cell)
3636  for (const unsigned int face : cell->face_indices())
3637  if (cell->face(face)->has_children() &&
3638  !cell->face(face)->at_boundary())
3639  {
3640  Assert(cell->reference_cell() ==
3641  ReferenceCells::get_hypercube<dim>(),
3642  ExcNotImplemented());
3643 
3644  // this line has children
3645  cell->face(face)->child(0)->vertex(1) =
3646  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3647  2;
3648  }
3649  }
3650  else if (dim == 3)
3651  {
3653  cell = triangulation.begin_active(),
3654  endc = triangulation.end();
3655  for (; cell != endc; ++cell)
3656  for (const unsigned int face : cell->face_indices())
3657  if (cell->face(face)->has_children() &&
3658  !cell->face(face)->at_boundary())
3659  {
3660  Assert(cell->reference_cell() ==
3661  ReferenceCells::get_hypercube<dim>(),
3662  ExcNotImplemented());
3663 
3664  // this face has hanging nodes
3665  cell->face(face)->child(0)->vertex(1) =
3666  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3667  2.0;
3668  cell->face(face)->child(0)->vertex(2) =
3669  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3670  2.0;
3671  cell->face(face)->child(1)->vertex(3) =
3672  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3673  2.0;
3674  cell->face(face)->child(2)->vertex(3) =
3675  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3676  2.0;
3677 
3678  // center of the face
3679  cell->face(face)->child(0)->vertex(3) =
3680  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3681  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3682  4.0;
3683  }
3684  }
3685 
3686  // Make sure FEValues notices that the mesh has changed
3687  triangulation.signals.mesh_movement();
3688  }
3689 
3690 
3691 
3692  template <class MeshType>
3693  std::vector<typename MeshType::active_cell_iterator>
3694  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3695  {
3696  std::vector<typename MeshType::active_cell_iterator> child_cells;
3697 
3698  if (cell->has_children())
3699  {
3700  for (unsigned int child = 0; child < cell->n_children(); ++child)
3701  if (cell->child(child)->has_children())
3702  {
3703  const std::vector<typename MeshType::active_cell_iterator>
3704  children = get_active_child_cells<MeshType>(cell->child(child));
3705  child_cells.insert(child_cells.end(),
3706  children.begin(),
3707  children.end());
3708  }
3709  else
3710  child_cells.push_back(cell->child(child));
3711  }
3712 
3713  return child_cells;
3714  }
3715 
3716 
3717 
3718  template <class MeshType>
3719  void
3721  const typename MeshType::active_cell_iterator & cell,
3722  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3723  {
3724  active_neighbors.clear();
3725  for (const unsigned int n : cell->face_indices())
3726  if (!cell->at_boundary(n))
3727  {
3728  if (MeshType::dimension == 1)
3729  {
3730  // check children of neighbor. note
3731  // that in 1d children of the neighbor
3732  // may be further refined. In 1d the
3733  // case is simple since we know what
3734  // children bound to the present cell
3735  typename MeshType::cell_iterator neighbor_child =
3736  cell->neighbor(n);
3737  if (!neighbor_child->is_active())
3738  {
3739  while (neighbor_child->has_children())
3740  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3741 
3742  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3743  ExcInternalError());
3744  }
3745  active_neighbors.push_back(neighbor_child);
3746  }
3747  else
3748  {
3749  if (cell->face(n)->has_children())
3750  // this neighbor has children. find
3751  // out which border to the present
3752  // cell
3753  for (unsigned int c = 0;
3754  c < cell->face(n)->n_active_descendants();
3755  ++c)
3756  active_neighbors.push_back(
3757  cell->neighbor_child_on_subface(n, c));
3758  else
3759  {
3760  // the neighbor must be active
3761  // himself
3762  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3763  active_neighbors.push_back(cell->neighbor(n));
3764  }
3765  }
3766  }
3767  }
3768 
3769 
3770 
3771  namespace internal
3772  {
3773  namespace ProjectToObject
3774  {
3787  struct CrossDerivative
3788  {
3789  const unsigned int direction_0;
3790  const unsigned int direction_1;
3791 
3792  CrossDerivative(const unsigned int d0, const unsigned int d1);
3793  };
3794 
3795  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3796  const unsigned int d1)
3797  : direction_0(d0)
3798  , direction_1(d1)
3799  {}
3800 
3801 
3802 
3807  template <typename F>
3808  inline auto
3809  centered_first_difference(const double center,
3810  const double step,
3811  const F &f) -> decltype(f(center) - f(center))
3812  {
3813  return (f(center + step) - f(center - step)) / (2.0 * step);
3814  }
3815 
3816 
3817 
3822  template <typename F>
3823  inline auto
3824  centered_second_difference(const double center,
3825  const double step,
3826  const F &f) -> decltype(f(center) - f(center))
3827  {
3828  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3829  (step * step);
3830  }
3831 
3832 
3833 
3843  template <int structdim, typename F>
3844  inline auto
3845  cross_stencil(
3846  const CrossDerivative cross_derivative,
3848  const double step,
3849  const F &f) -> decltype(f(center) - f(center))
3850  {
3852  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3853  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3854  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3855  1.0 / 3.0 * f(center - simplex_vector) +
3856  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3857  step;
3858  }
3859 
3860 
3861 
3868  template <int spacedim, int structdim, typename F>
3869  inline double
3870  gradient_entry(
3871  const unsigned int row_n,
3872  const unsigned int dependent_direction,
3873  const Point<spacedim> &p0,
3875  const double step,
3876  const F & f)
3877  {
3879  dependent_direction <
3881  ExcMessage("This function assumes that the last weight is a "
3882  "dependent variable (and hence we cannot take its "
3883  "derivative directly)."));
3884  Assert(row_n != dependent_direction,
3885  ExcMessage(
3886  "We cannot differentiate with respect to the variable "
3887  "that is assumed to be dependent."));
3888 
3889  const Point<spacedim> manifold_point = f(center);
3890  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3891  {row_n, dependent_direction}, center, step, f);
3892  double entry = 0.0;
3893  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3894  entry +=
3895  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3896  return entry;
3897  }
3898 
3904  template <typename Iterator, int spacedim, int structdim>
3906  project_to_d_linear_object(const Iterator & object,
3907  const Point<spacedim> &trial_point)
3908  {
3909  // let's look at this for simplicity for a quad (structdim==2) in a
3910  // space with spacedim>2 (notate trial_point by y): all points on the
3911  // surface are given by
3912  // x(\xi) = sum_i v_i phi_x(\xi)
3913  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3914  // reference coordinates of the quad. so what we are trying to do is
3915  // find a point x on the surface that is closest to the point y. there
3916  // are different ways to solve this problem, but in the end it's a
3917  // nonlinear problem and we have to find reference coordinates \xi so
3918  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3919  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3920  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3921  // have to use a Newton method to find the answer. This leads to the
3922  // following formulation of Newton steps:
3923  //
3924  // Given \xi_k, find \delta\xi_k so that
3925  // H_k \delta\xi_k = - F_k
3926  // where H_k is an approximation to the second derivatives of J at
3927  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3928  // number of times until the right hand side is small enough. As a
3929  // stopping criterion, we terminate if ||\delta\xi||<eps.
3930  //
3931  // As for the Hessian, the best choice would be
3932  // H_k = J''(\xi_k)
3933  // but we'll opt for the simpler Gauss-Newton form
3934  // H_k = A^T A
3935  // i.e.
3936  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3937  // \partial_n phi_i *
3938  // \partial_m phi_j
3939  // we start at xi=(0.5, 0.5).
3940  Point<structdim> xi;
3941  for (unsigned int d = 0; d < structdim; ++d)
3942  xi[d] = 0.5;
3943 
3944  Point<spacedim> x_k;
3945  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3946  x_k += object->vertex(i) *
3948 
3949  do
3950  {
3952  for (const unsigned int i :
3954  F_k +=
3955  (x_k - trial_point) * object->vertex(i) *
3957  i);
3958 
3960  for (const unsigned int i :
3962  for (const unsigned int j :
3964  {
3967  xi, i),
3969  xi, j));
3970  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3971  }
3972 
3973  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3974  xi += delta_xi;
3975 
3976  x_k = Point<spacedim>();
3977  for (const unsigned int i :
3979  x_k += object->vertex(i) *
3981 
3982  if (delta_xi.norm() < 1e-7)
3983  break;
3984  }
3985  while (true);
3986 
3987  return x_k;
3988  }
3989  } // namespace ProjectToObject
3990  } // namespace internal
3991 
3992 
3993 
3994  namespace internal
3995  {
3996  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3997  // inside the project_to_object function below.
3998  template <int structdim>
3999  inline bool
4000  weights_are_ok(
4002  {
4003  // clang has trouble figuring out structdim here, so define it
4004  // again:
4005  static const std::size_t n_vertices_per_cell =
4007  n_independent_components;
4008  std::array<double, n_vertices_per_cell> copied_weights;
4009  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4010  {
4011  copied_weights[i] = v[i];
4012  if (v[i] < 0.0 || v[i] > 1.0)
4013  return false;
4014  }
4015 
4016  // check the sum: try to avoid some roundoff errors by summing in order
4017  std::sort(copied_weights.begin(), copied_weights.end());
4018  const double sum =
4019  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4020  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4021  }
4022  } // namespace internal
4023 
4024  template <typename Iterator>
4027  const Iterator & object,
4029  {
4030  const int spacedim = Iterator::AccessorType::space_dimension;
4031  const int structdim = Iterator::AccessorType::structure_dimension;
4032 
4033  Point<spacedim> projected_point = trial_point;
4034 
4035  if (structdim >= spacedim)
4036  return projected_point;
4037  else if (structdim == 1 || structdim == 2)
4038  {
4039  using namespace internal::ProjectToObject;
4040  // Try to use the special flat algorithm for quads (this is better
4041  // than the general algorithm in 3D). This does not take into account
4042  // whether projected_point is outside the quad, but we optimize along
4043  // lines below anyway:
4044  const int dim = Iterator::AccessorType::dimension;
4045  const Manifold<dim, spacedim> &manifold = object->get_manifold();
4046  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4047  &manifold) != nullptr)
4048  {
4049  projected_point =
4050  project_to_d_linear_object<Iterator, spacedim, structdim>(
4051  object, trial_point);
4052  }
4053  else
4054  {
4055  // We want to find a point on the convex hull (defined by the
4056  // vertices of the object and the manifold description) that is
4057  // relatively close to the trial point. This has a few issues:
4058  //
4059  // 1. For a general convex hull we are not guaranteed that a unique
4060  // minimum exists.
4061  // 2. The independent variables in the optimization process are the
4062  // weights given to Manifold::get_new_point, which must sum to 1,
4063  // so we cannot use standard finite differences to approximate a
4064  // gradient.
4065  //
4066  // There is not much we can do about 1., but for 2. we can derive
4067  // finite difference stencils that work on a structdim-dimensional
4068  // simplex and rewrite the optimization problem to use those
4069  // instead. Consider the structdim 2 case and let
4070  //
4071  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4072  // c2, c3})
4073  //
4074  // where {c0, c1, c2, c3} are the weights for the four vertices on
4075  // the quadrilateral. We seek to minimize the Euclidean distance
4076  // between F(...) and trial_point. We can solve for c3 in terms of
4077  // the other weights and get, for one coordinate direction
4078  //
4079  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4080  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4081  //
4082  // where we substitute back in for c3 after taking the
4083  // derivative. We can compute a stencil for the cross derivative
4084  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4085  // (and gradient_entry computes the sum over the independent
4086  // variables). Below, we somewhat arbitrarily pick the last
4087  // component as the dependent one.
4088  //
4089  // Since we can now calculate derivatives of the objective
4090  // function we can use gradient descent to minimize it.
4091  //
4092  // Of course, this is much simpler in the structdim = 1 case (we
4093  // could rewrite the projection as a 1D optimization problem), but
4094  // to reduce the potential for bugs we use the same code in both
4095  // cases.
4096  const double step_size = object->diameter() / 64.0;
4097 
4098  constexpr unsigned int n_vertices_per_cell =
4100 
4101  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4102  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4103  ++vertex_n)
4104  vertices[vertex_n] = object->vertex(vertex_n);
4105 
4106  auto get_point_from_weights =
4107  [&](const Tensor<1, n_vertices_per_cell> &weights)
4108  -> Point<spacedim> {
4109  return object->get_manifold().get_new_point(
4110  make_array_view(vertices.begin(), vertices.end()),
4111  make_array_view(weights.begin_raw(), weights.end_raw()));
4112  };
4113 
4114  // pick the initial weights as (normalized) inverse distances from
4115  // the trial point:
4116  Tensor<1, n_vertices_per_cell> guess_weights;
4117  double guess_weights_sum = 0.0;
4118  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4119  ++vertex_n)
4120  {
4121  const double distance =
4122  vertices[vertex_n].distance(trial_point);
4123  if (distance == 0.0)
4124  {
4125  guess_weights = 0.0;
4126  guess_weights[vertex_n] = 1.0;
4127  guess_weights_sum = 1.0;
4128  break;
4129  }
4130  else
4131  {
4132  guess_weights[vertex_n] = 1.0 / distance;
4133  guess_weights_sum += guess_weights[vertex_n];
4134  }
4135  }
4136  guess_weights /= guess_weights_sum;
4137  Assert(internal::weights_are_ok<structdim>(guess_weights),
4138  ExcInternalError());
4139 
4140  // The optimization algorithm consists of two parts:
4141  //
4142  // 1. An outer loop where we apply the gradient descent algorithm.
4143  // 2. An inner loop where we do a line search to find the optimal
4144  // length of the step one should take in the gradient direction.
4145  //
4146  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4147  {
4148  const unsigned int dependent_direction =
4149  n_vertices_per_cell - 1;
4150  Tensor<1, n_vertices_per_cell> current_gradient;
4151  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4152  ++row_n)
4153  {
4154  if (row_n != dependent_direction)
4155  {
4156  current_gradient[row_n] =
4157  gradient_entry<spacedim, structdim>(
4158  row_n,
4159  dependent_direction,
4160  trial_point,
4161  guess_weights,
4162  step_size,
4163  get_point_from_weights);
4164 
4165  current_gradient[dependent_direction] -=
4166  current_gradient[row_n];
4167  }
4168  }
4169 
4170  // We need to travel in the -gradient direction, as noted
4171  // above, but we may not want to take a full step in that
4172  // direction; instead, guess that we will go -0.5*gradient and
4173  // do quasi-Newton iteration to pick the best multiplier. The
4174  // goal is to find a scalar alpha such that
4175  //
4176  // F(x - alpha g)
4177  //
4178  // is minimized, where g is the gradient and F is the
4179  // objective function. To find the optimal value we find roots
4180  // of the derivative of the objective function with respect to
4181  // alpha by Newton iteration, where we approximate the first
4182  // and second derivatives of F(x - alpha g) with centered
4183  // finite differences.
4184  double gradient_weight = -0.5;
4185  auto gradient_weight_objective_function =
4186  [&](const double gradient_weight_guess) -> double {
4187  return (trial_point -
4188  get_point_from_weights(guess_weights +
4189  gradient_weight_guess *
4190  current_gradient))
4191  .norm_square();
4192  };
4193 
4194  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4195  {
4196  const double update_numerator = centered_first_difference(
4197  gradient_weight,
4198  step_size,
4199  gradient_weight_objective_function);
4200  const double update_denominator =
4201  centered_second_difference(
4202  gradient_weight,
4203  step_size,
4204  gradient_weight_objective_function);
4205 
4206  // avoid division by zero. Note that we limit the gradient
4207  // weight below
4208  if (std::abs(update_denominator) == 0.0)
4209  break;
4210  gradient_weight =
4211  gradient_weight - update_numerator / update_denominator;
4212 
4213  // Put a fairly lenient bound on the largest possible
4214  // gradient (things tend to be locally flat, so the gradient
4215  // itself is usually small)
4216  if (std::abs(gradient_weight) > 10)
4217  {
4218  gradient_weight = -10.0;
4219  break;
4220  }
4221  }
4222 
4223  // It only makes sense to take convex combinations with weights
4224  // between zero and one. If the update takes us outside of this
4225  // region then rescale the update to stay within the region and
4226  // try again
4227  Tensor<1, n_vertices_per_cell> tentative_weights =
4228  guess_weights + gradient_weight * current_gradient;
4229 
4230  double new_gradient_weight = gradient_weight;
4231  for (unsigned int iteration_count = 0; iteration_count < 40;
4232  ++iteration_count)
4233  {
4234  if (internal::weights_are_ok<structdim>(tentative_weights))
4235  break;
4236 
4237  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4238  {
4239  if (tentative_weights[i] < 0.0)
4240  {
4241  tentative_weights -=
4242  (tentative_weights[i] / current_gradient[i]) *
4243  current_gradient;
4244  }
4245  if (tentative_weights[i] < 0.0 ||
4246  1.0 < tentative_weights[i])
4247  {
4248  new_gradient_weight /= 2.0;
4249  tentative_weights =
4250  guess_weights +
4251  new_gradient_weight * current_gradient;
4252  }
4253  }
4254  }
4255 
4256  // the update might still send us outside the valid region, so
4257  // check again and quit if the update is still not valid
4258  if (!internal::weights_are_ok<structdim>(tentative_weights))
4259  break;
4260 
4261  // if we cannot get closer by traveling in the gradient
4262  // direction then quit
4263  if (get_point_from_weights(tentative_weights)
4264  .distance(trial_point) <
4265  get_point_from_weights(guess_weights).distance(trial_point))
4266  guess_weights = tentative_weights;
4267  else
4268  break;
4269  Assert(internal::weights_are_ok<structdim>(guess_weights),
4270  ExcInternalError());
4271  }
4272  Assert(internal::weights_are_ok<structdim>(guess_weights),
4273  ExcInternalError());
4274  projected_point = get_point_from_weights(guess_weights);
4275  }
4276 
4277  // if structdim == 2 and the optimal point is not on the interior then
4278  // we may be able to get a more accurate result by projecting onto the
4279  // lines.
4280  if (structdim == 2)
4281  {
4282  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4283  line_projections;
4284  for (unsigned int line_n = 0;
4285  line_n < GeometryInfo<structdim>::lines_per_cell;
4286  ++line_n)
4287  {
4288  line_projections[line_n] =
4289  project_to_object(object->line(line_n), trial_point);
4290  }
4291  std::sort(line_projections.begin(),
4292  line_projections.end(),
4293  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4294  return a.distance(trial_point) <
4295  b.distance(trial_point);
4296  });
4297  if (line_projections[0].distance(trial_point) <
4298  projected_point.distance(trial_point))
4299  projected_point = line_projections[0];
4300  }
4301  }
4302  else
4303  {
4304  Assert(false, ExcNotImplemented());
4305  return projected_point;
4306  }
4307 
4308  return projected_point;
4309  }
4310 
4311 
4312 
4313  namespace internal
4314  {
4315  template <typename DataType,
4316  typename MeshType,
4317  typename MeshCellIteratorType>
4318  inline void
4319  exchange_cell_data(
4320  const MeshType &mesh,
4321  const std::function<
4322  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4323  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4324  & unpack,
4325  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4326  const std::function<void(
4327  const std::function<void(const MeshCellIteratorType &,
4328  const types::subdomain_id)> &)> &process_cells,
4329  const std::function<std::set<types::subdomain_id>(
4330  const parallel::TriangulationBase<MeshType::dimension,
4331  MeshType::space_dimension> &)>
4332  &compute_ghost_owners)
4333  {
4334 # ifndef DEAL_II_WITH_MPI
4335  (void)mesh;
4336  (void)pack;
4337  (void)unpack;
4338  (void)cell_filter;
4339  (void)process_cells;
4340  (void)compute_ghost_owners;
4341  Assert(false, ExcNeedsMPI());
4342 # else
4343  constexpr int dim = MeshType::dimension;
4344  constexpr int spacedim = MeshType::space_dimension;
4345  auto tria =
4346  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4347  &mesh.get_triangulation());
4348  Assert(
4349  tria != nullptr,
4350  ExcMessage(
4351  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4352 
4353  if (const auto tria = dynamic_cast<
4355  &mesh.get_triangulation()))
4356  {
4357  Assert(
4358  tria->with_artificial_cells(),
4359  ExcMessage(
4360  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4361  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4362  "operate on a single layer of ghost cells. However, you have "
4363  "given a Triangulation object of type "
4364  "parallel::shared::Triangulation without artificial cells "
4365  "resulting in arbitrary numbers of ghost layers."));
4366  }
4367 
4368  // build list of cells to request for each neighbor
4369  std::set<::types::subdomain_id> ghost_owners =
4370  compute_ghost_owners(*tria);
4371  std::map<::types::subdomain_id,
4372  std::vector<typename CellId::binary_type>>
4373  neighbor_cell_list;
4374 
4375  for (const auto ghost_owner : ghost_owners)
4376  neighbor_cell_list[ghost_owner] = {};
4377 
4378  process_cells([&](const auto &cell, const auto key) -> void {
4379  if (cell_filter(cell))
4380  {
4381  constexpr int spacedim = MeshType::space_dimension;
4382  neighbor_cell_list[key].emplace_back(
4383  cell->id().template to_binary<spacedim>());
4384  }
4385  });
4386 
4387  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4388  ExcInternalError());
4389 
4390 
4391  // Before sending & receiving, make sure we protect this section with
4392  // a mutex:
4393  static Utilities::MPI::CollectiveMutex mutex;
4395  mutex, tria->get_communicator());
4396 
4397  const int mpi_tag =
4399  const int mpi_tag_reply =
4401 
4402  // send our requests
4403  std::vector<MPI_Request> requests(ghost_owners.size());
4404  {
4405  unsigned int idx = 0;
4406  for (const auto &it : neighbor_cell_list)
4407  {
4408  // send the data about the relevant cells
4409  const int ierr = MPI_Isend(it.second.data(),
4410  it.second.size() * sizeof(it.second[0]),
4411  MPI_BYTE,
4412  it.first,
4413  mpi_tag,
4415  &requests[idx]);
4416  AssertThrowMPI(ierr);
4417  ++idx;
4418  }
4419  }
4420 
4421  // receive requests and reply with the results
4422  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4423  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4424 
4425  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4426  {
4427  MPI_Status status;
4428  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4429  mpi_tag,
4431  &status);
4432  AssertThrowMPI(ierr);
4433 
4434  int len;
4435  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4436  AssertThrowMPI(ierr);
4437  Assert(len % sizeof(typename CellId::binary_type) == 0,
4438  ExcInternalError());
4439 
4440  const unsigned int n_cells =
4441  len / sizeof(typename CellId::binary_type);
4442  std::vector<typename CellId::binary_type> cells_with_requests(
4443  n_cells);
4444  std::vector<DataType> data_to_send;
4445  data_to_send.reserve(n_cells);
4446  std::vector<bool> cell_carries_data(n_cells, false);
4447 
4448  ierr = MPI_Recv(cells_with_requests.data(),
4449  len,
4450  MPI_BYTE,
4451  status.MPI_SOURCE,
4452  status.MPI_TAG,
4454  &status);
4455  AssertThrowMPI(ierr);
4456 
4457  // store data for each cell
4458  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4459  {
4460  const auto cell =
4461  tria->create_cell_iterator(CellId(cells_with_requests[c]));
4462 
4463  MeshCellIteratorType mesh_it(tria,
4464  cell->level(),
4465  cell->index(),
4466  &mesh);
4467 
4468  const std_cxx17::optional<DataType> data = pack(mesh_it);
4469  if (data)
4470  {
4471  data_to_send.emplace_back(*data);
4472  cell_carries_data[c] = true;
4473  }
4474  }
4475 
4476  // collect data for sending the reply in a buffer
4477 
4478  // (a) make room for storing the local offsets in case we receive
4479  // other data
4480  sendbuffers[idx].resize(sizeof(std::size_t));
4481 
4482  // (b) append the actual data and store how much memory it
4483  // corresponds to, which we then insert into the leading position of
4484  // the sendbuffer
4485  std::size_t size_of_send =
4486  Utilities::pack(data_to_send,
4487  sendbuffers[idx],
4488  /*enable_compression*/ false);
4489  std::memcpy(sendbuffers[idx].data(),
4490  &size_of_send,
4491  sizeof(std::size_t));
4492 
4493  // (c) append information of certain cells that got left out in case
4494  // we need it
4495  if (data_to_send.size() < n_cells)
4496  Utilities::pack(cell_carries_data,
4497  sendbuffers[idx],
4498  /*enable_compression*/ false);
4499 
4500  // send data
4501  ierr = MPI_Isend(sendbuffers[idx].data(),
4502  sendbuffers[idx].size(),
4503  MPI_BYTE,
4504  status.MPI_SOURCE,
4505  mpi_tag_reply,
4507  &reply_requests[idx]);
4508  AssertThrowMPI(ierr);
4509  }
4510 
4511  // finally receive the replies
4512  std::vector<char> receive;
4513  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4514  {
4515  MPI_Status status;
4516  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4517  mpi_tag_reply,
4519  &status);
4520  AssertThrowMPI(ierr);
4521 
4522  int len;
4523  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4524  AssertThrowMPI(ierr);
4525 
4526  receive.resize(len);
4527 
4528  ierr = MPI_Recv(receive.data(),
4529  len,
4530  MPI_BYTE,
4531  status.MPI_SOURCE,
4532  status.MPI_TAG,
4534  &status);
4535  AssertThrowMPI(ierr);
4536 
4537  // (a) first determine the length of the data section in the
4538  // received buffer
4539  auto data_iterator = receive.begin();
4540  std::size_t size_of_received_data =
4541  Utilities::unpack<std::size_t>(data_iterator,
4542  data_iterator + sizeof(std::size_t));
4543  data_iterator += sizeof(std::size_t);
4544 
4545  // (b) unpack the data section in the indicated region
4546  auto received_data = Utilities::unpack<std::vector<DataType>>(
4547  data_iterator,
4548  data_iterator + size_of_received_data,
4549  /*enable_compression*/ false);
4550  data_iterator += size_of_received_data;
4551 
4552  // (c) check if the received data contained fewer entries than the
4553  // number of cells we identified in the beginning, in which case we
4554  // need to extract the boolean vector with the relevant information
4555  const std::vector<typename CellId::binary_type> &this_cell_list =
4556  neighbor_cell_list[status.MPI_SOURCE];
4557  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4558  std::vector<bool> cells_with_data;
4559  if (received_data.size() < this_cell_list.size())
4560  {
4561  cells_with_data = Utilities::unpack<std::vector<bool>>(
4562  data_iterator, receive.end(), /*enable_compression*/ false);
4563  AssertDimension(cells_with_data.size(), this_cell_list.size());
4564  }
4565 
4566  // (d) go through the received data and call the user-provided
4567  // unpack function
4568  auto received_data_iterator = received_data.begin();
4569  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4570  if (cells_with_data.empty() || cells_with_data[c])
4571  {
4573  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4574 
4575  MeshCellIteratorType cell(tria,
4576  tria_cell->level(),
4577  tria_cell->index(),
4578  &mesh);
4579 
4580  unpack(cell, *received_data_iterator);
4581  ++received_data_iterator;
4582  }
4583  }
4584 
4585  // make sure that all communication is finished
4586  // when we leave this function.
4587  if (requests.size() > 0)
4588  {
4589  const int ierr =
4590  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4591  AssertThrowMPI(ierr);
4592  }
4593  if (reply_requests.size() > 0)
4594  {
4595  const int ierr = MPI_Waitall(reply_requests.size(),
4596  reply_requests.data(),
4597  MPI_STATUSES_IGNORE);
4598  AssertThrowMPI(ierr);
4599  }
4600 
4601 
4602 # endif // DEAL_II_WITH_MPI
4603  }
4604 
4605  } // namespace internal
4606 
4607  template <typename DataType, typename MeshType>
4608  inline void
4610  const MeshType & mesh,
4611  const std::function<std_cxx17::optional<DataType>(
4612  const typename MeshType::active_cell_iterator &)> &pack,
4613  const std::function<void(const typename MeshType::active_cell_iterator &,
4614  const DataType &)> & unpack,
4615  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4616  &cell_filter)
4617  {
4618 # ifndef DEAL_II_WITH_MPI
4619  (void)mesh;
4620  (void)pack;
4621  (void)unpack;
4622  (void)cell_filter;
4623  Assert(false, ExcNeedsMPI());
4624 # else
4625  internal::exchange_cell_data<DataType,
4626  MeshType,
4627  typename MeshType::active_cell_iterator>(
4628  mesh,
4629  pack,
4630  unpack,
4631  cell_filter,
4632  [&](const auto &process) {
4633  for (const auto &cell : mesh.active_cell_iterators())
4634  if (cell->is_ghost())
4635  process(cell, cell->subdomain_id());
4636  },
4637  [](const auto &tria) { return tria.ghost_owners(); });
4638 # endif
4639  }
4640 
4641 
4642 
4643  template <typename DataType, typename MeshType>
4644  inline void
4646  const MeshType & mesh,
4647  const std::function<std_cxx17::optional<DataType>(
4648  const typename MeshType::level_cell_iterator &)> &pack,
4649  const std::function<void(const typename MeshType::level_cell_iterator &,
4650  const DataType &)> & unpack,
4651  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4652  &cell_filter)
4653  {
4654 # ifndef DEAL_II_WITH_MPI
4655  (void)mesh;
4656  (void)pack;
4657  (void)unpack;
4658  (void)cell_filter;
4659  Assert(false, ExcNeedsMPI());
4660 # else
4661  internal::exchange_cell_data<DataType,
4662  MeshType,
4663  typename MeshType::level_cell_iterator>(
4664  mesh,
4665  pack,
4666  unpack,
4667  cell_filter,
4668  [&](const auto &process) {
4669  for (const auto &cell : mesh.cell_iterators())
4670  if (cell->is_ghost_on_level())
4671  process(cell, cell->level_subdomain_id());
4672  },
4673  [](const auto &tria) { return tria.level_ghost_owners(); });
4674 # endif
4675  }
4676 } // namespace GridTools
4677 
4678 #endif // DOXYGEN
4679 
4681 
4682 #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:6642
typename VectorType::value_type value_type
Definition: grid_tools.h:3412
const unsigned int n_subdivisions
Definition: grid_tools.h:3493
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:6579
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:6596
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:6626
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:6807
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:1293
__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:6398
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:4925
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:5017
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:4950
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6521
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:6239
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:6335
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:5045
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:5930
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3732
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:3085
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2084
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4072
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2227
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:4883
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:6194
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:936
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:3003
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:6175
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:4363
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:5227
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:5198
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:4390
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2050
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:5572
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:5165
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:4178
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:2738
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:869
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3232
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:4347
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1959
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:5543
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:3831
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:3330
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:3766
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:4277
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:2259
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:3795
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:2610
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:3379
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:4417
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5133
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:5740
@ 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:3352
void serialize(Archive &archive, const unsigned int version)
void save(Archive &ar, const unsigned int) const
Definition: grid_tools.h:3326
std::vector< CellId > cell_ids
Definition: grid_tools.h:3309
std::bitset< 3 > orientation
Definition: grid_tools.h:2736
FullMatrix< double > matrix
Definition: grid_tools.h:2750
unsigned int face_idx[2]
Definition: grid_tools.h:2729
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1181
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1156
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