Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_grid_tools_h
17 #define dealii_grid_tools_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 #include <deal.II/base/point.h>
25 
27 
29 
31 
32 #include <deal.II/fe/fe_values.h>
33 #include <deal.II/fe/mapping.h>
34 
37 #include <deal.II/grid/manifold.h>
38 #include <deal.II/grid/tria.h>
41 
46 
47 #include <deal.II/numerics/rtree.h>
48 
49 #include <boost/archive/binary_iarchive.hpp>
50 #include <boost/archive/binary_oarchive.hpp>
51 #include <boost/random/mersenne_twister.hpp>
52 #include <boost/serialization/array.hpp>
53 #include <boost/serialization/vector.hpp>
54 
55 #ifdef DEAL_II_WITH_ZLIB
56 # include <boost/iostreams/device/back_inserter.hpp>
57 # include <boost/iostreams/filter/gzip.hpp>
58 # include <boost/iostreams/filtering_stream.hpp>
59 # include <boost/iostreams/stream.hpp>
60 #endif
61 
62 #include <bitset>
63 #include <optional>
64 #include <set>
65 
66 #ifdef DEAL_II_HAVE_CXX20
67 # include <concepts>
68 #endif
69 
70 
72 
73 // Forward declarations
74 #ifndef DOXYGEN
75 namespace parallel
76 {
77  namespace distributed
78  {
79  template <int dim, int spacedim>
80  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
81  class Triangulation;
82  }
83 } // namespace parallel
84 
85 namespace hp
86 {
87  template <int, int>
88  class MappingCollection;
89 }
90 
91 class SparsityPattern;
92 
93 namespace GridTools
94 {
95  template <int dim, int spacedim>
96  class Cache;
97 }
98 #endif
99 
100 namespace internal
101 {
102  template <int dim, int spacedim, typename MeshType>
103  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
105  {
106  public:
107 #ifndef _MSC_VER
108  using type = typename MeshType::active_cell_iterator;
109 #else
111 #endif
112  };
113 
114 #ifdef _MSC_VER
115  template <int dim, int spacedim>
116  class ActiveCellIterator<dim, spacedim, DoFHandler<dim, spacedim>>
117  {
118  public:
119  using type =
121  };
122 #endif
123 } // namespace internal
124 
133 namespace GridTools
134 {
221  template <int dim, typename Transformation, int spacedim>
223  (std::invocable<Transformation, Point<spacedim>> &&
224  std::assignable_from<
225  Point<spacedim> &,
226  std::invoke_result_t<Transformation, Point<spacedim>>>))
227  void transform(const Transformation &transformation,
228  Triangulation<dim, spacedim> &triangulation);
229 
236  template <int dim, int spacedim>
237  void
238  shift(const Tensor<1, spacedim> &shift_vector,
239  Triangulation<dim, spacedim> &triangulation);
240 
241 
252  template <int dim, int spacedim>
253  void
254  rotate(const double angle, Triangulation<dim, spacedim> &triangulation);
255 
268  template <int dim>
269  void
270  rotate(const Tensor<1, 3, double> &axis,
271  const double angle,
272  Triangulation<dim, 3> &triangulation);
273 
289  template <int dim>
290  DEAL_II_DEPRECATED void
291  rotate(const double angle,
292  const unsigned int axis,
293  Triangulation<dim, 3> &triangulation);
294 
352  template <int dim>
353  void
354  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
355  Triangulation<dim> &tria,
356  const Function<dim, double> *coefficient = nullptr,
357  const bool solve_for_absolute_positions = false);
358 
366  template <int dim, int spacedim>
367  void
368  scale(const double scaling_factor,
369  Triangulation<dim, spacedim> &triangulation);
370 
390  template <int dim, int spacedim>
391  void
393  const double factor,
394  Triangulation<dim, spacedim> &triangulation,
395  const bool keep_boundary = true,
396  const unsigned int seed = boost::random::mt19937::default_seed);
397 
487  template <int dim, int spacedim>
488  void
490  const double limit_angle_fraction = .75);
491 
555  template <int dim, int spacedim>
556 #ifndef DOXYGEN
557  std::tuple<
558  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
559  std::vector<std::vector<Point<dim>>>,
560  std::vector<std::vector<unsigned int>>>
561 #else
562  return_type
563 #endif
565  const Cache<dim, spacedim> &cache,
566  const std::vector<Point<spacedim>> &points,
568  &cell_hint =
570 
604  template <int dim, int spacedim>
605 #ifndef DOXYGEN
606  std::tuple<
607  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
608  std::vector<std::vector<Point<dim>>>,
609  std::vector<std::vector<unsigned int>>,
610  std::vector<unsigned int>>
611 #else
612  return_type
613 #endif
615  const Cache<dim, spacedim> &cache,
616  const std::vector<Point<spacedim>> &points,
618  &cell_hint =
620 
700  template <int dim, int spacedim>
701 #ifndef DOXYGEN
702  std::tuple<
703  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
704  std::vector<std::vector<Point<dim>>>,
705  std::vector<std::vector<unsigned int>>,
706  std::vector<std::vector<Point<spacedim>>>,
707  std::vector<std::vector<unsigned int>>>
708 #else
709  return_type
710 #endif
712  const GridTools::Cache<dim, spacedim> &cache,
713  const std::vector<Point<spacedim>> &local_points,
714  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
715  const double tolerance = 1e-10,
716  const std::vector<bool> &marked_vertices = {},
717  const bool enforce_unique_mapping = true);
718 
719  namespace internal
720  {
735  template <int dim, int spacedim>
737  {
739 
746  void
747  finalize_setup();
748 
752  unsigned int n_searched_points;
753 
760  std::vector<std::tuple<std::pair<int, int>,
761  unsigned int,
762  unsigned int,
763  Point<dim>,
765  unsigned int>>
767 
771  std::vector<unsigned int> send_ranks;
772 
778  std::vector<unsigned int> send_ptrs;
779 
790  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
792 
796  std::vector<unsigned int> recv_ranks;
797 
803  std::vector<unsigned int> recv_ptrs;
804  };
805 
815  template <int dim, int spacedim>
818  const GridTools::Cache<dim, spacedim> &cache,
819  const std::vector<Point<spacedim>> &points,
820  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
821  const std::vector<bool> &marked_vertices,
822  const double tolerance,
823  const bool perform_handshake,
824  const bool enforce_unique_mapping = false);
825 
826 
834  template <int structdim, int spacedim>
836  {
841  std::array<::Point<spacedim>, structdim + 1>;
842 
851  std::vector<std::tuple<std::pair<int, int>,
852  unsigned int,
853  unsigned int,
856 
865  std::vector<std::tuple<unsigned int, unsigned int, IntersectionType>>
867 
871  std::vector<unsigned int> recv_ptrs;
872 
888  template <int dim>
890  spacedim>
892  const unsigned int n_points_1D,
894  const Mapping<dim, spacedim> &mapping,
895  const bool consistent_numbering_of_sender_and_receiver = false) const;
896 
897  private:
905  std::map<unsigned int, std::vector<unsigned int>>
907  const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
908  &point_recv_components,
909  const MPI_Comm comm) const;
910  };
911 
919  template <int structdim, int dim, int spacedim>
922  const Cache<dim, spacedim> &cache,
923  const std::vector<std::vector<Point<spacedim>>> &intersection_requests,
924  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
925  const std::vector<bool> &marked_vertices,
926  const double tolerance);
927 
928  } // namespace internal
929 
939  template <int spacedim>
940  unsigned int
941  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
942  const Point<spacedim> &p);
943 
970  template <int dim, template <int, int> class MeshType, int spacedim>
972  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
973  unsigned int find_closest_vertex(
974  const MeshType<dim, spacedim> &mesh,
975  const Point<spacedim> &p,
976  const std::vector<bool> &marked_vertices = {});
977 
1004  template <int dim, template <int, int> class MeshType, int spacedim>
1006  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1007  unsigned int find_closest_vertex(
1008  const Mapping<dim, spacedim> &mapping,
1009  const MeshType<dim, spacedim> &mesh,
1010  const Point<spacedim> &p,
1011  const std::vector<bool> &marked_vertices = {});
1012 
1013 
1036  template <typename MeshType>
1037  DEAL_II_CXX20_REQUIRES((concepts::is_triangulation_or_dof_handler<MeshType>))
1038 #ifndef _MSC_VER
1039  std::vector<typename MeshType::active_cell_iterator>
1040 #else
1041  std::vector<
1042  typename ::internal::ActiveCellIterator<MeshType::dimension,
1043  MeshType::space_dimension,
1044  MeshType>::type>
1045 #endif
1046  find_cells_adjacent_to_vertex(const MeshType &container,
1047  const unsigned int vertex_index);
1048 
1114  template <int dim, template <int, int> class MeshType, int spacedim>
1116  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1117 #ifndef _MSC_VER
1118  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1119 #else
1120  std::pair<typename ::internal::
1121  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1122  Point<dim>>
1123 #endif
1125  const MeshType<dim, spacedim> &mesh,
1126  const Point<spacedim> &p,
1127  const std::vector<bool> &marked_vertices = {},
1128  const double tolerance = 1.e-10);
1129 
1140  template <int dim, template <int, int> class MeshType, int spacedim>
1142  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1143 #ifndef _MSC_VER
1144  typename MeshType<dim, spacedim>::active_cell_iterator
1145 #else
1146  typename ::internal::
1147  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1148 #endif
1149  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1150  const Point<spacedim> &p,
1151  const std::vector<bool> &marked_vertices = {},
1152  const double tolerance = 1.e-10);
1153 
1160  template <int dim, int spacedim>
1161  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1162  Point<dim>>
1164  const hp::MappingCollection<dim, spacedim> &mapping,
1165  const DoFHandler<dim, spacedim> &mesh,
1166  const Point<spacedim> &p,
1167  const double tolerance = 1.e-10);
1168 
1220  template <int dim, int spacedim>
1221  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1222  Point<dim>>
1224  const Cache<dim, spacedim> &cache,
1225  const Point<spacedim> &p,
1228  const std::vector<bool> &marked_vertices = {},
1229  const double tolerance = 1.e-10);
1230 
1247  template <int dim, template <int, int> class MeshType, int spacedim>
1249  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1250 #ifndef _MSC_VER
1251  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1252 #else
1253  std::pair<typename ::internal::
1254  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1255  Point<dim>>
1256 #endif
1258  const Mapping<dim, spacedim> &mapping,
1259  const MeshType<dim, spacedim> &mesh,
1260  const Point<spacedim> &p,
1261  const std::vector<
1262  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1264  const std::vector<std::vector<Tensor<1, spacedim>>>
1265  &vertex_to_cell_centers,
1266  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1267  typename MeshType<dim, spacedim>::active_cell_iterator(),
1268  const std::vector<bool> &marked_vertices = {},
1269  const RTree<std::pair<Point<spacedim>, unsigned int>> &
1270  used_vertices_rtree = RTree<std::pair<Point<spacedim>, unsigned int>>{},
1271  const double tolerance = 1.e-10,
1272  const RTree<
1273  std::pair<BoundingBox<spacedim>,
1275  *relevant_cell_bounding_boxes_rtree = nullptr);
1276 
1306  template <int dim, template <int, int> class MeshType, int spacedim>
1308  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1309 #ifndef _MSC_VER
1310  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1311  Point<dim>>>
1312 #else
1313  std::vector<std::pair<
1314  typename ::internal::
1315  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1316  Point<dim>>>
1317 #endif
1319  const Mapping<dim, spacedim> &mapping,
1320  const MeshType<dim, spacedim> &mesh,
1321  const Point<spacedim> &p,
1322  const double tolerance,
1323  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1324  Point<dim>> &first_cell,
1325  const std::vector<
1326  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1327  *vertex_to_cells = nullptr);
1328 
1338  template <int dim, template <int, int> class MeshType, int spacedim>
1340  (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1341 #ifndef _MSC_VER
1342  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1343  Point<dim>>>
1344 #else
1345  std::vector<std::pair<
1346  typename ::internal::
1347  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1348  Point<dim>>>
1349 #endif
1351  const Mapping<dim, spacedim> &mapping,
1352  const MeshType<dim, spacedim> &mesh,
1353  const Point<spacedim> &p,
1354  const double tolerance = 1e-10,
1355  const std::vector<bool> &marked_vertices = {});
1356 
1380  template <typename MeshType>
1381  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1382  std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
1383  const typename MeshType::cell_iterator &cell);
1384 
1411  template <typename MeshType>
1414  const typename MeshType::active_cell_iterator &cell,
1415  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1416 
1468  template <typename MeshType>
1470  std::
1471  vector<typename MeshType::active_cell_iterator> compute_active_cell_halo_layer(
1472  const MeshType &mesh,
1473  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1474  &predicate);
1475 
1476 
1486  template <typename MeshType>
1488  std::
1489  vector<typename MeshType::cell_iterator> compute_cell_halo_layer_on_level(
1490  const MeshType &mesh,
1491  const std::function<bool(const typename MeshType::cell_iterator &)>
1492  &predicate,
1493  const unsigned int level);
1494 
1495 
1510  template <typename MeshType>
1512  std::vector<
1513  typename MeshType::
1514  active_cell_iterator> compute_ghost_cell_halo_layer(const MeshType &mesh);
1515 
1566  template <typename MeshType>
1568  std::
1569  vector<typename MeshType::active_cell_iterator> compute_active_cell_layer_within_distance(
1570  const MeshType &mesh,
1571  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1572  &predicate,
1573  const double layer_thickness);
1574 
1599  template <typename MeshType>
1601  std::
1602  vector<typename MeshType::active_cell_iterator> compute_ghost_cell_layer_within_distance(
1603  const MeshType &mesh,
1604  const double layer_thickness);
1605 
1678  template <typename MeshType>
1680  std::
1681  vector<BoundingBox<MeshType::space_dimension>> compute_mesh_predicate_bounding_box(
1682  const MeshType &mesh,
1683  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1684  &predicate,
1685  const unsigned int refinement_level = 0,
1686  const bool allow_merge = false,
1687  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1688 
1716  template <int spacedim>
1717 #ifndef DOXYGEN
1718  std::tuple<std::vector<std::vector<unsigned int>>,
1719  std::map<unsigned int, unsigned int>,
1720  std::map<unsigned int, std::vector<unsigned int>>>
1721 #else
1722  return_type
1723 #endif
1725  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1726  const std::vector<Point<spacedim>> &points);
1727 
1728 
1763  template <int spacedim>
1764 #ifndef DOXYGEN
1765  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1766  std::map<unsigned int, unsigned int>,
1767  std::map<unsigned int, std::vector<unsigned int>>>
1768 #else
1769  return_type
1770 #endif
1772  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1773  const std::vector<Point<spacedim>> &points);
1774 
1787  template <int dim, int spacedim>
1788  std::vector<std::vector<Tensor<1, spacedim>>>
1790  const Triangulation<dim, spacedim> &mesh,
1791  const std::vector<
1793  &vertex_to_cells);
1794 
1795 
1803  template <int dim, int spacedim>
1804  unsigned int
1807  const Point<spacedim> &position,
1808  const Mapping<dim, spacedim> &mapping =
1809  (ReferenceCells::get_hypercube<dim>()
1810 #ifndef _MSC_VER
1811  .template get_default_linear_mapping<dim, spacedim>()
1812 #else
1813  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1814 #endif
1815  ));
1816 
1828  template <int dim, int spacedim>
1829  std::map<unsigned int, types::global_vertex_index>
1832 
1859  template <int dim, int spacedim>
1860  void
1861  partition_triangulation(const unsigned int n_partitions,
1863  const SparsityTools::Partitioner partitioner =
1865 
1876  template <int dim, int spacedim>
1877  void
1878  partition_triangulation(const unsigned int n_partitions,
1879  const std::vector<unsigned int> &cell_weights,
1881  const SparsityTools::Partitioner partitioner =
1883 
1929  template <int dim, int spacedim>
1930  void
1931  partition_triangulation(const unsigned int n_partitions,
1932  const SparsityPattern &cell_connection_graph,
1934  const SparsityTools::Partitioner partitioner =
1936 
1947  template <int dim, int spacedim>
1948  void
1949  partition_triangulation(const unsigned int n_partitions,
1950  const std::vector<unsigned int> &cell_weights,
1951  const SparsityPattern &cell_connection_graph,
1953  const SparsityTools::Partitioner partitioner =
1955 
1970  template <int dim, int spacedim>
1971  void
1972  partition_triangulation_zorder(const unsigned int n_partitions,
1974  const bool group_siblings = true);
1975 
1987  template <int dim, int spacedim>
1988  void
1990 
1998  template <int dim, int spacedim>
1999  std::vector<types::subdomain_id>
2001  const std::vector<CellId> &cell_ids);
2002 
2013  template <int dim, int spacedim>
2014  void
2016  std::vector<types::subdomain_id> &subdomain);
2017 
2032  template <int dim, int spacedim>
2033  unsigned int
2036  const types::subdomain_id subdomain);
2037 
2067  template <int dim, int spacedim>
2068  std::vector<bool>
2070 
2097  template <int dim, int spacedim>
2101  &distorted_cells,
2103 
2104 
2105 
2155  template <typename MeshType>
2156  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2157  std::vector<typename MeshType::active_cell_iterator> get_patch_around_cell(
2158  const typename MeshType::active_cell_iterator &cell);
2159 
2160 
2182  template <class Container>
2183  std::vector<typename Container::cell_iterator>
2185  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2186 
2253  template <class Container>
2254  void
2256  const std::vector<typename Container::active_cell_iterator> &patch,
2258  &local_triangulation,
2259  std::map<
2260  typename Triangulation<Container::dimension,
2261  Container::space_dimension>::active_cell_iterator,
2262  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2263 
2295  template <int dim, int spacedim>
2296  std::map<
2298  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2300 
2301 
2316  template <typename CellIterator>
2318  {
2322  CellIterator cell[2];
2323 
2328  unsigned int face_idx[2];
2329 
2335  std::bitset<3> orientation;
2336 
2350 
2354  std::size_t
2356  };
2357 
2358 
2425  template <typename FaceIterator>
2426  std::optional<std::bitset<3>>
2428  const FaceIterator &face1,
2429  const FaceIterator &face2,
2430  const unsigned int direction,
2434 
2493  template <typename MeshType>
2494  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2496  const MeshType &mesh,
2497  const types::boundary_id b_id1,
2498  const types::boundary_id b_id2,
2499  const unsigned int direction,
2501  &matched_pairs,
2502  const Tensor<1, MeshType::space_dimension> &offset =
2505 
2506 
2531  template <typename MeshType>
2532  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2534  const MeshType &mesh,
2535  const types::boundary_id b_id,
2536  const unsigned int direction,
2537  std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
2538  &matched_pairs,
2539  const ::Tensor<1, MeshType::space_dimension> &offset =
2540  ::Tensor<1, MeshType::space_dimension>(),
2541  const FullMatrix<double> &matrix = FullMatrix<double>());
2542 
2569  template <int dim, int spacedim>
2570  void
2572  const bool reset_boundary_ids = false);
2573 
2595  template <int dim, int spacedim>
2596  void
2598  const std::vector<types::boundary_id> &src_boundary_ids,
2599  const std::vector<types::manifold_id> &dst_manifold_ids,
2600  Triangulation<dim, spacedim> &tria,
2601  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2602 
2632  template <int dim, int spacedim>
2633  void
2635  const bool compute_face_ids = false);
2636 
2661  template <int dim, int spacedim>
2662  void
2665  const std::function<types::manifold_id(
2666  const std::set<types::manifold_id> &)> &disambiguation_function =
2667  [](const std::set<types::manifold_id> &manifold_ids) {
2668  if (manifold_ids.size() == 1)
2669  return *manifold_ids.begin();
2670  else
2672  },
2673  bool overwrite_only_flat_manifold_ids = true);
2762  template <typename DataType, typename MeshType>
2763  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2765  const MeshType &mesh,
2766  const std::function<std::optional<DataType>(
2767  const typename MeshType::active_cell_iterator &)> &pack,
2768  const std::function<void(const typename MeshType::active_cell_iterator &,
2769  const DataType &)> &unpack,
2770  const std::function<bool(const typename MeshType::active_cell_iterator &)>
2771  &cell_filter =
2772  always_return<typename MeshType::active_cell_iterator, bool>{true});
2773 
2786  template <typename DataType, typename MeshType>
2787  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2789  const MeshType &mesh,
2790  const std::function<std::optional<DataType>(
2791  const typename MeshType::level_cell_iterator &)> &pack,
2792  const std::function<void(const typename MeshType::level_cell_iterator &,
2793  const DataType &)> &unpack,
2794  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
2795  cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
2796  true});
2797 
2798  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2799  * boxes @p local_bboxes.
2800  *
2801  * This function is meant to exchange bounding boxes describing the locally
2802  * owned cells in a distributed triangulation obtained with the function
2803  * GridTools::compute_mesh_predicate_bounding_box .
2804  *
2805  * The output vector's size is the number of processes of the MPI
2806  * communicator:
2807  * its i-th entry contains the vector @p local_bboxes of the i-th process.
2808  */
2809  template <int spacedim>
2810  std::vector<std::vector<BoundingBox<spacedim>>>
2812  const std::vector<BoundingBox<spacedim>> &local_bboxes,
2813  const MPI_Comm mpi_communicator);
2814 
2847  template <int spacedim>
2850  const std::vector<BoundingBox<spacedim>> &local_description,
2851  const MPI_Comm mpi_communicator);
2852 
2870  template <int dim, int spacedim>
2871  void
2874  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
2875  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
2876 
2896  template <int dim, int spacedim>
2897  std::map<unsigned int, std::set<::types::subdomain_id>>
2900 
2921  template <int dim, typename VectorType>
2923  {
2924  public:
2928  using value_type = typename VectorType::value_type;
2929 
2934  const FiniteElement<dim, dim> &fe,
2935  const unsigned int n_subdivisions = 1,
2936  const double tolerance = 1e-10);
2937 
2948  void
2949  process(const DoFHandler<dim> &background_dof_handler,
2950  const VectorType &ls_vector,
2951  const double iso_level,
2952  std::vector<Point<dim>> &vertices,
2953  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
2954 
2959  void
2960  process(const DoFHandler<dim> &background_dof_handler,
2961  const VectorType &ls_vector,
2962  const double iso_level,
2963  std::vector<Point<dim>> &vertices) const;
2964 
2974  void
2976  const VectorType &ls_vector,
2977  const double iso_level,
2978  std::vector<Point<dim>> &vertices,
2979  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
2985  void
2987  const VectorType &ls_vector,
2988  const double iso_level,
2989  std::vector<Point<dim>> &vertices) const;
2990 
2991  private:
2996  static Quadrature<dim>
2997  create_quadrature_rule(const unsigned int n_subdivisions);
2998 
3002  void
3003  process_cell(std::vector<value_type> &ls_values,
3004  const std::vector<Point<dim>> &points,
3005  const double iso_level,
3006  std::vector<Point<dim>> &vertices,
3007  std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
3008  const bool write_back_cell_data = true) const;
3009 
3013  void
3014  process_sub_cell(const std::vector<value_type> &,
3015  const std::vector<Point<1>> &,
3016  const std::vector<unsigned int> &,
3017  const double,
3018  std::vector<Point<1>> &,
3019  std::vector<CellData<1>> &,
3020  const bool) const
3021  {
3022  AssertThrow(false, ExcNotImplemented());
3023  }
3024 
3031  void
3032  process_sub_cell(const std::vector<value_type> &ls_values,
3033  const std::vector<Point<2>> &points,
3034  const std::vector<unsigned int> &mask,
3035  const double iso_level,
3036  std::vector<Point<2>> &vertices,
3037  std::vector<CellData<1>> &cells,
3038  const bool write_back_cell_data) const;
3039 
3043  void
3044  process_sub_cell(const std::vector<value_type> &ls_values,
3045  const std::vector<Point<3>> &points,
3046  const std::vector<unsigned int> &mask,
3047  const double iso_level,
3048  std::vector<Point<3>> &vertices,
3049  std::vector<CellData<2>> &cells,
3050  const bool write_back_cell_data) const;
3051 
3056  const unsigned int n_subdivisions;
3057 
3062  const double tolerance;
3063 
3069  };
3070 
3071 
3072 
3082  int,
3083  << "The number of partitions you gave is " << arg1
3084  << ", but must be greater than zero.");
3089  int,
3090  << "The subdomain id " << arg1
3091  << " has no cells associated with it.");
3096 
3101  double,
3102  << "The scaling factor must be positive, but it is " << arg1
3103  << '.');
3104 
3109  unsigned int,
3110  << "The given vertex with index " << arg1
3111  << " is not used in the given triangulation.");
3112 
3115 } /*namespace GridTools*/
3116 
3117 /* ----------------- Template function --------------- */
3118 
3119 #ifndef DOXYGEN
3120 
3121 namespace GridTools
3122 {
3123  template <int dim>
3124  double
3125  cell_measure(
3126  const std::vector<Point<dim>> &all_vertices,
3127  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3128  {
3129  // We forward call to the ArrayView version:
3130  const ArrayView<const unsigned int> view(
3132  return cell_measure(all_vertices, view);
3133  }
3134 
3135 
3136 
3137  // This specialization is defined here so that the general template in the
3138  // source file doesn't need to have further 1d overloads for the internal
3139  // functions it calls.
3140  template <>
3144  {
3145  return {};
3146  }
3147 
3148 
3149 
3150  template <int dim, typename Transformation, int spacedim>
3152  (std::invocable<Transformation, Point<spacedim>> &&
3153  std::assignable_from<
3154  Point<spacedim> &,
3155  std::invoke_result_t<Transformation, Point<spacedim>>>))
3156  void transform(const Transformation &transformation,
3157  Triangulation<dim, spacedim> &triangulation)
3158  {
3159  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3160 
3161  // loop over all active cells, and
3162  // transform those vertices that
3163  // have not yet been touched. note
3164  // that we get to all vertices in
3165  // the triangulation by only
3166  // visiting the active cells.
3168  cell = triangulation.begin_active(),
3169  endc = triangulation.end();
3170  for (; cell != endc; ++cell)
3171  for (const unsigned int v : cell->vertex_indices())
3172  if (treated_vertices[cell->vertex_index(v)] == false)
3173  {
3174  // transform this vertex
3175  cell->vertex(v) = transformation(cell->vertex(v));
3176  // and mark it as treated
3177  treated_vertices[cell->vertex_index(v)] = true;
3178  };
3179 
3180 
3181  // now fix any vertices on hanging nodes so that we don't create any holes
3182  if (dim == 2)
3183  {
3185  cell = triangulation.begin_active(),
3186  endc = triangulation.end();
3187  for (; cell != endc; ++cell)
3188  for (const unsigned int face : cell->face_indices())
3189  if (cell->face(face)->has_children() &&
3190  !cell->face(face)->at_boundary())
3191  {
3192  Assert(cell->reference_cell() ==
3193  ReferenceCells::get_hypercube<dim>(),
3194  ExcNotImplemented());
3195 
3196  // this line has children
3197  cell->face(face)->child(0)->vertex(1) =
3198  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3199  2;
3200  }
3201  }
3202  else if (dim == 3)
3203  {
3205  cell = triangulation.begin_active(),
3206  endc = triangulation.end();
3207  for (; cell != endc; ++cell)
3208  for (const unsigned int face : cell->face_indices())
3209  if (cell->face(face)->has_children() &&
3210  !cell->face(face)->at_boundary())
3211  {
3212  Assert(cell->reference_cell() ==
3213  ReferenceCells::get_hypercube<dim>(),
3214  ExcNotImplemented());
3215 
3216  // this face has hanging nodes
3217  cell->face(face)->child(0)->vertex(1) =
3218  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3219  2.0;
3220  cell->face(face)->child(0)->vertex(2) =
3221  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3222  2.0;
3223  cell->face(face)->child(1)->vertex(3) =
3224  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3225  2.0;
3226  cell->face(face)->child(2)->vertex(3) =
3227  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3228  2.0;
3229 
3230  // center of the face
3231  cell->face(face)->child(0)->vertex(3) =
3232  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3233  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3234  4.0;
3235  }
3236  }
3237 
3238  // Make sure FEValues notices that the mesh has changed
3239  triangulation.signals.mesh_movement();
3240  }
3241 
3242 
3243 
3244  template <typename MeshType>
3245  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3246  std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
3247  const typename MeshType::cell_iterator &cell)
3248  {
3249  std::vector<typename MeshType::active_cell_iterator> child_cells;
3250 
3251  if (cell->has_children())
3252  {
3253  for (unsigned int child = 0; child < cell->n_children(); ++child)
3254  if (cell->child(child)->has_children())
3255  {
3256  const std::vector<typename MeshType::active_cell_iterator>
3257  children = get_active_child_cells<MeshType>(cell->child(child));
3258  child_cells.insert(child_cells.end(),
3259  children.begin(),
3260  children.end());
3261  }
3262  else
3263  child_cells.push_back(cell->child(child));
3264  }
3265 
3266  return child_cells;
3267  }
3268 
3269 
3270 
3271  template <typename MeshType>
3272  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3273  void get_active_neighbors(
3274  const typename MeshType::active_cell_iterator &cell,
3275  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3276  {
3277  active_neighbors.clear();
3278  for (const unsigned int n : cell->face_indices())
3279  if (!cell->at_boundary(n))
3280  {
3281  if (MeshType::dimension == 1)
3282  {
3283  // check children of neighbor. note
3284  // that in 1d children of the neighbor
3285  // may be further refined. In 1d the
3286  // case is simple since we know what
3287  // children bound to the present cell
3288  typename MeshType::cell_iterator neighbor_child =
3289  cell->neighbor(n);
3290  if (!neighbor_child->is_active())
3291  {
3292  while (neighbor_child->has_children())
3293  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3294 
3295  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3296  ExcInternalError());
3297  }
3298  active_neighbors.push_back(neighbor_child);
3299  }
3300  else
3301  {
3302  if (cell->face(n)->has_children())
3303  // this neighbor has children. find
3304  // out which border to the present
3305  // cell
3306  for (unsigned int c = 0;
3307  c < cell->face(n)->n_active_descendants();
3308  ++c)
3309  active_neighbors.push_back(
3310  cell->neighbor_child_on_subface(n, c));
3311  else
3312  {
3313  // the neighbor must be active
3314  // himself
3315  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3316  active_neighbors.push_back(cell->neighbor(n));
3317  }
3318  }
3319  }
3320  }
3321 
3322 
3323 
3324  template <typename CellIterator>
3325  std::size_t
3327  {
3328  return sizeof(*this) + matrix.memory_consumption();
3329  }
3330 
3331 
3332 
3333  namespace internal
3334  {
3335  template <typename DataType,
3336  typename MeshType,
3337  typename MeshCellIteratorType>
3338  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3339  inline void exchange_cell_data(
3340  const MeshType &mesh,
3341  const std::function<std::optional<DataType>(const MeshCellIteratorType &)>
3342  &pack,
3343  const std::function<void(const MeshCellIteratorType &, const DataType &)>
3344  &unpack,
3345  const std::function<bool(const MeshCellIteratorType &)> &cell_filter,
3346  const std::function<void(
3347  const std::function<void(const MeshCellIteratorType &,
3348  const types::subdomain_id)> &)> &process_cells,
3349  const std::function<std::set<types::subdomain_id>(
3350  const parallel::TriangulationBase<MeshType::dimension,
3351  MeshType::space_dimension> &)>
3352  &compute_ghost_owners)
3353  {
3354 # ifndef DEAL_II_WITH_MPI
3355  (void)mesh;
3356  (void)pack;
3357  (void)unpack;
3358  (void)cell_filter;
3359  (void)process_cells;
3360  (void)compute_ghost_owners;
3361  Assert(false, ExcNeedsMPI());
3362 # else
3363  constexpr int dim = MeshType::dimension;
3364  constexpr int spacedim = MeshType::space_dimension;
3365  auto tria =
3366  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3367  &mesh.get_triangulation());
3368  Assert(
3369  tria != nullptr,
3370  ExcMessage(
3371  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3372 
3373  if (const auto tria = dynamic_cast<
3375  &mesh.get_triangulation()))
3376  {
3377  Assert(
3378  tria->with_artificial_cells(),
3379  ExcMessage(
3380  "The functions GridTools::exchange_cell_data_to_ghosts() and "
3381  "GridTools::exchange_cell_data_to_level_ghosts() can only "
3382  "operate on a single layer of ghost cells. However, you have "
3383  "given a Triangulation object of type "
3384  "parallel::shared::Triangulation without artificial cells "
3385  "resulting in an arbitrary number of ghost layers. "
3386  "To use this function for a Triangulation object of type "
3387  "parallel::shared::Triangulation, make sure to create the "
3388  "Triangulation object with allow_artificial_cells set to true. "
3389  "This results in a parallel::shared::Triangulation with only "
3390  "a single layer of ghost cells."));
3391  }
3392 
3393  // build list of cells to request for each neighbor
3394  std::set<::types::subdomain_id> ghost_owners =
3395  compute_ghost_owners(*tria);
3396  std::map<::types::subdomain_id,
3397  std::vector<typename CellId::binary_type>>
3398  neighbor_cell_list;
3399 
3400  for (const auto ghost_owner : ghost_owners)
3401  neighbor_cell_list[ghost_owner] = {};
3402 
3403  process_cells([&](const auto &cell, const auto key) -> void {
3404  if (cell_filter(cell))
3405  {
3406  constexpr int spacedim = MeshType::space_dimension;
3407  neighbor_cell_list[key].emplace_back(
3408  cell->id().template to_binary<spacedim>());
3409  }
3410  });
3411 
3412  Assert(ghost_owners.size() == neighbor_cell_list.size(),
3413  ExcInternalError());
3414 
3415 
3416  // Before sending & receiving, make sure we protect this section with
3417  // a mutex:
3418  static Utilities::MPI::CollectiveMutex mutex;
3420  mutex, tria->get_communicator());
3421 
3422  const int mpi_tag =
3424  const int mpi_tag_reply =
3426 
3427  // send our requests
3428  std::vector<MPI_Request> requests(ghost_owners.size());
3429  {
3430  unsigned int idx = 0;
3431  for (const auto &it : neighbor_cell_list)
3432  {
3433  // send the data about the relevant cells
3434  const int ierr = MPI_Isend(it.second.data(),
3435  it.second.size() * sizeof(it.second[0]),
3436  MPI_BYTE,
3437  it.first,
3438  mpi_tag,
3440  &requests[idx]);
3441  AssertThrowMPI(ierr);
3442  ++idx;
3443  }
3444  }
3445 
3446  // receive requests and reply with the results
3447  std::vector<MPI_Request> reply_requests(ghost_owners.size());
3448  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
3449 
3450  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
3451  {
3452  MPI_Status status;
3453  int ierr = MPI_Probe(MPI_ANY_SOURCE,
3454  mpi_tag,
3456  &status);
3457  AssertThrowMPI(ierr);
3458 
3459  int len;
3460  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3461  AssertThrowMPI(ierr);
3462  Assert(len % sizeof(typename CellId::binary_type) == 0,
3463  ExcInternalError());
3464 
3465  const unsigned int n_cells =
3466  len / sizeof(typename CellId::binary_type);
3467  std::vector<typename CellId::binary_type> cells_with_requests(
3468  n_cells);
3469  std::vector<DataType> data_to_send;
3470  data_to_send.reserve(n_cells);
3471  std::vector<bool> cell_carries_data(n_cells, false);
3472 
3473  ierr = MPI_Recv(cells_with_requests.data(),
3474  len,
3475  MPI_BYTE,
3476  status.MPI_SOURCE,
3477  status.MPI_TAG,
3479  &status);
3480  AssertThrowMPI(ierr);
3481 
3482  // store data for each cell
3483  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
3484  {
3485  const auto cell =
3486  tria->create_cell_iterator(CellId(cells_with_requests[c]));
3487 
3488  MeshCellIteratorType mesh_it(tria,
3489  cell->level(),
3490  cell->index(),
3491  &mesh);
3492 
3493  const std::optional<DataType> data = pack(mesh_it);
3494  if (data)
3495  {
3496  data_to_send.emplace_back(*data);
3497  cell_carries_data[c] = true;
3498  }
3499  }
3500 
3501  // collect data for sending the reply in a buffer
3502 
3503  // (a) make room for storing the local offsets in case we receive
3504  // other data
3505  sendbuffers[idx].resize(sizeof(std::size_t));
3506 
3507  // (b) append the actual data and store how much memory it
3508  // corresponds to, which we then insert into the leading position of
3509  // the sendbuffer
3510  std::size_t size_of_send =
3511  Utilities::pack(data_to_send,
3512  sendbuffers[idx],
3513  /*enable_compression*/ false);
3514  std::memcpy(sendbuffers[idx].data(),
3515  &size_of_send,
3516  sizeof(std::size_t));
3517 
3518  // (c) append information of certain cells that got left out in case
3519  // we need it
3520  if (data_to_send.size() < n_cells)
3521  Utilities::pack(cell_carries_data,
3522  sendbuffers[idx],
3523  /*enable_compression*/ false);
3524 
3525  // send data
3526  ierr = MPI_Isend(sendbuffers[idx].data(),
3527  sendbuffers[idx].size(),
3528  MPI_BYTE,
3529  status.MPI_SOURCE,
3530  mpi_tag_reply,
3532  &reply_requests[idx]);
3533  AssertThrowMPI(ierr);
3534  }
3535 
3536  // finally receive the replies
3537  std::vector<char> receive;
3538  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
3539  {
3540  MPI_Status status;
3541  int ierr = MPI_Probe(MPI_ANY_SOURCE,
3542  mpi_tag_reply,
3544  &status);
3545  AssertThrowMPI(ierr);
3546 
3547  int len;
3548  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3549  AssertThrowMPI(ierr);
3550 
3551  receive.resize(len);
3552 
3553  ierr = MPI_Recv(receive.data(),
3554  len,
3555  MPI_BYTE,
3556  status.MPI_SOURCE,
3557  status.MPI_TAG,
3559  &status);
3560  AssertThrowMPI(ierr);
3561 
3562  // (a) first determine the length of the data section in the
3563  // received buffer
3564  auto data_iterator = receive.begin();
3565  std::size_t size_of_received_data =
3566  Utilities::unpack<std::size_t>(data_iterator,
3567  data_iterator + sizeof(std::size_t));
3568  data_iterator += sizeof(std::size_t);
3569 
3570  // (b) unpack the data section in the indicated region
3571  auto received_data = Utilities::unpack<std::vector<DataType>>(
3572  data_iterator,
3573  data_iterator + size_of_received_data,
3574  /*enable_compression*/ false);
3575  data_iterator += size_of_received_data;
3576 
3577  // (c) check if the received data contained fewer entries than the
3578  // number of cells we identified in the beginning, in which case we
3579  // need to extract the boolean vector with the relevant information
3580  const std::vector<typename CellId::binary_type> &this_cell_list =
3581  neighbor_cell_list[status.MPI_SOURCE];
3582  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
3583  std::vector<bool> cells_with_data;
3584  if (received_data.size() < this_cell_list.size())
3585  {
3586  cells_with_data = Utilities::unpack<std::vector<bool>>(
3587  data_iterator, receive.end(), /*enable_compression*/ false);
3588  AssertDimension(cells_with_data.size(), this_cell_list.size());
3589  }
3590 
3591  // (d) go through the received data and call the user-provided
3592  // unpack function
3593  auto received_data_iterator = received_data.begin();
3594  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
3595  if (cells_with_data.empty() || cells_with_data[c])
3596  {
3598  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
3599 
3600  MeshCellIteratorType cell(tria,
3601  tria_cell->level(),
3602  tria_cell->index(),
3603  &mesh);
3604 
3605  unpack(cell, *received_data_iterator);
3606  ++received_data_iterator;
3607  }
3608  }
3609 
3610  // make sure that all communication is finished
3611  // when we leave this function.
3612  if (requests.size() > 0)
3613  {
3614  const int ierr =
3615  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
3616  AssertThrowMPI(ierr);
3617  }
3618  if (reply_requests.size() > 0)
3619  {
3620  const int ierr = MPI_Waitall(reply_requests.size(),
3621  reply_requests.data(),
3622  MPI_STATUSES_IGNORE);
3623  AssertThrowMPI(ierr);
3624  }
3625 
3626 
3627 # endif // DEAL_II_WITH_MPI
3628  }
3629 
3630  } // namespace internal
3631 
3632  template <typename DataType, typename MeshType>
3633  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3634  inline void exchange_cell_data_to_ghosts(
3635  const MeshType &mesh,
3636  const std::function<std::optional<DataType>(
3637  const typename MeshType::active_cell_iterator &)> &pack,
3638  const std::function<void(const typename MeshType::active_cell_iterator &,
3639  const DataType &)> &unpack,
3640  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3641  &cell_filter)
3642  {
3643 # ifndef DEAL_II_WITH_MPI
3644  (void)mesh;
3645  (void)pack;
3646  (void)unpack;
3647  (void)cell_filter;
3648  Assert(false, ExcNeedsMPI());
3649 # else
3650  internal::exchange_cell_data<DataType,
3651  MeshType,
3652  typename MeshType::active_cell_iterator>(
3653  mesh,
3654  pack,
3655  unpack,
3656  cell_filter,
3657  [&](const auto &process) {
3658  for (const auto &cell : mesh.active_cell_iterators())
3659  if (cell->is_ghost())
3660  process(cell, cell->subdomain_id());
3661  },
3662  [](const auto &tria) { return tria.ghost_owners(); });
3663 # endif
3664  }
3665 
3666 
3667 
3668  template <typename DataType, typename MeshType>
3669  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3671  const MeshType &mesh,
3672  const std::function<std::optional<DataType>(
3673  const typename MeshType::level_cell_iterator &)> &pack,
3674  const std::function<void(const typename MeshType::level_cell_iterator &,
3675  const DataType &)> &unpack,
3676  const std::function<bool(const typename MeshType::level_cell_iterator &)>
3677  &cell_filter)
3678  {
3679 # ifndef DEAL_II_WITH_MPI
3680  (void)mesh;
3681  (void)pack;
3682  (void)unpack;
3683  (void)cell_filter;
3684  Assert(false, ExcNeedsMPI());
3685 # else
3686  internal::exchange_cell_data<DataType,
3687  MeshType,
3688  typename MeshType::level_cell_iterator>(
3689  mesh,
3690  pack,
3691  unpack,
3692  cell_filter,
3693  [&](const auto &process) {
3694  for (const auto &cell : mesh.cell_iterators())
3695  if (cell->is_ghost_on_level())
3696  process(cell, cell->level_subdomain_id());
3697  },
3698  [](const auto &tria) { return tria.level_ghost_owners(); });
3699 # endif
3700  }
3701 } // namespace GridTools
3702 
3703 #endif // DOXYGEN
3704 
3706 
3707 #endif
Definition: cell_id.h:72
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:81
typename VectorType::value_type value_type
Definition: grid_tools.h:2928
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:5109
void process_cell(const typename DoFHandler< dim >::active_cell_iterator &cell, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 >> &cells) const
Definition: grid_tools.cc:5145
void process_sub_cell(const std::vector< value_type > &, const std::vector< Point< 1 >> &, const std::vector< unsigned int > &, const double, std::vector< Point< 1 >> &, std::vector< CellData< 1 >> &, const bool) const
Definition: grid_tools.h:3014
const unsigned int n_subdivisions
Definition: grid_tools.h:3056
MarchingCubeAlgorithm(const Mapping< dim, dim > &mapping, const FiniteElement< dim, dim > &fe, const unsigned int n_subdivisions=1, const double tolerance=1e-10)
Definition: grid_tools.cc:5057
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:5074
Definition: tensor.h:516
virtual MPI_Comm get_communicator() const
cell_iterator create_cell_iterator(const CellId &cell_id) const
cell_iterator begin(const unsigned int level=0) const
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:108
#define DEAL_II_DEPRECATED
Definition: config.h:178
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
unsigned int level
Definition: grid_out.cc:4617
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:472
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1947
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:517
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:441
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:2806
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:2898
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:2831
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:2926
void random(DoFHandler< dim, spacedim > &dof_handler)
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition: grid_tools.cc:3942
DistributedComputeIntersectionLocationsInternal< structdim, spacedim > distributed_compute_intersection_locations(const Cache< dim, spacedim > &cache, const std::vector< std::vector< Point< spacedim >>> &intersection_requests, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance)
Definition: grid_tools.cc:4323
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:1168
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:257
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:2007
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:4776
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:2764
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:201
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:4571
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &local_points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const double tolerance=1e-10, const std::vector< bool > &marked_vertices={}, const bool enforce_unique_mapping=true)
Definition: grid_tools.cc:3529
double cell_measure(const std::vector< Point< dim >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:1080
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2298
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:3014
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:191
std::vector< typename MeshType::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType &container, const unsigned int vertex_index)
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_halo_layer(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const unsigned int direction, std::vector< PeriodicFacePair< typename MeshType::cell_iterator >> &matched_pairs, const Tensor< 1, MeshType::space_dimension > &offset=::Tensor< 1, MeshType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
std::optional< std::bitset< 3 > > orthogonal_equality(const FaceIterator &face1, const FaceIterator &face2, const unsigned int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_layer_within_distance(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const double layer_thickness)
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:3359
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm mpi_communicator)
Definition: grid_tools.cc:4616
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2113
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:761
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:4899
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:1315
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:2282
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:3330
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)
std::vector< typename MeshType::cell_iterator > compute_cell_halo_layer_on_level(const MeshType &mesh, const std::function< bool(const typename MeshType::cell_iterator &)> &predicate, const unsigned int level)
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:1768
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
std::vector< types::subdomain_id > get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const std::vector< CellId > &cell_ids)
Definition: grid_tools.cc:2212
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition: grid_tools.cc:415
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)
std::vector< std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > > find_all_active_cells_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const double tolerance, const std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim >> &first_cell, const std::vector< std::set< typename MeshType< dim, spacedim >::active_cell_iterator >> *vertex_to_cells=nullptr)
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:1413
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm mpi_communicator)
Definition: grid_tools.cc:4713
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
@ exchange_cell_data_request
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
@ exchange_cell_data_reply
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1343
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1500
constexpr bool is_triangulation_or_dof_handler
concept is_triangulation_or_dof_handler
Definition: hp.h:118
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:14894
static const unsigned int invalid_unsigned_int
Definition: types.h:221
const types::manifold_id flat_manifold_id
Definition: types.h:326
Definition: types.h:33
unsigned int manifold_id
Definition: types.h:157
unsigned int global_dof_index
Definition: types.h:82
unsigned int subdomain_id
Definition: types.h:44
unsigned int boundary_id
Definition: types.h:145
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:144
std::bitset< 3 > orientation
Definition: grid_tools.h:2335
std::size_t memory_consumption() const
FullMatrix< double > matrix
Definition: grid_tools.h:2349
unsigned int face_idx[2]
Definition: grid_tools.h:2328
GridTools::internal::DistributedComputePointLocationsInternal< dim, spacedim > convert_to_distributed_compute_point_locations_internal(const unsigned int n_points_1D, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping, const bool consistent_numbering_of_sender_and_receiver=false) const
std::map< unsigned int, std::vector< unsigned int > > communicate_indices(const std::vector< std::tuple< unsigned int, unsigned int, unsigned int >> &point_recv_components, const MPI_Comm comm) const
Definition: grid_tools.cc:4211
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, IntersectionType > > send_components
Definition: grid_tools.h:855
std::vector< std::tuple< unsigned int, unsigned int, IntersectionType > > recv_components
Definition: grid_tools.h:866
std::array<::Point< spacedim >, structdim+1 > IntersectionType
Definition: grid_tools.h:841
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:791
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:766
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
const MPI_Comm comm
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
const ::Triangulation< dim, spacedim > & tria