deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50:00+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\}}\)
Loading...
Searching...
No Matches
grid_tools.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_grid_tools_h
16#define dealii_grid_tools_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/point.h>
24
26
28
30
32#include <deal.II/fe/mapping.h>
33
37#include <deal.II/grid/tria.h>
40
45
47
48#include <boost/archive/binary_iarchive.hpp>
49#include <boost/archive/binary_oarchive.hpp>
50#include <boost/random/mersenne_twister.hpp>
51#include <boost/serialization/array.hpp>
52#include <boost/serialization/vector.hpp>
53
54#ifdef DEAL_II_WITH_ZLIB
55# include <boost/iostreams/device/back_inserter.hpp>
56# include <boost/iostreams/filter/gzip.hpp>
57# include <boost/iostreams/filtering_stream.hpp>
58# include <boost/iostreams/stream.hpp>
59#endif
60
61#include <optional>
62#include <set>
63
64#ifdef DEAL_II_HAVE_CXX20
65# include <concepts>
66#endif
67
68
70
71// Forward declarations
72#ifndef DOXYGEN
73namespace parallel
74{
75 namespace distributed
76 {
77 template <int dim, int spacedim>
79 class Triangulation;
80 }
81} // namespace parallel
82
83namespace hp
84{
85 template <int, int>
86 class MappingCollection;
87}
88
89class SparsityPattern;
90
91namespace GridTools
92{
93 template <int dim, int spacedim>
94 class Cache;
95}
96#endif
97
98namespace internal
99{
100 template <int dim, int spacedim, typename MeshType>
103 {
104 public:
105#ifndef _MSC_VER
106 using type = typename MeshType::active_cell_iterator;
107#else
109#endif
110 };
111
112#ifdef _MSC_VER
113 template <int dim, int spacedim>
114 class ActiveCellIterator<dim, spacedim, DoFHandler<dim, spacedim>>
115 {
116 public:
117 using type =
119 };
120#endif
121} // namespace internal
122
131namespace GridTools
132{
219 template <int dim, typename Transformation, int spacedim>
221 (std::invocable<Transformation, Point<spacedim>> &&
222 std::assignable_from<
224 std::invoke_result_t<Transformation, Point<spacedim>>>))
225 void transform(const Transformation &transformation,
226 Triangulation<dim, spacedim> &triangulation);
227
234 template <int dim, int spacedim>
235 void
236 shift(const Tensor<1, spacedim> &shift_vector,
237 Triangulation<dim, spacedim> &triangulation);
238
239
250 template <int dim, int spacedim>
251 void
252 rotate(const double angle, Triangulation<dim, spacedim> &triangulation);
253
266 template <int dim>
267 void
268 rotate(const Tensor<1, 3, double> &axis,
269 const double angle,
271
329 template <int dim>
330 void
331 laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
332 Triangulation<dim> &tria,
333 const Function<dim, double> *coefficient = nullptr,
334 const bool solve_for_absolute_positions = false);
335
343 template <int dim, int spacedim>
344 void
345 scale(const double scaling_factor,
346 Triangulation<dim, spacedim> &triangulation);
347
377 template <int dim, int spacedim>
378 void
380 const double factor,
381 Triangulation<dim, spacedim> &triangulation,
382 const bool keep_boundary = true,
383 const unsigned int seed = boost::random::mt19937::default_seed);
384
474 template <int dim, int spacedim>
475 void
476 regularize_corner_cells(Triangulation<dim, spacedim> &tria,
477 const double limit_angle_fraction = .75);
478
542 template <int dim, int spacedim>
543#ifndef DOXYGEN
544 std::tuple<
545 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
546 std::vector<std::vector<Point<dim>>>,
547 std::vector<std::vector<unsigned int>>>
548#else
549 return_type
550#endif
552 const Cache<dim, spacedim> &cache,
553 const std::vector<Point<spacedim>> &points,
555 &cell_hint =
557
591 template <int dim, int spacedim>
592#ifndef DOXYGEN
593 std::tuple<
594 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
595 std::vector<std::vector<Point<dim>>>,
596 std::vector<std::vector<unsigned int>>,
597 std::vector<unsigned int>>
598#else
599 return_type
600#endif
602 const Cache<dim, spacedim> &cache,
603 const std::vector<Point<spacedim>> &points,
605 &cell_hint =
607
687 template <int dim, int spacedim>
688#ifndef DOXYGEN
689 std::tuple<
690 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
691 std::vector<std::vector<Point<dim>>>,
692 std::vector<std::vector<unsigned int>>,
693 std::vector<std::vector<Point<spacedim>>>,
694 std::vector<std::vector<unsigned int>>>
695#else
696 return_type
697#endif
700 const std::vector<Point<spacedim>> &local_points,
701 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
702 const double tolerance = 1e-10,
703 const std::vector<bool> &marked_vertices = {},
704 const bool enforce_unique_mapping = true);
705
706 namespace internal
707 {
722 template <int dim, int spacedim>
724 {
726
733 void
735
739 unsigned int n_searched_points;
740
747 std::vector<std::tuple<std::pair<int, int>,
748 unsigned int,
749 unsigned int,
752 unsigned int>>
754
758 std::vector<unsigned int> send_ranks;
759
765 std::vector<unsigned int> send_ptrs;
766
777 std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
779
783 std::vector<unsigned int> recv_ranks;
784
790 std::vector<unsigned int> recv_ptrs;
791 };
792
802 template <int dim, int spacedim>
806 const std::vector<Point<spacedim>> &points,
807 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
808 const std::vector<bool> &marked_vertices,
809 const double tolerance,
810 const bool perform_handshake,
811 const bool enforce_unique_mapping = false);
812
813
821 template <int structdim, int spacedim>
823 {
828 std::array<::Point<spacedim>, structdim + 1>;
829
838 std::vector<std::tuple<std::pair<int, int>,
839 unsigned int,
840 unsigned int,
843
852 std::vector<std::tuple<unsigned int, unsigned int, IntersectionType>>
854
858 std::vector<unsigned int> recv_ptrs;
859
880 template <int dim>
882 spacedim>
884 const unsigned int n_points_1D,
886 const Mapping<dim, spacedim> &mapping,
887 std::vector<Quadrature<spacedim>> *mapped_quadratures_recv_comp =
888 nullptr,
889 const bool consistent_numbering_of_sender_and_receiver = false) const;
890
891 private:
899 std::map<unsigned int, std::vector<unsigned int>>
901 const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
902 &point_recv_components,
903 const MPI_Comm comm) const;
904 };
905
913 template <int structdim, int dim, int spacedim>
916 const Cache<dim, spacedim> &cache,
917 const std::vector<std::vector<Point<spacedim>>> &intersection_requests,
918 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
919 const std::vector<bool> &marked_vertices,
920 const double tolerance);
921
922 } // namespace internal
923
933 template <int spacedim>
934 unsigned int
935 find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
936 const Point<spacedim> &p);
937
964 template <int dim, template <int, int> class MeshType, int spacedim>
966 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
967 unsigned int find_closest_vertex(
968 const MeshType<dim, spacedim> &mesh,
969 const Point<spacedim> &p,
970 const std::vector<bool> &marked_vertices = {});
971
998 template <int dim, template <int, int> class MeshType, int spacedim>
1000 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1001 unsigned int find_closest_vertex(
1002 const Mapping<dim, spacedim> &mapping,
1003 const MeshType<dim, spacedim> &mesh,
1004 const Point<spacedim> &p,
1005 const std::vector<bool> &marked_vertices = {});
1006
1007
1030 template <typename MeshType>
1032#ifndef _MSC_VER
1033 std::vector<typename MeshType::active_cell_iterator>
1034#else
1035 std::vector<
1036 typename ::internal::ActiveCellIterator<MeshType::dimension,
1037 MeshType::space_dimension,
1038 MeshType>::type>
1039#endif
1040 find_cells_adjacent_to_vertex(const MeshType &container,
1041 const unsigned int vertex_index);
1042
1108 template <int dim, template <int, int> class MeshType, int spacedim>
1110 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1111#ifndef _MSC_VER
1112 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1113#else
1114 std::pair<typename ::internal::
1115 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1116 Point<dim>>
1117#endif
1119 const MeshType<dim, spacedim> &mesh,
1120 const Point<spacedim> &p,
1121 const std::vector<bool> &marked_vertices = {},
1122 const double tolerance = 1.e-10);
1123
1134 template <int dim, template <int, int> class MeshType, int spacedim>
1136 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1137#ifndef _MSC_VER
1138 typename MeshType<dim, spacedim>::active_cell_iterator
1139#else
1140 typename ::internal::
1141 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1142#endif
1143 find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1144 const Point<spacedim> &p,
1145 const std::vector<bool> &marked_vertices = {},
1146 const double tolerance = 1.e-10);
1147
1154 template <int dim, int spacedim>
1155 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1156 Point<dim>>
1159 const DoFHandler<dim, spacedim> &mesh,
1160 const Point<spacedim> &p,
1161 const double tolerance = 1.e-10);
1162
1214 template <int dim, int spacedim>
1215 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1216 Point<dim>>
1218 const Cache<dim, spacedim> &cache,
1219 const Point<spacedim> &p,
1222 const std::vector<bool> &marked_vertices = {},
1223 const double tolerance = 1.e-10);
1224
1241 template <int dim, template <int, int> class MeshType, int spacedim>
1243 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1244#ifndef _MSC_VER
1245 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1246#else
1247 std::pair<typename ::internal::
1248 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1249 Point<dim>>
1250#endif
1252 const Mapping<dim, spacedim> &mapping,
1253 const MeshType<dim, spacedim> &mesh,
1254 const Point<spacedim> &p,
1255 const std::vector<
1256 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1258 const std::vector<std::vector<Tensor<1, spacedim>>>
1259 &vertex_to_cell_centers,
1260 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1261 typename MeshType<dim, spacedim>::active_cell_iterator(),
1262 const std::vector<bool> &marked_vertices = {},
1263 const RTree<std::pair<Point<spacedim>, unsigned int>> &
1264 used_vertices_rtree = RTree<std::pair<Point<spacedim>, unsigned int>>{},
1265 const double tolerance = 1.e-10,
1266 const RTree<
1267 std::pair<BoundingBox<spacedim>,
1269 *relevant_cell_bounding_boxes_rtree = nullptr);
1270
1300 template <int dim, template <int, int> class MeshType, int spacedim>
1302 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1303#ifndef _MSC_VER
1304 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1305 Point<dim>>>
1306#else
1307 std::vector<std::pair<
1308 typename ::internal::
1309 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1310 Point<dim>>>
1311#endif
1313 const Mapping<dim, spacedim> &mapping,
1314 const MeshType<dim, spacedim> &mesh,
1315 const Point<spacedim> &p,
1316 const double tolerance,
1317 const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1318 Point<dim>> &first_cell,
1319 const std::vector<
1320 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1321 *vertex_to_cells = nullptr);
1322
1332 template <int dim, template <int, int> class MeshType, int spacedim>
1334 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1335#ifndef _MSC_VER
1336 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1337 Point<dim>>>
1338#else
1339 std::vector<std::pair<
1340 typename ::internal::
1341 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1342 Point<dim>>>
1343#endif
1345 const Mapping<dim, spacedim> &mapping,
1346 const MeshType<dim, spacedim> &mesh,
1347 const Point<spacedim> &p,
1348 const double tolerance = 1e-10,
1349 const std::vector<bool> &marked_vertices = {});
1350
1374 template <typename MeshType>
1376 std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
1377 const typename MeshType::cell_iterator &cell);
1378
1405 template <typename MeshType>
1406 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1408 const typename MeshType::active_cell_iterator &cell,
1409 std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1410
1462 template <typename MeshType>
1463 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1464 std::
1465 vector<typename MeshType::active_cell_iterator> compute_active_cell_halo_layer(
1466 const MeshType &mesh,
1467 const std::function<bool(const typename MeshType::active_cell_iterator &)>
1468 &predicate);
1469
1470
1480 template <typename MeshType>
1481 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1482 std::
1483 vector<typename MeshType::cell_iterator> compute_cell_halo_layer_on_level(
1484 const MeshType &mesh,
1485 const std::function<bool(const typename MeshType::cell_iterator &)>
1486 &predicate,
1487 const unsigned int level);
1488
1489
1504 template <typename MeshType>
1505 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1506 std::vector<
1507 typename MeshType::
1508 active_cell_iterator> compute_ghost_cell_halo_layer(const MeshType &mesh);
1509
1560 template <typename MeshType>
1561 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1562 std::
1563 vector<typename MeshType::active_cell_iterator> compute_active_cell_layer_within_distance(
1564 const MeshType &mesh,
1565 const std::function<bool(const typename MeshType::active_cell_iterator &)>
1566 &predicate,
1567 const double layer_thickness);
1568
1593 template <typename MeshType>
1594 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1595 std::
1596 vector<typename MeshType::active_cell_iterator> compute_ghost_cell_layer_within_distance(
1597 const MeshType &mesh,
1598 const double layer_thickness);
1599
1672 template <typename MeshType>
1673 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1674 std::
1675 vector<BoundingBox<MeshType::space_dimension>> compute_mesh_predicate_bounding_box(
1676 const MeshType &mesh,
1677 const std::function<bool(const typename MeshType::active_cell_iterator &)>
1678 &predicate,
1679 const unsigned int refinement_level = 0,
1680 const bool allow_merge = false,
1681 const unsigned int max_boxes = numbers::invalid_unsigned_int);
1682
1710 template <int spacedim>
1711#ifndef DOXYGEN
1712 std::tuple<std::vector<std::vector<unsigned int>>,
1713 std::map<unsigned int, unsigned int>,
1714 std::map<unsigned int, std::vector<unsigned int>>>
1715#else
1716 return_type
1717#endif
1719 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1720 const std::vector<Point<spacedim>> &points);
1721
1722
1757 template <int spacedim>
1758#ifndef DOXYGEN
1759 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1760 std::map<unsigned int, unsigned int>,
1761 std::map<unsigned int, std::vector<unsigned int>>>
1762#else
1763 return_type
1764#endif
1766 const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1767 const std::vector<Point<spacedim>> &points);
1768
1781 template <int dim, int spacedim>
1782 std::vector<std::vector<Tensor<1, spacedim>>>
1784 const Triangulation<dim, spacedim> &mesh,
1785 const std::vector<
1787 &vertex_to_cells);
1788
1789
1797 template <int dim, int spacedim>
1798 unsigned int
1801 const Point<spacedim> &position,
1802 const Mapping<dim, spacedim> &mapping =
1803 (ReferenceCells::get_hypercube<dim>()
1804#ifndef _MSC_VER
1805 .template get_default_linear_mapping<dim, spacedim>()
1806#else
1807 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1808#endif
1809 ));
1810
1822 template <int dim, int spacedim>
1823 std::map<unsigned int, types::global_vertex_index>
1826
1853 template <int dim, int spacedim>
1854 void
1855 partition_triangulation(const unsigned int n_partitions,
1857 const SparsityTools::Partitioner partitioner =
1859
1870 template <int dim, int spacedim>
1871 void
1872 partition_triangulation(const unsigned int n_partitions,
1873 const std::vector<unsigned int> &cell_weights,
1875 const SparsityTools::Partitioner partitioner =
1877
1923 template <int dim, int spacedim>
1924 void
1925 partition_triangulation(const unsigned int n_partitions,
1926 const SparsityPattern &cell_connection_graph,
1928 const SparsityTools::Partitioner partitioner =
1930
1941 template <int dim, int spacedim>
1942 void
1943 partition_triangulation(const unsigned int n_partitions,
1944 const std::vector<unsigned int> &cell_weights,
1945 const SparsityPattern &cell_connection_graph,
1947 const SparsityTools::Partitioner partitioner =
1949
1964 template <int dim, int spacedim>
1965 void
1966 partition_triangulation_zorder(const unsigned int n_partitions,
1968 const bool group_siblings = true);
1969
1981 template <int dim, int spacedim>
1982 void
1984
1992 template <int dim, int spacedim>
1993 std::vector<types::subdomain_id>
1995 const std::vector<CellId> &cell_ids);
1996
2007 template <int dim, int spacedim>
2008 void
2010 std::vector<types::subdomain_id> &subdomain);
2011
2026 template <int dim, int spacedim>
2027 unsigned int
2030 const types::subdomain_id subdomain);
2031
2061 template <int dim, int spacedim>
2062 std::vector<bool>
2064
2091 template <int dim, int spacedim>
2095 &distorted_cells,
2097
2098
2099
2149 template <typename MeshType>
2151 std::vector<typename MeshType::active_cell_iterator> get_patch_around_cell(
2152 const typename MeshType::active_cell_iterator &cell);
2153
2154
2176 template <class Container>
2177 std::vector<typename Container::cell_iterator>
2179 const std::vector<typename Container::active_cell_iterator> &patch_cells);
2180
2247 template <class Container>
2248 void
2250 const std::vector<typename Container::active_cell_iterator> &patch,
2252 &local_triangulation,
2253 std::map<
2254 typename Triangulation<Container::dimension,
2255 Container::space_dimension>::active_cell_iterator,
2256 typename Container::active_cell_iterator> &patch_to_global_tria_map);
2257
2289 template <int dim, int spacedim>
2290 std::map<
2292 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2294
2295
2310 template <typename CellIterator>
2312 {
2316 CellIterator cell[2];
2317
2322 unsigned int face_idx[2];
2323
2330
2344
2348 std::size_t
2350 };
2351
2352
2376 template <typename FaceIterator>
2377 std::optional<types::geometric_orientation>
2379 const FaceIterator &face1,
2380 const FaceIterator &face2,
2381 const unsigned int direction,
2384 const FullMatrix<double> &matrix = FullMatrix<double>(),
2385 const double abs_tol = 1e-10);
2386
2446 template <typename MeshType>
2449 const MeshType &mesh,
2450 const types::boundary_id b_id1,
2451 const types::boundary_id b_id2,
2452 const unsigned int direction,
2454 &matched_pairs,
2457 const FullMatrix<double> &matrix = FullMatrix<double>(),
2458 const double abs_tol = 1e-10);
2459
2460
2481 template <typename MeshType>
2484 const MeshType &mesh,
2485 const types::boundary_id b_id,
2486 const unsigned int direction,
2487 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
2488 &matched_pairs,
2489 const ::Tensor<1, MeshType::space_dimension> &offset =
2490 ::Tensor<1, MeshType::space_dimension>(),
2491 const FullMatrix<double> &matrix = FullMatrix<double>(),
2492 const double abs_tol = 1e-10);
2493
2520 template <int dim, int spacedim>
2521 void
2522 copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
2523 const bool reset_boundary_ids = false);
2524
2546 template <int dim, int spacedim>
2547 void
2549 const std::vector<types::boundary_id> &src_boundary_ids,
2550 const std::vector<types::manifold_id> &dst_manifold_ids,
2551 Triangulation<dim, spacedim> &tria,
2552 const std::vector<types::boundary_id> &reset_boundary_ids = {});
2553
2583 template <int dim, int spacedim>
2584 void
2586 const bool compute_face_ids = false);
2587
2612 template <int dim, int spacedim>
2613 void
2616 const std::function<types::manifold_id(
2617 const std::set<types::manifold_id> &)> &disambiguation_function =
2618 [](const std::set<types::manifold_id> &manifold_ids) {
2619 if (manifold_ids.size() == 1)
2620 return *manifold_ids.begin();
2621 else
2623 },
2624 bool overwrite_only_flat_manifold_ids = true);
2713 template <typename DataType, typename MeshType>
2716 const MeshType &mesh,
2717 const std::function<std::optional<DataType>(
2718 const typename MeshType::active_cell_iterator &)> &pack,
2719 const std::function<void(const typename MeshType::active_cell_iterator &,
2720 const DataType &)> &unpack,
2721 const std::function<bool(const typename MeshType::active_cell_iterator &)>
2722 &cell_filter =
2723 always_return<typename MeshType::active_cell_iterator, bool>{true});
2724
2737 template <typename DataType, typename MeshType>
2740 const MeshType &mesh,
2741 const std::function<std::optional<DataType>(
2742 const typename MeshType::level_cell_iterator &)> &pack,
2743 const std::function<void(const typename MeshType::level_cell_iterator &,
2744 const DataType &)> &unpack,
2745 const std::function<bool(const typename MeshType::level_cell_iterator &)> &
2746 cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
2747 true});
2748
2749 /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2750 * boxes @p local_bboxes.
2751 *
2752 * This function is meant to exchange bounding boxes describing the locally
2753 * owned cells in a distributed triangulation obtained with the function
2754 * GridTools::compute_mesh_predicate_bounding_box .
2755 *
2756 * The output vector's size is the number of processes of the MPI
2757 * communicator:
2758 * its i-th entry contains the vector @p local_bboxes of the i-th process.
2759 */
2760 template <int spacedim>
2761 std::vector<std::vector<BoundingBox<spacedim>>>
2763 const std::vector<BoundingBox<spacedim>> &local_bboxes,
2764 const MPI_Comm mpi_communicator);
2765
2798 template <int spacedim>
2801 const std::vector<BoundingBox<spacedim>> &local_description,
2802 const MPI_Comm mpi_communicator);
2803
2821 template <int dim, int spacedim>
2822 void
2824 const Triangulation<dim, spacedim> &tria,
2825 std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
2826 std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
2827
2847 template <int dim, int spacedim>
2848 std::map<unsigned int, std::set<::types::subdomain_id>>
2850 const Triangulation<dim, spacedim> &tria);
2851
2872 template <int dim, typename VectorType>
2874 {
2875 public:
2879 using value_type = typename VectorType::value_type;
2880
2885 const FiniteElement<dim, dim> &fe,
2886 const unsigned int n_subdivisions = 1,
2887 const double tolerance = 1e-10);
2888
2899 void
2900 process(const DoFHandler<dim> &background_dof_handler,
2901 const VectorType &ls_vector,
2902 const double iso_level,
2903 std::vector<Point<dim>> &vertices,
2904 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
2905
2910 void
2911 process(const DoFHandler<dim> &background_dof_handler,
2912 const VectorType &ls_vector,
2913 const double iso_level,
2914 std::vector<Point<dim>> &vertices) const;
2915
2925 void
2927 const VectorType &ls_vector,
2928 const double iso_level,
2929 std::vector<Point<dim>> &vertices,
2930 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
2936 void
2938 const VectorType &ls_vector,
2939 const double iso_level,
2940 std::vector<Point<dim>> &vertices) const;
2941
2942 private:
2947 static Quadrature<dim>
2948 create_quadrature_rule(const unsigned int n_subdivisions);
2949
2953 void
2954 process_cell(std::vector<value_type> &ls_values,
2955 const std::vector<Point<dim>> &points,
2956 const double iso_level,
2957 std::vector<Point<dim>> &vertices,
2958 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
2959 const bool write_back_cell_data = true) const;
2960
2964 void
2965 process_sub_cell(const std::vector<value_type> &,
2966 const std::vector<Point<1>> &,
2967 const std::vector<unsigned int> &,
2968 const double,
2969 std::vector<Point<1>> &,
2970 std::vector<CellData<1>> &,
2971 const bool) const
2972 {
2974 }
2975
2982 void
2983 process_sub_cell(const std::vector<value_type> &ls_values,
2984 const std::vector<Point<2>> &points,
2985 const std::vector<unsigned int> &mask,
2986 const double iso_level,
2987 std::vector<Point<2>> &vertices,
2988 std::vector<CellData<1>> &cells,
2989 const bool write_back_cell_data) const;
2990
2994 void
2995 process_sub_cell(const std::vector<value_type> &ls_values,
2996 const std::vector<Point<3>> &points,
2997 const std::vector<unsigned int> &mask,
2998 const double iso_level,
2999 std::vector<Point<3>> &vertices,
3000 std::vector<CellData<2>> &cells,
3001 const bool write_back_cell_data) const;
3002
3007 const unsigned int n_subdivisions;
3008
3013 const double tolerance;
3014
3020 };
3021
3022
3023
3033 int,
3034 << "The number of partitions you gave is " << arg1
3035 << ", but must be greater than zero.");
3040 int,
3041 << "The subdomain id " << arg1
3042 << " has no cells associated with it.");
3047
3052 double,
3053 << "The scaling factor must be positive, but it is " << arg1
3054 << '.');
3055
3060 unsigned int,
3061 << "The given vertex with index " << arg1
3062 << " is not used in the given triangulation.");
3063
3066} /*namespace GridTools*/
3067
3068/* ----------------- Template function --------------- */
3069
3070#ifndef DOXYGEN
3071
3072namespace GridTools
3073{
3074 template <int dim>
3075 double
3077 const std::vector<Point<dim>> &all_vertices,
3078 const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3079 {
3080 // We forward call to the ArrayView version:
3083 return cell_measure(all_vertices, view);
3084 }
3085
3086
3087
3088 // This specialization is defined here so that the general template in the
3089 // source file doesn't need to have further 1d overloads for the internal
3090 // functions it calls.
3091 template <>
3095 {
3096 return {};
3097 }
3098
3099
3100
3101 template <int dim, typename Transformation, int spacedim>
3103 (std::invocable<Transformation, Point<spacedim>> &&
3104 std::assignable_from<
3106 std::invoke_result_t<Transformation, Point<spacedim>>>))
3107 void transform(const Transformation &transformation,
3108 Triangulation<dim, spacedim> &triangulation)
3109 {
3110 std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3111
3112 // loop over all active cells, and
3113 // transform those vertices that
3114 // have not yet been touched. note
3115 // that we get to all vertices in
3116 // the triangulation by only
3117 // visiting the active cells.
3119 cell = triangulation.begin_active(),
3120 endc = triangulation.end();
3121 for (; cell != endc; ++cell)
3122 for (const unsigned int v : cell->vertex_indices())
3123 if (treated_vertices[cell->vertex_index(v)] == false)
3124 {
3125 // transform this vertex
3126 cell->vertex(v) = transformation(cell->vertex(v));
3127 // and mark it as treated
3128 treated_vertices[cell->vertex_index(v)] = true;
3129 };
3130
3131
3132 // now fix any vertices on hanging nodes so that we don't create any holes
3133 if (dim == 2)
3134 {
3136 cell = triangulation.begin_active(),
3137 endc = triangulation.end();
3138 for (; cell != endc; ++cell)
3139 for (const unsigned int face : cell->face_indices())
3140 if (cell->face(face)->has_children() &&
3141 !cell->face(face)->at_boundary())
3142 {
3143 Assert(cell->reference_cell() ==
3144 ReferenceCells::get_hypercube<dim>(),
3146
3147 // this line has children
3148 cell->face(face)->child(0)->vertex(1) =
3149 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3150 2.0;
3151 }
3152 }
3153 else if (dim == 3)
3154 {
3156 cell = triangulation.begin_active(),
3157 endc = triangulation.end();
3158 for (; cell != endc; ++cell)
3159 for (const unsigned int face : cell->face_indices())
3160 if (cell->face(face)->has_children() &&
3161 !cell->face(face)->at_boundary())
3162 {
3163 if (static_cast<std::uint8_t>(
3164 cell->face(face)->refinement_case()) ==
3166 {
3167 Assert(cell->reference_cell() ==
3168 ReferenceCells::get_hypercube<dim>(),
3170
3171 // this face has hanging nodes
3172 cell->face(face)->child(0)->vertex(1) =
3173 (cell->face(face)->vertex(0) +
3174 cell->face(face)->vertex(1)) /
3175 2.0;
3176 cell->face(face)->child(0)->vertex(2) =
3177 (cell->face(face)->vertex(0) +
3178 cell->face(face)->vertex(2)) /
3179 2.0;
3180 cell->face(face)->child(1)->vertex(3) =
3181 (cell->face(face)->vertex(1) +
3182 cell->face(face)->vertex(3)) /
3183 2.0;
3184 cell->face(face)->child(2)->vertex(3) =
3185 (cell->face(face)->vertex(2) +
3186 cell->face(face)->vertex(3)) /
3187 2.0;
3188
3189 // center of the face
3190 cell->face(face)->child(0)->vertex(3) =
3191 (cell->face(face)->vertex(0) +
3192 cell->face(face)->vertex(1) +
3193 cell->face(face)->vertex(2) +
3194 cell->face(face)->vertex(3)) /
3195 4.0;
3196 }
3197 else
3198 {
3199 // Special case for anisotropic refinement
3200 for (unsigned int line = 0;
3201 line < GeometryInfo<dim - 1>::faces_per_cell;
3202 line++)
3203 if (cell->face(face)->line(line)->has_children())
3204 cell->face(face)->line(line)->child(0)->vertex(1) =
3205 (cell->face(face)->line(line)->vertex(0) +
3206 cell->face(face)->line(line)->vertex(1)) /
3207 2.0;
3208 }
3209 }
3210 }
3211
3212 // Make sure FEValues notices that the mesh has changed
3213 triangulation.signals.mesh_movement();
3214 }
3215
3216
3217
3218 template <typename MeshType>
3220 std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
3221 const typename MeshType::cell_iterator &cell)
3222 {
3223 std::vector<typename MeshType::active_cell_iterator> child_cells;
3224
3225 if (cell->has_children())
3226 {
3227 for (unsigned int child = 0; child < cell->n_children(); ++child)
3228 if (cell->child(child)->has_children())
3229 {
3230 const std::vector<typename MeshType::active_cell_iterator>
3231 children = get_active_child_cells<MeshType>(cell->child(child));
3232 child_cells.insert(child_cells.end(),
3233 children.begin(),
3234 children.end());
3235 }
3236 else
3237 child_cells.push_back(cell->child(child));
3238 }
3239
3240 return child_cells;
3241 }
3242
3243
3244
3245 template <typename MeshType>
3248 const typename MeshType::active_cell_iterator &cell,
3249 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3250 {
3251 active_neighbors.clear();
3252 for (const unsigned int n : cell->face_indices())
3253 if (!cell->at_boundary(n))
3254 {
3255 if (MeshType::dimension == 1)
3256 {
3257 // check children of neighbor. note
3258 // that in 1d children of the neighbor
3259 // may be further refined. In 1d the
3260 // case is simple since we know what
3261 // children bound to the present cell
3262 typename MeshType::cell_iterator neighbor_child =
3263 cell->neighbor(n);
3264 if (!neighbor_child->is_active())
3265 {
3266 while (neighbor_child->has_children())
3267 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3268
3269 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3271 }
3272 active_neighbors.push_back(neighbor_child);
3273 }
3274 else
3275 {
3276 if (cell->face(n)->has_children())
3277 // this neighbor has children. find
3278 // out which border to the present
3279 // cell
3280 for (unsigned int c = 0;
3281 c < cell->face(n)->n_active_descendants();
3282 ++c)
3283 active_neighbors.push_back(
3284 cell->neighbor_child_on_subface(n, c));
3285 else
3286 {
3287 // the neighbor must be active
3288 // himself
3289 Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3290 active_neighbors.push_back(cell->neighbor(n));
3291 }
3292 }
3293 }
3294 }
3295
3296
3297
3298 template <typename CellIterator>
3299 std::size_t
3301 {
3302 return sizeof(*this) + matrix.memory_consumption();
3303 }
3304
3305
3306
3307 namespace internal
3308 {
3309 template <typename DataType,
3310 typename MeshType,
3311 typename MeshCellIteratorType>
3313 inline void exchange_cell_data(
3314 const MeshType &mesh,
3315 const std::function<std::optional<DataType>(const MeshCellIteratorType &)>
3316 &pack,
3317 const std::function<void(const MeshCellIteratorType &, const DataType &)>
3318 &unpack,
3319 const std::function<bool(const MeshCellIteratorType &)> &cell_filter,
3320 const std::function<void(
3321 const std::function<void(const MeshCellIteratorType &,
3322 const types::subdomain_id)> &)> &process_cells,
3323 const std::function<std::set<types::subdomain_id>(
3324 const parallel::TriangulationBase<MeshType::dimension,
3325 MeshType::space_dimension> &)>
3326 &compute_ghost_owners)
3327 {
3328# ifndef DEAL_II_WITH_MPI
3329 (void)mesh;
3330 (void)pack;
3331 (void)unpack;
3332 (void)cell_filter;
3333 (void)process_cells;
3334 (void)compute_ghost_owners;
3335 Assert(false, ExcNeedsMPI());
3336# else
3337 constexpr int dim = MeshType::dimension;
3338 constexpr int spacedim = MeshType::space_dimension;
3339 auto tria =
3340 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3341 &mesh.get_triangulation());
3342 Assert(
3343 tria != nullptr,
3344 ExcMessage(
3345 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3346
3347 if (const auto tria = dynamic_cast<
3349 &mesh.get_triangulation()))
3350 {
3351 Assert(
3352 tria->with_artificial_cells(),
3353 ExcMessage(
3354 "The functions GridTools::exchange_cell_data_to_ghosts() and "
3355 "GridTools::exchange_cell_data_to_level_ghosts() can only "
3356 "operate on a single layer of ghost cells. However, you have "
3357 "given a Triangulation object of type "
3358 "parallel::shared::Triangulation without artificial cells "
3359 "resulting in an arbitrary number of ghost layers. "
3360 "To use this function for a Triangulation object of type "
3361 "parallel::shared::Triangulation, make sure to create the "
3362 "Triangulation object with allow_artificial_cells set to true. "
3363 "This results in a parallel::shared::Triangulation with only "
3364 "a single layer of ghost cells."));
3365 }
3366
3367 // build list of cells to request for each neighbor
3368 std::set<::types::subdomain_id> ghost_owners =
3369 compute_ghost_owners(*tria);
3370 std::map<::types::subdomain_id,
3371 std::vector<typename CellId::binary_type>>
3372 neighbor_cell_list;
3373
3374 for (const auto ghost_owner : ghost_owners)
3375 neighbor_cell_list[ghost_owner] = {};
3376
3377 process_cells([&](const auto &cell, const auto key) -> void {
3378 if (cell_filter(cell))
3379 {
3380 constexpr int spacedim = MeshType::space_dimension;
3381 neighbor_cell_list[key].emplace_back(
3382 cell->id().template to_binary<spacedim>());
3383 }
3384 });
3385
3386 Assert(ghost_owners.size() == neighbor_cell_list.size(),
3388
3389
3390 // Before sending & receiving, make sure we protect this section with
3391 // a mutex:
3394 mutex, tria->get_mpi_communicator());
3395
3396 const int mpi_tag =
3398 const int mpi_tag_reply =
3400
3401 // send our requests
3402 std::vector<MPI_Request> requests(ghost_owners.size());
3403 {
3404 unsigned int idx = 0;
3405 for (const auto &it : neighbor_cell_list)
3406 {
3407 // send the data about the relevant cells
3408 const int ierr = MPI_Isend(it.second.data(),
3409 it.second.size() * sizeof(it.second[0]),
3410 MPI_BYTE,
3411 it.first,
3412 mpi_tag,
3413 tria->get_mpi_communicator(),
3414 &requests[idx]);
3415 AssertThrowMPI(ierr);
3416 ++idx;
3417 }
3418 }
3419
3420 // receive requests and reply with the results
3421 std::vector<MPI_Request> reply_requests(ghost_owners.size());
3422 std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
3423
3424 for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
3425 {
3426 MPI_Status status;
3427 int ierr = MPI_Probe(MPI_ANY_SOURCE,
3428 mpi_tag,
3429 tria->get_mpi_communicator(),
3430 &status);
3431 AssertThrowMPI(ierr);
3432
3433 int len;
3434 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3435 AssertThrowMPI(ierr);
3436 Assert(len % sizeof(typename CellId::binary_type) == 0,
3438
3439 const unsigned int n_cells =
3440 len / sizeof(typename CellId::binary_type);
3441 std::vector<typename CellId::binary_type> cells_with_requests(
3442 n_cells);
3443 std::vector<DataType> data_to_send;
3444 data_to_send.reserve(n_cells);
3445 std::vector<bool> cell_carries_data(n_cells, false);
3446
3447 ierr = MPI_Recv(cells_with_requests.data(),
3448 len,
3449 MPI_BYTE,
3450 status.MPI_SOURCE,
3451 status.MPI_TAG,
3452 tria->get_mpi_communicator(),
3453 &status);
3454 AssertThrowMPI(ierr);
3455
3456 // store data for each cell
3457 for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
3458 {
3459 const auto cell =
3460 tria->create_cell_iterator(CellId(cells_with_requests[c]));
3461
3462 MeshCellIteratorType mesh_it(tria,
3463 cell->level(),
3464 cell->index(),
3465 &mesh);
3466
3467 const std::optional<DataType> data = pack(mesh_it);
3468 if (data)
3469 {
3470 data_to_send.emplace_back(*data);
3471 cell_carries_data[c] = true;
3472 }
3473 }
3474
3475 // collect data for sending the reply in a buffer
3476
3477 // (a) make room for storing the local offsets in case we receive
3478 // other data
3479 sendbuffers[idx].resize(sizeof(std::size_t));
3480
3481 // (b) append the actual data and store how much memory it
3482 // corresponds to, which we then insert into the leading position of
3483 // the sendbuffer
3484 std::size_t size_of_send =
3485 Utilities::pack(data_to_send,
3486 sendbuffers[idx],
3487 /*enable_compression*/ false);
3488 std::memcpy(sendbuffers[idx].data(),
3489 &size_of_send,
3490 sizeof(std::size_t));
3491
3492 // (c) append information of certain cells that got left out in case
3493 // we need it
3494 if (data_to_send.size() < n_cells)
3495 Utilities::pack(cell_carries_data,
3496 sendbuffers[idx],
3497 /*enable_compression*/ false);
3498
3499 // send data
3500 ierr = MPI_Isend(sendbuffers[idx].data(),
3501 sendbuffers[idx].size(),
3502 MPI_BYTE,
3503 status.MPI_SOURCE,
3504 mpi_tag_reply,
3505 tria->get_mpi_communicator(),
3506 &reply_requests[idx]);
3507 AssertThrowMPI(ierr);
3508 }
3509
3510 // finally receive the replies
3511 std::vector<char> receive;
3512 for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
3513 {
3514 MPI_Status status;
3515 int ierr = MPI_Probe(MPI_ANY_SOURCE,
3516 mpi_tag_reply,
3517 tria->get_mpi_communicator(),
3518 &status);
3519 AssertThrowMPI(ierr);
3520
3521 int len;
3522 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3523 AssertThrowMPI(ierr);
3524
3525 receive.resize(len);
3526
3527 ierr = MPI_Recv(receive.data(),
3528 len,
3529 MPI_BYTE,
3530 status.MPI_SOURCE,
3531 status.MPI_TAG,
3532 tria->get_mpi_communicator(),
3533 &status);
3534 AssertThrowMPI(ierr);
3535
3536 // (a) first determine the length of the data section in the
3537 // received buffer
3538 auto data_iterator = receive.begin();
3539 std::size_t size_of_received_data =
3540 Utilities::unpack<std::size_t>(data_iterator,
3541 data_iterator + sizeof(std::size_t));
3542 data_iterator += sizeof(std::size_t);
3543
3544 // (b) unpack the data section in the indicated region
3545 auto received_data = Utilities::unpack<std::vector<DataType>>(
3546 data_iterator,
3547 data_iterator + size_of_received_data,
3548 /*enable_compression*/ false);
3549 data_iterator += size_of_received_data;
3550
3551 // (c) check if the received data contained fewer entries than the
3552 // number of cells we identified in the beginning, in which case we
3553 // need to extract the boolean vector with the relevant information
3554 const std::vector<typename CellId::binary_type> &this_cell_list =
3555 neighbor_cell_list[status.MPI_SOURCE];
3556 AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
3557 std::vector<bool> cells_with_data;
3558 if (received_data.size() < this_cell_list.size())
3559 {
3560 cells_with_data = Utilities::unpack<std::vector<bool>>(
3561 data_iterator, receive.end(), /*enable_compression*/ false);
3562 AssertDimension(cells_with_data.size(), this_cell_list.size());
3563 }
3564
3565 // (d) go through the received data and call the user-provided
3566 // unpack function
3567 auto received_data_iterator = received_data.begin();
3568 for (unsigned int c = 0; c < this_cell_list.size(); ++c)
3569 if (cells_with_data.empty() || cells_with_data[c])
3570 {
3572 tria_cell = tria->create_cell_iterator(this_cell_list[c]);
3573
3574 MeshCellIteratorType cell(tria,
3575 tria_cell->level(),
3576 tria_cell->index(),
3577 &mesh);
3578
3579 unpack(cell, *received_data_iterator);
3580 ++received_data_iterator;
3581 }
3582 }
3583
3584 // make sure that all communication is finished
3585 // when we leave this function.
3586 if (requests.size() > 0)
3587 {
3588 const int ierr =
3589 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
3590 AssertThrowMPI(ierr);
3591 }
3592 if (reply_requests.size() > 0)
3593 {
3594 const int ierr = MPI_Waitall(reply_requests.size(),
3595 reply_requests.data(),
3596 MPI_STATUSES_IGNORE);
3597 AssertThrowMPI(ierr);
3598 }
3599
3600
3601# endif // DEAL_II_WITH_MPI
3602 }
3603
3604 } // namespace internal
3605
3606 template <typename DataType, typename MeshType>
3608 inline void exchange_cell_data_to_ghosts(
3609 const MeshType &mesh,
3610 const std::function<std::optional<DataType>(
3611 const typename MeshType::active_cell_iterator &)> &pack,
3612 const std::function<void(const typename MeshType::active_cell_iterator &,
3613 const DataType &)> &unpack,
3614 const std::function<bool(const typename MeshType::active_cell_iterator &)>
3615 &cell_filter)
3616 {
3617# ifndef DEAL_II_WITH_MPI
3618 (void)mesh;
3619 (void)pack;
3620 (void)unpack;
3621 (void)cell_filter;
3622 Assert(false, ExcNeedsMPI());
3623# else
3624 internal::exchange_cell_data<DataType,
3625 MeshType,
3626 typename MeshType::active_cell_iterator>(
3627 mesh,
3628 pack,
3629 unpack,
3630 cell_filter,
3631 [&](const auto &process) {
3632 for (const auto &cell : mesh.active_cell_iterators())
3633 if (cell->is_ghost())
3634 process(cell, cell->subdomain_id());
3635 },
3636 [](const auto &tria) { return tria.ghost_owners(); });
3637# endif
3638 }
3639
3640
3641
3642 template <typename DataType, typename MeshType>
3645 const MeshType &mesh,
3646 const std::function<std::optional<DataType>(
3647 const typename MeshType::level_cell_iterator &)> &pack,
3648 const std::function<void(const typename MeshType::level_cell_iterator &,
3649 const DataType &)> &unpack,
3650 const std::function<bool(const typename MeshType::level_cell_iterator &)>
3651 &cell_filter)
3652 {
3653# ifndef DEAL_II_WITH_MPI
3654 (void)mesh;
3655 (void)pack;
3656 (void)unpack;
3657 (void)cell_filter;
3658 Assert(false, ExcNeedsMPI());
3659# else
3660 internal::exchange_cell_data<DataType,
3661 MeshType,
3662 typename MeshType::level_cell_iterator>(
3663 mesh,
3664 pack,
3665 unpack,
3666 cell_filter,
3667 [&](const auto &process) {
3668 for (const auto &cell : mesh.cell_iterators())
3669 if (cell->is_ghost_on_level())
3670 process(cell, cell->level_subdomain_id());
3671 },
3672 [](const auto &tria) { return tria.level_ghost_owners(); });
3673# endif
3674 }
3675} // namespace GridTools
3676
3677#endif // DOXYGEN
3678
3680
3681#endif
std::array< unsigned int, 4 > binary_type
Definition cell_id.h:80
typename VectorType::value_type value_type
const unsigned int n_subdivisions
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
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
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
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
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
cell_iterator create_cell_iterator(const CellId &cell_id) const
virtual MPI_Comm get_mpi_communicator() const
typename MeshType::active_cell_iterator type
Definition grid_tools.h:106
#define DEAL_II_DEPRECATED
Definition config.h:281
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:243
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
unsigned int level
Definition grid_out.cc:4635
unsigned int vertex_indices[2]
#define DeclException0(Exception0)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNeedsMPI()
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
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={})
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)
std::vector< index_type > data
Definition mpi.cc:750
std::size_t size
Definition mpi.cc:749
const MPI_Comm comm
Definition mpi.cc:928
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)
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)
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)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim > > &local_description, const MPI_Comm mpi_communicator)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Triangulation< dim, spacedim >::DistortedCellList fix_up_distorted_child_cells(const typename Triangulation< dim, spacedim >::DistortedCellList &distorted_cells, Triangulation< dim, spacedim > &triangulation)
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)
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())
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 >(), const double abs_tol=1e-10)
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim > > &vertices, const Point< spacedim > &p)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
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 >()))
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)
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
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 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::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)
std::map< unsigned int, types::global_vertex_index > compute_local_to_global_vertex_index_map(const Triangulation< dim, spacedim > &triangulation)
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 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)
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
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)
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const std::vector< Point< spacedim > > &points)
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)
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)
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())
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
std::vector< types::subdomain_id > get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const std::vector< CellId > &cell_ids)
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
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)
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)
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::optional< types::geometric_orientation > 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 >(), const double abs_tol=1e-10)
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim > > &local_bboxes, const MPI_Comm mpi_communicator)
double cell_measure(const std::vector< Point< dim > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
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)
@ matrix
Contents is actually a matrix.
@ 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
std::size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1382
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition utilities.h:1539
Definition hp.h:117
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14906
constexpr types::manifold_id flat_manifold_id
Definition types.h:342
STL namespace.
Definition types.h:32
unsigned int global_dof_index
Definition types.h:94
unsigned int subdomain_id
Definition types.h:52
unsigned char geometric_orientation
Definition types.h:40
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition rtree.h:145
std::size_t memory_consumption() const
FullMatrix< double > matrix
types::geometric_orientation orientation
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
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, std::vector< Quadrature< spacedim > > *mapped_quadratures_recv_comp=nullptr, const bool consistent_numbering_of_sender_and_receiver=false) const
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, IntersectionType > > send_components
Definition grid_tools.h:842
std::vector< std::tuple< unsigned int, unsigned int, IntersectionType > > recv_components
Definition grid_tools.h:853
std::array<::Point< spacedim >, structdim+1 > IntersectionType
Definition grid_tools.h:828
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition grid_tools.h:778
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition grid_tools.h:753