Processing math: 0%
 deal.II version GIT relicensing-2846-g6fb608615f 2025-03-15 04:10: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
2375 template <typename FaceIterator>
2376 std::optional<types::geometric_orientation>
2378 const FaceIterator &face1,
2379 const FaceIterator &face2,
2380 const unsigned int direction,
2383 const FullMatrix<double> &matrix = FullMatrix<double>());
2384
2443 template <typename MeshType>
2446 const MeshType &mesh,
2447 const types::boundary_id b_id1,
2448 const types::boundary_id b_id2,
2449 const unsigned int direction,
2451 &matched_pairs,
2454 const FullMatrix<double> &matrix = FullMatrix<double>());
2455
2456
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
2519 template <int dim, int spacedim>
2520 void
2521 copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
2522 const bool reset_boundary_ids = false);
2523
2545 template <int dim, int spacedim>
2546 void
2548 const std::vector<types::boundary_id> &src_boundary_ids,
2549 const std::vector<types::manifold_id> &dst_manifold_ids,
2550 Triangulation<dim, spacedim> &tria,
2551 const std::vector<types::boundary_id> &reset_boundary_ids = {});
2552
2582 template <int dim, int spacedim>
2583 void
2585 const bool compute_face_ids = false);
2586
2611 template <int dim, int spacedim>
2612 void
2615 const std::function<types::manifold_id(
2616 const std::set<types::manifold_id> &)> &disambiguation_function =
2617 [](const std::set<types::manifold_id> &manifold_ids) {
2618 if (manifold_ids.size() == 1)
2619 return *manifold_ids.begin();
2620 else
2622 },
2623 bool overwrite_only_flat_manifold_ids = true);
2712 template <typename DataType, typename MeshType>
2715 const MeshType &mesh,
2716 const std::function<std::optional<DataType>(
2717 const typename MeshType::active_cell_iterator &)> &pack,
2718 const std::function<void(const typename MeshType::active_cell_iterator &,
2719 const DataType &)> &unpack,
2720 const std::function<bool(const typename MeshType::active_cell_iterator &)>
2721 &cell_filter =
2722 always_return<typename MeshType::active_cell_iterator, bool>{true});
2723
2736 template <typename DataType, typename MeshType>
2739 const MeshType &mesh,
2740 const std::function<std::optional<DataType>(
2741 const typename MeshType::level_cell_iterator &)> &pack,
2742 const std::function<void(const typename MeshType::level_cell_iterator &,
2743 const DataType &)> &unpack,
2744 const std::function<bool(const typename MeshType::level_cell_iterator &)> &
2745 cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
2746 true});
2747
2748 /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2749 * boxes @p local_bboxes.
2750 *
2751 * This function is meant to exchange bounding boxes describing the locally
2752 * owned cells in a distributed triangulation obtained with the function
2753 * GridTools::compute_mesh_predicate_bounding_box .
2754 *
2755 * The output vector's size is the number of processes of the MPI
2756 * communicator:
2757 * its i-th entry contains the vector @p local_bboxes of the i-th process.
2758 */
2759 template <int spacedim>
2760 std::vector<std::vector<BoundingBox<spacedim>>>
2762 const std::vector<BoundingBox<spacedim>> &local_bboxes,
2763 const MPI_Comm mpi_communicator);
2764
2797 template <int spacedim>
2800 const std::vector<BoundingBox<spacedim>> &local_description,
2801 const MPI_Comm mpi_communicator);
2802
2820 template <int dim, int spacedim>
2821 void
2823 const Triangulation<dim, spacedim> &tria,
2824 std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
2825 std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
2826
2846 template <int dim, int spacedim>
2847 std::map<unsigned int, std::set<::types::subdomain_id>>
2849 const Triangulation<dim, spacedim> &tria);
2850
2871 template <int dim, typename VectorType>
2873 {
2874 public:
2878 using value_type = typename VectorType::value_type;
2879
2884 const FiniteElement<dim, dim> &fe,
2885 const unsigned int n_subdivisions = 1,
2886 const double tolerance = 1e-10);
2887
2898 void
2899 process(const DoFHandler<dim> &background_dof_handler,
2900 const VectorType &ls_vector,
2901 const double iso_level,
2902 std::vector<Point<dim>> &vertices,
2903 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
2904
2909 void
2910 process(const DoFHandler<dim> &background_dof_handler,
2911 const VectorType &ls_vector,
2912 const double iso_level,
2913 std::vector<Point<dim>> &vertices) const;
2914
2924 void
2926 const VectorType &ls_vector,
2927 const double iso_level,
2928 std::vector<Point<dim>> &vertices,
2929 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
2935 void
2937 const VectorType &ls_vector,
2938 const double iso_level,
2939 std::vector<Point<dim>> &vertices) const;
2940
2941 private:
2946 static Quadrature<dim>
2947 create_quadrature_rule(const unsigned int n_subdivisions);
2948
2952 void
2953 process_cell(std::vector<value_type> &ls_values,
2954 const std::vector<Point<dim>> &points,
2955 const double iso_level,
2956 std::vector<Point<dim>> &vertices,
2957 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
2958 const bool write_back_cell_data = true) const;
2959
2963 void
2964 process_sub_cell(const std::vector<value_type> &,
2965 const std::vector<Point<1>> &,
2966 const std::vector<unsigned int> &,
2967 const double,
2968 std::vector<Point<1>> &,
2969 std::vector<CellData<1>> &,
2970 const bool) const
2971 {
2973 }
2974
2981 void
2982 process_sub_cell(const std::vector<value_type> &ls_values,
2983 const std::vector<Point<2>> &points,
2984 const std::vector<unsigned int> &mask,
2985 const double iso_level,
2986 std::vector<Point<2>> &vertices,
2987 std::vector<CellData<1>> &cells,
2988 const bool write_back_cell_data) const;
2989
2993 void
2994 process_sub_cell(const std::vector<value_type> &ls_values,
2995 const std::vector<Point<3>> &points,
2996 const std::vector<unsigned int> &mask,
2997 const double iso_level,
2998 std::vector<Point<3>> &vertices,
2999 std::vector<CellData<2>> &cells,
3000 const bool write_back_cell_data) const;
3001
3006 const unsigned int n_subdivisions;
3007
3012 const double tolerance;
3013
3019 };
3020
3021
3022
3032 int,
3033 << "The number of partitions you gave is " << arg1
3034 << ", but must be greater than zero.");
3039 int,
3040 << "The subdomain id " << arg1
3041 << " has no cells associated with it.");
3046
3051 double,
3052 << "The scaling factor must be positive, but it is " << arg1
3053 << '.');
3054
3059 unsigned int,
3060 << "The given vertex with index " << arg1
3061 << " is not used in the given triangulation.");
3062
3065} /*namespace GridTools*/
3066
3067/* ----------------- Template function --------------- */
3068
3069#ifndef DOXYGEN
3070
3071namespace GridTools
3072{
3073 template <int dim>
3074 double
3076 const std::vector<Point<dim>> &all_vertices,
3077 const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3078 {
3079 // We forward call to the ArrayView version:
3082 return cell_measure(all_vertices, view);
3083 }
3084
3085
3086
3087 // This specialization is defined here so that the general template in the
3088 // source file doesn't need to have further 1d overloads for the internal
3089 // functions it calls.
3090 template <>
3094 {
3095 return {};
3096 }
3097
3098
3099
3100 template <int dim, typename Transformation, int spacedim>
3102 (std::invocable<Transformation, Point<spacedim>> &&
3103 std::assignable_from<
3105 std::invoke_result_t<Transformation, Point<spacedim>>>))
3106 void transform(const Transformation &transformation,
3107 Triangulation<dim, spacedim> &triangulation)
3108 {
3109 std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3110
3111 // loop over all active cells, and
3112 // transform those vertices that
3113 // have not yet been touched. note
3114 // that we get to all vertices in
3115 // the triangulation by only
3116 // visiting the active cells.
3118 cell = triangulation.begin_active(),
3119 endc = triangulation.end();
3120 for (; cell != endc; ++cell)
3121 for (const unsigned int v : cell->vertex_indices())
3122 if (treated_vertices[cell->vertex_index(v)] == false)
3123 {
3124 // transform this vertex
3125 cell->vertex(v) = transformation(cell->vertex(v));
3126 // and mark it as treated
3127 treated_vertices[cell->vertex_index(v)] = true;
3128 };
3129
3130
3131 // now fix any vertices on hanging nodes so that we don't create any holes
3132 if (dim == 2)
3133 {
3135 cell = triangulation.begin_active(),
3136 endc = triangulation.end();
3137 for (; cell != endc; ++cell)
3138 for (const unsigned int face : cell->face_indices())
3139 if (cell->face(face)->has_children() &&
3140 !cell->face(face)->at_boundary())
3141 {
3142 Assert(cell->reference_cell() ==
3143 ReferenceCells::get_hypercube<dim>(),
3145
3146 // this line has children
3147 cell->face(face)->child(0)->vertex(1) =
3148 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3149 2.0;
3150 }
3151 }
3152 else if (dim == 3)
3153 {
3155 cell = triangulation.begin_active(),
3156 endc = triangulation.end();
3157 for (; cell != endc; ++cell)
3158 for (const unsigned int face : cell->face_indices())
3159 if (cell->face(face)->has_children() &&
3160 !cell->face(face)->at_boundary())
3161 {
3162 if (static_cast<std::uint8_t>(
3163 cell->face(face)->refinement_case()) ==
3165 {
3166 Assert(cell->reference_cell() ==
3167 ReferenceCells::get_hypercube<dim>(),
3169
3170 // this face has hanging nodes
3171 cell->face(face)->child(0)->vertex(1) =
3172 (cell->face(face)->vertex(0) +
3173 cell->face(face)->vertex(1)) /
3174 2.0;
3175 cell->face(face)->child(0)->vertex(2) =
3176 (cell->face(face)->vertex(0) +
3177 cell->face(face)->vertex(2)) /
3178 2.0;
3179 cell->face(face)->child(1)->vertex(3) =
3180 (cell->face(face)->vertex(1) +
3181 cell->face(face)->vertex(3)) /
3182 2.0;
3183 cell->face(face)->child(2)->vertex(3) =
3184 (cell->face(face)->vertex(2) +
3185 cell->face(face)->vertex(3)) /
3186 2.0;
3187
3188 // center of the face
3189 cell->face(face)->child(0)->vertex(3) =
3190 (cell->face(face)->vertex(0) +
3191 cell->face(face)->vertex(1) +
3192 cell->face(face)->vertex(2) +
3193 cell->face(face)->vertex(3)) /
3194 4.0;
3195 }
3196 else
3197 {
3198 // Special case for anisotropic refinement
3199 for (unsigned int line = 0;
3200 line < GeometryInfo<dim - 1>::faces_per_cell;
3201 line++)
3202 if (cell->face(face)->line(line)->has_children())
3203 cell->face(face)->line(line)->child(0)->vertex(1) =
3204 (cell->face(face)->line(line)->vertex(0) +
3205 cell->face(face)->line(line)->vertex(1)) /
3206 2.0;
3207 }
3208 }
3209 }
3210
3211 // Make sure FEValues notices that the mesh has changed
3212 triangulation.signals.mesh_movement();
3213 }
3214
3215
3216
3217 template <typename MeshType>
3219 std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
3220 const typename MeshType::cell_iterator &cell)
3221 {
3222 std::vector<typename MeshType::active_cell_iterator> child_cells;
3223
3224 if (cell->has_children())
3225 {
3226 for (unsigned int child = 0; child < cell->n_children(); ++child)
3227 if (cell->child(child)->has_children())
3228 {
3229 const std::vector<typename MeshType::active_cell_iterator>
3230 children = get_active_child_cells<MeshType>(cell->child(child));
3231 child_cells.insert(child_cells.end(),
3232 children.begin(),
3233 children.end());
3234 }
3235 else
3236 child_cells.push_back(cell->child(child));
3237 }
3238
3239 return child_cells;
3240 }
3241
3242
3243
3244 template <typename MeshType>
3247 const typename MeshType::active_cell_iterator &cell,
3248 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3249 {
3250 active_neighbors.clear();
3251 for (const unsigned int n : cell->face_indices())
3252 if (!cell->at_boundary(n))
3253 {
3254 if (MeshType::dimension == 1)
3255 {
3256 // check children of neighbor. note
3257 // that in 1d children of the neighbor
3258 // may be further refined. In 1d the
3259 // case is simple since we know what
3260 // children bound to the present cell
3261 typename MeshType::cell_iterator neighbor_child =
3262 cell->neighbor(n);
3263 if (!neighbor_child->is_active())
3264 {
3265 while (neighbor_child->has_children())
3266 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3267
3268 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3270 }
3271 active_neighbors.push_back(neighbor_child);
3272 }
3273 else
3274 {
3275 if (cell->face(n)->has_children())
3276 // this neighbor has children. find
3277 // out which border to the present
3278 // cell
3279 for (unsigned int c = 0;
3280 c < cell->face(n)->n_active_descendants();
3281 ++c)
3282 active_neighbors.push_back(
3283 cell->neighbor_child_on_subface(n, c));
3284 else
3285 {
3286 // the neighbor must be active
3287 // himself
3288 Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3289 active_neighbors.push_back(cell->neighbor(n));
3290 }
3291 }
3292 }
3293 }
3294
3295
3296
3297 template <typename CellIterator>
3298 std::size_t
3300 {
3301 return sizeof(*this) + matrix.memory_consumption();
3302 }
3303
3304
3305
3306 namespace internal
3307 {
3308 template <typename DataType,
3309 typename MeshType,
3310 typename MeshCellIteratorType>
3312 inline void exchange_cell_data(
3313 const MeshType &mesh,
3314 const std::function<std::optional<DataType>(const MeshCellIteratorType &)>
3315 &pack,
3316 const std::function<void(const MeshCellIteratorType &, const DataType &)>
3317 &unpack,
3318 const std::function<bool(const MeshCellIteratorType &)> &cell_filter,
3319 const std::function<void(
3320 const std::function<void(const MeshCellIteratorType &,
3321 const types::subdomain_id)> &)> &process_cells,
3322 const std::function<std::set<types::subdomain_id>(
3323 const parallel::TriangulationBase<MeshType::dimension,
3324 MeshType::space_dimension> &)>
3325 &compute_ghost_owners)
3326 {
3327# ifndef DEAL_II_WITH_MPI
3328 (void)mesh;
3329 (void)pack;
3330 (void)unpack;
3331 (void)cell_filter;
3332 (void)process_cells;
3333 (void)compute_ghost_owners;
3334 Assert(false, ExcNeedsMPI());
3335# else
3336 constexpr int dim = MeshType::dimension;
3337 constexpr int spacedim = MeshType::space_dimension;
3338 auto tria =
3339 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3340 &mesh.get_triangulation());
3341 Assert(
3342 tria != nullptr,
3343 ExcMessage(
3344 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3345
3346 if (const auto tria = dynamic_cast<
3348 &mesh.get_triangulation()))
3349 {
3350 Assert(
3351 tria->with_artificial_cells(),
3352 ExcMessage(
3353 "The functions GridTools::exchange_cell_data_to_ghosts() and "
3354 "GridTools::exchange_cell_data_to_level_ghosts() can only "
3355 "operate on a single layer of ghost cells. However, you have "
3356 "given a Triangulation object of type "
3357 "parallel::shared::Triangulation without artificial cells "
3358 "resulting in an arbitrary number of ghost layers. "
3359 "To use this function for a Triangulation object of type "
3360 "parallel::shared::Triangulation, make sure to create the "
3361 "Triangulation object with allow_artificial_cells set to true. "
3362 "This results in a parallel::shared::Triangulation with only "
3363 "a single layer of ghost cells."));
3364 }
3365
3366 // build list of cells to request for each neighbor
3367 std::set<::types::subdomain_id> ghost_owners =
3368 compute_ghost_owners(*tria);
3369 std::map<::types::subdomain_id,
3370 std::vector<typename CellId::binary_type>>
3371 neighbor_cell_list;
3372
3373 for (const auto ghost_owner : ghost_owners)
3374 neighbor_cell_list[ghost_owner] = {};
3375
3376 process_cells([&](const auto &cell, const auto key) -> void {
3377 if (cell_filter(cell))
3378 {
3379 constexpr int spacedim = MeshType::space_dimension;
3380 neighbor_cell_list[key].emplace_back(
3381 cell->id().template to_binary<spacedim>());
3382 }
3383 });
3384
3385 Assert(ghost_owners.size() == neighbor_cell_list.size(),
3387
3388
3389 // Before sending & receiving, make sure we protect this section with
3390 // a mutex:
3393 mutex, tria->get_mpi_communicator());
3394
3395 const int mpi_tag =
3397 const int mpi_tag_reply =
3399
3400 // send our requests
3401 std::vector<MPI_Request> requests(ghost_owners.size());
3402 {
3403 unsigned int idx = 0;
3404 for (const auto &it : neighbor_cell_list)
3405 {
3406 // send the data about the relevant cells
3407 const int ierr = MPI_Isend(it.second.data(),
3408 it.second.size() * sizeof(it.second[0]),
3409 MPI_BYTE,
3410 it.first,
3411 mpi_tag,
3412 tria->get_mpi_communicator(),
3413 &requests[idx]);
3414 AssertThrowMPI(ierr);
3415 ++idx;
3416 }
3417 }
3418
3419 // receive requests and reply with the results
3420 std::vector<MPI_Request> reply_requests(ghost_owners.size());
3421 std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
3422
3423 for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
3424 {
3425 MPI_Status status;
3426 int ierr = MPI_Probe(MPI_ANY_SOURCE,
3427 mpi_tag,
3428 tria->get_mpi_communicator(),
3429 &status);
3430 AssertThrowMPI(ierr);
3431
3432 int len;
3433 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3434 AssertThrowMPI(ierr);
3435 Assert(len % sizeof(typename CellId::binary_type) == 0,
3437
3438 const unsigned int n_cells =
3439 len / sizeof(typename CellId::binary_type);
3440 std::vector<typename CellId::binary_type> cells_with_requests(
3441 n_cells);
3442 std::vector<DataType> data_to_send;
3443 data_to_send.reserve(n_cells);
3444 std::vector<bool> cell_carries_data(n_cells, false);
3445
3446 ierr = MPI_Recv(cells_with_requests.data(),
3447 len,
3448 MPI_BYTE,
3449 status.MPI_SOURCE,
3450 status.MPI_TAG,
3451 tria->get_mpi_communicator(),
3452 &status);
3453 AssertThrowMPI(ierr);
3454
3455 // store data for each cell
3456 for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
3457 {
3458 const auto cell =
3459 tria->create_cell_iterator(CellId(cells_with_requests[c]));
3460
3461 MeshCellIteratorType mesh_it(tria,
3462 cell->level(),
3463 cell->index(),
3464 &mesh);
3465
3466 const std::optional<DataType> data = pack(mesh_it);
3467 if (data)
3468 {
3469 data_to_send.emplace_back(*data);
3470 cell_carries_data[c] = true;
3471 }
3472 }
3473
3474 // collect data for sending the reply in a buffer
3475
3476 // (a) make room for storing the local offsets in case we receive
3477 // other data
3478 sendbuffers[idx].resize(sizeof(std::size_t));
3479
3480 // (b) append the actual data and store how much memory it
3481 // corresponds to, which we then insert into the leading position of
3482 // the sendbuffer
3483 std::size_t size_of_send =
3484 Utilities::pack(data_to_send,
3485 sendbuffers[idx],
3486 /*enable_compression*/ false);
3487 std::memcpy(sendbuffers[idx].data(),
3488 &size_of_send,
3489 sizeof(std::size_t));
3490
3491 // (c) append information of certain cells that got left out in case
3492 // we need it
3493 if (data_to_send.size() < n_cells)
3494 Utilities::pack(cell_carries_data,
3495 sendbuffers[idx],
3496 /*enable_compression*/ false);
3497
3498 // send data
3499 ierr = MPI_Isend(sendbuffers[idx].data(),
3500 sendbuffers[idx].size(),
3501 MPI_BYTE,
3502 status.MPI_SOURCE,
3503 mpi_tag_reply,
3504 tria->get_mpi_communicator(),
3505 &reply_requests[idx]);
3506 AssertThrowMPI(ierr);
3507 }
3508
3509 // finally receive the replies
3510 std::vector<char> receive;
3511 for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
3512 {
3513 MPI_Status status;
3514 int ierr = MPI_Probe(MPI_ANY_SOURCE,
3515 mpi_tag_reply,
3516 tria->get_mpi_communicator(),
3517 &status);
3518 AssertThrowMPI(ierr);
3519
3520 int len;
3521 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3522 AssertThrowMPI(ierr);
3523
3524 receive.resize(len);
3525
3526 ierr = MPI_Recv(receive.data(),
3527 len,
3528 MPI_BYTE,
3529 status.MPI_SOURCE,
3530 status.MPI_TAG,
3531 tria->get_mpi_communicator(),
3532 &status);
3533 AssertThrowMPI(ierr);
3534
3535 // (a) first determine the length of the data section in the
3536 // received buffer
3537 auto data_iterator = receive.begin();
3538 std::size_t size_of_received_data =
3539 Utilities::unpack<std::size_t>(data_iterator,
3540 data_iterator + sizeof(std::size_t));
3541 data_iterator += sizeof(std::size_t);
3542
3543 // (b) unpack the data section in the indicated region
3544 auto received_data = Utilities::unpack<std::vector<DataType>>(
3545 data_iterator,
3546 data_iterator + size_of_received_data,
3547 /*enable_compression*/ false);
3548 data_iterator += size_of_received_data;
3549
3550 // (c) check if the received data contained fewer entries than the
3551 // number of cells we identified in the beginning, in which case we
3552 // need to extract the boolean vector with the relevant information
3553 const std::vector<typename CellId::binary_type> &this_cell_list =
3554 neighbor_cell_list[status.MPI_SOURCE];
3555 AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
3556 std::vector<bool> cells_with_data;
3557 if (received_data.size() < this_cell_list.size())
3558 {
3559 cells_with_data = Utilities::unpack<std::vector<bool>>(
3560 data_iterator, receive.end(), /*enable_compression*/ false);
3561 AssertDimension(cells_with_data.size(), this_cell_list.size());
3562 }
3563
3564 // (d) go through the received data and call the user-provided
3565 // unpack function
3566 auto received_data_iterator = received_data.begin();
3567 for (unsigned int c = 0; c < this_cell_list.size(); ++c)
3568 if (cells_with_data.empty() || cells_with_data[c])
3569 {
3571 tria_cell = tria->create_cell_iterator(this_cell_list[c]);
3572
3573 MeshCellIteratorType cell(tria,
3574 tria_cell->level(),
3575 tria_cell->index(),
3576 &mesh);
3577
3578 unpack(cell, *received_data_iterator);
3579 ++received_data_iterator;
3580 }
3581 }
3582
3583 // make sure that all communication is finished
3584 // when we leave this function.
3585 if (requests.size() > 0)
3586 {
3587 const int ierr =
3588 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
3589 AssertThrowMPI(ierr);
3590 }
3591 if (reply_requests.size() > 0)
3592 {
3593 const int ierr = MPI_Waitall(reply_requests.size(),
3594 reply_requests.data(),
3595 MPI_STATUSES_IGNORE);
3596 AssertThrowMPI(ierr);
3597 }
3598
3599
3600# endif // DEAL_II_WITH_MPI
3601 }
3602
3603 } // namespace internal
3604
3605 template <typename DataType, typename MeshType>
3607 inline void exchange_cell_data_to_ghosts(
3608 const MeshType &mesh,
3609 const std::function<std::optional<DataType>(
3610 const typename MeshType::active_cell_iterator &)> &pack,
3611 const std::function<void(const typename MeshType::active_cell_iterator &,
3612 const DataType &)> &unpack,
3613 const std::function<bool(const typename MeshType::active_cell_iterator &)>
3614 &cell_filter)
3615 {
3616# ifndef DEAL_II_WITH_MPI
3617 (void)mesh;
3618 (void)pack;
3619 (void)unpack;
3620 (void)cell_filter;
3621 Assert(false, ExcNeedsMPI());
3622# else
3623 internal::exchange_cell_data<DataType,
3624 MeshType,
3625 typename MeshType::active_cell_iterator>(
3626 mesh,
3627 pack,
3628 unpack,
3629 cell_filter,
3630 [&](const auto &process) {
3631 for (const auto &cell : mesh.active_cell_iterators())
3632 if (cell->is_ghost())
3633 process(cell, cell->subdomain_id());
3634 },
3635 [](const auto &tria) { return tria.ghost_owners(); });
3636# endif
3637 }
3638
3639
3640
3641 template <typename DataType, typename MeshType>
3644 const MeshType &mesh,
3645 const std::function<std::optional<DataType>(
3646 const typename MeshType::level_cell_iterator &)> &pack,
3647 const std::function<void(const typename MeshType::level_cell_iterator &,
3648 const DataType &)> &unpack,
3649 const std::function<bool(const typename MeshType::level_cell_iterator &)>
3650 &cell_filter)
3651 {
3652# ifndef DEAL_II_WITH_MPI
3653 (void)mesh;
3654 (void)pack;
3655 (void)unpack;
3656 (void)cell_filter;
3657 Assert(false, ExcNeedsMPI());
3658# else
3659 internal::exchange_cell_data<DataType,
3660 MeshType,
3661 typename MeshType::level_cell_iterator>(
3662 mesh,
3663 pack,
3664 unpack,
3665 cell_filter,
3666 [&](const auto &process) {
3667 for (const auto &cell : mesh.cell_iterators())
3668 if (cell->is_ghost_on_level())
3669 process(cell, cell->level_subdomain_id());
3670 },
3671 [](const auto &tria) { return tria.level_ghost_owners(); });
3672# endif
3673 }
3674} // namespace GridTools
3675
3676#endif // DOXYGEN
3677
3679
3680#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:283
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:245
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
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:746
std::size_t size
Definition mpi.cc:745
const MPI_Comm comm
Definition mpi.cc:924
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())
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 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_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 >())
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:14902
constexpr types::manifold_id flat_manifold_id
Definition types.h:336
STL namespace.
Definition types.h:32
unsigned int global_dof_index
Definition types.h:90
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