Reference documentation for deal.II version Git abc1da1e07 2020-07-06 12:09:19 -0400
\(\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 - 2020 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 
25 
27 
29 
30 # include <deal.II/fe/mapping.h>
31 # include <deal.II/fe/mapping_q1.h>
32 
33 # include <deal.II/grid/manifold.h>
34 # include <deal.II/grid/tria.h>
37 
38 # include <deal.II/hp/dof_handler.h>
39 
41 
42 # include <deal.II/numerics/rtree.h>
43 
45 # include <boost/archive/binary_iarchive.hpp>
46 # include <boost/archive/binary_oarchive.hpp>
47 # include <boost/geometry/index/rtree.hpp>
48 # include <boost/serialization/array.hpp>
49 # include <boost/serialization/vector.hpp>
50 
51 # ifdef DEAL_II_WITH_ZLIB
52 # include <boost/iostreams/device/back_inserter.hpp>
53 # include <boost/iostreams/filter/gzip.hpp>
54 # include <boost/iostreams/filtering_stream.hpp>
55 # include <boost/iostreams/stream.hpp>
56 # endif
58 
59 # include <bitset>
60 # include <list>
61 # include <set>
62 
64 
65 // Forward declarations
66 # ifndef DOXYGEN
67 namespace parallel
68 {
69  namespace distributed
70  {
71  template <int, int>
72  class Triangulation;
73  }
74 } // namespace parallel
75 
76 namespace hp
77 {
78  template <int, int>
79  class MappingCollection;
80 }
81 
82 class SparsityPattern;
83 # endif
84 
85 namespace internal
86 {
87  template <int dim, int spacedim, class MeshType>
89  {
90  public:
91 # ifndef _MSC_VER
92  using type = typename MeshType::active_cell_iterator;
93 # else
95 # endif
96  };
97 
98 # ifdef _MSC_VER
99  template <int dim, int spacedim>
100  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
101  {
102  public:
103  using type =
105  };
106 # endif
107 
108 
109 # ifdef _MSC_VER
110  template <int dim, int spacedim>
111  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim>>
112  {
113  public:
114  using type =
116  };
117 # endif
118 } // namespace internal
119 
128 namespace GridTools
129 {
130  template <int dim, int spacedim>
131  class Cache;
132 
137 
144  template <int dim, int spacedim>
145  double
147 
174  template <int dim, int spacedim>
175  double
177  const Mapping<dim, spacedim> & mapping =
179 
190  template <int dim, int spacedim>
191  double
193  const Mapping<dim, spacedim> & mapping =
195 
206  template <int dim, int spacedim>
207  double
209  const Mapping<dim, spacedim> & mapping =
211 
221  template <int dim>
222  double
223  cell_measure(
224  const std::vector<Point<dim>> &all_vertices,
226 
236  template <int dim>
237  double
238  cell_measure(const std::vector<Point<dim>> & all_vertices,
240 
246  template <int dim, typename T>
247  double
248  cell_measure(const T &, ...);
249 
276  template <int dim>
279  const Triangulation<dim> &triangulation,
280  const Quadrature<dim> & quadrature);
281 
289  template <int dim>
290  double
292  const Triangulation<dim> &triangulation,
293  const Quadrature<dim> & quadrature);
294 
308  template <int dim, int spacedim>
311 
329  template <typename Iterator>
332  const Iterator & object,
334 
347  template <int dim, int spacedim>
348  std::
349  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
351 
357 
374  template <int dim, int spacedim>
375  void
377  std::vector<CellData<dim>> & cells,
378  SubCellData & subcelldata);
379 
397  template <int dim, int spacedim>
398  void
399  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
400  std::vector<CellData<dim>> & cells,
401  SubCellData & subcelldata,
402  std::vector<unsigned int> & considered_vertices,
403  const double tol = 1e-12);
404 
410 
464  template <int dim, typename Transformation, int spacedim>
465  void
466  transform(const Transformation & transformation,
467  Triangulation<dim, spacedim> &triangulation);
468 
474  template <int dim, int spacedim>
475  void
476  shift(const Tensor<1, spacedim> & shift_vector,
477  Triangulation<dim, spacedim> &triangulation);
478 
479 
489  template <int dim>
490  void
491  rotate(const double angle, Triangulation<dim> &triangulation);
492 
505  template <int dim>
506  void
507  rotate(const double angle,
508  const unsigned int axis,
509  Triangulation<dim, 3> &triangulation);
510 
568  template <int dim>
569  void
570  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
571  Triangulation<dim> & tria,
572  const Function<dim, double> *coefficient = nullptr,
573  const bool solve_for_absolute_positions = false);
574 
580  template <int dim, int spacedim>
581  std::map<unsigned int, Point<spacedim>>
583 
591  template <int dim, int spacedim>
592  void
593  scale(const double scaling_factor,
594  Triangulation<dim, spacedim> &triangulation);
595 
606  template <int dim, int spacedim>
607  void
608  distort_random(const double factor,
609  Triangulation<dim, spacedim> &triangulation,
610  const bool keep_boundary = true);
611 
645  template <int dim, int spacedim>
646  void
648  const bool isotropic = false,
649  const unsigned int max_iterations = 100);
650 
675  template <int dim, int spacedim>
676  void
678  const double max_ratio = 1.6180339887,
679  const unsigned int max_iterations = 5);
680 
770  template <int dim, int spacedim>
771  void
773  const double limit_angle_fraction = .75);
774 
780 
830  template <int dim, int spacedim>
831 # ifndef DOXYGEN
832  std::tuple<
833  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
834  std::vector<std::vector<Point<dim>>>,
835  std::vector<std::vector<unsigned int>>>
836 # else
837  return_type
838 # endif
840  const Cache<dim, spacedim> & cache,
841  const std::vector<Point<spacedim>> &points,
843  &cell_hint =
845 
874  template <int dim, int spacedim>
875 # ifndef DOXYGEN
876  std::tuple<
877  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
878  std::vector<std::vector<Point<dim>>>,
879  std::vector<std::vector<unsigned int>>,
880  std::vector<unsigned int>>
881 # else
882  return_type
883 # endif
885  const Cache<dim, spacedim> & cache,
886  const std::vector<Point<spacedim>> &points,
888  &cell_hint =
890 
952  template <int dim, int spacedim>
953 # ifndef DOXYGEN
954  std::tuple<
955  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
956  std::vector<std::vector<Point<dim>>>,
957  std::vector<std::vector<unsigned int>>,
958  std::vector<std::vector<Point<spacedim>>>,
959  std::vector<std::vector<unsigned int>>>
960 # else
961  return_type
962 # endif
964  const GridTools::Cache<dim, spacedim> & cache,
965  const std::vector<Point<spacedim>> & local_points,
966  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes);
967 
1004  template <int dim, int spacedim>
1005  std::map<unsigned int, Point<spacedim>>
1007  const Mapping<dim, spacedim> & mapping =
1009 
1019  template <int spacedim>
1020  unsigned int
1021  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1022  const Point<spacedim> & p);
1023 
1047  template <int dim, template <int, int> class MeshType, int spacedim>
1048  unsigned int
1049  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1050  const Point<spacedim> & p,
1051  const std::vector<bool> & marked_vertices = {});
1052 
1076  template <int dim, template <int, int> class MeshType, int spacedim>
1077  unsigned int
1079  const MeshType<dim, spacedim> &mesh,
1080  const Point<spacedim> & p,
1081  const std::vector<bool> & marked_vertices = {});
1082 
1083 
1108  template <int dim, template <int, int> class MeshType, int spacedim>
1109 # ifndef _MSC_VER
1110  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1111 # else
1112  std::vector<
1113  typename ::internal::
1114  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1115 # endif
1116  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1117  const unsigned int vertex_index);
1118 
1186  template <int dim, template <int, int> class MeshType, int spacedim>
1187 # ifndef _MSC_VER
1188  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1189 # else
1190  std::pair<typename ::internal::
1191  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1192  Point<dim>>
1193 # endif
1195  const MeshType<dim, spacedim> &mesh,
1196  const Point<spacedim> & p,
1197  const std::vector<bool> &marked_vertices = {},
1198  const double tolerance = 1.e-10);
1199 
1207  template <int dim, template <int, int> class MeshType, int spacedim>
1208 # ifndef _MSC_VER
1209  typename MeshType<dim, spacedim>::active_cell_iterator
1210 # else
1211  typename ::internal::
1212  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1213 # endif
1214  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1215  const Point<spacedim> & p,
1216  const std::vector<bool> &marked_vertices = {},
1217  const double tolerance = 1.e-10);
1218 
1225  template <int dim, int spacedim>
1226  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1227  Point<dim>>
1229  const hp::MappingCollection<dim, spacedim> &mapping,
1230  const DoFHandler<dim, spacedim> & mesh,
1231  const Point<spacedim> & p,
1232  const double tolerance = 1.e-10);
1233 
1282  template <int dim, int spacedim>
1283  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1284  Point<dim>>
1286  const Cache<dim, spacedim> &cache,
1287  const Point<spacedim> & p,
1290  const std::vector<bool> &marked_vertices = {},
1291  const double tolerance = 1.e-10);
1292 
1306  template <int dim, template <int, int> class MeshType, int spacedim>
1307 # ifndef _MSC_VER
1308  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1309 # else
1310  std::pair<typename ::internal::
1311  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1312  Point<dim>>
1313 # endif
1315  const Mapping<dim, spacedim> & mapping,
1316  const MeshType<dim, spacedim> &mesh,
1317  const Point<spacedim> & p,
1318  const std::vector<
1319  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1321  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1322  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1323  typename MeshType<dim, spacedim>::active_cell_iterator(),
1324  const std::vector<bool> & marked_vertices = {},
1325  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1326  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1327  const double tolerance = 1.e-10);
1328 
1349  template <int dim, template <int, int> class MeshType, int spacedim>
1350 # ifndef _MSC_VER
1351  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1352  Point<dim>>>
1353 # else
1354  std::vector<std::pair<
1355  typename ::internal::
1356  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1357  Point<dim>>>
1358 # endif
1360  const Mapping<dim, spacedim> & mapping,
1361  const MeshType<dim, spacedim> &mesh,
1362  const Point<spacedim> & p,
1363  const double tolerance,
1364  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1365  Point<dim>> & first_cell);
1366 
1373  template <int dim, template <int, int> class MeshType, int spacedim>
1374 # ifndef _MSC_VER
1375  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1376  Point<dim>>>
1377 # else
1378  std::vector<std::pair<
1379  typename ::internal::
1380  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1381  Point<dim>>>
1382 # endif
1384  const Mapping<dim, spacedim> & mapping,
1385  const MeshType<dim, spacedim> &mesh,
1386  const Point<spacedim> & p,
1387  const double tolerance = 1e-10,
1388  const std::vector<bool> & marked_vertices = {});
1389 
1411  template <class MeshType>
1412  std::vector<typename MeshType::active_cell_iterator>
1413  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1414 
1439  template <class MeshType>
1440  void
1442  const typename MeshType::active_cell_iterator & cell,
1443  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1444 
1495  template <class MeshType>
1496  std::vector<typename MeshType::active_cell_iterator>
1498  const MeshType &mesh,
1499  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1500  &predicate);
1501 
1502 
1510  template <class MeshType>
1511  std::vector<typename MeshType::cell_iterator>
1513  const MeshType &mesh,
1514  const std::function<bool(const typename MeshType::cell_iterator &)>
1515  & predicate,
1516  const unsigned int level);
1517 
1518 
1532  template <class MeshType>
1533  std::vector<typename MeshType::active_cell_iterator>
1534  compute_ghost_cell_halo_layer(const MeshType &mesh);
1535 
1585  template <class MeshType>
1586  std::vector<typename MeshType::active_cell_iterator>
1588  const MeshType &mesh,
1589  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1590  & predicate,
1591  const double layer_thickness);
1592 
1616  template <class MeshType>
1617  std::vector<typename MeshType::active_cell_iterator>
1618  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1619  const double layer_thickness);
1620 
1636  template <class MeshType>
1637  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1639  const MeshType &mesh,
1640  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1641  &predicate);
1642 
1695  template <class MeshType>
1696  std::vector<BoundingBox<MeshType::space_dimension>>
1698  const MeshType &mesh,
1699  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1700  & predicate,
1701  const unsigned int refinement_level = 0,
1702  const bool allow_merge = false,
1703  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1704 
1732  template <int spacedim>
1733 # ifndef DOXYGEN
1734  std::tuple<std::vector<std::vector<unsigned int>>,
1735  std::map<unsigned int, unsigned int>,
1736  std::map<unsigned int, std::vector<unsigned int>>>
1737 # else
1738  return_type
1739 # endif
1741  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1742  const std::vector<Point<spacedim>> & points);
1743 
1744 
1779  template <int spacedim>
1780 # ifndef DOXYGEN
1781  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1782  std::map<unsigned int, unsigned int>,
1783  std::map<unsigned int, std::vector<unsigned int>>>
1784 # else
1785  return_type
1786 # endif
1788  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1789  const std::vector<Point<spacedim>> & points);
1790 
1791 
1800  template <int dim, int spacedim>
1801  std::vector<
1802  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1803  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1804 
1817  template <int dim, int spacedim>
1818  std::vector<std::vector<Tensor<1, spacedim>>>
1820  const Triangulation<dim, spacedim> &mesh,
1821  const std::vector<
1823  &vertex_to_cells);
1824 
1825 
1833  template <int dim, int spacedim>
1834  unsigned int
1837  const Point<spacedim> & position,
1838  const Mapping<dim, spacedim> & mapping =
1840 
1852  template <int dim, int spacedim>
1853  std::map<unsigned int, types::global_vertex_index>
1856 
1868  template <int dim, int spacedim>
1869  std::pair<unsigned int, double>
1872 
1878 
1887  template <int dim, int spacedim>
1888  void
1890  const Triangulation<dim, spacedim> &triangulation,
1891  DynamicSparsityPattern & connectivity);
1892 
1901  template <int dim, int spacedim>
1902  void
1904  const Triangulation<dim, spacedim> &triangulation,
1905  DynamicSparsityPattern & connectivity);
1906 
1915  template <int dim, int spacedim>
1916  void
1918  const Triangulation<dim, spacedim> &triangulation,
1919  const unsigned int level,
1920  DynamicSparsityPattern & connectivity);
1921 
1942  template <int dim, int spacedim>
1943  void
1944  partition_triangulation(const unsigned int n_partitions,
1945  Triangulation<dim, spacedim> & triangulation,
1946  const SparsityTools::Partitioner partitioner =
1948 
1959  template <int dim, int spacedim>
1960  void
1961  partition_triangulation(const unsigned int n_partitions,
1962  const std::vector<unsigned int> &cell_weights,
1963  Triangulation<dim, spacedim> & triangulation,
1964  const SparsityTools::Partitioner partitioner =
1966 
2012  template <int dim, int spacedim>
2013  void
2014  partition_triangulation(const unsigned int n_partitions,
2015  const SparsityPattern & cell_connection_graph,
2016  Triangulation<dim, spacedim> &triangulation,
2017  const SparsityTools::Partitioner partitioner =
2019 
2030  template <int dim, int spacedim>
2031  void
2032  partition_triangulation(const unsigned int n_partitions,
2033  const std::vector<unsigned int> &cell_weights,
2034  const SparsityPattern & cell_connection_graph,
2035  Triangulation<dim, spacedim> &triangulation,
2036  const SparsityTools::Partitioner partitioner =
2038 
2053  template <int dim, int spacedim>
2054  void
2055  partition_triangulation_zorder(const unsigned int n_partitions,
2056  Triangulation<dim, spacedim> &triangulation,
2057  const bool group_siblings = true);
2058 
2070  template <int dim, int spacedim>
2071  void
2073 
2084  template <int dim, int spacedim>
2085  void
2087  std::vector<types::subdomain_id> & subdomain);
2088 
2103  template <int dim, int spacedim>
2104  unsigned int
2106  const Triangulation<dim, spacedim> &triangulation,
2107  const types::subdomain_id subdomain);
2108 
2109 
2139  template <int dim, int spacedim>
2140  std::vector<bool>
2142 
2148 
2182  template <typename MeshType>
2183  std::list<std::pair<typename MeshType::cell_iterator,
2184  typename MeshType::cell_iterator>>
2185  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2186 
2196  template <int dim, int spacedim>
2197  bool
2199  const Triangulation<dim, spacedim> &mesh_2);
2200 
2210  template <typename MeshType>
2211  bool
2212  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2213 
2219 
2235  template <int dim, int spacedim>
2239  & distorted_cells,
2240  Triangulation<dim, spacedim> &triangulation);
2241 
2242 
2243 
2252 
2253 
2291  template <class MeshType>
2292  std::vector<typename MeshType::active_cell_iterator>
2293  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2294 
2295 
2317  template <class Container>
2318  std::vector<typename Container::cell_iterator>
2320  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2321 
2388  template <class Container>
2389  void
2391  const std::vector<typename Container::active_cell_iterator> &patch,
2393  &local_triangulation,
2394  std::map<
2395  typename Triangulation<Container::dimension,
2396  Container::space_dimension>::active_cell_iterator,
2397  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2398 
2430  template <int dim, int spacedim>
2431  std::map<
2433  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2435 
2436 
2443 
2449  template <typename CellIterator>
2451  {
2455  CellIterator cell[2];
2456 
2461  unsigned int face_idx[2];
2462 
2468  std::bitset<3> orientation;
2469 
2483  };
2484 
2485 
2549  template <typename FaceIterator>
2550  bool orthogonal_equality(
2551  std::bitset<3> & orientation,
2552  const FaceIterator & face1,
2553  const FaceIterator & face2,
2554  const int direction,
2558 
2559 
2563  template <typename FaceIterator>
2564  bool
2566  const FaceIterator & face1,
2567  const FaceIterator & face2,
2568  const int direction,
2572 
2573 
2630  template <typename MeshType>
2631  void
2633  const MeshType & mesh,
2634  const types::boundary_id b_id1,
2635  const types::boundary_id b_id2,
2636  const int direction,
2638  & matched_pairs,
2639  const Tensor<1, MeshType::space_dimension> &offset =
2642 
2643 
2666  template <typename MeshType>
2667  void
2669  const MeshType & mesh,
2670  const types::boundary_id b_id,
2671  const int direction,
2673  & matched_pairs,
2674  const ::Tensor<1, MeshType::space_dimension> &offset =
2677 
2683 
2704  template <int dim, int spacedim>
2705  void
2707  const bool reset_boundary_ids = false);
2708 
2730  template <int dim, int spacedim>
2731  void
2733  const std::vector<types::boundary_id> &src_boundary_ids,
2734  const std::vector<types::manifold_id> &dst_manifold_ids,
2736  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2737 
2767  template <int dim, int spacedim>
2768  void
2770  const bool compute_face_ids = false);
2771 
2796  template <int dim, int spacedim>
2797  void
2800  const std::function<types::manifold_id(
2801  const std::set<types::manifold_id> &)> &disambiguation_function =
2802  [](const std::set<types::manifold_id> &manifold_ids) {
2803  if (manifold_ids.size() == 1)
2804  return *manifold_ids.begin();
2805  else
2807  },
2808  bool overwrite_only_flat_manifold_ids = true);
2895  template <typename DataType, typename MeshType>
2896  void
2898  const MeshType & mesh,
2899  const std::function<std_cxx17::optional<DataType>(
2900  const typename MeshType::active_cell_iterator &)> &pack,
2901  const std::function<void(const typename MeshType::active_cell_iterator &,
2902  const DataType &)> & unpack,
2903  const std::function<bool(const typename MeshType::active_cell_iterator &)>
2904  &cell_filter =
2905  [](const typename MeshType::active_cell_iterator &) { return true; });
2906 
2917  template <typename DataType, typename MeshType>
2918  void
2920  const MeshType & mesh,
2921  const std::function<std_cxx17::optional<DataType>(
2922  const typename MeshType::level_cell_iterator &)> &pack,
2923  const std::function<void(const typename MeshType::level_cell_iterator &,
2924  const DataType &)> & unpack,
2925  const std::function<bool(const typename MeshType::level_cell_iterator &)>
2926  &cell_filter =
2927  [](const typename MeshType::level_cell_iterator &) { return true; });
2928 
2929  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2930  * boxes @p local_bboxes.
2931  *
2932  * This function is meant to exchange bounding boxes describing the locally
2933  * owned cells in a distributed triangulation obtained with the function
2934  * GridTools::compute_mesh_predicate_bounding_box .
2935  *
2936  * The output vector's size is the number of processes of the MPI
2937  * communicator:
2938  * its i-th entry contains the vector @p local_bboxes of the i-th process.
2939  */
2940  template <int spacedim>
2941  std::vector<std::vector<BoundingBox<spacedim>>>
2943  const std::vector<BoundingBox<spacedim>> &local_bboxes,
2944  MPI_Comm mpi_communicator);
2945 
2978  template <int spacedim>
2981  const std::vector<BoundingBox<spacedim>> &local_description,
2982  MPI_Comm mpi_communicator);
2983 
3001  template <int dim, int spacedim>
3002  void
3004  const Triangulation<dim, spacedim> & tria,
3005  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3006  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3007 
3014  template <int dim, int spacedim>
3015  std::map<unsigned int, std::set<::types::subdomain_id>>
3017  const Triangulation<dim, spacedim> &tria);
3018 
3031  template <int dim, typename T>
3033  {
3037  std::vector<CellId> cell_ids;
3038 
3042  std::vector<T> data;
3043 
3051  template <class Archive>
3052  void
3053  save(Archive &ar, const unsigned int version) const;
3054 
3059  template <class Archive>
3060  void
3061  load(Archive &ar, const unsigned int version);
3062 
3063 # ifdef DOXYGEN
3064 
3068  template <class Archive>
3069  void
3070  serialize(Archive &archive, const unsigned int version);
3071 # else
3072  // This macro defines the serialize() method that is compatible with
3073  // the templated save() and load() method that have been implemented.
3074  BOOST_SERIALIZATION_SPLIT_MEMBER()
3075 # endif
3076  };
3077 
3082 
3087  int,
3088  << "The number of partitions you gave is " << arg1
3089  << ", but must be greater than zero.");
3094  int,
3095  << "The subdomain id " << arg1
3096  << " has no cells associated with it.");
3101 
3106  double,
3107  << "The scaling factor must be positive, but it is " << arg1
3108  << ".");
3112  template <int N>
3114  Point<N>,
3115  << "The point <" << arg1
3116  << "> could not be found inside any of the "
3117  << "coarse grid cells.");
3121  template <int N>
3123  Point<N>,
3124  << "The point <" << arg1
3125  << "> could not be found inside any of the "
3126  << "subcells of a coarse grid cell.");
3127 
3132  unsigned int,
3133  << "The given vertex with index " << arg1
3134  << " is not used in the given triangulation.");
3135 
3136 
3139 } /*namespace GridTools*/
3140 
3141 
3142 
3143 /* ----------------- Template function --------------- */
3144 
3145 # ifndef DOXYGEN
3146 
3147 namespace GridTools
3148 {
3149  template <int dim>
3150  double
3151  cell_measure(
3152  const std::vector<Point<dim>> &all_vertices,
3153  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3154  {
3155  // We forward call to the ArrayView version:
3156  const ArrayView<const unsigned int> view(
3157  indices, GeometryInfo<dim>::vertices_per_cell);
3158  return cell_measure(all_vertices, view);
3159  }
3160 
3161  template <int dim, typename T>
3162  double
3163  cell_measure(const T &, ...)
3164  {
3165  Assert(false, ExcNotImplemented());
3166  return std::numeric_limits<double>::quiet_NaN();
3167  }
3168 
3169 
3170 
3171  // This specialization is defined here so that the general template in the
3172  // source file doesn't need to have further 1D overloads for the internal
3173  // functions it calls.
3174  template <>
3178  {
3179  return {};
3180  }
3181 
3182 
3183 
3184  template <int dim, typename Predicate, int spacedim>
3185  void
3186  transform(const Predicate & predicate,
3187  Triangulation<dim, spacedim> &triangulation)
3188  {
3189  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3190 
3191  // loop over all active cells, and
3192  // transform those vertices that
3193  // have not yet been touched. note
3194  // that we get to all vertices in
3195  // the triangulation by only
3196  // visiting the active cells.
3198  cell = triangulation.begin_active(),
3199  endc = triangulation.end();
3200  for (; cell != endc; ++cell)
3201  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
3202  if (treated_vertices[cell->vertex_index(v)] == false)
3203  {
3204  // transform this vertex
3205  cell->vertex(v) = predicate(cell->vertex(v));
3206  // and mark it as treated
3207  treated_vertices[cell->vertex_index(v)] = true;
3208  };
3209 
3210 
3211  // now fix any vertices on hanging nodes so that we don't create any holes
3212  if (dim == 2)
3213  {
3215  cell = triangulation.begin_active(),
3216  endc = triangulation.end();
3217  for (; cell != endc; ++cell)
3218  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3219  if (cell->face(face)->has_children() &&
3220  !cell->face(face)->at_boundary())
3221  {
3222  // this line has children
3223  cell->face(face)->child(0)->vertex(1) =
3224  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3225  2;
3226  }
3227  }
3228  else if (dim == 3)
3229  {
3231  cell = triangulation.begin_active(),
3232  endc = triangulation.end();
3233  for (; cell != endc; ++cell)
3234  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3235  if (cell->face(face)->has_children() &&
3236  !cell->face(face)->at_boundary())
3237  {
3238  // this face has hanging nodes
3239  cell->face(face)->child(0)->vertex(1) =
3240  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3241  2.0;
3242  cell->face(face)->child(0)->vertex(2) =
3243  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3244  2.0;
3245  cell->face(face)->child(1)->vertex(3) =
3246  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3247  2.0;
3248  cell->face(face)->child(2)->vertex(3) =
3249  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3250  2.0;
3251 
3252  // center of the face
3253  cell->face(face)->child(0)->vertex(3) =
3254  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3255  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3256  4.0;
3257  }
3258  }
3259 
3260  // Make sure FEValues notices that the mesh has changed
3261  triangulation.signals.mesh_movement();
3262  }
3263 
3264 
3265 
3266  template <class MeshType>
3267  std::vector<typename MeshType::active_cell_iterator>
3268  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3269  {
3270  std::vector<typename MeshType::active_cell_iterator> child_cells;
3271 
3272  if (cell->has_children())
3273  {
3274  for (unsigned int child = 0; child < cell->n_children(); ++child)
3275  if (cell->child(child)->has_children())
3276  {
3277  const std::vector<typename MeshType::active_cell_iterator>
3278  children = get_active_child_cells<MeshType>(cell->child(child));
3279  child_cells.insert(child_cells.end(),
3280  children.begin(),
3281  children.end());
3282  }
3283  else
3284  child_cells.push_back(cell->child(child));
3285  }
3286 
3287  return child_cells;
3288  }
3289 
3290 
3291 
3292  template <class MeshType>
3293  void
3295  const typename MeshType::active_cell_iterator & cell,
3296  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3297  {
3298  active_neighbors.clear();
3299  for (const unsigned int n :
3301  if (!cell->at_boundary(n))
3302  {
3303  if (MeshType::dimension == 1)
3304  {
3305  // check children of neighbor. note
3306  // that in 1d children of the neighbor
3307  // may be further refined. In 1d the
3308  // case is simple since we know what
3309  // children bound to the present cell
3310  typename MeshType::cell_iterator neighbor_child =
3311  cell->neighbor(n);
3312  if (!neighbor_child->is_active())
3313  {
3314  while (neighbor_child->has_children())
3315  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3316 
3317  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3318  ExcInternalError());
3319  }
3320  active_neighbors.push_back(neighbor_child);
3321  }
3322  else
3323  {
3324  if (cell->face(n)->has_children())
3325  // this neighbor has children. find
3326  // out which border to the present
3327  // cell
3328  for (unsigned int c = 0;
3329  c < cell->face(n)->number_of_children();
3330  ++c)
3331  active_neighbors.push_back(
3332  cell->neighbor_child_on_subface(n, c));
3333  else
3334  {
3335  // the neighbor must be active
3336  // himself
3337  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3338  active_neighbors.push_back(cell->neighbor(n));
3339  }
3340  }
3341  }
3342  }
3343 
3344 
3345 
3346  namespace internal
3347  {
3348  namespace ProjectToObject
3349  {
3362  struct CrossDerivative
3363  {
3364  const unsigned int direction_0;
3365  const unsigned int direction_1;
3366 
3367  CrossDerivative(const unsigned int d0, const unsigned int d1);
3368  };
3369 
3370  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3371  const unsigned int d1)
3372  : direction_0(d0)
3373  , direction_1(d1)
3374  {}
3375 
3376 
3377 
3382  template <typename F>
3383  inline auto
3384  centered_first_difference(const double center,
3385  const double step,
3386  const F &f) -> decltype(f(center) - f(center))
3387  {
3388  return (f(center + step) - f(center - step)) / (2.0 * step);
3389  }
3390 
3391 
3392 
3397  template <typename F>
3398  inline auto
3399  centered_second_difference(const double center,
3400  const double step,
3401  const F &f) -> decltype(f(center) - f(center))
3402  {
3403  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3404  (step * step);
3405  }
3406 
3407 
3408 
3418  template <int structdim, typename F>
3419  inline auto
3420  cross_stencil(
3421  const CrossDerivative cross_derivative,
3423  const double step,
3424  const F &f) -> decltype(f(center) - f(center))
3425  {
3427  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3428  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3429  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3430  1.0 / 3.0 * f(center - simplex_vector) +
3431  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3432  step;
3433  }
3434 
3435 
3436 
3443  template <int spacedim, int structdim, typename F>
3444  inline double
3445  gradient_entry(
3446  const unsigned int row_n,
3447  const unsigned int dependent_direction,
3448  const Point<spacedim> &p0,
3450  const double step,
3451  const F & f)
3452  {
3454  dependent_direction <
3456  ExcMessage("This function assumes that the last weight is a "
3457  "dependent variable (and hence we cannot take its "
3458  "derivative directly)."));
3459  Assert(row_n != dependent_direction,
3460  ExcMessage(
3461  "We cannot differentiate with respect to the variable "
3462  "that is assumed to be dependent."));
3463 
3464  const Point<spacedim> manifold_point = f(center);
3465  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3466  {row_n, dependent_direction}, center, step, f);
3467  double entry = 0.0;
3468  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3469  entry +=
3470  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3471  return entry;
3472  }
3473 
3479  template <typename Iterator, int spacedim, int structdim>
3481  project_to_d_linear_object(const Iterator & object,
3482  const Point<spacedim> &trial_point)
3483  {
3484  // let's look at this for simplicity for a quad (structdim==2) in a
3485  // space with spacedim>2 (notate trial_point by y): all points on the
3486  // surface are given by
3487  // x(\xi) = sum_i v_i phi_x(\xi)
3488  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3489  // reference coordinates of the quad. so what we are trying to do is
3490  // find a point x on the surface that is closest to the point y. there
3491  // are different ways to solve this problem, but in the end it's a
3492  // nonlinear problem and we have to find reference coordinates \xi so
3493  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3494  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3495  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3496  // have to use a Newton method to find the answer. This leads to the
3497  // following formulation of Newton steps:
3498  //
3499  // Given \xi_k, find \delta\xi_k so that
3500  // H_k \delta\xi_k = - F_k
3501  // where H_k is an approximation to the second derivatives of J at
3502  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3503  // number of times until the right hand side is small enough. As a
3504  // stopping criterion, we terminate if ||\delta\xi||<eps.
3505  //
3506  // As for the Hessian, the best choice would be
3507  // H_k = J''(\xi_k)
3508  // but we'll opt for the simpler Gauss-Newton form
3509  // H_k = A^T A
3510  // i.e.
3511  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3512  // \partial_n phi_i *
3513  // \partial_m phi_j
3514  // we start at xi=(0.5, 0.5).
3515  Point<structdim> xi;
3516  for (unsigned int d = 0; d < structdim; ++d)
3517  xi[d] = 0.5;
3518 
3519  Point<spacedim> x_k;
3520  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3521  x_k += object->vertex(i) *
3523 
3524  do
3525  {
3527  for (const unsigned int i :
3529  F_k +=
3530  (x_k - trial_point) * object->vertex(i) *
3532  i);
3533 
3535  for (const unsigned int i :
3537  for (const unsigned int j :
3539  {
3542  xi, i),
3544  xi, j));
3545  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3546  }
3547 
3548  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3549  xi += delta_xi;
3550 
3551  x_k = Point<spacedim>();
3552  for (const unsigned int i :
3554  x_k += object->vertex(i) *
3556 
3557  if (delta_xi.norm() < 1e-7)
3558  break;
3559  }
3560  while (true);
3561 
3562  return x_k;
3563  }
3564  } // namespace ProjectToObject
3565  } // namespace internal
3566 
3567 
3568 
3569  namespace internal
3570  {
3571  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3572  // inside the project_to_object function below.
3573  template <int structdim>
3574  inline bool
3575  weights_are_ok(
3577  {
3578  // clang has trouble figuring out structdim here, so define it
3579  // again:
3580  static const std::size_t n_vertices_per_cell =
3582  n_independent_components;
3583  std::array<double, n_vertices_per_cell> copied_weights;
3584  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3585  {
3586  copied_weights[i] = v[i];
3587  if (v[i] < 0.0 || v[i] > 1.0)
3588  return false;
3589  }
3590 
3591  // check the sum: try to avoid some roundoff errors by summing in order
3592  std::sort(copied_weights.begin(), copied_weights.end());
3593  const double sum =
3594  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3595  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3596  }
3597  } // namespace internal
3598 
3599  template <typename Iterator>
3602  const Iterator & object,
3604  {
3605  const int spacedim = Iterator::AccessorType::space_dimension;
3606  const int structdim = Iterator::AccessorType::structure_dimension;
3607 
3608  Point<spacedim> projected_point = trial_point;
3609 
3610  if (structdim >= spacedim)
3611  return projected_point;
3612  else if (structdim == 1 || structdim == 2)
3613  {
3614  using namespace internal::ProjectToObject;
3615  // Try to use the special flat algorithm for quads (this is better
3616  // than the general algorithm in 3D). This does not take into account
3617  // whether projected_point is outside the quad, but we optimize along
3618  // lines below anyway:
3619  const int dim = Iterator::AccessorType::dimension;
3620  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3621  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3622  &manifold) != nullptr)
3623  {
3624  projected_point =
3625  project_to_d_linear_object<Iterator, spacedim, structdim>(
3626  object, trial_point);
3627  }
3628  else
3629  {
3630  // We want to find a point on the convex hull (defined by the
3631  // vertices of the object and the manifold description) that is
3632  // relatively close to the trial point. This has a few issues:
3633  //
3634  // 1. For a general convex hull we are not guaranteed that a unique
3635  // minimum exists.
3636  // 2. The independent variables in the optimization process are the
3637  // weights given to Manifold::get_new_point, which must sum to 1,
3638  // so we cannot use standard finite differences to approximate a
3639  // gradient.
3640  //
3641  // There is not much we can do about 1., but for 2. we can derive
3642  // finite difference stencils that work on a structdim-dimensional
3643  // simplex and rewrite the optimization problem to use those
3644  // instead. Consider the structdim 2 case and let
3645  //
3646  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3647  // c2, c3})
3648  //
3649  // where {c0, c1, c2, c3} are the weights for the four vertices on
3650  // the quadrilateral. We seek to minimize the Euclidean distance
3651  // between F(...) and trial_point. We can solve for c3 in terms of
3652  // the other weights and get, for one coordinate direction
3653  //
3654  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3655  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3656  //
3657  // where we substitute back in for c3 after taking the
3658  // derivative. We can compute a stencil for the cross derivative
3659  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3660  // (and gradient_entry computes the sum over the independent
3661  // variables). Below, we somewhat arbitrarily pick the last
3662  // component as the dependent one.
3663  //
3664  // Since we can now calculate derivatives of the objective
3665  // function we can use gradient descent to minimize it.
3666  //
3667  // Of course, this is much simpler in the structdim = 1 case (we
3668  // could rewrite the projection as a 1D optimization problem), but
3669  // to reduce the potential for bugs we use the same code in both
3670  // cases.
3671  const double step_size = object->diameter() / 64.0;
3672 
3673  constexpr unsigned int n_vertices_per_cell =
3675 
3676  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3677  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3678  ++vertex_n)
3679  vertices[vertex_n] = object->vertex(vertex_n);
3680 
3681  auto get_point_from_weights =
3682  [&](const Tensor<1, n_vertices_per_cell> &weights)
3683  -> Point<spacedim> {
3684  return object->get_manifold().get_new_point(
3685  make_array_view(vertices.begin(), vertices.end()),
3686  make_array_view(weights.begin_raw(), weights.end_raw()));
3687  };
3688 
3689  // pick the initial weights as (normalized) inverse distances from
3690  // the trial point:
3691  Tensor<1, n_vertices_per_cell> guess_weights;
3692  double guess_weights_sum = 0.0;
3693  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3694  ++vertex_n)
3695  {
3696  const double distance =
3697  vertices[vertex_n].distance(trial_point);
3698  if (distance == 0.0)
3699  {
3700  guess_weights = 0.0;
3701  guess_weights[vertex_n] = 1.0;
3702  guess_weights_sum = 1.0;
3703  break;
3704  }
3705  else
3706  {
3707  guess_weights[vertex_n] = 1.0 / distance;
3708  guess_weights_sum += guess_weights[vertex_n];
3709  }
3710  }
3711  guess_weights /= guess_weights_sum;
3712  Assert(internal::weights_are_ok<structdim>(guess_weights),
3713  ExcInternalError());
3714 
3715  // The optimization algorithm consists of two parts:
3716  //
3717  // 1. An outer loop where we apply the gradient descent algorithm.
3718  // 2. An inner loop where we do a line search to find the optimal
3719  // length of the step one should take in the gradient direction.
3720  //
3721  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3722  {
3723  const unsigned int dependent_direction =
3724  n_vertices_per_cell - 1;
3725  Tensor<1, n_vertices_per_cell> current_gradient;
3726  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3727  ++row_n)
3728  {
3729  if (row_n != dependent_direction)
3730  {
3731  current_gradient[row_n] =
3732  gradient_entry<spacedim, structdim>(
3733  row_n,
3734  dependent_direction,
3735  trial_point,
3736  guess_weights,
3737  step_size,
3738  get_point_from_weights);
3739 
3740  current_gradient[dependent_direction] -=
3741  current_gradient[row_n];
3742  }
3743  }
3744 
3745  // We need to travel in the -gradient direction, as noted
3746  // above, but we may not want to take a full step in that
3747  // direction; instead, guess that we will go -0.5*gradient and
3748  // do quasi-Newton iteration to pick the best multiplier. The
3749  // goal is to find a scalar alpha such that
3750  //
3751  // F(x - alpha g)
3752  //
3753  // is minimized, where g is the gradient and F is the
3754  // objective function. To find the optimal value we find roots
3755  // of the derivative of the objective function with respect to
3756  // alpha by Newton iteration, where we approximate the first
3757  // and second derivatives of F(x - alpha g) with centered
3758  // finite differences.
3759  double gradient_weight = -0.5;
3760  auto gradient_weight_objective_function =
3761  [&](const double gradient_weight_guess) -> double {
3762  return (trial_point -
3763  get_point_from_weights(guess_weights +
3764  gradient_weight_guess *
3765  current_gradient))
3766  .norm_square();
3767  };
3768 
3769  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3770  {
3771  const double update_numerator = centered_first_difference(
3772  gradient_weight,
3773  step_size,
3774  gradient_weight_objective_function);
3775  const double update_denominator =
3776  centered_second_difference(
3777  gradient_weight,
3778  step_size,
3779  gradient_weight_objective_function);
3780 
3781  // avoid division by zero. Note that we limit the gradient
3782  // weight below
3783  if (std::abs(update_denominator) == 0.0)
3784  break;
3785  gradient_weight =
3786  gradient_weight - update_numerator / update_denominator;
3787 
3788  // Put a fairly lenient bound on the largest possible
3789  // gradient (things tend to be locally flat, so the gradient
3790  // itself is usually small)
3791  if (std::abs(gradient_weight) > 10)
3792  {
3793  gradient_weight = -10.0;
3794  break;
3795  }
3796  }
3797 
3798  // It only makes sense to take convex combinations with weights
3799  // between zero and one. If the update takes us outside of this
3800  // region then rescale the update to stay within the region and
3801  // try again
3802  Tensor<1, n_vertices_per_cell> tentative_weights =
3803  guess_weights + gradient_weight * current_gradient;
3804 
3805  double new_gradient_weight = gradient_weight;
3806  for (unsigned int iteration_count = 0; iteration_count < 40;
3807  ++iteration_count)
3808  {
3809  if (internal::weights_are_ok<structdim>(tentative_weights))
3810  break;
3811 
3812  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3813  {
3814  if (tentative_weights[i] < 0.0)
3815  {
3816  tentative_weights -=
3817  (tentative_weights[i] / current_gradient[i]) *
3818  current_gradient;
3819  }
3820  if (tentative_weights[i] < 0.0 ||
3821  1.0 < tentative_weights[i])
3822  {
3823  new_gradient_weight /= 2.0;
3824  tentative_weights =
3825  guess_weights +
3826  new_gradient_weight * current_gradient;
3827  }
3828  }
3829  }
3830 
3831  // the update might still send us outside the valid region, so
3832  // check again and quit if the update is still not valid
3833  if (!internal::weights_are_ok<structdim>(tentative_weights))
3834  break;
3835 
3836  // if we cannot get closer by traveling in the gradient
3837  // direction then quit
3838  if (get_point_from_weights(tentative_weights)
3839  .distance(trial_point) <
3840  get_point_from_weights(guess_weights).distance(trial_point))
3841  guess_weights = tentative_weights;
3842  else
3843  break;
3844  Assert(internal::weights_are_ok<structdim>(guess_weights),
3845  ExcInternalError());
3846  }
3847  Assert(internal::weights_are_ok<structdim>(guess_weights),
3848  ExcInternalError());
3849  projected_point = get_point_from_weights(guess_weights);
3850  }
3851 
3852  // if structdim == 2 and the optimal point is not on the interior then
3853  // we may be able to get a more accurate result by projecting onto the
3854  // lines.
3855  if (structdim == 2)
3856  {
3857  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
3858  line_projections;
3859  for (unsigned int line_n = 0;
3860  line_n < GeometryInfo<structdim>::lines_per_cell;
3861  ++line_n)
3862  {
3863  line_projections[line_n] =
3864  project_to_object(object->line(line_n), trial_point);
3865  }
3866  std::sort(line_projections.begin(),
3867  line_projections.end(),
3868  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
3869  return a.distance(trial_point) <
3870  b.distance(trial_point);
3871  });
3872  if (line_projections[0].distance(trial_point) <
3873  projected_point.distance(trial_point))
3874  projected_point = line_projections[0];
3875  }
3876  }
3877  else
3878  {
3879  Assert(false, ExcNotImplemented());
3880  return projected_point;
3881  }
3882 
3883  return projected_point;
3884  }
3885 
3886 
3887 
3888  template <int dim, typename T>
3889  template <class Archive>
3890  void
3892  const unsigned int /*version*/) const
3893  {
3894  Assert(cell_ids.size() == data.size(),
3895  ExcDimensionMismatch(cell_ids.size(), data.size()));
3896  // archive the cellids in an efficient binary format
3897  const std::size_t n_cells = cell_ids.size();
3898  ar & n_cells;
3899  for (const auto &id : cell_ids)
3900  {
3901  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3902  ar & binary_cell_id;
3903  }
3904 
3905  ar &data;
3906  }
3907 
3908 
3909 
3910  template <int dim, typename T>
3911  template <class Archive>
3912  void
3914  const unsigned int /*version*/)
3915  {
3916  std::size_t n_cells;
3917  ar & n_cells;
3918  cell_ids.clear();
3919  cell_ids.reserve(n_cells);
3920  for (unsigned int c = 0; c < n_cells; ++c)
3921  {
3922  CellId::binary_type value;
3923  ar & value;
3924  cell_ids.emplace_back(value);
3925  }
3926  ar &data;
3927  }
3928 
3929 
3930  namespace internal
3931  {
3932  template <typename DataType,
3933  typename MeshType,
3934  typename MeshCellIteratorType>
3935  void
3936  exchange_cell_data(
3937  const MeshType &mesh,
3938  const std::function<
3939  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
3940  const std::function<void(const MeshCellIteratorType &, const DataType &)>
3941  & unpack,
3942  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
3943  const std::function<void(
3944  const std::function<void(const MeshCellIteratorType &,
3945  const types::subdomain_id)> &)> &process_cells,
3946  const std::function<std::set<types::subdomain_id>(
3947  const parallel::TriangulationBase<MeshType::dimension,
3948  MeshType::space_dimension> &)>
3949  &compute_ghost_owners)
3950  {
3951 # ifndef DEAL_II_WITH_MPI
3952  (void)mesh;
3953  (void)pack;
3954  (void)unpack;
3955  (void)cell_filter;
3956  (void)process_cells;
3957  (void)compute_ghost_owners;
3958  Assert(false,
3959  ExcMessage("GridTools::exchange_cell_data() requires MPI."));
3960 # else
3961  constexpr int dim = MeshType::dimension;
3962  constexpr int spacedim = MeshType::space_dimension;
3963  auto tria =
3964  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3965  &mesh.get_triangulation());
3966  Assert(
3967  tria != nullptr,
3968  ExcMessage(
3969  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3970 
3971  // build list of cells to request for each neighbor
3972  std::set<::types::subdomain_id> ghost_owners =
3973  compute_ghost_owners(*tria);
3974  std::map<::types::subdomain_id,
3975  std::vector<typename CellId::binary_type>>
3976  neighbor_cell_list;
3977 
3978  for (const auto ghost_owner : ghost_owners)
3979  neighbor_cell_list[ghost_owner] = {};
3980 
3981  process_cells([&](const auto &cell, const auto key) {
3982  if (cell_filter(cell))
3983  neighbor_cell_list[key].emplace_back(
3984  cell->id().template to_binary<spacedim>());
3985  });
3986 
3987  Assert(ghost_owners.size() == neighbor_cell_list.size(),
3988  ExcInternalError());
3989 
3990 
3991  // Before sending & receiving, make sure we protect this section with
3992  // a mutex:
3993  static Utilities::MPI::CollectiveMutex mutex;
3995  mutex, tria->get_communicator());
3996 
3997  const int mpi_tag =
3999  const int mpi_tag_reply =
4001 
4002  // send our requests:
4003  std::vector<MPI_Request> requests(ghost_owners.size());
4004  {
4005  unsigned int idx = 0;
4006  for (const auto &it : neighbor_cell_list)
4007  {
4008  // send the data about the relevant cells
4009  const int ierr = MPI_Isend(it.second.data(),
4010  it.second.size() * sizeof(it.second[0]),
4011  MPI_BYTE,
4012  it.first,
4013  mpi_tag,
4014  tria->get_communicator(),
4015  &requests[idx]);
4016  AssertThrowMPI(ierr);
4017  ++idx;
4018  }
4019  }
4020 
4021  using DestinationToBufferMap =
4022  std::map<::types::subdomain_id,
4024  DestinationToBufferMap destination_to_data_buffer_map;
4025 
4026  // receive requests and reply with the ghost indices
4027  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4028  ghost_owners.size());
4029  std::vector<std::vector<types::global_dof_index>>
4030  send_dof_numbers_and_indices(ghost_owners.size());
4031  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4032  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4033 
4034  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4035  {
4036  MPI_Status status;
4037  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4038  mpi_tag,
4039  tria->get_communicator(),
4040  &status);
4041  AssertThrowMPI(ierr);
4042 
4043  int len;
4044  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4045  AssertThrowMPI(ierr);
4046  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4047  ExcInternalError());
4048 
4049  const unsigned int n_cells =
4050  len / sizeof(typename CellId::binary_type);
4051  cell_data_to_send[idx].resize(n_cells);
4052 
4053  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4054  len,
4055  MPI_BYTE,
4056  status.MPI_SOURCE,
4057  status.MPI_TAG,
4058  tria->get_communicator(),
4059  &status);
4060  AssertThrowMPI(ierr);
4061 
4062  // store data for each cell
4063  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4064  {
4065  const auto cell =
4066  CellId(cell_data_to_send[idx][c]).to_cell(*tria);
4067 
4068  MeshCellIteratorType mesh_it(tria,
4069  cell->level(),
4070  cell->index(),
4071  &mesh);
4072  const std_cxx17::optional<DataType> data = pack(mesh_it);
4073 
4074  if (data)
4075  {
4076  typename DestinationToBufferMap::iterator p =
4077  destination_to_data_buffer_map
4078  .insert(std::make_pair(
4079  idx,
4080  GridTools::CellDataTransferBuffer<dim, DataType>()))
4081  .first;
4082 
4083  p->second.cell_ids.emplace_back(cell->id());
4084  p->second.data.emplace_back(*data);
4085  }
4086  }
4087 
4088  // send reply
4089  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4090  destination_to_data_buffer_map[idx];
4091 
4092  sendbuffers[idx] =
4093  Utilities::pack(data, /*enable_compression*/ false);
4094  ierr = MPI_Isend(sendbuffers[idx].data(),
4095  sendbuffers[idx].size(),
4096  MPI_BYTE,
4097  status.MPI_SOURCE,
4098  mpi_tag_reply,
4099  tria->get_communicator(),
4100  &reply_requests[idx]);
4101  AssertThrowMPI(ierr);
4102  }
4103 
4104  // finally receive the replies
4105  std::vector<char> receive;
4106  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4107  {
4108  MPI_Status status;
4109  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4110  mpi_tag_reply,
4111  tria->get_communicator(),
4112  &status);
4113  AssertThrowMPI(ierr);
4114 
4115  int len;
4116  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4117  AssertThrowMPI(ierr);
4118 
4119  receive.resize(len);
4120 
4121  char *ptr = receive.data();
4122  ierr = MPI_Recv(ptr,
4123  len,
4124  MPI_BYTE,
4125  status.MPI_SOURCE,
4126  status.MPI_TAG,
4127  tria->get_communicator(),
4128  &status);
4129  AssertThrowMPI(ierr);
4130 
4131  auto cellinfo =
4132  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4133  receive, /*enable_compression*/ false);
4134 
4135  DataType *data = cellinfo.data.data();
4136  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4137  {
4139  tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
4140 
4141  MeshCellIteratorType cell(tria,
4142  tria_cell->level(),
4143  tria_cell->index(),
4144  &mesh);
4145 
4146  unpack(cell, *data);
4147  }
4148  }
4149 
4150  // make sure that all communication is finished
4151  // when we leave this function.
4152  if (requests.size() > 0)
4153  {
4154  const int ierr =
4155  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4156  AssertThrowMPI(ierr);
4157  }
4158  if (reply_requests.size() > 0)
4159  {
4160  const int ierr = MPI_Waitall(reply_requests.size(),
4161  reply_requests.data(),
4162  MPI_STATUSES_IGNORE);
4163  AssertThrowMPI(ierr);
4164  }
4165 
4166 
4167 # endif // DEAL_II_WITH_MPI
4168  }
4169 
4170  } // namespace internal
4171 
4172  template <typename DataType, typename MeshType>
4173  void
4175  const MeshType & mesh,
4176  const std::function<std_cxx17::optional<DataType>(
4177  const typename MeshType::active_cell_iterator &)> &pack,
4178  const std::function<void(const typename MeshType::active_cell_iterator &,
4179  const DataType &)> & unpack,
4180  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4181  &cell_filter)
4182  {
4183 # ifndef DEAL_II_WITH_MPI
4184  (void)mesh;
4185  (void)pack;
4186  (void)unpack;
4187  (void)cell_filter;
4188  Assert(false,
4189  ExcMessage(
4190  "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
4191 # else
4192  internal::exchange_cell_data<DataType,
4193  MeshType,
4194  typename MeshType::active_cell_iterator>(
4195  mesh,
4196  pack,
4197  unpack,
4198  cell_filter,
4199  [&](const auto &process) {
4200  for (const auto &cell : mesh.active_cell_iterators())
4201  if (cell->is_ghost())
4202  process(cell, cell->subdomain_id());
4203  },
4204  [](const auto &tria) { return tria.ghost_owners(); });
4205 # endif
4206  }
4207 
4208 
4209 
4210  template <typename DataType, typename MeshType>
4211  void
4213  const MeshType & mesh,
4214  const std::function<std_cxx17::optional<DataType>(
4215  const typename MeshType::level_cell_iterator &)> &pack,
4216  const std::function<void(const typename MeshType::level_cell_iterator &,
4217  const DataType &)> & unpack,
4218  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4219  &cell_filter)
4220  {
4221 # ifndef DEAL_II_WITH_MPI
4222  (void)mesh;
4223  (void)pack;
4224  (void)unpack;
4225  (void)cell_filter;
4226  Assert(false,
4227  ExcMessage(
4228  "GridTools::exchange_cell_data_to_level_ghosts() requires MPI."));
4229 # else
4230  internal::exchange_cell_data<DataType,
4231  MeshType,
4232  typename MeshType::level_cell_iterator>(
4233  mesh,
4234  pack,
4235  unpack,
4236  cell_filter,
4237  [&](const auto &process) {
4238  for (const auto &cell : mesh.cell_iterators())
4239  if (cell->level_subdomain_id() !=
4241  !cell->is_locally_owned_on_level())
4242  process(cell, cell->level_subdomain_id());
4243  },
4244  [](const auto &tria) { return tria.level_ghost_owners(); });
4245 # endif
4246  }
4247 } // namespace GridTools
4248 
4249 # endif
4250 
4252 
4253 /*---------------------------- grid_tools.h ---------------------------*/
4254 /* end of #ifndef dealii_grid_tools_h */
4255 #endif
4256 /*---------------------------- grid_tools.h ---------------------------*/
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:3816
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:3600
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:213
unsigned int n_vertices() const
const types::manifold_id flat_manifold_id
Definition: types.h:273
static const unsigned int invalid_unsigned_int
Definition: types.h:191
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:982
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3575
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:78
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true)
Definition: grid_tools.cc:1015
static ::ExceptionBase & ExcPointNotFoundInCoarseGrid(Point< N > arg1)
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:194
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: grid_tools.cc:5187
Contents is actually a matrix.
Triangulation< dim, spacedim >::cell_iterator to_cell(const Triangulation< dim, spacedim > &tria) const
Definition: cell_id.cc:159
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:3784
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< std::set< typename MeshType< dim, spacedim >::active_cell_iterator >> &vertex_to_cell_map, const std::vector< std::vector< Tensor< 1, spacedim >>> &vertex_to_cell_centers, const typename MeshType< dim, spacedim >::active_cell_iterator &cell_hint=typename MeshType< dim, spacedim >::active_cell_iterator(), const std::vector< bool > &marked_vertices={}, const RTree< std::pair< Point< spacedim >, unsigned int >> &used_vertices_rtree=RTree< std::pair< Point< spacedim >, unsigned int >>{}, const double tolerance=1.e-10)
Definition: grid_tools.cc:1575
void load(Archive &ar, const unsigned int version)
return_type guess_point_owner(const RTree< std::pair< BoundingBox< spacedim >, unsigned int >> &covering_rtree, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:2022
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:3878
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:839
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:135
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2072
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
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)
Definition: grid_tools.cc:4954
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12228
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:2123
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, MPI_Comm mpi_communicator)
Definition: grid_tools.cc:5347
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:79
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
std::pair< Point< MeshType::space_dimension >, Point< MeshType::space_dimension > > compute_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3062
std::bitset< 3 > orientation
Definition: grid_tools.h:2468
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2921
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3033
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:145
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
const double angle
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:376
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:515
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2524
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
void partition_triangulation(const unsigned int n_partitions, const std::vector< unsigned int > &cell_weights, const SparsityPattern &cell_connection_graph, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:2735
cell_iterator end() const
Definition: tria.cc:12294
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:417
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:92
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
std::map< types::global_dof_index, std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandler< dim, spacedim > &dof_handler)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
static const char T
T sum(const T &t, const MPI_Comm &mpi_communicator)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2553
bool orthogonal_equality(const FaceIterator &face1, const FaceIterator &face2, const int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
#define Assert(cond, exc)
Definition: exceptions.h:1403
Signals signals
Definition: tria.h:2216
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3037
Abstract base class for mapping classes.
Definition: mapping.h:301
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
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:664
#define DeclException0(Exception0)
Definition: exceptions.h:470
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:1827
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, MPI_Comm mpi_communicator)
Definition: grid_tools.cc:5251
unsigned int level
Definition: grid_out.cc:4355
void save(Archive &ar, const unsigned int version) const
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:3849
Point< 3 > vertices[4]
void build_triangulation_from_patch(const std::vector< typename Container::active_cell_iterator > &patch, Triangulation< Container::dimension, Container::space_dimension > &local_triangulation, std::map< typename Triangulation< Container::dimension, Container::space_dimension >::active_cell_iterator, typename Container::active_cell_iterator > &patch_to_global_tria_map)
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:5410
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:3668
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:3533
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1182
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: hp.h:117
std::vector< std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_active_cells_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const double tolerance=1e-10, const std::vector< bool > &marked_vertices={})
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12833
Definition: cell_id.h:69
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator >> &vertex_to_cells)
Definition: grid_tools.cc:1512
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=[](const typename MeshType::active_cell_iterator &) { return true;})
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
unsigned int global_dof_index
Definition: types.h:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:302
__global__ void set(Number *val, const Number s, const size_type N)
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:4229
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
Definition: table.h:1907
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1684
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:413
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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:1380
double cell_measure(const T &,...)
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:260
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:2948
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2476
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1326
Point< 3 > center
FullMatrix< double > matrix
Definition: grid_tools.h:2482
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:3696
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:828
static ::ExceptionBase & ExcNotImplemented()
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:4193
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2978
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:2816
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
numbers::NumberTraits< Number >::real_type norm() const
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: grid_tools.cc:1744
unsigned int find_closest_vertex(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={})
Definition: grid_tools.cc:1318
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
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=[](const typename MeshType::level_cell_iterator &) { return true;})
static ::ExceptionBase & ExcPointNotFound(Point< N > arg1)
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id, const int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator >> &matched_pairs, const ::Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:819
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_layer_within_distance(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const double layer_thickness)
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:2962
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
static ::ExceptionBase & ExcInternalError()
void delete_duplicated_vertices(std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
Definition: grid_tools.cc:616
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:5532