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