Reference documentation for deal.II version Git cd8b20072f 2020-09-24 16:28:23 -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 : cell->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 : cell->face_indices())
3219  if (cell->face(face)->has_children() &&
3220  !cell->face(face)->at_boundary())
3221  {
3222  Assert(cell->reference_cell_type() ==
3224  ExcNotImplemented());
3225 
3226  // this line has children
3227  cell->face(face)->child(0)->vertex(1) =
3228  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3229  2;
3230  }
3231  }
3232  else if (dim == 3)
3233  {
3235  cell = triangulation.begin_active(),
3236  endc = triangulation.end();
3237  for (; cell != endc; ++cell)
3238  for (const unsigned int face : cell->face_indices())
3239  if (cell->face(face)->has_children() &&
3240  !cell->face(face)->at_boundary())
3241  {
3242  Assert(cell->reference_cell_type() ==
3244  ExcNotImplemented());
3245 
3246  // this face has hanging nodes
3247  cell->face(face)->child(0)->vertex(1) =
3248  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3249  2.0;
3250  cell->face(face)->child(0)->vertex(2) =
3251  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3252  2.0;
3253  cell->face(face)->child(1)->vertex(3) =
3254  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3255  2.0;
3256  cell->face(face)->child(2)->vertex(3) =
3257  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3258  2.0;
3259 
3260  // center of the face
3261  cell->face(face)->child(0)->vertex(3) =
3262  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3263  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3264  4.0;
3265  }
3266  }
3267 
3268  // Make sure FEValues notices that the mesh has changed
3269  triangulation.signals.mesh_movement();
3270  }
3271 
3272 
3273 
3274  template <class MeshType>
3275  std::vector<typename MeshType::active_cell_iterator>
3276  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3277  {
3278  std::vector<typename MeshType::active_cell_iterator> child_cells;
3279 
3280  if (cell->has_children())
3281  {
3282  for (unsigned int child = 0; child < cell->n_children(); ++child)
3283  if (cell->child(child)->has_children())
3284  {
3285  const std::vector<typename MeshType::active_cell_iterator>
3286  children = get_active_child_cells<MeshType>(cell->child(child));
3287  child_cells.insert(child_cells.end(),
3288  children.begin(),
3289  children.end());
3290  }
3291  else
3292  child_cells.push_back(cell->child(child));
3293  }
3294 
3295  return child_cells;
3296  }
3297 
3298 
3299 
3300  template <class MeshType>
3301  void
3303  const typename MeshType::active_cell_iterator & cell,
3304  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3305  {
3306  active_neighbors.clear();
3307  for (const unsigned int n : cell->face_indices())
3308  if (!cell->at_boundary(n))
3309  {
3310  if (MeshType::dimension == 1)
3311  {
3312  // check children of neighbor. note
3313  // that in 1d children of the neighbor
3314  // may be further refined. In 1d the
3315  // case is simple since we know what
3316  // children bound to the present cell
3317  typename MeshType::cell_iterator neighbor_child =
3318  cell->neighbor(n);
3319  if (!neighbor_child->is_active())
3320  {
3321  while (neighbor_child->has_children())
3322  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3323 
3324  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3325  ExcInternalError());
3326  }
3327  active_neighbors.push_back(neighbor_child);
3328  }
3329  else
3330  {
3331  if (cell->face(n)->has_children())
3332  // this neighbor has children. find
3333  // out which border to the present
3334  // cell
3335  for (unsigned int c = 0;
3336  c < cell->face(n)->number_of_children();
3337  ++c)
3338  active_neighbors.push_back(
3339  cell->neighbor_child_on_subface(n, c));
3340  else
3341  {
3342  // the neighbor must be active
3343  // himself
3344  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3345  active_neighbors.push_back(cell->neighbor(n));
3346  }
3347  }
3348  }
3349  }
3350 
3351 
3352 
3353  namespace internal
3354  {
3355  namespace ProjectToObject
3356  {
3369  struct CrossDerivative
3370  {
3371  const unsigned int direction_0;
3372  const unsigned int direction_1;
3373 
3374  CrossDerivative(const unsigned int d0, const unsigned int d1);
3375  };
3376 
3377  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3378  const unsigned int d1)
3379  : direction_0(d0)
3380  , direction_1(d1)
3381  {}
3382 
3383 
3384 
3389  template <typename F>
3390  inline auto
3391  centered_first_difference(const double center,
3392  const double step,
3393  const F &f) -> decltype(f(center) - f(center))
3394  {
3395  return (f(center + step) - f(center - step)) / (2.0 * step);
3396  }
3397 
3398 
3399 
3404  template <typename F>
3405  inline auto
3406  centered_second_difference(const double center,
3407  const double step,
3408  const F &f) -> decltype(f(center) - f(center))
3409  {
3410  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3411  (step * step);
3412  }
3413 
3414 
3415 
3425  template <int structdim, typename F>
3426  inline auto
3427  cross_stencil(
3428  const CrossDerivative cross_derivative,
3430  const double step,
3431  const F &f) -> decltype(f(center) - f(center))
3432  {
3434  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3435  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3436  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3437  1.0 / 3.0 * f(center - simplex_vector) +
3438  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3439  step;
3440  }
3441 
3442 
3443 
3450  template <int spacedim, int structdim, typename F>
3451  inline double
3452  gradient_entry(
3453  const unsigned int row_n,
3454  const unsigned int dependent_direction,
3455  const Point<spacedim> &p0,
3457  const double step,
3458  const F & f)
3459  {
3461  dependent_direction <
3463  ExcMessage("This function assumes that the last weight is a "
3464  "dependent variable (and hence we cannot take its "
3465  "derivative directly)."));
3466  Assert(row_n != dependent_direction,
3467  ExcMessage(
3468  "We cannot differentiate with respect to the variable "
3469  "that is assumed to be dependent."));
3470 
3471  const Point<spacedim> manifold_point = f(center);
3472  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3473  {row_n, dependent_direction}, center, step, f);
3474  double entry = 0.0;
3475  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3476  entry +=
3477  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3478  return entry;
3479  }
3480 
3486  template <typename Iterator, int spacedim, int structdim>
3488  project_to_d_linear_object(const Iterator & object,
3489  const Point<spacedim> &trial_point)
3490  {
3491  // let's look at this for simplicity for a quad (structdim==2) in a
3492  // space with spacedim>2 (notate trial_point by y): all points on the
3493  // surface are given by
3494  // x(\xi) = sum_i v_i phi_x(\xi)
3495  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3496  // reference coordinates of the quad. so what we are trying to do is
3497  // find a point x on the surface that is closest to the point y. there
3498  // are different ways to solve this problem, but in the end it's a
3499  // nonlinear problem and we have to find reference coordinates \xi so
3500  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3501  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3502  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3503  // have to use a Newton method to find the answer. This leads to the
3504  // following formulation of Newton steps:
3505  //
3506  // Given \xi_k, find \delta\xi_k so that
3507  // H_k \delta\xi_k = - F_k
3508  // where H_k is an approximation to the second derivatives of J at
3509  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3510  // number of times until the right hand side is small enough. As a
3511  // stopping criterion, we terminate if ||\delta\xi||<eps.
3512  //
3513  // As for the Hessian, the best choice would be
3514  // H_k = J''(\xi_k)
3515  // but we'll opt for the simpler Gauss-Newton form
3516  // H_k = A^T A
3517  // i.e.
3518  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3519  // \partial_n phi_i *
3520  // \partial_m phi_j
3521  // we start at xi=(0.5, 0.5).
3522  Point<structdim> xi;
3523  for (unsigned int d = 0; d < structdim; ++d)
3524  xi[d] = 0.5;
3525 
3526  Point<spacedim> x_k;
3527  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3528  x_k += object->vertex(i) *
3530 
3531  do
3532  {
3534  for (const unsigned int i :
3536  F_k +=
3537  (x_k - trial_point) * object->vertex(i) *
3539  i);
3540 
3542  for (const unsigned int i :
3544  for (const unsigned int j :
3546  {
3549  xi, i),
3551  xi, j));
3552  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3553  }
3554 
3555  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3556  xi += delta_xi;
3557 
3558  x_k = Point<spacedim>();
3559  for (const unsigned int i :
3561  x_k += object->vertex(i) *
3563 
3564  if (delta_xi.norm() < 1e-7)
3565  break;
3566  }
3567  while (true);
3568 
3569  return x_k;
3570  }
3571  } // namespace ProjectToObject
3572  } // namespace internal
3573 
3574 
3575 
3576  namespace internal
3577  {
3578  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3579  // inside the project_to_object function below.
3580  template <int structdim>
3581  inline bool
3582  weights_are_ok(
3584  {
3585  // clang has trouble figuring out structdim here, so define it
3586  // again:
3587  static const std::size_t n_vertices_per_cell =
3589  n_independent_components;
3590  std::array<double, n_vertices_per_cell> copied_weights;
3591  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3592  {
3593  copied_weights[i] = v[i];
3594  if (v[i] < 0.0 || v[i] > 1.0)
3595  return false;
3596  }
3597 
3598  // check the sum: try to avoid some roundoff errors by summing in order
3599  std::sort(copied_weights.begin(), copied_weights.end());
3600  const double sum =
3601  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3602  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3603  }
3604  } // namespace internal
3605 
3606  template <typename Iterator>
3609  const Iterator & object,
3611  {
3612  const int spacedim = Iterator::AccessorType::space_dimension;
3613  const int structdim = Iterator::AccessorType::structure_dimension;
3614 
3615  Point<spacedim> projected_point = trial_point;
3616 
3617  if (structdim >= spacedim)
3618  return projected_point;
3619  else if (structdim == 1 || structdim == 2)
3620  {
3621  using namespace internal::ProjectToObject;
3622  // Try to use the special flat algorithm for quads (this is better
3623  // than the general algorithm in 3D). This does not take into account
3624  // whether projected_point is outside the quad, but we optimize along
3625  // lines below anyway:
3626  const int dim = Iterator::AccessorType::dimension;
3627  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3628  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3629  &manifold) != nullptr)
3630  {
3631  projected_point =
3632  project_to_d_linear_object<Iterator, spacedim, structdim>(
3633  object, trial_point);
3634  }
3635  else
3636  {
3637  // We want to find a point on the convex hull (defined by the
3638  // vertices of the object and the manifold description) that is
3639  // relatively close to the trial point. This has a few issues:
3640  //
3641  // 1. For a general convex hull we are not guaranteed that a unique
3642  // minimum exists.
3643  // 2. The independent variables in the optimization process are the
3644  // weights given to Manifold::get_new_point, which must sum to 1,
3645  // so we cannot use standard finite differences to approximate a
3646  // gradient.
3647  //
3648  // There is not much we can do about 1., but for 2. we can derive
3649  // finite difference stencils that work on a structdim-dimensional
3650  // simplex and rewrite the optimization problem to use those
3651  // instead. Consider the structdim 2 case and let
3652  //
3653  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3654  // c2, c3})
3655  //
3656  // where {c0, c1, c2, c3} are the weights for the four vertices on
3657  // the quadrilateral. We seek to minimize the Euclidean distance
3658  // between F(...) and trial_point. We can solve for c3 in terms of
3659  // the other weights and get, for one coordinate direction
3660  //
3661  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3662  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3663  //
3664  // where we substitute back in for c3 after taking the
3665  // derivative. We can compute a stencil for the cross derivative
3666  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3667  // (and gradient_entry computes the sum over the independent
3668  // variables). Below, we somewhat arbitrarily pick the last
3669  // component as the dependent one.
3670  //
3671  // Since we can now calculate derivatives of the objective
3672  // function we can use gradient descent to minimize it.
3673  //
3674  // Of course, this is much simpler in the structdim = 1 case (we
3675  // could rewrite the projection as a 1D optimization problem), but
3676  // to reduce the potential for bugs we use the same code in both
3677  // cases.
3678  const double step_size = object->diameter() / 64.0;
3679 
3680  constexpr unsigned int n_vertices_per_cell =
3682 
3683  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3684  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3685  ++vertex_n)
3686  vertices[vertex_n] = object->vertex(vertex_n);
3687 
3688  auto get_point_from_weights =
3689  [&](const Tensor<1, n_vertices_per_cell> &weights)
3690  -> Point<spacedim> {
3691  return object->get_manifold().get_new_point(
3692  make_array_view(vertices.begin(), vertices.end()),
3693  make_array_view(weights.begin_raw(), weights.end_raw()));
3694  };
3695 
3696  // pick the initial weights as (normalized) inverse distances from
3697  // the trial point:
3698  Tensor<1, n_vertices_per_cell> guess_weights;
3699  double guess_weights_sum = 0.0;
3700  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3701  ++vertex_n)
3702  {
3703  const double distance =
3704  vertices[vertex_n].distance(trial_point);
3705  if (distance == 0.0)
3706  {
3707  guess_weights = 0.0;
3708  guess_weights[vertex_n] = 1.0;
3709  guess_weights_sum = 1.0;
3710  break;
3711  }
3712  else
3713  {
3714  guess_weights[vertex_n] = 1.0 / distance;
3715  guess_weights_sum += guess_weights[vertex_n];
3716  }
3717  }
3718  guess_weights /= guess_weights_sum;
3719  Assert(internal::weights_are_ok<structdim>(guess_weights),
3720  ExcInternalError());
3721 
3722  // The optimization algorithm consists of two parts:
3723  //
3724  // 1. An outer loop where we apply the gradient descent algorithm.
3725  // 2. An inner loop where we do a line search to find the optimal
3726  // length of the step one should take in the gradient direction.
3727  //
3728  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3729  {
3730  const unsigned int dependent_direction =
3731  n_vertices_per_cell - 1;
3732  Tensor<1, n_vertices_per_cell> current_gradient;
3733  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3734  ++row_n)
3735  {
3736  if (row_n != dependent_direction)
3737  {
3738  current_gradient[row_n] =
3739  gradient_entry<spacedim, structdim>(
3740  row_n,
3741  dependent_direction,
3742  trial_point,
3743  guess_weights,
3744  step_size,
3745  get_point_from_weights);
3746 
3747  current_gradient[dependent_direction] -=
3748  current_gradient[row_n];
3749  }
3750  }
3751 
3752  // We need to travel in the -gradient direction, as noted
3753  // above, but we may not want to take a full step in that
3754  // direction; instead, guess that we will go -0.5*gradient and
3755  // do quasi-Newton iteration to pick the best multiplier. The
3756  // goal is to find a scalar alpha such that
3757  //
3758  // F(x - alpha g)
3759  //
3760  // is minimized, where g is the gradient and F is the
3761  // objective function. To find the optimal value we find roots
3762  // of the derivative of the objective function with respect to
3763  // alpha by Newton iteration, where we approximate the first
3764  // and second derivatives of F(x - alpha g) with centered
3765  // finite differences.
3766  double gradient_weight = -0.5;
3767  auto gradient_weight_objective_function =
3768  [&](const double gradient_weight_guess) -> double {
3769  return (trial_point -
3770  get_point_from_weights(guess_weights +
3771  gradient_weight_guess *
3772  current_gradient))
3773  .norm_square();
3774  };
3775 
3776  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3777  {
3778  const double update_numerator = centered_first_difference(
3779  gradient_weight,
3780  step_size,
3781  gradient_weight_objective_function);
3782  const double update_denominator =
3783  centered_second_difference(
3784  gradient_weight,
3785  step_size,
3786  gradient_weight_objective_function);
3787 
3788  // avoid division by zero. Note that we limit the gradient
3789  // weight below
3790  if (std::abs(update_denominator) == 0.0)
3791  break;
3792  gradient_weight =
3793  gradient_weight - update_numerator / update_denominator;
3794 
3795  // Put a fairly lenient bound on the largest possible
3796  // gradient (things tend to be locally flat, so the gradient
3797  // itself is usually small)
3798  if (std::abs(gradient_weight) > 10)
3799  {
3800  gradient_weight = -10.0;
3801  break;
3802  }
3803  }
3804 
3805  // It only makes sense to take convex combinations with weights
3806  // between zero and one. If the update takes us outside of this
3807  // region then rescale the update to stay within the region and
3808  // try again
3809  Tensor<1, n_vertices_per_cell> tentative_weights =
3810  guess_weights + gradient_weight * current_gradient;
3811 
3812  double new_gradient_weight = gradient_weight;
3813  for (unsigned int iteration_count = 0; iteration_count < 40;
3814  ++iteration_count)
3815  {
3816  if (internal::weights_are_ok<structdim>(tentative_weights))
3817  break;
3818 
3819  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3820  {
3821  if (tentative_weights[i] < 0.0)
3822  {
3823  tentative_weights -=
3824  (tentative_weights[i] / current_gradient[i]) *
3825  current_gradient;
3826  }
3827  if (tentative_weights[i] < 0.0 ||
3828  1.0 < tentative_weights[i])
3829  {
3830  new_gradient_weight /= 2.0;
3831  tentative_weights =
3832  guess_weights +
3833  new_gradient_weight * current_gradient;
3834  }
3835  }
3836  }
3837 
3838  // the update might still send us outside the valid region, so
3839  // check again and quit if the update is still not valid
3840  if (!internal::weights_are_ok<structdim>(tentative_weights))
3841  break;
3842 
3843  // if we cannot get closer by traveling in the gradient
3844  // direction then quit
3845  if (get_point_from_weights(tentative_weights)
3846  .distance(trial_point) <
3847  get_point_from_weights(guess_weights).distance(trial_point))
3848  guess_weights = tentative_weights;
3849  else
3850  break;
3851  Assert(internal::weights_are_ok<structdim>(guess_weights),
3852  ExcInternalError());
3853  }
3854  Assert(internal::weights_are_ok<structdim>(guess_weights),
3855  ExcInternalError());
3856  projected_point = get_point_from_weights(guess_weights);
3857  }
3858 
3859  // if structdim == 2 and the optimal point is not on the interior then
3860  // we may be able to get a more accurate result by projecting onto the
3861  // lines.
3862  if (structdim == 2)
3863  {
3864  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
3865  line_projections;
3866  for (unsigned int line_n = 0;
3867  line_n < GeometryInfo<structdim>::lines_per_cell;
3868  ++line_n)
3869  {
3870  line_projections[line_n] =
3871  project_to_object(object->line(line_n), trial_point);
3872  }
3873  std::sort(line_projections.begin(),
3874  line_projections.end(),
3875  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
3876  return a.distance(trial_point) <
3877  b.distance(trial_point);
3878  });
3879  if (line_projections[0].distance(trial_point) <
3880  projected_point.distance(trial_point))
3881  projected_point = line_projections[0];
3882  }
3883  }
3884  else
3885  {
3886  Assert(false, ExcNotImplemented());
3887  return projected_point;
3888  }
3889 
3890  return projected_point;
3891  }
3892 
3893 
3894 
3895  template <int dim, typename T>
3896  template <class Archive>
3897  void
3899  const unsigned int /*version*/) const
3900  {
3901  Assert(cell_ids.size() == data.size(),
3902  ExcDimensionMismatch(cell_ids.size(), data.size()));
3903  // archive the cellids in an efficient binary format
3904  const std::size_t n_cells = cell_ids.size();
3905  ar & n_cells;
3906  for (const auto &id : cell_ids)
3907  {
3908  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3909  ar & binary_cell_id;
3910  }
3911 
3912  ar &data;
3913  }
3914 
3915 
3916 
3917  template <int dim, typename T>
3918  template <class Archive>
3919  void
3921  const unsigned int /*version*/)
3922  {
3923  std::size_t n_cells;
3924  ar & n_cells;
3925  cell_ids.clear();
3926  cell_ids.reserve(n_cells);
3927  for (unsigned int c = 0; c < n_cells; ++c)
3928  {
3929  CellId::binary_type value;
3930  ar & value;
3931  cell_ids.emplace_back(value);
3932  }
3933  ar &data;
3934  }
3935 
3936 
3937  namespace internal
3938  {
3939  template <typename DataType,
3940  typename MeshType,
3941  typename MeshCellIteratorType>
3942  void
3943  exchange_cell_data(
3944  const MeshType &mesh,
3945  const std::function<
3946  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
3947  const std::function<void(const MeshCellIteratorType &, const DataType &)>
3948  & unpack,
3949  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
3950  const std::function<void(
3951  const std::function<void(const MeshCellIteratorType &,
3952  const types::subdomain_id)> &)> &process_cells,
3953  const std::function<std::set<types::subdomain_id>(
3954  const parallel::TriangulationBase<MeshType::dimension,
3955  MeshType::space_dimension> &)>
3956  &compute_ghost_owners)
3957  {
3958 # ifndef DEAL_II_WITH_MPI
3959  (void)mesh;
3960  (void)pack;
3961  (void)unpack;
3962  (void)cell_filter;
3963  (void)process_cells;
3964  (void)compute_ghost_owners;
3965  Assert(false,
3966  ExcMessage("GridTools::exchange_cell_data() requires MPI."));
3967 # else
3968  constexpr int dim = MeshType::dimension;
3969  constexpr int spacedim = MeshType::space_dimension;
3970  auto tria =
3971  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3972  &mesh.get_triangulation());
3973  Assert(
3974  tria != nullptr,
3975  ExcMessage(
3976  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3977 
3978  // build list of cells to request for each neighbor
3979  std::set<::types::subdomain_id> ghost_owners =
3980  compute_ghost_owners(*tria);
3981  std::map<::types::subdomain_id,
3982  std::vector<typename CellId::binary_type>>
3983  neighbor_cell_list;
3984 
3985  for (const auto ghost_owner : ghost_owners)
3986  neighbor_cell_list[ghost_owner] = {};
3987 
3988  process_cells([&](const auto &cell, const auto key) {
3989  if (cell_filter(cell))
3990  neighbor_cell_list[key].emplace_back(
3991  cell->id().template to_binary<spacedim>());
3992  });
3993 
3994  Assert(ghost_owners.size() == neighbor_cell_list.size(),
3995  ExcInternalError());
3996 
3997 
3998  // Before sending & receiving, make sure we protect this section with
3999  // a mutex:
4000  static Utilities::MPI::CollectiveMutex mutex;
4002  mutex, tria->get_communicator());
4003 
4004  const int mpi_tag =
4006  const int mpi_tag_reply =
4008 
4009  // send our requests:
4010  std::vector<MPI_Request> requests(ghost_owners.size());
4011  {
4012  unsigned int idx = 0;
4013  for (const auto &it : neighbor_cell_list)
4014  {
4015  // send the data about the relevant cells
4016  const int ierr = MPI_Isend(it.second.data(),
4017  it.second.size() * sizeof(it.second[0]),
4018  MPI_BYTE,
4019  it.first,
4020  mpi_tag,
4021  tria->get_communicator(),
4022  &requests[idx]);
4023  AssertThrowMPI(ierr);
4024  ++idx;
4025  }
4026  }
4027 
4028  using DestinationToBufferMap =
4029  std::map<::types::subdomain_id,
4031  DestinationToBufferMap destination_to_data_buffer_map;
4032 
4033  // receive requests and reply with the ghost indices
4034  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4035  ghost_owners.size());
4036  std::vector<std::vector<types::global_dof_index>>
4037  send_dof_numbers_and_indices(ghost_owners.size());
4038  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4039  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4040 
4041  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4042  {
4043  MPI_Status status;
4044  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4045  mpi_tag,
4046  tria->get_communicator(),
4047  &status);
4048  AssertThrowMPI(ierr);
4049 
4050  int len;
4051  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4052  AssertThrowMPI(ierr);
4053  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4054  ExcInternalError());
4055 
4056  const unsigned int n_cells =
4057  len / sizeof(typename CellId::binary_type);
4058  cell_data_to_send[idx].resize(n_cells);
4059 
4060  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4061  len,
4062  MPI_BYTE,
4063  status.MPI_SOURCE,
4064  status.MPI_TAG,
4065  tria->get_communicator(),
4066  &status);
4067  AssertThrowMPI(ierr);
4068 
4069  // store data for each cell
4070  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4071  {
4072  const auto cell =
4073  CellId(cell_data_to_send[idx][c]).to_cell(*tria);
4074 
4075  MeshCellIteratorType mesh_it(tria,
4076  cell->level(),
4077  cell->index(),
4078  &mesh);
4079  const std_cxx17::optional<DataType> data = pack(mesh_it);
4080 
4081  if (data)
4082  {
4083  typename DestinationToBufferMap::iterator p =
4084  destination_to_data_buffer_map
4085  .insert(std::make_pair(
4086  idx,
4087  GridTools::CellDataTransferBuffer<dim, DataType>()))
4088  .first;
4089 
4090  p->second.cell_ids.emplace_back(cell->id());
4091  p->second.data.emplace_back(*data);
4092  }
4093  }
4094 
4095  // send reply
4096  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4097  destination_to_data_buffer_map[idx];
4098 
4099  sendbuffers[idx] =
4100  Utilities::pack(data, /*enable_compression*/ false);
4101  ierr = MPI_Isend(sendbuffers[idx].data(),
4102  sendbuffers[idx].size(),
4103  MPI_BYTE,
4104  status.MPI_SOURCE,
4105  mpi_tag_reply,
4106  tria->get_communicator(),
4107  &reply_requests[idx]);
4108  AssertThrowMPI(ierr);
4109  }
4110 
4111  // finally receive the replies
4112  std::vector<char> receive;
4113  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4114  {
4115  MPI_Status status;
4116  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4117  mpi_tag_reply,
4118  tria->get_communicator(),
4119  &status);
4120  AssertThrowMPI(ierr);
4121 
4122  int len;
4123  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4124  AssertThrowMPI(ierr);
4125 
4126  receive.resize(len);
4127 
4128  char *ptr = receive.data();
4129  ierr = MPI_Recv(ptr,
4130  len,
4131  MPI_BYTE,
4132  status.MPI_SOURCE,
4133  status.MPI_TAG,
4134  tria->get_communicator(),
4135  &status);
4136  AssertThrowMPI(ierr);
4137 
4138  auto cellinfo =
4139  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4140  receive, /*enable_compression*/ false);
4141 
4142  DataType *data = cellinfo.data.data();
4143  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4144  {
4146  tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
4147 
4148  MeshCellIteratorType cell(tria,
4149  tria_cell->level(),
4150  tria_cell->index(),
4151  &mesh);
4152 
4153  unpack(cell, *data);
4154  }
4155  }
4156 
4157  // make sure that all communication is finished
4158  // when we leave this function.
4159  if (requests.size() > 0)
4160  {
4161  const int ierr =
4162  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4163  AssertThrowMPI(ierr);
4164  }
4165  if (reply_requests.size() > 0)
4166  {
4167  const int ierr = MPI_Waitall(reply_requests.size(),
4168  reply_requests.data(),
4169  MPI_STATUSES_IGNORE);
4170  AssertThrowMPI(ierr);
4171  }
4172 
4173 
4174 # endif // DEAL_II_WITH_MPI
4175  }
4176 
4177  } // namespace internal
4178 
4179  template <typename DataType, typename MeshType>
4180  void
4182  const MeshType & mesh,
4183  const std::function<std_cxx17::optional<DataType>(
4184  const typename MeshType::active_cell_iterator &)> &pack,
4185  const std::function<void(const typename MeshType::active_cell_iterator &,
4186  const DataType &)> & unpack,
4187  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4188  &cell_filter)
4189  {
4190 # ifndef DEAL_II_WITH_MPI
4191  (void)mesh;
4192  (void)pack;
4193  (void)unpack;
4194  (void)cell_filter;
4195  Assert(false,
4196  ExcMessage(
4197  "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
4198 # else
4199  internal::exchange_cell_data<DataType,
4200  MeshType,
4201  typename MeshType::active_cell_iterator>(
4202  mesh,
4203  pack,
4204  unpack,
4205  cell_filter,
4206  [&](const auto &process) {
4207  for (const auto &cell : mesh.active_cell_iterators())
4208  if (cell->is_ghost())
4209  process(cell, cell->subdomain_id());
4210  },
4211  [](const auto &tria) { return tria.ghost_owners(); });
4212 # endif
4213  }
4214 
4215 
4216 
4217  template <typename DataType, typename MeshType>
4218  void
4220  const MeshType & mesh,
4221  const std::function<std_cxx17::optional<DataType>(
4222  const typename MeshType::level_cell_iterator &)> &pack,
4223  const std::function<void(const typename MeshType::level_cell_iterator &,
4224  const DataType &)> & unpack,
4225  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4226  &cell_filter)
4227  {
4228 # ifndef DEAL_II_WITH_MPI
4229  (void)mesh;
4230  (void)pack;
4231  (void)unpack;
4232  (void)cell_filter;
4233  Assert(false,
4234  ExcMessage(
4235  "GridTools::exchange_cell_data_to_level_ghosts() requires MPI."));
4236 # else
4237  internal::exchange_cell_data<DataType,
4238  MeshType,
4239  typename MeshType::level_cell_iterator>(
4240  mesh,
4241  pack,
4242  unpack,
4243  cell_filter,
4244  [&](const auto &process) {
4245  for (const auto &cell : mesh.cell_iterators())
4246  if (cell->level_subdomain_id() !=
4248  !cell->is_locally_owned_on_level())
4249  process(cell, cell->level_subdomain_id());
4250  },
4251  [](const auto &tria) { return tria.level_ghost_owners(); });
4252 # endif
4253  }
4254 } // namespace GridTools
4255 
4256 # endif
4257 
4259 
4260 /*---------------------------- grid_tools.h ---------------------------*/
4261 /* end of #ifndef dealii_grid_tools_h */
4262 #endif
4263 /*---------------------------- 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:3792
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:3577
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:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:978
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3552
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:1010
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:193
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:5194
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:3760
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:1567
Type get_hypercube(const unsigned int dim)
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:2019
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:3854
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:835
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:134
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2069
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:4961
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11138
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:2118
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:5354
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:3043
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:2902
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3014
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:507
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2505
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:2716
cell_iterator end() const
Definition: tria.cc:11204
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:415
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:2534
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:1411
Signals signals
Definition: tria.h:2238
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:665
#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:1824
#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:5258
unsigned int level
Definition: grid_out.cc:4338
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:3825
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:5417
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:3644
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:3510
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:11743
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:1504
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:293
__global__ void set(Number *val, const Number s, const size_type N)
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:4205
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1692
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:1372
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:259
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:2929
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2471
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:3672
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:824
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:4169
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2959
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:2797
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:1742
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:1310
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:815
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:2943
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:612
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:5540