Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2001 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_grid_tools_h
17#define dealii_grid_tools_h
18
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/base/bounding_box.h>
23#include <deal.II/base/geometry_info.h>
24#include <deal.II/base/point.h>
25#include <deal.II/base/std_cxx17/optional.h>
26
27#include <deal.II/boost_adaptors/bounding_box.h>
28
29#include <deal.II/distributed/shared_tria.h>
30
31#include <deal.II/dofs/dof_handler.h>
32
33#include <deal.II/fe/fe_values.h>
34#include <deal.II/fe/mapping.h>
35
36#include <deal.II/grid/manifold.h>
37#include <deal.II/grid/tria.h>
38#include <deal.II/grid/tria_accessor.h>
39#include <deal.II/grid/tria_iterator.h>
40
41#include <deal.II/lac/la_parallel_vector.h>
42#include <deal.II/lac/la_vector.h>
43#include <deal.II/lac/petsc_vector.h>
44#include <deal.II/lac/sparsity_tools.h>
45#include <deal.II/lac/trilinos_vector.h>
46
47#include <deal.II/numerics/rtree.h>
48
49#include <boost/archive/binary_iarchive.hpp>
50#include <boost/archive/binary_oarchive.hpp>
51#include <boost/random/mersenne_twister.hpp>
52#include <boost/serialization/array.hpp>
53#include <boost/serialization/vector.hpp>
54
55#ifdef DEAL_II_WITH_ZLIB
56# include <boost/iostreams/device/back_inserter.hpp>
57# include <boost/iostreams/filter/gzip.hpp>
58# include <boost/iostreams/filtering_stream.hpp>
59# include <boost/iostreams/stream.hpp>
60#endif
61
62#include <bitset>
63#include <list>
64#include <set>
65
66#ifdef DEAL_II_HAVE_CXX20
67# include <concepts>
68#endif
69
70
71DEAL_II_NAMESPACE_OPEN
72
73// Forward declarations
74#ifndef DOXYGEN
75namespace parallel
76{
77 namespace distributed
78 {
79 template <int dim, int spacedim>
80 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
81 class Triangulation;
82 }
83} // namespace parallel
84
85namespace hp
86{
87 template <int, int>
88 class MappingCollection;
89}
90
91class SparsityPattern;
92
93namespace GridTools
94{
95 template <int dim, int spacedim>
96 class Cache;
97}
98#endif
99
100namespace internal
101{
102 template <int dim, int spacedim, class MeshType>
103 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
104 class ActiveCellIterator
105 {
106 public:
107#ifndef _MSC_VER
108 using type = typename MeshType::active_cell_iterator;
109#else
110 using type = TriaActiveIterator<::CellAccessor<dim, spacedim>>;
111#endif
112 };
113
114#ifdef _MSC_VER
115 template <int dim, int spacedim>
116 class ActiveCellIterator<dim, spacedim, DoFHandler<dim, spacedim>>
117 {
118 public:
119 using type =
120 TriaActiveIterator<::DoFCellAccessor<dim, spacedim, false>>;
121 };
122#endif
123} // namespace internal
124
133namespace GridTools
134{
146 template <int dim, int spacedim>
147 double
148 diameter(const Triangulation<dim, spacedim> &tria);
149
173 template <int dim, int spacedim>
174 double
175 volume(const Triangulation<dim, spacedim> &tria);
176
204 template <int dim, int spacedim>
205 double
206 volume(const Triangulation<dim, spacedim> &tria,
207 const Mapping<dim, spacedim> & mapping);
208
219 template <int dim, int spacedim>
220 double
221 minimal_cell_diameter(
222 const Triangulation<dim, spacedim> &triangulation,
223 const Mapping<dim, spacedim> & mapping =
224 (ReferenceCells::get_hypercube<dim>()
225#ifndef _MSC_VER
226 .template get_default_linear_mapping<dim, spacedim>()
227#else
228 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
229#endif
230 ));
231
242 template <int dim, int spacedim>
243 double
244 maximal_cell_diameter(
245 const Triangulation<dim, spacedim> &triangulation,
246 const Mapping<dim, spacedim> & mapping =
247 (ReferenceCells::get_hypercube<dim>()
248#ifndef _MSC_VER
249 .template get_default_linear_mapping<dim, spacedim>()
250#else
251 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
252#endif
253 ));
254
266 template <int dim>
267 DEAL_II_DEPRECATED double
268 cell_measure(
269 const std::vector<Point<dim>> &all_vertices,
270 const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]);
271
291 template <int dim>
292 double
293 cell_measure(const std::vector<Point<dim>> & all_vertices,
294 const ArrayView<const unsigned int> &vertex_indices);
295
318 template <int dim, int spacedim>
319 std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
320 affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices);
321
351 template <int dim>
352 Vector<double>
353 compute_aspect_ratio_of_cells(const Mapping<dim> & mapping,
354 const Triangulation<dim> &triangulation,
355 const Quadrature<dim> & quadrature);
356
364 template <int dim>
365 double
366 compute_maximum_aspect_ratio(const Mapping<dim> & mapping,
367 const Triangulation<dim> &triangulation,
368 const Quadrature<dim> & quadrature);
369
383 template <int dim, int spacedim>
384 BoundingBox<spacedim>
385 compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);
386
404 template <typename Iterator>
405 Point<Iterator::AccessorType::space_dimension>
406 project_to_object(
407 const Iterator & object,
408 const Point<Iterator::AccessorType::space_dimension> &trial_point);
409
422 template <int dim, int spacedim>
423 std::
424 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
425 get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
426
449 template <int dim, int spacedim>
450 void
451 delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
452 std::vector<CellData<dim>> & cells,
453 SubCellData & subcelldata);
454
473 template <int dim, int spacedim>
474 void
475 delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
476 std::vector<CellData<dim>> & cells,
477 SubCellData & subcelldata,
478 std::vector<unsigned int> & considered_vertices,
479 const double tol = 1e-12);
480
488 template <int dim>
489 void
490 delete_duplicated_vertices(std::vector<Point<dim>> &vertices,
491 const double tol = 1e-12);
492
511 template <int dim, int spacedim>
512 void
513 invert_all_negative_measure_cells(
514 const std::vector<Point<spacedim>> &all_vertices,
515 std::vector<CellData<dim>> & cells);
516
526 template <int dim, int spacedim>
527 std::size_t
528 invert_cells_with_negative_measure(
529 const std::vector<Point<spacedim>> &all_vertices,
530 std::vector<CellData<dim>> & cells);
531
542 template <int dim>
543 void
544 consistently_order_cells(std::vector<CellData<dim>> &cells);
545
633 template <int dim, typename Transformation, int spacedim>
634 DEAL_II_CXX20_REQUIRES(
635 (std::invocable<Transformation, Point<spacedim>> &&
636 std::assignable_from<
637 Point<spacedim> &,
638 std::invoke_result_t<Transformation, Point<spacedim>>>))
639 void transform(const Transformation & transformation,
640 Triangulation<dim, spacedim> &triangulation);
641
648 template <int dim, int spacedim>
649 void
650 shift(const Tensor<1, spacedim> & shift_vector,
651 Triangulation<dim, spacedim> &triangulation);
652
653
664 template <int dim, int spacedim>
665 void
666 rotate(const double angle, Triangulation<dim, spacedim> &triangulation);
667
680 template <int dim>
681 void
682 rotate(const Tensor<1, 3, double> &axis,
683 const double angle,
684 Triangulation<dim, 3> & triangulation);
685
701 template <int dim>
702 DEAL_II_DEPRECATED void
703 rotate(const double angle,
704 const unsigned int axis,
705 Triangulation<dim, 3> &triangulation);
706
764 template <int dim>
765 void
766 laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
767 Triangulation<dim> & tria,
768 const Function<dim, double> *coefficient = nullptr,
769 const bool solve_for_absolute_positions = false);
770
776 template <int dim, int spacedim>
777 std::map<unsigned int, Point<spacedim>>
778 get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria);
779
787 template <int dim, int spacedim>
788 void
789 scale(const double scaling_factor,
790 Triangulation<dim, spacedim> &triangulation);
791
811 template <int dim, int spacedim>
812 void
813 distort_random(
814 const double factor,
815 Triangulation<dim, spacedim> &triangulation,
816 const bool keep_boundary = true,
817 const unsigned int seed = boost::random::mt19937::default_seed);
818
852 template <int dim, int spacedim>
853 void
854 remove_hanging_nodes(Triangulation<dim, spacedim> &tria,
855 const bool isotropic = false,
856 const unsigned int max_iterations = 100);
857
882 template <int dim, int spacedim>
883 void
884 remove_anisotropy(Triangulation<dim, spacedim> &tria,
885 const double max_ratio = 1.6180339887,
886 const unsigned int max_iterations = 5);
887
977 template <int dim, int spacedim>
978 void
979 regularize_corner_cells(Triangulation<dim, spacedim> &tria,
980 const double limit_angle_fraction = .75);
981
1045 template <int dim, int spacedim>
1046#ifndef DOXYGEN
1047 std::tuple<
1048 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1049 std::vector<std::vector<Point<dim>>>,
1050 std::vector<std::vector<unsigned int>>>
1051#else
1052 return_type
1053#endif
1054 compute_point_locations(
1055 const Cache<dim, spacedim> & cache,
1056 const std::vector<Point<spacedim>> &points,
1057 const typename Triangulation<dim, spacedim>::active_cell_iterator
1058 &cell_hint =
1059 typename Triangulation<dim, spacedim>::active_cell_iterator());
1060
1094 template <int dim, int spacedim>
1095#ifndef DOXYGEN
1096 std::tuple<
1097 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1098 std::vector<std::vector<Point<dim>>>,
1099 std::vector<std::vector<unsigned int>>,
1100 std::vector<unsigned int>>
1101#else
1102 return_type
1103#endif
1104 compute_point_locations_try_all(
1105 const Cache<dim, spacedim> & cache,
1106 const std::vector<Point<spacedim>> &points,
1107 const typename Triangulation<dim, spacedim>::active_cell_iterator
1108 &cell_hint =
1109 typename Triangulation<dim, spacedim>::active_cell_iterator());
1110
1190 template <int dim, int spacedim>
1191#ifndef DOXYGEN
1192 std::tuple<
1193 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1194 std::vector<std::vector<Point<dim>>>,
1195 std::vector<std::vector<unsigned int>>,
1196 std::vector<std::vector<Point<spacedim>>>,
1197 std::vector<std::vector<unsigned int>>>
1198#else
1199 return_type
1200#endif
1201 distributed_compute_point_locations(
1202 const GridTools::Cache<dim, spacedim> & cache,
1203 const std::vector<Point<spacedim>> & local_points,
1204 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1205 const double tolerance = 1e-10,
1206 const std::vector<bool> & marked_vertices = {},
1207 const bool enforce_unique_mapping = true);
1208
1209 namespace internal
1210 {
1225 template <int dim, int spacedim>
1226 struct DistributedComputePointLocationsInternal
1227 {
1228 DistributedComputePointLocationsInternal();
1229
1236 void
1237 finalize_setup();
1238
1242 unsigned int n_searched_points;
1243
1250 std::vector<std::tuple<std::pair<int, int>,
1251 unsigned int,
1252 unsigned int,
1253 Point<dim>,
1254 Point<spacedim>,
1255 unsigned int>>
1256 send_components;
1257
1261 std::vector<unsigned int> send_ranks;
1262
1268 std::vector<unsigned int> send_ptrs;
1269
1280 std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1281 recv_components;
1282
1286 std::vector<unsigned int> recv_ranks;
1287
1293 std::vector<unsigned int> recv_ptrs;
1294 };
1295
1305 template <int dim, int spacedim>
1306 DistributedComputePointLocationsInternal<dim, spacedim>
1307 distributed_compute_point_locations(
1308 const GridTools::Cache<dim, spacedim> & cache,
1309 const std::vector<Point<spacedim>> & points,
1310 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1311 const std::vector<bool> & marked_vertices,
1312 const double tolerance,
1313 const bool perform_handshake,
1314 const bool enforce_unique_mapping = false);
1315
1316
1324 template <int structdim, int spacedim>
1325 struct DistributedComputeIntersectionLocationsInternal
1326 {
1330 using IntersectionType =
1331 std::array<::Point<spacedim>, structdim + 1>;
1332
1341 std::vector<std::tuple<std::pair<int, int>,
1342 unsigned int,
1343 unsigned int,
1344 IntersectionType>>
1345 send_components;
1346
1355 std::vector<std::tuple<unsigned int, unsigned int, IntersectionType>>
1356 recv_components;
1357
1361 std::vector<unsigned int> recv_ptrs;
1362
1378 template <int dim>
1379 GridTools::internal::DistributedComputePointLocationsInternal<dim,
1380 spacedim>
1381 convert_to_distributed_compute_point_locations_internal(
1382 const unsigned int n_points_1D,
1383 const Triangulation<dim, spacedim> &tria,
1384 const Mapping<dim, spacedim> & mapping,
1385 const bool consistent_numbering_of_sender_and_receiver = false) const;
1386
1387 private:
1395 std::map<unsigned int, std::vector<unsigned int>>
1396 communicate_indices(
1397 const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1398 & point_recv_components,
1399 const MPI_Comm comm) const;
1400 };
1401
1409 template <int structdim, int dim, int spacedim>
1410 DistributedComputeIntersectionLocationsInternal<structdim, spacedim>
1411 distributed_compute_intersection_locations(
1412 const Cache<dim, spacedim> & cache,
1413 const std::vector<std::vector<Point<spacedim>>> &intersection_requests,
1414 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1415 const std::vector<bool> & marked_vertices,
1416 const double tolerance);
1417
1418 } // namespace internal
1419
1456 template <int dim, int spacedim>
1457 std::map<unsigned int, Point<spacedim>>
1458 extract_used_vertices(
1459 const Triangulation<dim, spacedim> &container,
1460 const Mapping<dim, spacedim> & mapping =
1461 (ReferenceCells::get_hypercube<dim>()
1462#ifndef _MSC_VER
1463 .template get_default_linear_mapping<dim, spacedim>()
1464#else
1465 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1466#endif
1467 ));
1468
1478 template <int spacedim>
1479 unsigned int
1480 find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1481 const Point<spacedim> & p);
1482
1509 template <int dim, template <int, int> class MeshType, int spacedim>
1510 DEAL_II_CXX20_REQUIRES(
1511 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1512 unsigned int find_closest_vertex(
1513 const MeshType<dim, spacedim> &mesh,
1514 const Point<spacedim> & p,
1515 const std::vector<bool> & marked_vertices = {});
1516
1543 template <int dim, template <int, int> class MeshType, int spacedim>
1544 DEAL_II_CXX20_REQUIRES(
1545 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1546 unsigned int find_closest_vertex(
1547 const Mapping<dim, spacedim> & mapping,
1548 const MeshType<dim, spacedim> &mesh,
1549 const Point<spacedim> & p,
1550 const std::vector<bool> & marked_vertices = {});
1551
1552
1575 template <class MeshType>
1576 DEAL_II_CXX20_REQUIRES((concepts::is_triangulation_or_dof_handler<MeshType>))
1577#ifndef _MSC_VER
1578 std::vector<typename MeshType::active_cell_iterator>
1579#else
1580 std::vector<
1581 typename ::internal::ActiveCellIterator<MeshType::dimension,
1582 MeshType::space_dimension,
1583 MeshType>::type>
1584#endif
1585 find_cells_adjacent_to_vertex(const MeshType & container,
1586 const unsigned int vertex_index);
1587
1653 template <int dim, template <int, int> class MeshType, int spacedim>
1654 DEAL_II_CXX20_REQUIRES(
1655 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1656#ifndef _MSC_VER
1657 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1658#else
1659 std::pair<typename ::internal::
1660 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1661 Point<dim>>
1662#endif
1663 find_active_cell_around_point(const Mapping<dim, spacedim> & mapping,
1664 const MeshType<dim, spacedim> &mesh,
1665 const Point<spacedim> & p,
1666 const std::vector<bool> &marked_vertices = {},
1667 const double tolerance = 1.e-10);
1668
1679 template <int dim, template <int, int> class MeshType, int spacedim>
1680 DEAL_II_CXX20_REQUIRES(
1681 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1682#ifndef _MSC_VER
1683 typename MeshType<dim, spacedim>::active_cell_iterator
1684#else
1685 typename ::internal::
1686 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1687#endif
1688 find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1689 const Point<spacedim> & p,
1690 const std::vector<bool> &marked_vertices = {},
1691 const double tolerance = 1.e-10);
1692
1699 template <int dim, int spacedim>
1700 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1701 Point<dim>>
1702 find_active_cell_around_point(
1703 const hp::MappingCollection<dim, spacedim> &mapping,
1704 const DoFHandler<dim, spacedim> & mesh,
1705 const Point<spacedim> & p,
1706 const double tolerance = 1.e-10);
1707
1759 template <int dim, int spacedim>
1760 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1761 Point<dim>>
1762 find_active_cell_around_point(
1763 const Cache<dim, spacedim> &cache,
1764 const Point<spacedim> & p,
1765 const typename Triangulation<dim, spacedim>::active_cell_iterator &
1766 cell_hint = typename Triangulation<dim, spacedim>::active_cell_iterator(),
1767 const std::vector<bool> &marked_vertices = {},
1768 const double tolerance = 1.e-10);
1769
1786 template <int dim, template <int, int> class MeshType, int spacedim>
1787 DEAL_II_CXX20_REQUIRES(
1788 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1789#ifndef _MSC_VER
1790 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1791#else
1792 std::pair<typename ::internal::
1793 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1794 Point<dim>>
1795#endif
1796 find_active_cell_around_point(
1797 const Mapping<dim, spacedim> & mapping,
1798 const MeshType<dim, spacedim> &mesh,
1799 const Point<spacedim> & p,
1800 const std::vector<
1801 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1802 &vertex_to_cell_map,
1803 const std::vector<std::vector<Tensor<1, spacedim>>>
1804 &vertex_to_cell_centers,
1805 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1806 typename MeshType<dim, spacedim>::active_cell_iterator(),
1807 const std::vector<bool> &marked_vertices = {},
1808 const RTree<std::pair<Point<spacedim>, unsigned int>> &
1809 used_vertices_rtree = RTree<std::pair<Point<spacedim>, unsigned int>>{},
1810 const double tolerance = 1.e-10,
1811 const RTree<
1812 std::pair<BoundingBox<spacedim>,
1813 typename Triangulation<dim, spacedim>::active_cell_iterator>>
1814 *relevant_cell_bounding_boxes_rtree = nullptr);
1815
1840 template <int dim, template <int, int> class MeshType, int spacedim>
1841 DEAL_II_CXX20_REQUIRES(
1842 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1843#ifndef _MSC_VER
1844 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1845 Point<dim>>>
1846#else
1847 std::vector<std::pair<
1848 typename ::internal::
1849 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1850 Point<dim>>>
1851#endif
1852 find_all_active_cells_around_point(
1853 const Mapping<dim, spacedim> & mapping,
1854 const MeshType<dim, spacedim> &mesh,
1855 const Point<spacedim> & p,
1856 const double tolerance,
1857 const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1858 Point<dim>> & first_cell);
1859
1869 template <int dim, template <int, int> class MeshType, int spacedim>
1870 DEAL_II_CXX20_REQUIRES(
1871 (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
1872#ifndef _MSC_VER
1873 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1874 Point<dim>>>
1875#else
1876 std::vector<std::pair<
1877 typename ::internal::
1878 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1879 Point<dim>>>
1880#endif
1881 find_all_active_cells_around_point(
1882 const Mapping<dim, spacedim> & mapping,
1883 const MeshType<dim, spacedim> &mesh,
1884 const Point<spacedim> & p,
1885 const double tolerance = 1e-10,
1886 const std::vector<bool> & marked_vertices = {});
1887
1911 template <class MeshType>
1912 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1913 std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
1914 const typename MeshType::cell_iterator &cell);
1915
1942 template <class MeshType>
1943 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
1944 void get_active_neighbors(
1945 const typename MeshType::active_cell_iterator & cell,
1946 std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1947
1999 template <class MeshType>
2000 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2001 std::
2002 vector<typename MeshType::active_cell_iterator> compute_active_cell_halo_layer(
2003 const MeshType &mesh,
2004 const std::function<bool(const typename MeshType::active_cell_iterator &)>
2005 &predicate);
2006
2007
2017 template <class MeshType>
2018 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2019 std::
2020 vector<typename MeshType::cell_iterator> compute_cell_halo_layer_on_level(
2021 const MeshType &mesh,
2022 const std::function<bool(const typename MeshType::cell_iterator &)>
2023 & predicate,
2024 const unsigned int level);
2025
2026
2041 template <class MeshType>
2042 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2043 std::vector<
2044 typename MeshType::
2045 active_cell_iterator> compute_ghost_cell_halo_layer(const MeshType &mesh);
2046
2097 template <class MeshType>
2098 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2099 std::
2100 vector<typename MeshType::active_cell_iterator> compute_active_cell_layer_within_distance(
2101 const MeshType &mesh,
2102 const std::function<bool(const typename MeshType::active_cell_iterator &)>
2103 & predicate,
2104 const double layer_thickness);
2105
2130 template <class MeshType>
2131 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2132 std::
2133 vector<typename MeshType::active_cell_iterator> compute_ghost_cell_layer_within_distance(
2134 const MeshType &mesh,
2135 const double layer_thickness);
2136
2154 template <class MeshType>
2155 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2156 std::pair<
2157 Point<MeshType::space_dimension>,
2158 Point<MeshType::
2159 space_dimension>> compute_bounding_box(const MeshType &mesh,
2160 const std::function<bool(
2161 const typename MeshType::
2162 active_cell_iterator &)>
2163 &predicate);
2164
2237 template <class MeshType>
2238 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2239 std::
2240 vector<BoundingBox<MeshType::space_dimension>> compute_mesh_predicate_bounding_box(
2241 const MeshType &mesh,
2242 const std::function<bool(const typename MeshType::active_cell_iterator &)>
2243 & predicate,
2244 const unsigned int refinement_level = 0,
2245 const bool allow_merge = false,
2246 const unsigned int max_boxes = numbers::invalid_unsigned_int);
2247
2275 template <int spacedim>
2276#ifndef DOXYGEN
2277 std::tuple<std::vector<std::vector<unsigned int>>,
2278 std::map<unsigned int, unsigned int>,
2279 std::map<unsigned int, std::vector<unsigned int>>>
2280#else
2281 return_type
2282#endif
2283 guess_point_owner(
2284 const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
2285 const std::vector<Point<spacedim>> & points);
2286
2287
2322 template <int spacedim>
2323#ifndef DOXYGEN
2324 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2325 std::map<unsigned int, unsigned int>,
2326 std::map<unsigned int, std::vector<unsigned int>>>
2327#else
2328 return_type
2329#endif
2330 guess_point_owner(
2331 const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2332 const std::vector<Point<spacedim>> & points);
2333
2334
2343 template <int dim, int spacedim>
2344 std::vector<
2345 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2346 vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
2347
2360 template <int dim, int spacedim>
2361 std::vector<std::vector<Tensor<1, spacedim>>>
2362 vertex_to_cell_centers_directions(
2363 const Triangulation<dim, spacedim> &mesh,
2364 const std::vector<
2365 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2366 &vertex_to_cells);
2367
2368
2376 template <int dim, int spacedim>
2377 unsigned int
2378 find_closest_vertex_of_cell(
2379 const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
2380 const Point<spacedim> & position,
2381 const Mapping<dim, spacedim> & mapping =
2382 (ReferenceCells::get_hypercube<dim>()
2383#ifndef _MSC_VER
2384 .template get_default_linear_mapping<dim, spacedim>()
2385#else
2386 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2387#endif
2388 ));
2389
2401 template <int dim, int spacedim>
2402 std::map<unsigned int, types::global_vertex_index>
2403 compute_local_to_global_vertex_index_map(
2404 const parallel::distributed::Triangulation<dim, spacedim> &triangulation);
2405
2417 template <int dim, int spacedim>
2418 std::pair<unsigned int, double>
2419 get_longest_direction(
2420 typename Triangulation<dim, spacedim>::active_cell_iterator cell);
2421
2436 template <int dim, int spacedim>
2437 void
2438 get_face_connectivity_of_cells(
2439 const Triangulation<dim, spacedim> &triangulation,
2440 DynamicSparsityPattern & connectivity);
2441
2450 template <int dim, int spacedim>
2451 void
2452 get_vertex_connectivity_of_cells(
2453 const Triangulation<dim, spacedim> &triangulation,
2454 DynamicSparsityPattern & connectivity);
2455
2464 template <int dim, int spacedim>
2465 void
2466 get_vertex_connectivity_of_cells_on_level(
2467 const Triangulation<dim, spacedim> &triangulation,
2468 const unsigned int level,
2469 DynamicSparsityPattern & connectivity);
2470
2491 template <int dim, int spacedim>
2492 void
2493 partition_triangulation(const unsigned int n_partitions,
2494 Triangulation<dim, spacedim> & triangulation,
2495 const SparsityTools::Partitioner partitioner =
2496 SparsityTools::Partitioner::metis);
2497
2508 template <int dim, int spacedim>
2509 void
2510 partition_triangulation(const unsigned int n_partitions,
2511 const std::vector<unsigned int> &cell_weights,
2512 Triangulation<dim, spacedim> & triangulation,
2513 const SparsityTools::Partitioner partitioner =
2514 SparsityTools::Partitioner::metis);
2515
2561 template <int dim, int spacedim>
2562 void
2563 partition_triangulation(const unsigned int n_partitions,
2564 const SparsityPattern & cell_connection_graph,
2565 Triangulation<dim, spacedim> &triangulation,
2566 const SparsityTools::Partitioner partitioner =
2567 SparsityTools::Partitioner::metis);
2568
2579 template <int dim, int spacedim>
2580 void
2581 partition_triangulation(const unsigned int n_partitions,
2582 const std::vector<unsigned int> &cell_weights,
2583 const SparsityPattern & cell_connection_graph,
2584 Triangulation<dim, spacedim> &triangulation,
2585 const SparsityTools::Partitioner partitioner =
2586 SparsityTools::Partitioner::metis);
2587
2602 template <int dim, int spacedim>
2603 void
2604 partition_triangulation_zorder(const unsigned int n_partitions,
2605 Triangulation<dim, spacedim> &triangulation,
2606 const bool group_siblings = true);
2607
2619 template <int dim, int spacedim>
2620 void
2621 partition_multigrid_levels(Triangulation<dim, spacedim> &triangulation);
2622
2630 template <int dim, int spacedim>
2631 std::vector<types::subdomain_id>
2632 get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
2633 const std::vector<CellId> & cell_ids);
2634
2645 template <int dim, int spacedim>
2646 void
2647 get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
2648 std::vector<types::subdomain_id> & subdomain);
2649
2664 template <int dim, int spacedim>
2665 unsigned int
2666 count_cells_with_subdomain_association(
2667 const Triangulation<dim, spacedim> &triangulation,
2668 const types::subdomain_id subdomain);
2669
2699 template <int dim, int spacedim>
2700 std::vector<bool>
2701 get_locally_owned_vertices(const Triangulation<dim, spacedim> &triangulation);
2702
2744 template <typename MeshType>
2745 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2746 std::list<std::pair<
2747 typename MeshType::cell_iterator,
2748 typename MeshType::cell_iterator>> get_finest_common_cells(const MeshType
2749 &mesh_1,
2750 const MeshType
2751 &mesh_2);
2752
2762 template <int dim, int spacedim>
2763 bool
2764 have_same_coarse_mesh(const Triangulation<dim, spacedim> &mesh_1,
2765 const Triangulation<dim, spacedim> &mesh_2);
2766
2778 template <typename MeshType>
2779 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2780 bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2781
2803 template <int dim, int spacedim>
2804 typename Triangulation<dim, spacedim>::DistortedCellList
2805 fix_up_distorted_child_cells(
2806 const typename Triangulation<dim, spacedim>::DistortedCellList
2807 & distorted_cells,
2808 Triangulation<dim, spacedim> &triangulation);
2809
2810
2811
2861 template <class MeshType>
2862 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
2863 std::vector<typename MeshType::active_cell_iterator> get_patch_around_cell(
2864 const typename MeshType::active_cell_iterator &cell);
2865
2866
2888 template <class Container>
2889 std::vector<typename Container::cell_iterator>
2890 get_cells_at_coarsest_common_level(
2891 const std::vector<typename Container::active_cell_iterator> &patch_cells);
2892
2959 template <class Container>
2960 void
2961 build_triangulation_from_patch(
2962 const std::vector<typename Container::active_cell_iterator> &patch,
2963 Triangulation<Container::dimension, Container::space_dimension>
2964 &local_triangulation,
2965 std::map<
2966 typename Triangulation<Container::dimension,
2967 Container::space_dimension>::active_cell_iterator,
2968 typename Container::active_cell_iterator> &patch_to_global_tria_map);
2969
3001 template <int dim, int spacedim>
3002 std::map<
3003 types::global_dof_index,
3004 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
3005 get_dof_to_support_patch_map(DoFHandler<dim, spacedim> &dof_handler);
3006
3007
3020 template <typename CellIterator>
3021 struct PeriodicFacePair
3022 {
3026 CellIterator cell[2];
3027
3032 unsigned int face_idx[2];
3033
3039 std::bitset<3> orientation;
3040
3053 FullMatrix<double> matrix;
3054
3058 std::size_t
3059 memory_consumption() const;
3060 };
3061
3062
3126 template <typename FaceIterator>
3127 bool
3128 orthogonal_equality(
3129 std::bitset<3> & orientation,
3130 const FaceIterator & face1,
3131 const FaceIterator & face2,
3132 const unsigned int direction,
3133 const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset =
3134 Tensor<1, FaceIterator::AccessorType::space_dimension>(),
3135 const FullMatrix<double> &matrix = FullMatrix<double>());
3136
3137
3141 template <typename FaceIterator>
3142 bool
3143 orthogonal_equality(
3144 const FaceIterator & face1,
3145 const FaceIterator & face2,
3146 const unsigned int direction,
3147 const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset =
3148 Tensor<1, FaceIterator::AccessorType::space_dimension>(),
3149 const FullMatrix<double> &matrix = FullMatrix<double>());
3150
3151
3210 template <typename MeshType>
3211 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3212 void collect_periodic_faces(
3213 const MeshType & mesh,
3214 const types::boundary_id b_id1,
3215 const types::boundary_id b_id2,
3216 const unsigned int direction,
3217 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
3218 & matched_pairs,
3219 const Tensor<1, MeshType::space_dimension> &offset =
3220 ::Tensor<1, MeshType::space_dimension>(),
3221 const FullMatrix<double> &matrix = FullMatrix<double>());
3222
3223
3248 template <typename MeshType>
3249 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3250 void collect_periodic_faces(
3251 const MeshType & mesh,
3252 const types::boundary_id b_id,
3253 const unsigned int direction,
3254 std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
3255 & matched_pairs,
3256 const ::Tensor<1, MeshType::space_dimension> &offset =
3257 ::Tensor<1, MeshType::space_dimension>(),
3258 const FullMatrix<double> &matrix = FullMatrix<double>());
3259
3286 template <int dim, int spacedim>
3287 void
3288 copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
3289 const bool reset_boundary_ids = false);
3290
3312 template <int dim, int spacedim>
3313 void
3314 map_boundary_to_manifold_ids(
3315 const std::vector<types::boundary_id> &src_boundary_ids,
3316 const std::vector<types::manifold_id> &dst_manifold_ids,
3317 Triangulation<dim, spacedim> & tria,
3318 const std::vector<types::boundary_id> &reset_boundary_ids = {});
3319
3349 template <int dim, int spacedim>
3350 void
3351 copy_material_to_manifold_id(Triangulation<dim, spacedim> &tria,
3352 const bool compute_face_ids = false);
3353
3378 template <int dim, int spacedim>
3379 void
3380 assign_co_dimensional_manifold_indicators(
3381 Triangulation<dim, spacedim> & tria,
3382 const std::function<types::manifold_id(
3383 const std::set<types::manifold_id> &)> &disambiguation_function =
3384 [](const std::set<types::manifold_id> &manifold_ids) {
3385 if (manifold_ids.size() == 1)
3386 return *manifold_ids.begin();
3387 else
3388 return numbers::flat_manifold_id;
3389 },
3390 bool overwrite_only_flat_manifold_ids = true);
3479 template <typename DataType, typename MeshType>
3480 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3481 void exchange_cell_data_to_ghosts(
3482 const MeshType & mesh,
3483 const std::function<std_cxx17::optional<DataType>(
3484 const typename MeshType::active_cell_iterator &)> &pack,
3485 const std::function<void(const typename MeshType::active_cell_iterator &,
3486 const DataType &)> & unpack,
3487 const std::function<bool(const typename MeshType::active_cell_iterator &)>
3488 &cell_filter =
3489 always_return<typename MeshType::active_cell_iterator, bool>{true});
3490
3503 template <typename DataType, typename MeshType>
3504 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
3505 void exchange_cell_data_to_level_ghosts(
3506 const MeshType & mesh,
3507 const std::function<std_cxx17::optional<DataType>(
3508 const typename MeshType::level_cell_iterator &)> &pack,
3509 const std::function<void(const typename MeshType::level_cell_iterator &,
3510 const DataType &)> & unpack,
3511 const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3512 cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
3513 true});
3514
3515 /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3516 * boxes @p local_bboxes.
3517 *
3518 * This function is meant to exchange bounding boxes describing the locally
3519 * owned cells in a distributed triangulation obtained with the function
3520 * GridTools::compute_mesh_predicate_bounding_box .
3521 *
3522 * The output vector's size is the number of processes of the MPI
3523 * communicator:
3524 * its i-th entry contains the vector @p local_bboxes of the i-th process.
3525 */
3526 template <int spacedim>
3527 std::vector<std::vector<BoundingBox<spacedim>>>
3528 exchange_local_bounding_boxes(
3529 const std::vector<BoundingBox<spacedim>> &local_bboxes,
3530 const MPI_Comm mpi_communicator);
3531
3564 template <int spacedim>
3565 RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
3566 build_global_description_tree(
3567 const std::vector<BoundingBox<spacedim>> &local_description,
3568 const MPI_Comm mpi_communicator);
3569
3587 template <int dim, int spacedim>
3588 void
3589 collect_coinciding_vertices(
3590 const Triangulation<dim, spacedim> & tria,
3591 std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3592 std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3593
3613 template <int dim, int spacedim>
3614 std::map<unsigned int, std::set<::types::subdomain_id>>
3615 compute_vertices_with_ghost_neighbors(
3616 const Triangulation<dim, spacedim> &tria);
3617
3633 template <int dim, typename T>
3634 struct DEAL_II_DEPRECATED CellDataTransferBuffer
3635 {
3639 std::vector<CellId> cell_ids;
3640
3644 std::vector<T> data;
3645
3654 template <class Archive>
3655 void
3656 save(Archive &ar, const unsigned int) const
3657 {
3658 // Implement the code inline as some compilers do otherwise complain
3659 // about the use of a deprecated interface.
3660 Assert(cell_ids.size() == data.size(),
3661 ExcDimensionMismatch(cell_ids.size(), data.size()));
3662 // archive the cellids in an efficient binary format
3663 const std::size_t n_cells = cell_ids.size();
3664 ar & n_cells;
3665 for (const auto &id : cell_ids)
3666 {
3667 CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3668 ar & binary_cell_id;
3669 }
3670
3671 ar &data;
3672 }
3673
3680 template <class Archive>
3681 void
3682 load(Archive &ar, const unsigned int)
3683 {
3684 // Implement the code inline as some compilers do otherwise complain
3685 // about the use of a deprecated interface.
3686 std::size_t n_cells;
3687 ar & n_cells;
3688 cell_ids.clear();
3689 cell_ids.reserve(n_cells);
3690 for (unsigned int c = 0; c < n_cells; ++c)
3691 {
3692 CellId::binary_type value;
3693 ar & value;
3694 cell_ids.emplace_back(value);
3695 }
3696 ar &data;
3697 }
3698
3699
3700#ifdef DOXYGEN
3706 template <class Archive>
3707 void
3708 serialize(Archive &archive, const unsigned int version);
3709#else
3710 // This macro defines the serialize() method that is compatible with
3711 // the templated save() and load() method that have been implemented.
3712 BOOST_SERIALIZATION_SPLIT_MEMBER()
3713#endif
3714 };
3715
3716
3717
3738 template <int dim, typename VectorType>
3739 class MarchingCubeAlgorithm
3740 {
3741 public:
3745 using value_type = typename VectorType::value_type;
3746
3750 MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3751 const FiniteElement<dim, dim> &fe,
3752 const unsigned int n_subdivisions = 1,
3753 const double tolerance = 1e-10);
3754
3765 void
3766 process(const DoFHandler<dim> & background_dof_handler,
3767 const VectorType & ls_vector,
3768 const double iso_level,
3769 std::vector<Point<dim>> &vertices,
3770 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3771
3776 void
3777 process(const DoFHandler<dim> & background_dof_handler,
3778 const VectorType & ls_vector,
3779 const double iso_level,
3780 std::vector<Point<dim>> &vertices) const;
3781
3791 void
3792 process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
3793 const VectorType & ls_vector,
3794 const double iso_level,
3795 std::vector<Point<dim>> & vertices,
3796 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
3802 void
3803 process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
3804 const VectorType & ls_vector,
3805 const double iso_level,
3806 std::vector<Point<dim>> &vertices) const;
3807
3808 private:
3813 static Quadrature<dim>
3814 create_quadrature_rule(const unsigned int n_subdivisions);
3815
3819 void
3820 process_cell(std::vector<value_type> & ls_values,
3821 const std::vector<Point<dim>> & points,
3822 const double iso_level,
3823 std::vector<Point<dim>> & vertices,
3824 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
3825 const bool write_back_cell_data = true) const;
3826
3830 void
3831 process_sub_cell(const std::vector<value_type> &,
3832 const std::vector<Point<1>> &,
3833 const std::vector<unsigned int> &,
3834 const double,
3835 std::vector<Point<1>> &,
3836 std::vector<CellData<1>> &,
3837 const bool) const
3838 {
3839 AssertThrow(false, ExcNotImplemented());
3840 }
3841
3848 void
3849 process_sub_cell(const std::vector<value_type> & ls_values,
3850 const std::vector<Point<2>> & points,
3851 const std::vector<unsigned int> &mask,
3852 const double iso_level,
3853 std::vector<Point<2>> & vertices,
3854 std::vector<CellData<1>> & cells,
3855 const bool write_back_cell_data) const;
3856
3860 void
3861 process_sub_cell(const std::vector<value_type> & ls_values,
3862 const std::vector<Point<3>> & points,
3863 const std::vector<unsigned int> &mask,
3864 const double iso_level,
3865 std::vector<Point<3>> & vertices,
3866 std::vector<CellData<2>> & cells,
3867 const bool write_back_cell_data) const;
3868
3873 const unsigned int n_subdivisions;
3874
3879 const double tolerance;
3880
3885 mutable FEValues<dim> fe_values;
3886 };
3887
3888
3889
3898 DeclException1(ExcInvalidNumberOfPartitions,
3899 int,
3900 << "The number of partitions you gave is " << arg1
3901 << ", but must be greater than zero.");
3905 DeclException1(ExcNonExistentSubdomain,
3906 int,
3907 << "The subdomain id " << arg1
3908 << " has no cells associated with it.");
3912 DeclException0(ExcTriangulationHasBeenRefined);
3913
3917 DeclException1(ExcScalingFactorNotPositive,
3918 double,
3919 << "The scaling factor must be positive, but it is " << arg1
3920 << '.');
3921
3925 DeclException1(ExcVertexNotUsed,
3926 unsigned int,
3927 << "The given vertex with index " << arg1
3928 << " is not used in the given triangulation.");
3929
3932} /*namespace GridTools*/
3933
3934
3944DeclExceptionMsg(ExcMeshNotOrientable,
3945 "The edges of the mesh are not consistently orientable.");
3946
3947
3948
3949/* ----------------- Template function --------------- */
3950
3951#ifndef DOXYGEN
3952
3953namespace GridTools
3954{
3955 template <int dim>
3956 double
3957 cell_measure(
3958 const std::vector<Point<dim>> &all_vertices,
3959 const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3960 {
3961 // We forward call to the ArrayView version:
3962 const ArrayView<const unsigned int> view(
3963 indices, GeometryInfo<dim>::vertices_per_cell);
3964 return cell_measure(all_vertices, view);
3965 }
3966
3967
3968
3969 // This specialization is defined here so that the general template in the
3970 // source file doesn't need to have further 1d overloads for the internal
3971 // functions it calls.
3972 template <>
3973 inline Triangulation<1, 1>::DistortedCellList
3974 fix_up_distorted_child_cells(const Triangulation<1, 1>::DistortedCellList &,
3975 Triangulation<1, 1> &)
3976 {
3977 return {};
3978 }
3979
3980
3981
3982 template <int dim, typename Transformation, int spacedim>
3983 DEAL_II_CXX20_REQUIRES(
3984 (std::invocable<Transformation, Point<spacedim>> &&
3985 std::assignable_from<
3986 Point<spacedim> &,
3987 std::invoke_result_t<Transformation, Point<spacedim>>>))
3988 void transform(const Transformation & transformation,
3989 Triangulation<dim, spacedim> &triangulation)
3990 {
3991 std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3992
3993 // loop over all active cells, and
3994 // transform those vertices that
3995 // have not yet been touched. note
3996 // that we get to all vertices in
3997 // the triangulation by only
3998 // visiting the active cells.
3999 typename Triangulation<dim, spacedim>::active_cell_iterator
4000 cell = triangulation.begin_active(),
4001 endc = triangulation.end();
4002 for (; cell != endc; ++cell)
4003 for (const unsigned int v : cell->vertex_indices())
4004 if (treated_vertices[cell->vertex_index(v)] == false)
4005 {
4006 // transform this vertex
4007 cell->vertex(v) = transformation(cell->vertex(v));
4008 // and mark it as treated
4009 treated_vertices[cell->vertex_index(v)] = true;
4010 };
4011
4012
4013 // now fix any vertices on hanging nodes so that we don't create any holes
4014 if (dim == 2)
4015 {
4016 typename Triangulation<dim, spacedim>::active_cell_iterator
4017 cell = triangulation.begin_active(),
4018 endc = triangulation.end();
4019 for (; cell != endc; ++cell)
4020 for (const unsigned int face : cell->face_indices())
4021 if (cell->face(face)->has_children() &&
4022 !cell->face(face)->at_boundary())
4023 {
4024 Assert(cell->reference_cell() ==
4025 ReferenceCells::get_hypercube<dim>(),
4026 ExcNotImplemented());
4027
4028 // this line has children
4029 cell->face(face)->child(0)->vertex(1) =
4030 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
4031 2;
4032 }
4033 }
4034 else if (dim == 3)
4035 {
4036 typename Triangulation<dim, spacedim>::active_cell_iterator
4037 cell = triangulation.begin_active(),
4038 endc = triangulation.end();
4039 for (; cell != endc; ++cell)
4040 for (const unsigned int face : cell->face_indices())
4041 if (cell->face(face)->has_children() &&
4042 !cell->face(face)->at_boundary())
4043 {
4044 Assert(cell->reference_cell() ==
4045 ReferenceCells::get_hypercube<dim>(),
4046 ExcNotImplemented());
4047
4048 // this face has hanging nodes
4049 cell->face(face)->child(0)->vertex(1) =
4050 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
4051 2.0;
4052 cell->face(face)->child(0)->vertex(2) =
4053 (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
4054 2.0;
4055 cell->face(face)->child(1)->vertex(3) =
4056 (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
4057 2.0;
4058 cell->face(face)->child(2)->vertex(3) =
4059 (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
4060 2.0;
4061
4062 // center of the face
4063 cell->face(face)->child(0)->vertex(3) =
4064 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
4065 cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
4066 4.0;
4067 }
4068 }
4069
4070 // Make sure FEValues notices that the mesh has changed
4071 triangulation.signals.mesh_movement();
4072 }
4073
4074
4075
4076 template <class MeshType>
4077 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
4078 std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
4079 const typename MeshType::cell_iterator &cell)
4080 {
4081 std::vector<typename MeshType::active_cell_iterator> child_cells;
4082
4083 if (cell->has_children())
4084 {
4085 for (unsigned int child = 0; child < cell->n_children(); ++child)
4086 if (cell->child(child)->has_children())
4087 {
4088 const std::vector<typename MeshType::active_cell_iterator>
4089 children = get_active_child_cells<MeshType>(cell->child(child));
4090 child_cells.insert(child_cells.end(),
4091 children.begin(),
4092 children.end());
4093 }
4094 else
4095 child_cells.push_back(cell->child(child));
4096 }
4097
4098 return child_cells;
4099 }
4100
4101
4102
4103 template <class MeshType>
4104 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
4105 void get_active_neighbors(
4106 const typename MeshType::active_cell_iterator & cell,
4107 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
4108 {
4109 active_neighbors.clear();
4110 for (const unsigned int n : cell->face_indices())
4111 if (!cell->at_boundary(n))
4112 {
4113 if (MeshType::dimension == 1)
4114 {
4115 // check children of neighbor. note
4116 // that in 1d children of the neighbor
4117 // may be further refined. In 1d the
4118 // case is simple since we know what
4119 // children bound to the present cell
4120 typename MeshType::cell_iterator neighbor_child =
4121 cell->neighbor(n);
4122 if (!neighbor_child->is_active())
4123 {
4124 while (neighbor_child->has_children())
4125 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
4126
4127 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
4128 ExcInternalError());
4129 }
4130 active_neighbors.push_back(neighbor_child);
4131 }
4132 else
4133 {
4134 if (cell->face(n)->has_children())
4135 // this neighbor has children. find
4136 // out which border to the present
4137 // cell
4138 for (unsigned int c = 0;
4139 c < cell->face(n)->n_active_descendants();
4140 ++c)
4141 active_neighbors.push_back(
4142 cell->neighbor_child_on_subface(n, c));
4143 else
4144 {
4145 // the neighbor must be active
4146 // himself
4147 Assert(cell->neighbor(n)->is_active(), ExcInternalError());
4148 active_neighbors.push_back(cell->neighbor(n));
4149 }
4150 }
4151 }
4152 }
4153
4154
4155
4156 template <typename CellIterator>
4157 std::size_t
4158 PeriodicFacePair<CellIterator>::memory_consumption() const
4159 {
4160 return sizeof(*this) + matrix.memory_consumption();
4161 }
4162
4163
4164
4165 namespace internal
4166 {
4167 namespace ProjectToObject
4168 {
4181 struct CrossDerivative
4182 {
4183 const unsigned int direction_0;
4184 const unsigned int direction_1;
4185
4186 CrossDerivative(const unsigned int d0, const unsigned int d1);
4187 };
4188
4189 inline CrossDerivative::CrossDerivative(const unsigned int d0,
4190 const unsigned int d1)
4191 : direction_0(d0)
4192 , direction_1(d1)
4193 {}
4194
4195
4196
4201 template <typename F>
4202 inline auto
4203 centered_first_difference(const double center,
4204 const double step,
4205 const F &f) -> decltype(f(center) - f(center))
4206 {
4207 return (f(center + step) - f(center - step)) / (2.0 * step);
4208 }
4209
4210
4211
4216 template <typename F>
4217 inline auto
4218 centered_second_difference(const double center,
4219 const double step,
4220 const F &f) -> decltype(f(center) - f(center))
4221 {
4222 return (f(center + step) - 2.0 * f(center) + f(center - step)) /
4223 (step * step);
4224 }
4225
4226
4227
4237 template <int structdim, typename F>
4238 inline auto
4239 cross_stencil(
4240 const CrossDerivative cross_derivative,
4241 const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
4242 const double step,
4243 const F &f) -> decltype(f(center) - f(center))
4244 {
4245 Tensor<1, GeometryInfo<structdim>::vertices_per_cell> simplex_vector;
4246 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
4247 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
4248 return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
4249 1.0 / 3.0 * f(center - simplex_vector) +
4250 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
4251 step;
4252 }
4253
4254
4255
4262 template <int spacedim, int structdim, typename F>
4263 inline double
4264 gradient_entry(
4265 const unsigned int row_n,
4266 const unsigned int dependent_direction,
4267 const Point<spacedim> &p0,
4268 const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
4269 const double step,
4270 const F & f)
4271 {
4272 Assert(row_n < GeometryInfo<structdim>::vertices_per_cell &&
4273 dependent_direction <
4274 GeometryInfo<structdim>::vertices_per_cell,
4275 ExcMessage("This function assumes that the last weight is a "
4276 "dependent variable (and hence we cannot take its "
4277 "derivative directly)."));
4278 Assert(row_n != dependent_direction,
4279 ExcMessage(
4280 "We cannot differentiate with respect to the variable "
4281 "that is assumed to be dependent."));
4282
4283 const Point<spacedim> manifold_point = f(center);
4284 const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
4285 {row_n, dependent_direction}, center, step, f);
4286 double entry = 0.0;
4287 for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
4288 entry +=
4289 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
4290 return entry;
4291 }
4292
4298 template <typename Iterator, int spacedim, int structdim>
4299 Point<spacedim>
4300 project_to_d_linear_object(const Iterator & object,
4301 const Point<spacedim> &trial_point)
4302 {
4303 // let's look at this for simplicity for a quadrilateral
4304 // (structdim==2) in a space with spacedim>2 (notate trial_point by
4305 // y): all points on the surface are given by
4306 // x(\xi) = sum_i v_i phi_x(\xi)
4307 // where v_i are the vertices of the quadrilateral, and
4308 // \xi=(\xi_1,\xi_2) are the reference coordinates of the
4309 // quadrilateral. so what we are trying to do is find a point x on the
4310 // surface that is closest to the point y. there are different ways to
4311 // solve this problem, but in the end it's a nonlinear problem and we
4312 // have to find reference coordinates \xi so that J(\xi) = 1/2 ||
4313 // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is
4314 // structdim-linear in \xi, so J(\xi) is a polynomial of degree
4315 // 2*structdim that we'd like to minimize. unless structdim==1, we'll
4316 // have to use a Newton method to find the answer. This leads to the
4317 // following formulation of Newton steps:
4318 //
4319 // Given \xi_k, find \delta\xi_k so that
4320 // H_k \delta\xi_k = - F_k
4321 // where H_k is an approximation to the second derivatives of J at
4322 // \xi_k, and F_k is the first derivative of J. We'll iterate this a
4323 // number of times until the right hand side is small enough. As a
4324 // stopping criterion, we terminate if ||\delta\xi||<eps.
4325 //
4326 // As for the Hessian, the best choice would be
4327 // H_k = J''(\xi_k)
4328 // but we'll opt for the simpler Gauss-Newton form
4329 // H_k = A^T A
4330 // i.e.
4331 // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
4332 // \partial_n phi_i *
4333 // \partial_m phi_j
4334 // we start at xi=(0.5, 0.5).
4335 Point<structdim> xi;
4336 for (unsigned int d = 0; d < structdim; ++d)
4337 xi[d] = 0.5;
4338
4339 Point<spacedim> x_k;
4340 for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
4341 x_k += object->vertex(i) *
4342 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
4343
4344 do
4345 {
4346 Tensor<1, structdim> F_k;
4347 for (const unsigned int i :
4348 GeometryInfo<structdim>::vertex_indices())
4349 F_k +=
4350 (x_k - trial_point) * object->vertex(i) *
4351 GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
4352 i);
4353
4354 Tensor<2, structdim> H_k;
4355 for (const unsigned int i :
4356 GeometryInfo<structdim>::vertex_indices())
4357 for (const unsigned int j :
4358 GeometryInfo<structdim>::vertex_indices())
4359 {
4360 Tensor<2, structdim> tmp = outer_product(
4361 GeometryInfo<structdim>::d_linear_shape_function_gradient(
4362 xi, i),
4363 GeometryInfo<structdim>::d_linear_shape_function_gradient(
4364 xi, j));
4365 H_k += (object->vertex(i) * object->vertex(j)) * tmp;
4366 }
4367
4368 const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
4369 xi += delta_xi;
4370
4371 x_k = Point<spacedim>();
4372 for (const unsigned int i :
4373 GeometryInfo<structdim>::vertex_indices())
4374 x_k += object->vertex(i) *
4375 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
4376
4377 if (delta_xi.norm() < 1e-7)
4378 break;
4379 }
4380 while (true);
4381
4382 return x_k;
4383 }
4384 } // namespace ProjectToObject
4385 } // namespace internal
4386
4387
4388
4389 namespace internal
4390 {
4391 // We hit an internal compiler error in ICC 15 if we define this as a lambda
4392 // inside the project_to_object function below.
4393 template <int structdim>
4394 inline bool
4395 weights_are_ok(
4396 const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v)
4397 {
4398 // clang has trouble figuring out structdim here, so define it
4399 // again:
4400 static const std::size_t n_vertices_per_cell =
4401 Tensor<1, GeometryInfo<structdim>::vertices_per_cell>::
4402 n_independent_components;
4403 std::array<double, n_vertices_per_cell> copied_weights;
4404 for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4405 {
4406 copied_weights[i] = v[i];
4407 if (v[i] < 0.0 || v[i] > 1.0)
4408 return false;
4409 }
4410
4411 // check the sum: try to avoid some roundoff errors by summing in order
4412 std::sort(copied_weights.begin(), copied_weights.end());
4413 const double sum =
4414 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4415 return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4416 }
4417 } // namespace internal
4418
4419 template <typename Iterator>
4420 Point<Iterator::AccessorType::space_dimension>
4421 project_to_object(
4422 const Iterator & object,
4423 const Point<Iterator::AccessorType::space_dimension> &trial_point)
4424 {
4425 const int spacedim = Iterator::AccessorType::space_dimension;
4426 const int structdim = Iterator::AccessorType::structure_dimension;
4427
4428 Point<spacedim> projected_point = trial_point;
4429
4430 if (structdim >= spacedim)
4431 return projected_point;
4432 else if (structdim == 1 || structdim == 2)
4433 {
4434 using namespace internal::ProjectToObject;
4435 // Try to use the special flat algorithm for quads (this is better
4436 // than the general algorithm in 3d). This does not take into account
4437 // whether projected_point is outside the quad, but we optimize along
4438 // lines below anyway:
4439 const int dim = Iterator::AccessorType::dimension;
4440 const Manifold<dim, spacedim> &manifold = object->get_manifold();
4441 if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4442 &manifold) != nullptr)
4443 {
4444 projected_point =
4445 project_to_d_linear_object<Iterator, spacedim, structdim>(
4446 object, trial_point);
4447 }
4448 else
4449 {
4450 // We want to find a point on the convex hull (defined by the
4451 // vertices of the object and the manifold description) that is
4452 // relatively close to the trial point. This has a few issues:
4453 //
4454 // 1. For a general convex hull we are not guaranteed that a unique
4455 // minimum exists.
4456 // 2. The independent variables in the optimization process are the
4457 // weights given to Manifold::get_new_point, which must sum to 1,
4458 // so we cannot use standard finite differences to approximate a
4459 // gradient.
4460 //
4461 // There is not much we can do about 1., but for 2. we can derive
4462 // finite difference stencils that work on a structdim-dimensional
4463 // simplex and rewrite the optimization problem to use those
4464 // instead. Consider the structdim 2 case and let
4465 //
4466 // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4467 // c2, c3})
4468 //
4469 // where {c0, c1, c2, c3} are the weights for the four vertices on
4470 // the quadrilateral. We seek to minimize the Euclidean distance
4471 // between F(...) and trial_point. We can solve for c3 in terms of
4472 // the other weights and get, for one coordinate direction
4473 //
4474 // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4475 // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4476 //
4477 // where we substitute back in for c3 after taking the
4478 // derivative. We can compute a stencil for the cross derivative
4479 // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4480 // (and gradient_entry computes the sum over the independent
4481 // variables). Below, we somewhat arbitrarily pick the last
4482 // component as the dependent one.
4483 //
4484 // Since we can now calculate derivatives of the objective
4485 // function we can use gradient descent to minimize it.
4486 //
4487 // Of course, this is much simpler in the structdim = 1 case (we
4488 // could rewrite the projection as a 1d optimization problem), but
4489 // to reduce the potential for bugs we use the same code in both
4490 // cases.
4491 const double step_size = object->diameter() / 64.0;
4492
4493 constexpr unsigned int n_vertices_per_cell =
4494 GeometryInfo<structdim>::vertices_per_cell;
4495
4496 std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4497 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4498 ++vertex_n)
4499 vertices[vertex_n] = object->vertex(vertex_n);
4500
4501 auto get_point_from_weights =
4502 [&](const Tensor<1, n_vertices_per_cell> &weights)
4503 -> Point<spacedim> {
4504 return object->get_manifold().get_new_point(
4505 make_array_view(vertices.begin(), vertices.end()),
4506 make_array_view(weights.begin_raw(), weights.end_raw()));
4507 };
4508
4509 // pick the initial weights as (normalized) inverse distances from
4510 // the trial point:
4511 Tensor<1, n_vertices_per_cell> guess_weights;
4512 double guess_weights_sum = 0.0;
4513 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4514 ++vertex_n)
4515 {
4516 const double distance =
4517 vertices[vertex_n].distance(trial_point);
4518 if (distance == 0.0)
4519 {
4520 guess_weights = 0.0;
4521 guess_weights[vertex_n] = 1.0;
4522 guess_weights_sum = 1.0;
4523 break;
4524 }
4525 else
4526 {
4527 guess_weights[vertex_n] = 1.0 / distance;
4528 guess_weights_sum += guess_weights[vertex_n];
4529 }
4530 }
4531 guess_weights /= guess_weights_sum;
4532 Assert(internal::weights_are_ok<structdim>(guess_weights),
4533 ExcInternalError());
4534
4535 // The optimization algorithm consists of two parts:
4536 //
4537 // 1. An outer loop where we apply the gradient descent algorithm.
4538 // 2. An inner loop where we do a line search to find the optimal
4539 // length of the step one should take in the gradient direction.
4540 //
4541 for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4542 {
4543 const unsigned int dependent_direction =
4544 n_vertices_per_cell - 1;
4545 Tensor<1, n_vertices_per_cell> current_gradient;
4546 for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4547 ++row_n)
4548 {
4549 if (row_n != dependent_direction)
4550 {
4551 current_gradient[row_n] =
4552 gradient_entry<spacedim, structdim>(
4553 row_n,
4554 dependent_direction,
4555 trial_point,
4556 guess_weights,
4557 step_size,
4558 get_point_from_weights);
4559
4560 current_gradient[dependent_direction] -=
4561 current_gradient[row_n];
4562 }
4563 }
4564
4565 // We need to travel in the -gradient direction, as noted
4566 // above, but we may not want to take a full step in that
4567 // direction; instead, guess that we will go -0.5*gradient and
4568 // do quasi-Newton iteration to pick the best multiplier. The
4569 // goal is to find a scalar alpha such that
4570 //
4571 // F(x - alpha g)
4572 //
4573 // is minimized, where g is the gradient and F is the
4574 // objective function. To find the optimal value we find roots
4575 // of the derivative of the objective function with respect to
4576 // alpha by Newton iteration, where we approximate the first
4577 // and second derivatives of F(x - alpha g) with centered
4578 // finite differences.
4579 double gradient_weight = -0.5;
4580 auto gradient_weight_objective_function =
4581 [&](const double gradient_weight_guess) -> double {
4582 return (trial_point -
4583 get_point_from_weights(guess_weights +
4584 gradient_weight_guess *
4585 current_gradient))
4586 .norm_square();
4587 };
4588
4589 for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4590 {
4591 const double update_numerator = centered_first_difference(
4592 gradient_weight,
4593 step_size,
4594 gradient_weight_objective_function);
4595 const double update_denominator =
4596 centered_second_difference(
4597 gradient_weight,
4598 step_size,
4599 gradient_weight_objective_function);
4600
4601 // avoid division by zero. Note that we limit the gradient
4602 // weight below
4603 if (std::abs(update_denominator) == 0.0)
4604 break;
4605 gradient_weight =
4606 gradient_weight - update_numerator / update_denominator;
4607
4608 // Put a fairly lenient bound on the largest possible
4609 // gradient (things tend to be locally flat, so the gradient
4610 // itself is usually small)
4611 if (std::abs(gradient_weight) > 10)
4612 {
4613 gradient_weight = -10.0;
4614 break;
4615 }
4616 }
4617
4618 // It only makes sense to take convex combinations with weights
4619 // between zero and one. If the update takes us outside of this
4620 // region then rescale the update to stay within the region and
4621 // try again
4622 Tensor<1, n_vertices_per_cell> tentative_weights =
4623 guess_weights + gradient_weight * current_gradient;
4624
4625 double new_gradient_weight = gradient_weight;
4626 for (unsigned int iteration_count = 0; iteration_count < 40;
4627 ++iteration_count)
4628 {
4629 if (internal::weights_are_ok<structdim>(tentative_weights))
4630 break;
4631
4632 for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4633 {
4634 if (tentative_weights[i] < 0.0)
4635 {
4636 tentative_weights -=
4637 (tentative_weights[i] / current_gradient[i]) *
4638 current_gradient;
4639 }
4640 if (tentative_weights[i] < 0.0 ||
4641 1.0 < tentative_weights[i])
4642 {
4643 new_gradient_weight /= 2.0;
4644 tentative_weights =
4645 guess_weights +
4646 new_gradient_weight * current_gradient;
4647 }
4648 }
4649 }
4650
4651 // the update might still send us outside the valid region, so
4652 // check again and quit if the update is still not valid
4653 if (!internal::weights_are_ok<structdim>(tentative_weights))
4654 break;
4655
4656 // if we cannot get closer by traveling in the gradient
4657 // direction then quit
4658 if (get_point_from_weights(tentative_weights)
4659 .distance(trial_point) <
4660 get_point_from_weights(guess_weights).distance(trial_point))
4661 guess_weights = tentative_weights;
4662 else
4663 break;
4664 Assert(internal::weights_are_ok<structdim>(guess_weights),
4665 ExcInternalError());
4666 }
4667 Assert(internal::weights_are_ok<structdim>(guess_weights),
4668 ExcInternalError());
4669 projected_point = get_point_from_weights(guess_weights);
4670 }
4671
4672 // if structdim == 2 and the optimal point is not on the interior then
4673 // we may be able to get a more accurate result by projecting onto the
4674 // lines.
4675 if (structdim == 2)
4676 {
4677 std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4678 line_projections;
4679 for (unsigned int line_n = 0;
4680 line_n < GeometryInfo<structdim>::lines_per_cell;
4681 ++line_n)
4682 {
4683 line_projections[line_n] =
4684 project_to_object(object->line(line_n), trial_point);
4685 }
4686 std::sort(line_projections.begin(),
4687 line_projections.end(),
4688 [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4689 return a.distance(trial_point) <
4690 b.distance(trial_point);
4691 });
4692 if (line_projections[0].distance(trial_point) <
4693 projected_point.distance(trial_point))
4694 projected_point = line_projections[0];
4695 }
4696 }
4697 else
4698 {
4699 Assert(false, ExcNotImplemented());
4700 return projected_point;
4701 }
4702
4703 return projected_point;
4704 }
4705
4706
4707
4708 namespace internal
4709 {
4710 template <typename DataType,
4711 typename MeshType,
4712 typename MeshCellIteratorType>
4713 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
4714 inline void exchange_cell_data(
4715 const MeshType &mesh,
4716 const std::function<
4717 std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4718 const std::function<void(const MeshCellIteratorType &, const DataType &)>
4719 & unpack,
4720 const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4721 const std::function<void(
4722 const std::function<void(const MeshCellIteratorType &,
4723 const types::subdomain_id)> &)> &process_cells,
4724 const std::function<std::set<types::subdomain_id>(
4725 const parallel::TriangulationBase<MeshType::dimension,
4726 MeshType::space_dimension> &)>
4727 &compute_ghost_owners)
4728 {
4729# ifndef DEAL_II_WITH_MPI
4730 (void)mesh;
4731 (void)pack;
4732 (void)unpack;
4733 (void)cell_filter;
4734 (void)process_cells;
4735 (void)compute_ghost_owners;
4736 Assert(false, ExcNeedsMPI());
4737# else
4738 constexpr int dim = MeshType::dimension;
4739 constexpr int spacedim = MeshType::space_dimension;
4740 auto tria =
4741 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4742 &mesh.get_triangulation());
4743 Assert(
4744 tria != nullptr,
4745 ExcMessage(
4746 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4747
4748 if (const auto tria = dynamic_cast<
4749 const parallel::shared::Triangulation<dim, spacedim> *>(
4750 &mesh.get_triangulation()))
4751 {
4752 Assert(
4753 tria->with_artificial_cells(),
4754 ExcMessage(
4755 "The functions GridTools::exchange_cell_data_to_ghosts() and "
4756 "GridTools::exchange_cell_data_to_level_ghosts() can only "
4757 "operate on a single layer of ghost cells. However, you have "
4758 "given a Triangulation object of type "
4759 "parallel::shared::Triangulation without artificial cells "
4760 "resulting in arbitrary numbers of ghost layers."));
4761 }
4762
4763 // build list of cells to request for each neighbor
4764 std::set<::types::subdomain_id> ghost_owners =
4765 compute_ghost_owners(*tria);
4766 std::map<::types::subdomain_id,
4767 std::vector<typename CellId::binary_type>>
4768 neighbor_cell_list;
4769
4770 for (const auto ghost_owner : ghost_owners)
4771 neighbor_cell_list[ghost_owner] = {};
4772
4773 process_cells([&](const auto &cell, const auto key) -> void {
4774 if (cell_filter(cell))
4775 {
4776 constexpr int spacedim = MeshType::space_dimension;
4777 neighbor_cell_list[key].emplace_back(
4778 cell->id().template to_binary<spacedim>());
4779 }
4780 });
4781
4782 Assert(ghost_owners.size() == neighbor_cell_list.size(),
4783 ExcInternalError());
4784
4785
4786 // Before sending & receiving, make sure we protect this section with
4787 // a mutex:
4788 static Utilities::MPI::CollectiveMutex mutex;
4789 Utilities::MPI::CollectiveMutex::ScopedLock lock(
4790 mutex, tria->get_communicator());
4791
4792 const int mpi_tag =
4793 Utilities::MPI::internal::Tags::exchange_cell_data_request;
4794 const int mpi_tag_reply =
4795 Utilities::MPI::internal::Tags::exchange_cell_data_reply;
4796
4797 // send our requests
4798 std::vector<MPI_Request> requests(ghost_owners.size());
4799 {
4800 unsigned int idx = 0;
4801 for (const auto &it : neighbor_cell_list)
4802 {
4803 // send the data about the relevant cells
4804 const int ierr = MPI_Isend(it.second.data(),
4805 it.second.size() * sizeof(it.second[0]),
4806 MPI_BYTE,
4807 it.first,
4808 mpi_tag,
4809 tria->get_communicator(),
4810 &requests[idx]);
4811 AssertThrowMPI(ierr);
4812 ++idx;
4813 }
4814 }
4815
4816 // receive requests and reply with the results
4817 std::vector<MPI_Request> reply_requests(ghost_owners.size());
4818 std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4819
4820 for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4821 {
4822 MPI_Status status;
4823 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4824 mpi_tag,
4825 tria->get_communicator(),
4826 &status);
4827 AssertThrowMPI(ierr);
4828
4829 int len;
4830 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4831 AssertThrowMPI(ierr);
4832 Assert(len % sizeof(typename CellId::binary_type) == 0,
4833 ExcInternalError());
4834
4835 const unsigned int n_cells =
4836 len / sizeof(typename CellId::binary_type);
4837 std::vector<typename CellId::binary_type> cells_with_requests(
4838 n_cells);
4839 std::vector<DataType> data_to_send;
4840 data_to_send.reserve(n_cells);
4841 std::vector<bool> cell_carries_data(n_cells, false);
4842
4843 ierr = MPI_Recv(cells_with_requests.data(),
4844 len,
4845 MPI_BYTE,
4846 status.MPI_SOURCE,
4847 status.MPI_TAG,
4848 tria->get_communicator(),
4849 &status);
4850 AssertThrowMPI(ierr);
4851
4852 // store data for each cell
4853 for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4854 {
4855 const auto cell =
4856 tria->create_cell_iterator(CellId(cells_with_requests[c]));
4857
4858 MeshCellIteratorType mesh_it(tria,
4859 cell->level(),
4860 cell->index(),
4861 &mesh);
4862
4863 const std_cxx17::optional<DataType> data = pack(mesh_it);
4864 if (data)
4865 {
4866 data_to_send.emplace_back(*data);
4867 cell_carries_data[c] = true;
4868 }
4869 }
4870
4871 // collect data for sending the reply in a buffer
4872
4873 // (a) make room for storing the local offsets in case we receive
4874 // other data
4875 sendbuffers[idx].resize(sizeof(std::size_t));
4876
4877 // (b) append the actual data and store how much memory it
4878 // corresponds to, which we then insert into the leading position of
4879 // the sendbuffer
4880 std::size_t size_of_send =
4881 Utilities::pack(data_to_send,
4882 sendbuffers[idx],
4883 /*enable_compression*/ false);
4884 std::memcpy(sendbuffers[idx].data(),
4885 &size_of_send,
4886 sizeof(std::size_t));
4887
4888 // (c) append information of certain cells that got left out in case
4889 // we need it
4890 if (data_to_send.size() < n_cells)
4891 Utilities::pack(cell_carries_data,
4892 sendbuffers[idx],
4893 /*enable_compression*/ false);
4894
4895 // send data
4896 ierr = MPI_Isend(sendbuffers[idx].data(),
4897 sendbuffers[idx].size(),
4898 MPI_BYTE,
4899 status.MPI_SOURCE,
4900 mpi_tag_reply,
4901 tria->get_communicator(),
4902 &reply_requests[idx]);
4903 AssertThrowMPI(ierr);
4904 }
4905
4906 // finally receive the replies
4907 std::vector<char> receive;
4908 for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4909 {
4910 MPI_Status status;
4911 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4912 mpi_tag_reply,
4913 tria->get_communicator(),
4914 &status);
4915 AssertThrowMPI(ierr);
4916
4917 int len;
4918 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4919 AssertThrowMPI(ierr);
4920
4921 receive.resize(len);
4922
4923 ierr = MPI_Recv(receive.data(),
4924 len,
4925 MPI_BYTE,
4926 status.MPI_SOURCE,
4927 status.MPI_TAG,
4928 tria->get_communicator(),
4929 &status);
4930 AssertThrowMPI(ierr);
4931
4932 // (a) first determine the length of the data section in the
4933 // received buffer
4934 auto data_iterator = receive.begin();
4935 std::size_t size_of_received_data =
4936 Utilities::unpack<std::size_t>(data_iterator,
4937 data_iterator + sizeof(std::size_t));
4938 data_iterator += sizeof(std::size_t);
4939
4940 // (b) unpack the data section in the indicated region
4941 auto received_data = Utilities::unpack<std::vector<DataType>>(
4942 data_iterator,
4943 data_iterator + size_of_received_data,
4944 /*enable_compression*/ false);
4945 data_iterator += size_of_received_data;
4946
4947 // (c) check if the received data contained fewer entries than the
4948 // number of cells we identified in the beginning, in which case we
4949 // need to extract the boolean vector with the relevant information
4950 const std::vector<typename CellId::binary_type> &this_cell_list =
4951 neighbor_cell_list[status.MPI_SOURCE];
4952 AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4953 std::vector<bool> cells_with_data;
4954 if (received_data.size() < this_cell_list.size())
4955 {
4956 cells_with_data = Utilities::unpack<std::vector<bool>>(
4957 data_iterator, receive.end(), /*enable_compression*/ false);
4958 AssertDimension(cells_with_data.size(), this_cell_list.size());
4959 }
4960
4961 // (d) go through the received data and call the user-provided
4962 // unpack function
4963 auto received_data_iterator = received_data.begin();
4964 for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4965 if (cells_with_data.empty() || cells_with_data[c])
4966 {
4967 const typename Triangulation<dim, spacedim>::cell_iterator
4968 tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4969
4970 MeshCellIteratorType cell(tria,
4971 tria_cell->level(),
4972 tria_cell->index(),
4973 &mesh);
4974
4975 unpack(cell, *received_data_iterator);
4976 ++received_data_iterator;
4977 }
4978 }
4979
4980 // make sure that all communication is finished
4981 // when we leave this function.
4982 if (requests.size() > 0)
4983 {
4984 const int ierr =
4985 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4986 AssertThrowMPI(ierr);
4987 }
4988 if (reply_requests.size() > 0)
4989 {
4990 const int ierr = MPI_Waitall(reply_requests.size(),
4991 reply_requests.data(),
4992 MPI_STATUSES_IGNORE);
4993 AssertThrowMPI(ierr);
4994 }
4995
4996
4997# endif // DEAL_II_WITH_MPI
4998 }
4999
5000 } // namespace internal
5001
5002 template <typename DataType, typename MeshType>
5003 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
5004 inline void exchange_cell_data_to_ghosts(
5005 const MeshType & mesh,
5006 const std::function<std_cxx17::optional<DataType>(
5007 const typename MeshType::active_cell_iterator &)> &pack,
5008 const std::function<void(const typename MeshType::active_cell_iterator &,
5009 const DataType &)> & unpack,
5010 const std::function<bool(const typename MeshType::active_cell_iterator &)>
5011 &cell_filter)
5012 {
5013# ifndef DEAL_II_WITH_MPI
5014 (void)mesh;
5015 (void)pack;
5016 (void)unpack;
5017 (void)cell_filter;
5018 Assert(false, ExcNeedsMPI());
5019# else
5020 internal::exchange_cell_data<DataType,
5021 MeshType,
5022 typename MeshType::active_cell_iterator>(
5023 mesh,
5024 pack,
5025 unpack,
5026 cell_filter,
5027 [&](const auto &process) {
5028 for (const auto &cell : mesh.active_cell_iterators())
5029 if (cell->is_ghost())
5030 process(cell, cell->subdomain_id());
5031 },
5032 [](const auto &tria) { return tria.ghost_owners(); });
5033# endif
5034 }
5035
5036
5037
5038 template <typename DataType, typename MeshType>
5039 DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
5040 inline void exchange_cell_data_to_level_ghosts(
5041 const MeshType & mesh,
5042 const std::function<std_cxx17::optional<DataType>(
5043 const typename MeshType::level_cell_iterator &)> &pack,
5044 const std::function<void(const typename MeshType::level_cell_iterator &,
5045 const DataType &)> & unpack,
5046 const std::function<bool(const typename MeshType::level_cell_iterator &)>
5047 &cell_filter)
5048 {
5049# ifndef DEAL_II_WITH_MPI
5050 (void)mesh;
5051 (void)pack;
5052 (void)unpack;
5053 (void)cell_filter;
5054 Assert(false, ExcNeedsMPI());
5055# else
5056 internal::exchange_cell_data<DataType,
5057 MeshType,
5058 typename MeshType::level_cell_iterator>(
5059 mesh,
5060 pack,
5061 unpack,
5062 cell_filter,
5063 [&](const auto &process) {
5064 for (const auto &cell : mesh.cell_iterators())
5065 if (cell->is_ghost_on_level())
5066 process(cell, cell->level_subdomain_id());
5067 },
5068 [](const auto &tria) { return tria.level_ghost_owners(); });
5069# endif
5070 }
5071} // namespace GridTools
5072
5073#endif // DOXYGEN
5074
5075DEAL_II_NAMESPACE_CLOSE
5076
5077#endif
make_array_view
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:704
bounding_box.h
bounding_box.h
CellId::binary_type
std::array< unsigned int, 4 > binary_type
Definition cell_id.h:81
GridTools::MarchingCubeAlgorithm::value_type
typename VectorType::value_type value_type
Definition grid_tools.h:3745
GridTools::MarchingCubeAlgorithm::n_subdivisions
const unsigned int n_subdivisions
Definition grid_tools.h:3873
GridTools::MarchingCubeAlgorithm::create_quadrature_rule
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition grid_tools.cc:7332
GridTools::MarchingCubeAlgorithm::process_cell
void process_cell(const typename DoFHandler< dim >::active_cell_iterator &cell, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim > > &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 > > &cells) const
Definition grid_tools.cc:7403
GridTools::MarchingCubeAlgorithm::process_sub_cell
void process_sub_cell(const std::vector< value_type > &, const std::vector< Point< 1 > > &, const std::vector< unsigned int > &, const double, std::vector< Point< 1 > > &, std::vector< CellData< 1 > > &, const bool) const
Definition grid_tools.h:3831
GridTools::MarchingCubeAlgorithm::process
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim > > &vertices, std::vector< CellData< dim==1 ? 1 :dim - 1 > > &cells) const
Definition grid_tools.cc:7367
MPI_Comm
Mapping
Abstract base class for mapping classes.
Definition mapping.h:317
Point
Definition point.h:112
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Tensor::norm
numbers::NumberTraits< Number >::real_type norm() const
Triangulation::get_communicator
virtual MPI_Comm get_communicator() const
Triangulation::create_cell_iterator
cell_iterator create_cell_iterator(const CellId &cell_id) const
internal::ActiveCellIterator::type
typename MeshType::active_cell_iterator type
Definition grid_tools.h:108
concepts::is_triangulation_or_dof_handler
Definition template_constraints.h:978
concepts::is_valid_dim_spacedim
Definition template_constraints.h:720
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition config.h:172
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_CXX20_REQUIRES
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:160
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
center
Point< 3 > center
Definition data_out_base.cc:268
vertices
Point< 3 > vertices[4]
Definition data_out_base.cc:273
dof_handler.h
fe_values.h
geometry_info.h
level
unsigned int level
Definition grid_out.cc:4618
vertex_indices
unsigned int vertex_indices[2]
Definition grid_tools.cc:1360
DeclException0
#define DeclException0(Exception0)
Definition exceptions.h:465
ExcMeshNotOrientable
static ::ExceptionBase & ExcMeshNotOrientable()
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
StandardExceptions::ExcNeedsMPI
static ::ExceptionBase & ExcNeedsMPI()
GridTools::ExcNonExistentSubdomain
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
GridTools::ExcVertexNotUsed
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition exceptions.h:1914
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
DeclExceptionMsg
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
GridTools::ExcTriangulationHasBeenRefined
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
GridTools::ExcInvalidNumberOfPartitions
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
GridTools::ExcScalingFactorNotPositive
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:441
GridTools::copy_boundary_to_manifold_id
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition grid_tools.cc:4989
GridTools::copy_material_to_manifold_id
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition grid_tools.cc:5081
GridTools::map_boundary_to_manifold_ids
void map_boundary_to_manifold_ids(const std::vector< types::boundary_id > &src_boundary_ids, const std::vector< types::manifold_id > &dst_manifold_ids, Triangulation< dim, spacedim > &tria, const std::vector< types::boundary_id > &reset_boundary_ids={})
Definition grid_tools.cc:5014
GridTools::assign_co_dimensional_manifold_indicators
void assign_co_dimensional_manifold_indicators(Triangulation< dim, spacedim > &tria, const std::function< types::manifold_id(const std::set< types::manifold_id > &)> &disambiguation_function=[](const std::set< types::manifold_id > &manifold_ids) { if(manifold_ids.size()==1) return *manifold_ids.begin();else return numbers::flat_manifold_id;}, bool overwrite_only_flat_manifold_ids=true)
Definition grid_tools.cc:5109
la_parallel_vector.h
la_vector.h
GridTools::internal::distributed_compute_point_locations
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim > > &points, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition grid_tools.cc:6178
GridTools::internal::distributed_compute_intersection_locations
DistributedComputeIntersectionLocationsInternal< structdim, spacedim > distributed_compute_intersection_locations(const Cache< dim, spacedim > &cache, const std::vector< std::vector< Point< spacedim > > > &intersection_requests, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const std::vector< bool > &marked_vertices, const double tolerance)
Definition grid_tools.cc:6563
GridTools::get_face_connectivity_of_cells
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition grid_tools.cc:3796
GridTools::delete_unused_vertices
void delete_unused_vertices(std::vector< Point< spacedim > > &vertices, std::vector< CellData< dim > > &cells, SubCellData &subcelldata)
Definition grid_tools.cc:656
GridTools::compute_mesh_predicate_bounding_box
std::vector< BoundingBox< MeshType::space_dimension > > compute_mesh_predicate_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const unsigned int refinement_level=0, const bool allow_merge=false, const unsigned int max_boxes=numbers::invalid_unsigned_int)
Definition grid_tools.cc:3147
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:2202
GridTools::build_global_description_tree
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim > > &local_description, const MPI_Comm mpi_communicator)
Definition grid_tools.cc:6971
GridTools::get_active_child_cells
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
GridTools::partition_triangulation_zorder
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition grid_tools.cc:4136
GridTools::get_all_vertices_at_boundary
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:2352
GridTools::fix_up_distorted_child_cells
Triangulation< dim, spacedim >::DistortedCellList fix_up_distorted_child_cells(const typename Triangulation< dim, spacedim >::DistortedCellList &distorted_cells, Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:4947
GridTools::get_cells_at_coarsest_common_level
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
Definition grid_tools_dof_handlers.cc:1479
GridTools::rotate
void rotate(const double angle, Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:2146
GridTools::compute_point_locations
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim > > &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition grid_tools.cc:5607
GridTools::find_closest_vertex
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim > > &vertices, const Point< spacedim > &p)
Definition grid_tools.cc:6830
GridTools::transform
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
GridTools::find_closest_vertex_of_cell
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition grid_tools.cc:3059
GridTools::find_active_cell_around_point
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)
Definition grid_tools_dof_handlers.cc:440
GridTools::extract_used_vertices
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition grid_tools.cc:6811
GridTools::get_active_neighbors
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
GridTools::get_finest_common_cells
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
Definition grid_tools_dof_handlers.cc:1148
GridTools::get_locally_owned_vertices
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:4427
GridTools::affine_cell_approximation
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim > > &vertices)
Definition grid_tools.cc:309
GridTools::regularize_corner_cells
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition grid_tools.cc:5291
GridTools::remove_anisotropy
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition grid_tools.cc:5262
GridTools::consistently_order_cells
void consistently_order_cells(std::vector< CellData< dim > > &cells)
Definition grid_tools.cc:2026
GridTools::minimal_cell_diameter
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition grid_tools.cc:4454
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:2136
GridTools::find_cells_adjacent_to_vertex
std::vector< typename MeshType::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType &container, const unsigned int vertex_index)
Definition grid_tools_dof_handlers.cc:202
GridTools::compute_active_cell_halo_layer
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)
Definition grid_tools_dof_handlers.cc:767
GridTools::cell_measure
double cell_measure(const std::vector< Point< dim > > &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< dim >::vertices_per_cell])
GridTools::compute_active_cell_layer_within_distance
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)
Definition grid_tools_dof_handlers.cc:911
GridTools::remove_hanging_nodes
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition grid_tools.cc:5229
GridTools::compute_ghost_cell_halo_layer
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
Definition grid_tools_dof_handlers.cc:888
GridTools::collect_periodic_faces
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 >())
Definition grid_tools_dof_handlers.cc:2232
GridTools::project_to_object
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
GridTools::laplace_transform
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)
GridTools::invert_cells_with_negative_measure
std::size_t invert_cells_with_negative_measure(const std::vector< Point< spacedim > > &all_vertices, std::vector< CellData< dim > > &cells)
Definition grid_tools.cc:923
GridTools::partition_multigrid_levels
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:4242
GridTools::delete_duplicated_vertices
void delete_duplicated_vertices(std::vector< Point< spacedim > > &all_vertices, std::vector< CellData< dim > > &cells, SubCellData &subcelldata, std::vector< unsigned int > &considered_vertices, const double tol=1e-12)
Definition grid_tools.cc:761
GridTools::compute_vertices_with_ghost_neighbors
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:7157
GridTools::collect_coinciding_vertices
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int > > &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition grid_tools.cc:7034
GridTools::count_cells_with_subdomain_association
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition grid_tools.cc:4411
GridTools::guess_point_owner
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const std::vector< Point< spacedim > > &points)
Definition grid_tools.cc:3294
GridTools::get_dof_to_support_patch_map
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)
Definition grid_tools_dof_handlers.cc:1742
GridTools::exchange_cell_data_to_ghosts
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
GridTools::compute_cell_halo_layer_on_level
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)
Definition grid_tools_dof_handlers.cc:805
GridTools::partition_triangulation
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition grid_tools.cc:3895
GridTools::compute_aspect_ratio_of_cells
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition grid_tools.cc:333
GridTools::vertex_to_cell_map
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:3392
GridTools::compute_ghost_cell_layer_within_distance
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
Definition grid_tools_dof_handlers.cc:1057
GridTools::compute_maximum_aspect_ratio
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition grid_tools.cc:399
GridTools::compute_point_locations_try_all
return_type compute_point_locations_try_all(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim > > &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition grid_tools.cc:5636
GridTools::have_same_coarse_mesh
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
Definition grid_tools_dof_handlers.cc:1256
GridTools::get_patch_around_cell
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
Definition grid_tools_dof_handlers.cc:1431
GridTools::exchange_cell_data_to_level_ghosts
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
GridTools::get_vertex_connectivity_of_cells
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition grid_tools.cc:3830
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:144
GridTools::get_subdomain_association
std::vector< types::subdomain_id > get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const std::vector< CellId > &cell_ids)
Definition grid_tools.cc:4341
GridTools::distort_random
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition grid_tools.cc:2384
GridTools::vertex_to_cell_centers_directions
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > &vertex_to_cells)
Definition grid_tools.cc:2731
GridTools::build_triangulation_from_patch
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)
Definition grid_tools_dof_handlers.cc:1525
GridTools::diameter
double diameter(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:88
GridTools::get_vertex_connectivity_of_cells_on_level
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition grid_tools.cc:3859
GridTools::invert_all_negative_measure_cells
void invert_all_negative_measure_cells(const std::vector< Point< spacedim > > &all_vertices, std::vector< CellData< dim > > &cells)
Definition grid_tools.cc:1003
GridTools::compute_local_to_global_vertex_index_map
std::map< unsigned int, types::global_vertex_index > compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:3441
GridTools::find_all_active_cells_around_point
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)
Definition grid_tools_dof_handlers.cc:613
GridTools::compute_bounding_box
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition grid_tools.cc:415
GridTools::maximal_cell_diameter
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition grid_tools.cc:4481
GridTools::get_longest_direction
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition grid_tools.cc:5197
GridTools::get_coarse_mesh_description
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition grid_tools.cc:555
GridTools::exchange_local_bounding_boxes
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim > > &local_bboxes, const MPI_Comm mpi_communicator)
Definition grid_tools.cc:6875
GridTools::distributed_compute_point_locations
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim > > &local_points, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bboxes, const double tolerance=1e-10, const std::vector< bool > &marked_vertices={}, const bool enforce_unique_mapping=true)
Definition grid_tools.cc:5806
GridTools::orthogonal_equality
bool orthogonal_equality(std::bitset< 3 > &orientation, 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 >())
Definition grid_tools_dof_handlers.cc:2458
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:58
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Utilities::MPI::internal::Tags::exchange_cell_data_request
@ exchange_cell_data_request
grid_tools.h: exchange_cell_ghosts()
Definition mpi_tags.h:66
Utilities::MPI::internal::Tags::exchange_cell_data_reply
@ exchange_cell_data_reply
grid_tools.h: exchange_cell_ghosts()
Definition mpi_tags.h:69
Utilities::MPI::internal::Tags::exchange_cell_data_to_ghosts
@ exchange_cell_data_to_ghosts
GridTools::exchange_cell_data_to_ghosts():
Definition mpi_tags.h:60
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm mpi_communicator)
Utilities::pack
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1351
Utilities::unpack
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition utilities.h:1510
hp
Definition hp.h:118
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:13826
numbers::flat_manifold_id
const types::manifold_id flat_manifold_id
Definition types.h:286
std_cxx17
Definition c++.h:141
std
STL namespace.
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6008
types
Definition types.h:33
types::global_dof_index
unsigned int global_dof_index
Definition types.h:82
types::subdomain_id
unsigned int subdomain_id
Definition types.h:44
types::boundary_id
unsigned int boundary_id
Definition types.h:141
triangulation
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition p4est_wrappers.cc:69
petsc_vector.h
RTree
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition rtree.h:144
shared_tria.h
sparsity_tools.h
GridTools::CellDataTransferBuffer::load
void load(Archive &ar, const unsigned int)
Definition grid_tools.h:3682
GridTools::CellDataTransferBuffer::serialize
void serialize(Archive &archive, const unsigned int version)
GridTools::CellDataTransferBuffer::save
void save(Archive &ar, const unsigned int) const
Definition grid_tools.h:3656
GridTools::CellDataTransferBuffer::cell_ids
std::vector< CellId > cell_ids
Definition grid_tools.h:3639
GridTools::PeriodicFacePair::orientation
std::bitset< 3 > orientation
Definition grid_tools.h:3039
GridTools::PeriodicFacePair::memory_consumption
std::size_t memory_consumption() const
GridTools::PeriodicFacePair::matrix
FullMatrix< double > matrix
Definition grid_tools.h:3053
GridTools::internal::DistributedComputeIntersectionLocationsInternal::communicate_indices
std::map< unsigned int, std::vector< unsigned int > > communicate_indices(const std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > &point_recv_components, const MPI_Comm comm) const
Definition grid_tools.cc:6451
GridTools::internal::DistributedComputeIntersectionLocationsInternal::convert_to_distributed_compute_point_locations_internal
GridTools::internal::DistributedComputePointLocationsInternal< dim, spacedim > convert_to_distributed_compute_point_locations_internal(const unsigned int n_points_1D, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping, const bool consistent_numbering_of_sender_and_receiver=false) const
GridTools::internal::DistributedComputeIntersectionLocationsInternal::send_components
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, IntersectionType > > send_components
Definition grid_tools.h:1345
GridTools::internal::DistributedComputeIntersectionLocationsInternal::recv_components
std::vector< std::tuple< unsigned int, unsigned int, IntersectionType > > recv_components
Definition grid_tools.h:1356
GridTools::internal::DistributedComputeIntersectionLocationsInternal::IntersectionType
std::array<::Point< spacedim >, structdim+1 > IntersectionType
Definition grid_tools.h:1331
GridTools::internal::DistributedComputePointLocationsInternal::recv_components
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition grid_tools.h:1281
GridTools::internal::DistributedComputePointLocationsInternal::send_components
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition grid_tools.h:1256
always_return
Definition template_constraints.h:355
invert
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Definition symmetric_tensor.h:3458
outer_product
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition symmetric_tensor.h:3509
tria_accessor.h
vertex_to_coinciding_vertex_group
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
Definition tria_description.cc:718
comm
const MPI_Comm comm
Definition tria_description.cc:712
coinciding_vertex_groups
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
Definition tria_description.cc:717
tria
const ::Triangulation< dim, spacedim > & tria
Definition tria_description.cc:704
tria_iterator.h
trilinos_vector.h
Generated by doxygen 1.9.7