Reference documentation for deal.II version 9.4.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\}}\)
grid_tools.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2022 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 #include <deal.II/fe/mapping_q1.h>
36 
37 #include <deal.II/grid/manifold.h>
38 #include <deal.II/grid/tria.h>
39 #include <deal.II/grid/tria_accessor.h>
40 #include <deal.II/grid/tria_iterator.h>
41 
42 #include <deal.II/hp/dof_handler.h>
43 
44 #include <deal.II/lac/la_parallel_vector.h>
45 #include <deal.II/lac/la_vector.h>
46 #include <deal.II/lac/petsc_vector.h>
47 #include <deal.II/lac/sparsity_tools.h>
48 #include <deal.II/lac/trilinos_vector.h>
49 
50 #include <deal.II/numerics/rtree.h>
51 
52 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
53 #include <boost/archive/binary_iarchive.hpp>
54 #include <boost/archive/binary_oarchive.hpp>
55 #include <boost/random/mersenne_twister.hpp>
56 #include <boost/serialization/array.hpp>
57 #include <boost/serialization/vector.hpp>
58 
59 #ifdef DEAL_II_WITH_ZLIB
60 # include <boost/iostreams/device/back_inserter.hpp>
61 # include <boost/iostreams/filter/gzip.hpp>
62 # include <boost/iostreams/filtering_stream.hpp>
63 # include <boost/iostreams/stream.hpp>
64 #endif
65 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
66 
67 #include <bitset>
68 #include <list>
69 #include <set>
70 
71 DEAL_II_NAMESPACE_OPEN
72 
73 // Forward declarations
74 #ifndef DOXYGEN
75 namespace parallel
76 {
77  namespace distributed
78  {
79  template <int, int>
80  class Triangulation;
81  }
82 } // namespace parallel
83 
84 namespace hp
85 {
86  template <int, int>
87  class MappingCollection;
88 }
89 
90 class SparsityPattern;
91 
92 namespace GridTools
93 {
94  template <int dim, int spacedim>
95  class Cache;
96 }
97 #endif
98 
99 namespace internal
100 {
101  template <int dim, int spacedim, class MeshType>
102  class ActiveCellIterator
103  {
104  public:
105 #ifndef _MSC_VER
106  using type = typename MeshType::active_cell_iterator;
107 #else
108  using type = TriaActiveIterator<::CellAccessor<dim, spacedim>>;
109 #endif
110  };
111 
112 #ifdef _MSC_VER
113  template <int dim, int spacedim>
114  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
115  {
116  public:
117  using type =
118  TriaActiveIterator<::DoFCellAccessor<dim, spacedim, false>>;
119  };
120 #endif
121 } // namespace internal
122 
131 namespace GridTools
132 {
137 
144  template <int dim, int spacedim>
145  double
146  diameter(const Triangulation<dim, spacedim> &tria);
147 
174  template <int dim, int spacedim>
175  double
176  volume(const Triangulation<dim, spacedim> &tria,
177  const Mapping<dim, spacedim> & mapping =
178  (ReferenceCells::get_hypercube<dim>()
179 #ifndef _MSC_VER
180  .template get_default_linear_mapping<dim, spacedim>()
181 #else
182  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
183 #endif
184  ));
185 
196  template <int dim, int spacedim>
197  double
198  minimal_cell_diameter(
199  const Triangulation<dim, spacedim> &triangulation,
200  const Mapping<dim, spacedim> & mapping =
201  (ReferenceCells::get_hypercube<dim>()
202 #ifndef _MSC_VER
203  .template get_default_linear_mapping<dim, spacedim>()
204 #else
205  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
206 #endif
207  ));
208 
219  template <int dim, int spacedim>
220  double
221  maximal_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 
243  template <int dim>
244  DEAL_II_DEPRECATED double
245  cell_measure(
246  const std::vector<Point<dim>> &all_vertices,
247  const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]);
248 
268  template <int dim>
269  double
270  cell_measure(const std::vector<Point<dim>> & all_vertices,
271  const ArrayView<const unsigned int> &vertex_indices);
272 
295  template <int dim, int spacedim>
296  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
297  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices);
298 
325  template <int dim>
326  Vector<double>
327  compute_aspect_ratio_of_cells(const Mapping<dim> & mapping,
328  const Triangulation<dim> &triangulation,
329  const Quadrature<dim> & quadrature);
330 
338  template <int dim>
339  double
340  compute_maximum_aspect_ratio(const Mapping<dim> & mapping,
341  const Triangulation<dim> &triangulation,
342  const Quadrature<dim> & quadrature);
343 
357  template <int dim, int spacedim>
358  BoundingBox<spacedim>
359  compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);
360 
378  template <typename Iterator>
379  Point<Iterator::AccessorType::space_dimension>
380  project_to_object(
381  const Iterator & object,
382  const Point<Iterator::AccessorType::space_dimension> &trial_point);
383 
396  template <int dim, int spacedim>
397  std::
398  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
399  get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
400 
406 
423  template <int dim, int spacedim>
424  void
425  delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
426  std::vector<CellData<dim>> & cells,
427  SubCellData & subcelldata);
428 
446  template <int dim, int spacedim>
447  void
448  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
449  std::vector<CellData<dim>> & cells,
450  SubCellData & subcelldata,
451  std::vector<unsigned int> & considered_vertices,
452  const double tol = 1e-12);
453 
472  template <int dim, int spacedim>
473  void
474  invert_all_negative_measure_cells(
475  const std::vector<Point<spacedim>> &all_vertices,
476  std::vector<CellData<dim>> & cells);
477 
487  template <int dim, int spacedim>
488  std::size_t
489  invert_cells_with_negative_measure(
490  const std::vector<Point<spacedim>> &all_vertices,
491  std::vector<CellData<dim>> & cells);
492 
503  template <int dim>
504  void
505  consistently_order_cells(std::vector<CellData<dim>> &cells);
506 
512 
575  template <int dim, typename Transformation, int spacedim>
576  void
577  transform(const Transformation & transformation,
578  Triangulation<dim, spacedim> &triangulation);
579 
586  template <int dim, int spacedim>
587  void
588  shift(const Tensor<1, spacedim> & shift_vector,
589  Triangulation<dim, spacedim> &triangulation);
590 
591 
602  template <int dim>
603  void
604  rotate(const double angle, Triangulation<dim> &triangulation);
605 
618  template <int dim>
619  void
620  rotate(const Tensor<1, 3, double> &axis,
621  const double angle,
622  Triangulation<dim, 3> & triangulation);
623 
639  template <int dim>
640  DEAL_II_DEPRECATED_EARLY void
641  rotate(const double angle,
642  const unsigned int axis,
643  Triangulation<dim, 3> &triangulation);
644 
702  template <int dim>
703  void
704  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
705  Triangulation<dim> & tria,
706  const Function<dim, double> *coefficient = nullptr,
707  const bool solve_for_absolute_positions = false);
708 
714  template <int dim, int spacedim>
715  std::map<unsigned int, Point<spacedim>>
716  get_all_vertices_at_boundary(const Triangulation<dim, spacedim> &tria);
717 
725  template <int dim, int spacedim>
726  void
727  scale(const double scaling_factor,
728  Triangulation<dim, spacedim> &triangulation);
729 
744  template <int dim, int spacedim>
745  void
746  distort_random(
747  const double factor,
748  Triangulation<dim, spacedim> &triangulation,
749  const bool keep_boundary = true,
750  const unsigned int seed = boost::random::mt19937::default_seed);
751 
785  template <int dim, int spacedim>
786  void
787  remove_hanging_nodes(Triangulation<dim, spacedim> &tria,
788  const bool isotropic = false,
789  const unsigned int max_iterations = 100);
790 
815  template <int dim, int spacedim>
816  void
817  remove_anisotropy(Triangulation<dim, spacedim> &tria,
818  const double max_ratio = 1.6180339887,
819  const unsigned int max_iterations = 5);
820 
910  template <int dim, int spacedim>
911  void
912  regularize_corner_cells(Triangulation<dim, spacedim> &tria,
913  const double limit_angle_fraction = .75);
914 
920 
975  template <int dim, int spacedim>
976 #ifndef DOXYGEN
977  std::tuple<
978  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
979  std::vector<std::vector<Point<dim>>>,
980  std::vector<std::vector<unsigned int>>>
981 #else
982  return_type
983 #endif
984  compute_point_locations(
985  const Cache<dim, spacedim> & cache,
986  const std::vector<Point<spacedim>> &points,
987  const typename Triangulation<dim, spacedim>::active_cell_iterator
988  &cell_hint =
989  typename Triangulation<dim, spacedim>::active_cell_iterator());
990 
1024  template <int dim, int spacedim>
1025 #ifndef DOXYGEN
1026  std::tuple<
1027  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1028  std::vector<std::vector<Point<dim>>>,
1029  std::vector<std::vector<unsigned int>>,
1030  std::vector<unsigned int>>
1031 #else
1032  return_type
1033 #endif
1034  compute_point_locations_try_all(
1035  const Cache<dim, spacedim> & cache,
1036  const std::vector<Point<spacedim>> &points,
1037  const typename Triangulation<dim, spacedim>::active_cell_iterator
1038  &cell_hint =
1039  typename Triangulation<dim, spacedim>::active_cell_iterator());
1040 
1110  template <int dim, int spacedim>
1111 #ifndef DOXYGEN
1112  std::tuple<
1113  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1114  std::vector<std::vector<Point<dim>>>,
1115  std::vector<std::vector<unsigned int>>,
1116  std::vector<std::vector<Point<spacedim>>>,
1117  std::vector<std::vector<unsigned int>>>
1118 #else
1119  return_type
1120 #endif
1121  distributed_compute_point_locations(
1122  const GridTools::Cache<dim, spacedim> & cache,
1123  const std::vector<Point<spacedim>> & local_points,
1124  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1125  const double tolerance = 1e-10);
1126 
1127  namespace internal
1128  {
1143  template <int dim, int spacedim>
1144  struct DistributedComputePointLocationsInternal
1145  {
1152  std::vector<std::tuple<std::pair<int, int>,
1153  unsigned int,
1154  unsigned int,
1155  Point<dim>,
1156  Point<spacedim>,
1157  unsigned int>>
1158  send_components;
1159 
1163  std::vector<unsigned int> send_ranks;
1164 
1170  std::vector<unsigned int> send_ptrs;
1171 
1182  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1183  recv_components;
1184 
1188  std::vector<unsigned int> recv_ranks;
1189 
1195  std::vector<unsigned int> recv_ptrs;
1196  };
1197 
1207  template <int dim, int spacedim>
1208  DistributedComputePointLocationsInternal<dim, spacedim>
1209  distributed_compute_point_locations(
1210  const GridTools::Cache<dim, spacedim> & cache,
1211  const std::vector<Point<spacedim>> & points,
1212  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1213  const std::vector<bool> & marked_vertices,
1214  const double tolerance,
1215  const bool perform_handshake,
1216  const bool enforce_unique_mapping = false);
1217 
1218  } // namespace internal
1219 
1256  template <int dim, int spacedim>
1257  std::map<unsigned int, Point<spacedim>>
1258  extract_used_vertices(
1259  const Triangulation<dim, spacedim> &container,
1260  const Mapping<dim, spacedim> & mapping =
1261  (ReferenceCells::get_hypercube<dim>()
1262 #ifndef _MSC_VER
1263  .template get_default_linear_mapping<dim, spacedim>()
1264 #else
1265  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1266 #endif
1267  ));
1268 
1278  template <int spacedim>
1279  unsigned int
1280  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1281  const Point<spacedim> & p);
1282 
1306  template <int dim, template <int, int> class MeshType, int spacedim>
1307  unsigned int
1308  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1309  const Point<spacedim> & p,
1310  const std::vector<bool> & marked_vertices = {});
1311 
1335  template <int dim, template <int, int> class MeshType, int spacedim>
1336  unsigned int
1337  find_closest_vertex(const Mapping<dim, spacedim> & mapping,
1338  const MeshType<dim, spacedim> &mesh,
1339  const Point<spacedim> & p,
1340  const std::vector<bool> & marked_vertices = {});
1341 
1342 
1362  template <int dim, template <int, int> class MeshType, int spacedim>
1363 #ifndef _MSC_VER
1364  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1365 #else
1366  std::vector<
1367  typename ::internal::
1368  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1369 #endif
1370  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1371  const unsigned int vertex_index);
1372 
1435  template <int dim, template <int, int> class MeshType, int spacedim>
1436 #ifndef _MSC_VER
1437  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1438 #else
1439  std::pair<typename ::internal::
1440  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1441  Point<dim>>
1442 #endif
1443  find_active_cell_around_point(const Mapping<dim, spacedim> & mapping,
1444  const MeshType<dim, spacedim> &mesh,
1445  const Point<spacedim> & p,
1446  const std::vector<bool> &marked_vertices = {},
1447  const double tolerance = 1.e-10);
1448 
1456  template <int dim, template <int, int> class MeshType, int spacedim>
1457 #ifndef _MSC_VER
1458  typename MeshType<dim, spacedim>::active_cell_iterator
1459 #else
1460  typename ::internal::
1461  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1462 #endif
1463  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1464  const Point<spacedim> & p,
1465  const std::vector<bool> &marked_vertices = {},
1466  const double tolerance = 1.e-10);
1467 
1474  template <int dim, int spacedim>
1475  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1476  Point<dim>>
1477  find_active_cell_around_point(
1478  const hp::MappingCollection<dim, spacedim> &mapping,
1479  const DoFHandler<dim, spacedim> & mesh,
1480  const Point<spacedim> & p,
1481  const double tolerance = 1.e-10);
1482 
1534  template <int dim, int spacedim>
1535  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1536  Point<dim>>
1537  find_active_cell_around_point(
1538  const Cache<dim, spacedim> &cache,
1539  const Point<spacedim> & p,
1540  const typename Triangulation<dim, spacedim>::active_cell_iterator &
1541  cell_hint = typename Triangulation<dim, spacedim>::active_cell_iterator(),
1542  const std::vector<bool> &marked_vertices = {},
1543  const double tolerance = 1.e-10);
1544 
1558  template <int dim, template <int, int> class MeshType, int spacedim>
1559 #ifndef _MSC_VER
1560  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1561 #else
1562  std::pair<typename ::internal::
1563  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1564  Point<dim>>
1565 #endif
1566  find_active_cell_around_point(
1567  const Mapping<dim, spacedim> & mapping,
1568  const MeshType<dim, spacedim> &mesh,
1569  const Point<spacedim> & p,
1570  const std::vector<
1571  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1572  & vertex_to_cell_map,
1573  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1574  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1575  typename MeshType<dim, spacedim>::active_cell_iterator(),
1576  const std::vector<bool> & marked_vertices = {},
1577  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1578  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1579  const double tolerance = 1.e-10,
1580  const RTree<
1581  std::pair<BoundingBox<spacedim>,
1582  typename Triangulation<dim, spacedim>::active_cell_iterator>>
1583  *relevant_cell_bounding_boxes_rtree = nullptr);
1584 
1606  template <int dim, template <int, int> class MeshType, int spacedim>
1607 #ifndef _MSC_VER
1608  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1609  Point<dim>>>
1610 #else
1611  std::vector<std::pair<
1612  typename ::internal::
1613  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1614  Point<dim>>>
1615 #endif
1616  find_all_active_cells_around_point(
1617  const Mapping<dim, spacedim> & mapping,
1618  const MeshType<dim, spacedim> &mesh,
1619  const Point<spacedim> & p,
1620  const double tolerance,
1621  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1622  Point<dim>> & first_cell);
1623 
1630  template <int dim, template <int, int> class MeshType, int spacedim>
1631 #ifndef _MSC_VER
1632  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1633  Point<dim>>>
1634 #else
1635  std::vector<std::pair<
1636  typename ::internal::
1637  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1638  Point<dim>>>
1639 #endif
1640  find_all_active_cells_around_point(
1641  const Mapping<dim, spacedim> & mapping,
1642  const MeshType<dim, spacedim> &mesh,
1643  const Point<spacedim> & p,
1644  const double tolerance = 1e-10,
1645  const std::vector<bool> & marked_vertices = {});
1646 
1668  template <class MeshType>
1669  std::vector<typename MeshType::active_cell_iterator>
1670  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1671 
1696  template <class MeshType>
1697  void
1698  get_active_neighbors(
1699  const typename MeshType::active_cell_iterator & cell,
1700  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1701 
1751  template <class MeshType>
1752  std::vector<typename MeshType::active_cell_iterator>
1753  compute_active_cell_halo_layer(
1754  const MeshType &mesh,
1755  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1756  &predicate);
1757 
1758 
1766  template <class MeshType>
1767  std::vector<typename MeshType::cell_iterator>
1768  compute_cell_halo_layer_on_level(
1769  const MeshType &mesh,
1770  const std::function<bool(const typename MeshType::cell_iterator &)>
1771  & predicate,
1772  const unsigned int level);
1773 
1774 
1787  template <class MeshType>
1788  std::vector<typename MeshType::active_cell_iterator>
1789  compute_ghost_cell_halo_layer(const MeshType &mesh);
1790 
1839  template <class MeshType>
1840  std::vector<typename MeshType::active_cell_iterator>
1841  compute_active_cell_layer_within_distance(
1842  const MeshType &mesh,
1843  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1844  & predicate,
1845  const double layer_thickness);
1846 
1869  template <class MeshType>
1870  std::vector<typename MeshType::active_cell_iterator>
1871  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1872  const double layer_thickness);
1873 
1889  template <class MeshType>
1890  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1891  compute_bounding_box(
1892  const MeshType &mesh,
1893  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1894  &predicate);
1895 
1948  template <class MeshType>
1949  std::vector<BoundingBox<MeshType::space_dimension>>
1950  compute_mesh_predicate_bounding_box(
1951  const MeshType &mesh,
1952  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1953  & predicate,
1954  const unsigned int refinement_level = 0,
1955  const bool allow_merge = false,
1956  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1957 
1985  template <int spacedim>
1986 #ifndef DOXYGEN
1987  std::tuple<std::vector<std::vector<unsigned int>>,
1988  std::map<unsigned int, unsigned int>,
1989  std::map<unsigned int, std::vector<unsigned int>>>
1990 #else
1991  return_type
1992 #endif
1993  guess_point_owner(
1994  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1995  const std::vector<Point<spacedim>> & points);
1996 
1997 
2032  template <int spacedim>
2033 #ifndef DOXYGEN
2034  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2035  std::map<unsigned int, unsigned int>,
2036  std::map<unsigned int, std::vector<unsigned int>>>
2037 #else
2038  return_type
2039 #endif
2040  guess_point_owner(
2041  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
2042  const std::vector<Point<spacedim>> & points);
2043 
2044 
2053  template <int dim, int spacedim>
2054  std::vector<
2055  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2056  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
2057 
2070  template <int dim, int spacedim>
2071  std::vector<std::vector<Tensor<1, spacedim>>>
2072  vertex_to_cell_centers_directions(
2073  const Triangulation<dim, spacedim> &mesh,
2074  const std::vector<
2075  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2076  &vertex_to_cells);
2077 
2078 
2086  template <int dim, int spacedim>
2087  unsigned int
2088  find_closest_vertex_of_cell(
2089  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
2090  const Point<spacedim> & position,
2091  const Mapping<dim, spacedim> & mapping =
2092  (ReferenceCells::get_hypercube<dim>()
2093 #ifndef _MSC_VER
2094  .template get_default_linear_mapping<dim, spacedim>()
2095 #else
2096  .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2097 #endif
2098  ));
2099 
2111  template <int dim, int spacedim>
2112  std::map<unsigned int, types::global_vertex_index>
2113  compute_local_to_global_vertex_index_map(
2114  const parallel::distributed::Triangulation<dim, spacedim> &triangulation);
2115 
2127  template <int dim, int spacedim>
2128  std::pair<unsigned int, double>
2129  get_longest_direction(
2130  typename Triangulation<dim, spacedim>::active_cell_iterator cell);
2131 
2137 
2146  template <int dim, int spacedim>
2147  void
2148  get_face_connectivity_of_cells(
2149  const Triangulation<dim, spacedim> &triangulation,
2150  DynamicSparsityPattern & connectivity);
2151 
2160  template <int dim, int spacedim>
2161  void
2162  get_vertex_connectivity_of_cells(
2163  const Triangulation<dim, spacedim> &triangulation,
2164  DynamicSparsityPattern & connectivity);
2165 
2174  template <int dim, int spacedim>
2175  void
2176  get_vertex_connectivity_of_cells_on_level(
2177  const Triangulation<dim, spacedim> &triangulation,
2178  const unsigned int level,
2179  DynamicSparsityPattern & connectivity);
2180 
2201  template <int dim, int spacedim>
2202  void
2203  partition_triangulation(const unsigned int n_partitions,
2204  Triangulation<dim, spacedim> & triangulation,
2205  const SparsityTools::Partitioner partitioner =
2206  SparsityTools::Partitioner::metis);
2207 
2218  template <int dim, int spacedim>
2219  void
2220  partition_triangulation(const unsigned int n_partitions,
2221  const std::vector<unsigned int> &cell_weights,
2222  Triangulation<dim, spacedim> & triangulation,
2223  const SparsityTools::Partitioner partitioner =
2224  SparsityTools::Partitioner::metis);
2225 
2271  template <int dim, int spacedim>
2272  void
2273  partition_triangulation(const unsigned int n_partitions,
2274  const SparsityPattern & cell_connection_graph,
2275  Triangulation<dim, spacedim> &triangulation,
2276  const SparsityTools::Partitioner partitioner =
2277  SparsityTools::Partitioner::metis);
2278 
2289  template <int dim, int spacedim>
2290  void
2291  partition_triangulation(const unsigned int n_partitions,
2292  const std::vector<unsigned int> &cell_weights,
2293  const SparsityPattern & cell_connection_graph,
2294  Triangulation<dim, spacedim> &triangulation,
2295  const SparsityTools::Partitioner partitioner =
2296  SparsityTools::Partitioner::metis);
2297 
2312  template <int dim, int spacedim>
2313  void
2314  partition_triangulation_zorder(const unsigned int n_partitions,
2315  Triangulation<dim, spacedim> &triangulation,
2316  const bool group_siblings = true);
2317 
2329  template <int dim, int spacedim>
2330  void
2331  partition_multigrid_levels(Triangulation<dim, spacedim> &triangulation);
2332 
2340  template <int dim, int spacedim>
2341  std::vector<types::subdomain_id>
2342  get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
2343  const std::vector<CellId> & cell_ids);
2344 
2355  template <int dim, int spacedim>
2356  void
2357  get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
2358  std::vector<types::subdomain_id> & subdomain);
2359 
2374  template <int dim, int spacedim>
2375  unsigned int
2376  count_cells_with_subdomain_association(
2377  const Triangulation<dim, spacedim> &triangulation,
2378  const types::subdomain_id subdomain);
2379 
2409  template <int dim, int spacedim>
2410  std::vector<bool>
2411  get_locally_owned_vertices(const Triangulation<dim, spacedim> &triangulation);
2412 
2418 
2452  template <typename MeshType>
2453  std::list<std::pair<typename MeshType::cell_iterator,
2454  typename MeshType::cell_iterator>>
2455  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2456 
2466  template <int dim, int spacedim>
2467  bool
2468  have_same_coarse_mesh(const Triangulation<dim, spacedim> &mesh_1,
2469  const Triangulation<dim, spacedim> &mesh_2);
2470 
2480  template <typename MeshType>
2481  bool
2482  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2483 
2489 
2505  template <int dim, int spacedim>
2506  typename Triangulation<dim, spacedim>::DistortedCellList
2507  fix_up_distorted_child_cells(
2508  const typename Triangulation<dim, spacedim>::DistortedCellList
2509  & distorted_cells,
2510  Triangulation<dim, spacedim> &triangulation);
2511 
2512 
2513 
2522 
2523 
2561  template <class MeshType>
2562  std::vector<typename MeshType::active_cell_iterator>
2563  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2564 
2565 
2587  template <class Container>
2588  std::vector<typename Container::cell_iterator>
2589  get_cells_at_coarsest_common_level(
2590  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2591 
2658  template <class Container>
2659  void
2660  build_triangulation_from_patch(
2661  const std::vector<typename Container::active_cell_iterator> &patch,
2662  Triangulation<Container::dimension, Container::space_dimension>
2663  &local_triangulation,
2664  std::map<
2665  typename Triangulation<Container::dimension,
2666  Container::space_dimension>::active_cell_iterator,
2667  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2668 
2700  template <int dim, int spacedim>
2701  std::map<
2702  types::global_dof_index,
2703  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2704  get_dof_to_support_patch_map(DoFHandler<dim, spacedim> &dof_handler);
2705 
2706 
2713 
2719  template <typename CellIterator>
2720  struct PeriodicFacePair
2721  {
2725  CellIterator cell[2];
2726 
2731  unsigned int face_idx[2];
2732 
2738  std::bitset<3> orientation;
2739 
2752  FullMatrix<double> matrix;
2753  };
2754 
2755 
2819  template <typename FaceIterator>
2820  bool
2821  orthogonal_equality(
2822  std::bitset<3> & orientation,
2823  const FaceIterator & face1,
2824  const FaceIterator & face2,
2825  const unsigned int direction,
2826  const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset =
2827  Tensor<1, FaceIterator::AccessorType::space_dimension>(),
2828  const FullMatrix<double> &matrix = FullMatrix<double>());
2829 
2830 
2834  template <typename FaceIterator>
2835  bool
2836  orthogonal_equality(
2837  const FaceIterator & face1,
2838  const FaceIterator & face2,
2839  const unsigned int direction,
2840  const Tensor<1, FaceIterator::AccessorType::space_dimension> &offset =
2841  Tensor<1, FaceIterator::AccessorType::space_dimension>(),
2842  const FullMatrix<double> &matrix = FullMatrix<double>());
2843 
2844 
2901  template <typename MeshType>
2902  void
2903  collect_periodic_faces(
2904  const MeshType & mesh,
2905  const types::boundary_id b_id1,
2906  const types::boundary_id b_id2,
2907  const unsigned int direction,
2908  std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
2909  & matched_pairs,
2910  const Tensor<1, MeshType::space_dimension> &offset =
2911  ::Tensor<1, MeshType::space_dimension>(),
2912  const FullMatrix<double> &matrix = FullMatrix<double>());
2913 
2914 
2937  template <typename MeshType>
2938  void
2939  collect_periodic_faces(
2940  const MeshType & mesh,
2941  const types::boundary_id b_id,
2942  const unsigned int direction,
2943  std::vector<PeriodicFacePair<typename MeshType::cell_iterator>>
2944  & matched_pairs,
2945  const ::Tensor<1, MeshType::space_dimension> &offset =
2946  ::Tensor<1, MeshType::space_dimension>(),
2947  const FullMatrix<double> &matrix = FullMatrix<double>());
2948 
2954 
2975  template <int dim, int spacedim>
2976  void
2977  copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
2978  const bool reset_boundary_ids = false);
2979 
3001  template <int dim, int spacedim>
3002  void
3003  map_boundary_to_manifold_ids(
3004  const std::vector<types::boundary_id> &src_boundary_ids,
3005  const std::vector<types::manifold_id> &dst_manifold_ids,
3006  Triangulation<dim, spacedim> & tria,
3007  const std::vector<types::boundary_id> &reset_boundary_ids = {});
3008 
3038  template <int dim, int spacedim>
3039  void
3040  copy_material_to_manifold_id(Triangulation<dim, spacedim> &tria,
3041  const bool compute_face_ids = false);
3042 
3067  template <int dim, int spacedim>
3068  void
3069  assign_co_dimensional_manifold_indicators(
3070  Triangulation<dim, spacedim> & tria,
3071  const std::function<types::manifold_id(
3072  const std::set<types::manifold_id> &)> &disambiguation_function =
3073  [](const std::set<types::manifold_id> &manifold_ids) {
3074  if (manifold_ids.size() == 1)
3075  return *manifold_ids.begin();
3076  else
3077  return numbers::flat_manifold_id;
3078  },
3079  bool overwrite_only_flat_manifold_ids = true);
3166  template <typename DataType, typename MeshType>
3167  void
3168  exchange_cell_data_to_ghosts(
3169  const MeshType & mesh,
3170  const std::function<std_cxx17::optional<DataType>(
3171  const typename MeshType::active_cell_iterator &)> &pack,
3172  const std::function<void(const typename MeshType::active_cell_iterator &,
3173  const DataType &)> & unpack,
3174  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3175  &cell_filter =
3176  always_return<typename MeshType::active_cell_iterator, bool>{true});
3177 
3188  template <typename DataType, typename MeshType>
3189  void
3190  exchange_cell_data_to_level_ghosts(
3191  const MeshType & mesh,
3192  const std::function<std_cxx17::optional<DataType>(
3193  const typename MeshType::level_cell_iterator &)> &pack,
3194  const std::function<void(const typename MeshType::level_cell_iterator &,
3195  const DataType &)> & unpack,
3196  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3197  cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
3198  true});
3199 
3200  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3201  * boxes @p local_bboxes.
3202  *
3203  * This function is meant to exchange bounding boxes describing the locally
3204  * owned cells in a distributed triangulation obtained with the function
3205  * GridTools::compute_mesh_predicate_bounding_box .
3206  *
3207  * The output vector's size is the number of processes of the MPI
3208  * communicator:
3209  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3210  */
3211  template <int spacedim>
3212  std::vector<std::vector<BoundingBox<spacedim>>>
3213  exchange_local_bounding_boxes(
3214  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3215  const MPI_Comm & mpi_communicator);
3216 
3249  template <int spacedim>
3250  RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
3251  build_global_description_tree(
3252  const std::vector<BoundingBox<spacedim>> &local_description,
3253  const MPI_Comm & mpi_communicator);
3254 
3272  template <int dim, int spacedim>
3273  void
3274  collect_coinciding_vertices(
3275  const Triangulation<dim, spacedim> & tria,
3276  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3277  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3278 
3285  template <int dim, int spacedim>
3286  std::map<unsigned int, std::set<::types::subdomain_id>>
3287  compute_vertices_with_ghost_neighbors(
3288  const Triangulation<dim, spacedim> &tria);
3289 
3305  template <int dim, typename T>
3306  struct DEAL_II_DEPRECATED CellDataTransferBuffer
3307  {
3311  std::vector<CellId> cell_ids;
3312 
3316  std::vector<T> data;
3317 
3326  template <class Archive>
3327  void
3328  save(Archive &ar, const unsigned int) const
3329  {
3330  // Implement the code inline as some compilers do otherwise complain
3331  // about the use of a deprecated interface.
3332  Assert(cell_ids.size() == data.size(),
3333  ExcDimensionMismatch(cell_ids.size(), data.size()));
3334  // archive the cellids in an efficient binary format
3335  const std::size_t n_cells = cell_ids.size();
3336  ar & n_cells;
3337  for (const auto &id : cell_ids)
3338  {
3339  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3340  ar & binary_cell_id;
3341  }
3342 
3343  ar &data;
3344  }
3345 
3352  template <class Archive>
3353  void
3354  load(Archive &ar, const unsigned int)
3355  {
3356  // Implement the code inline as some compilers do otherwise complain
3357  // about the use of a deprecated interface.
3358  std::size_t n_cells;
3359  ar & n_cells;
3360  cell_ids.clear();
3361  cell_ids.reserve(n_cells);
3362  for (unsigned int c = 0; c < n_cells; ++c)
3363  {
3364  CellId::binary_type value;
3365  ar & value;
3366  cell_ids.emplace_back(value);
3367  }
3368  ar &data;
3369  }
3370 
3371 
3372 #ifdef DOXYGEN
3378  template <class Archive>
3379  void
3380  serialize(Archive &archive, const unsigned int version);
3381 #else
3382  // This macro defines the serialize() method that is compatible with
3383  // the templated save() and load() method that have been implemented.
3384  BOOST_SERIALIZATION_SPLIT_MEMBER()
3385 #endif
3386  };
3387 
3388 
3389 
3407  template <int dim, typename VectorType>
3408  class MarchingCubeAlgorithm
3409  {
3410  public:
3414  using value_type = typename VectorType::value_type;
3415 
3419  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3420  const FiniteElement<dim, dim> &fe,
3421  const unsigned int n_subdivisions = 1,
3422  const double tolerance = 1e-10);
3423 
3428  void
3429  process(const DoFHandler<dim> & background_dof_handler,
3430  const VectorType & ls_vector,
3431  const double iso_level,
3432  std::vector<Point<dim>> & vertices,
3433  std::vector<CellData<dim - 1>> &cells) const;
3434 
3441  void
3442  process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
3443  const VectorType & ls_vector,
3444  const double iso_level,
3445  std::vector<Point<dim>> & vertices,
3446  std::vector<CellData<dim - 1>> &cells) const;
3447 
3448  private:
3453  static Quadrature<dim>
3454  create_quadrature_rule(const unsigned int n_subdivisions);
3455 
3459  void
3460  process_cell(std::vector<value_type> & ls_values,
3461  const std::vector<Point<dim>> & points,
3462  const double iso_level,
3463  std::vector<Point<dim>> & vertices,
3464  std::vector<CellData<dim - 1>> &cells) const;
3465 
3472  void
3473  process_sub_cell(const std::vector<value_type> & ls_values,
3474  const std::vector<Point<2>> & points,
3475  const std::vector<unsigned int> mask,
3476  const double iso_level,
3477  std::vector<Point<2>> & vertices,
3478  std::vector<CellData<1>> & cells) const;
3479 
3483  void
3484  process_sub_cell(const std::vector<value_type> & ls_values,
3485  const std::vector<Point<3>> & points,
3486  const std::vector<unsigned int> mask,
3487  const double iso_level,
3488  std::vector<Point<3>> & vertices,
3489  std::vector<CellData<2>> & cells) const;
3490 
3495  const unsigned int n_subdivisions;
3496 
3501  const double tolerance;
3502 
3507  mutable FEValues<dim> fe_values;
3508  };
3509 
3510 
3511 
3516 
3520  DeclException1(ExcInvalidNumberOfPartitions,
3521  int,
3522  << "The number of partitions you gave is " << arg1
3523  << ", but must be greater than zero.");
3527  DeclException1(ExcNonExistentSubdomain,
3528  int,
3529  << "The subdomain id " << arg1
3530  << " has no cells associated with it.");
3534  DeclException0(ExcTriangulationHasBeenRefined);
3535 
3539  DeclException1(ExcScalingFactorNotPositive,
3540  double,
3541  << "The scaling factor must be positive, but it is " << arg1
3542  << '.');
3543 
3547  DeclException1(ExcVertexNotUsed,
3548  unsigned int,
3549  << "The given vertex with index " << arg1
3550  << " is not used in the given triangulation.");
3551 
3554 } /*namespace GridTools*/
3555 
3556 
3566 DeclExceptionMsg(ExcMeshNotOrientable,
3567  "The edges of the mesh are not consistently orientable.");
3568 
3569 
3570 
3571 /* ----------------- Template function --------------- */
3572 
3573 #ifndef DOXYGEN
3574 
3575 namespace GridTools
3576 {
3577  template <int dim>
3578  double
3579  cell_measure(
3580  const std::vector<Point<dim>> &all_vertices,
3581  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3582  {
3583  // We forward call to the ArrayView version:
3584  const ArrayView<const unsigned int> view(
3585  indices, GeometryInfo<dim>::vertices_per_cell);
3586  return cell_measure(all_vertices, view);
3587  }
3588 
3589 
3590 
3591  // This specialization is defined here so that the general template in the
3592  // source file doesn't need to have further 1D overloads for the internal
3593  // functions it calls.
3594  template <>
3595  inline Triangulation<1, 1>::DistortedCellList
3596  fix_up_distorted_child_cells(const Triangulation<1, 1>::DistortedCellList &,
3597  Triangulation<1, 1> &)
3598  {
3599  return {};
3600  }
3601 
3602 
3603 
3604  template <int dim, typename Predicate, int spacedim>
3605  void
3606  transform(const Predicate & predicate,
3607  Triangulation<dim, spacedim> &triangulation)
3608  {
3609  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3610 
3611  // loop over all active cells, and
3612  // transform those vertices that
3613  // have not yet been touched. note
3614  // that we get to all vertices in
3615  // the triangulation by only
3616  // visiting the active cells.
3617  typename Triangulation<dim, spacedim>::active_cell_iterator
3618  cell = triangulation.begin_active(),
3619  endc = triangulation.end();
3620  for (; cell != endc; ++cell)
3621  for (const unsigned int v : cell->vertex_indices())
3622  if (treated_vertices[cell->vertex_index(v)] == false)
3623  {
3624  // transform this vertex
3625  cell->vertex(v) = predicate(cell->vertex(v));
3626  // and mark it as treated
3627  treated_vertices[cell->vertex_index(v)] = true;
3628  };
3629 
3630 
3631  // now fix any vertices on hanging nodes so that we don't create any holes
3632  if (dim == 2)
3633  {
3634  typename Triangulation<dim, spacedim>::active_cell_iterator
3635  cell = triangulation.begin_active(),
3636  endc = triangulation.end();
3637  for (; cell != endc; ++cell)
3638  for (const unsigned int face : cell->face_indices())
3639  if (cell->face(face)->has_children() &&
3640  !cell->face(face)->at_boundary())
3641  {
3642  Assert(cell->reference_cell() ==
3643  ReferenceCells::get_hypercube<dim>(),
3644  ExcNotImplemented());
3645 
3646  // this line has children
3647  cell->face(face)->child(0)->vertex(1) =
3648  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3649  2;
3650  }
3651  }
3652  else if (dim == 3)
3653  {
3654  typename Triangulation<dim, spacedim>::active_cell_iterator
3655  cell = triangulation.begin_active(),
3656  endc = triangulation.end();
3657  for (; cell != endc; ++cell)
3658  for (const unsigned int face : cell->face_indices())
3659  if (cell->face(face)->has_children() &&
3660  !cell->face(face)->at_boundary())
3661  {
3662  Assert(cell->reference_cell() ==
3663  ReferenceCells::get_hypercube<dim>(),
3664  ExcNotImplemented());
3665 
3666  // this face has hanging nodes
3667  cell->face(face)->child(0)->vertex(1) =
3668  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3669  2.0;
3670  cell->face(face)->child(0)->vertex(2) =
3671  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3672  2.0;
3673  cell->face(face)->child(1)->vertex(3) =
3674  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3675  2.0;
3676  cell->face(face)->child(2)->vertex(3) =
3677  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3678  2.0;
3679 
3680  // center of the face
3681  cell->face(face)->child(0)->vertex(3) =
3682  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3683  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3684  4.0;
3685  }
3686  }
3687 
3688  // Make sure FEValues notices that the mesh has changed
3689  triangulation.signals.mesh_movement();
3690  }
3691 
3692 
3693 
3694  template <class MeshType>
3695  std::vector<typename MeshType::active_cell_iterator>
3696  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3697  {
3698  std::vector<typename MeshType::active_cell_iterator> child_cells;
3699 
3700  if (cell->has_children())
3701  {
3702  for (unsigned int child = 0; child < cell->n_children(); ++child)
3703  if (cell->child(child)->has_children())
3704  {
3705  const std::vector<typename MeshType::active_cell_iterator>
3706  children = get_active_child_cells<MeshType>(cell->child(child));
3707  child_cells.insert(child_cells.end(),
3708  children.begin(),
3709  children.end());
3710  }
3711  else
3712  child_cells.push_back(cell->child(child));
3713  }
3714 
3715  return child_cells;
3716  }
3717 
3718 
3719 
3720  template <class MeshType>
3721  void
3722  get_active_neighbors(
3723  const typename MeshType::active_cell_iterator & cell,
3724  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3725  {
3726  active_neighbors.clear();
3727  for (const unsigned int n : cell->face_indices())
3728  if (!cell->at_boundary(n))
3729  {
3730  if (MeshType::dimension == 1)
3731  {
3732  // check children of neighbor. note
3733  // that in 1d children of the neighbor
3734  // may be further refined. In 1d the
3735  // case is simple since we know what
3736  // children bound to the present cell
3737  typename MeshType::cell_iterator neighbor_child =
3738  cell->neighbor(n);
3739  if (!neighbor_child->is_active())
3740  {
3741  while (neighbor_child->has_children())
3742  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3743 
3744  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3745  ExcInternalError());
3746  }
3747  active_neighbors.push_back(neighbor_child);
3748  }
3749  else
3750  {
3751  if (cell->face(n)->has_children())
3752  // this neighbor has children. find
3753  // out which border to the present
3754  // cell
3755  for (unsigned int c = 0;
3756  c < cell->face(n)->n_active_descendants();
3757  ++c)
3758  active_neighbors.push_back(
3759  cell->neighbor_child_on_subface(n, c));
3760  else
3761  {
3762  // the neighbor must be active
3763  // himself
3764  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3765  active_neighbors.push_back(cell->neighbor(n));
3766  }
3767  }
3768  }
3769  }
3770 
3771 
3772 
3773  namespace internal
3774  {
3775  namespace ProjectToObject
3776  {
3789  struct CrossDerivative
3790  {
3791  const unsigned int direction_0;
3792  const unsigned int direction_1;
3793 
3794  CrossDerivative(const unsigned int d0, const unsigned int d1);
3795  };
3796 
3797  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3798  const unsigned int d1)
3799  : direction_0(d0)
3800  , direction_1(d1)
3801  {}
3802 
3803 
3804 
3809  template <typename F>
3810  inline auto
3811  centered_first_difference(const double center,
3812  const double step,
3813  const F &f) -> decltype(f(center) - f(center))
3814  {
3815  return (f(center + step) - f(center - step)) / (2.0 * step);
3816  }
3817 
3818 
3819 
3824  template <typename F>
3825  inline auto
3826  centered_second_difference(const double center,
3827  const double step,
3828  const F &f) -> decltype(f(center) - f(center))
3829  {
3830  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3831  (step * step);
3832  }
3833 
3834 
3835 
3845  template <int structdim, typename F>
3846  inline auto
3847  cross_stencil(
3848  const CrossDerivative cross_derivative,
3849  const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
3850  const double step,
3851  const F &f) -> decltype(f(center) - f(center))
3852  {
3853  Tensor<1, GeometryInfo<structdim>::vertices_per_cell> simplex_vector;
3854  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3855  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3856  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3857  1.0 / 3.0 * f(center - simplex_vector) +
3858  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3859  step;
3860  }
3861 
3862 
3863 
3870  template <int spacedim, int structdim, typename F>
3871  inline double
3872  gradient_entry(
3873  const unsigned int row_n,
3874  const unsigned int dependent_direction,
3875  const Point<spacedim> &p0,
3876  const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
3877  const double step,
3878  const F & f)
3879  {
3880  Assert(row_n < GeometryInfo<structdim>::vertices_per_cell &&
3881  dependent_direction <
3882  GeometryInfo<structdim>::vertices_per_cell,
3883  ExcMessage("This function assumes that the last weight is a "
3884  "dependent variable (and hence we cannot take its "
3885  "derivative directly)."));
3886  Assert(row_n != dependent_direction,
3887  ExcMessage(
3888  "We cannot differentiate with respect to the variable "
3889  "that is assumed to be dependent."));
3890 
3891  const Point<spacedim> manifold_point = f(center);
3892  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3893  {row_n, dependent_direction}, center, step, f);
3894  double entry = 0.0;
3895  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3896  entry +=
3897  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3898  return entry;
3899  }
3900 
3906  template <typename Iterator, int spacedim, int structdim>
3907  Point<spacedim>
3908  project_to_d_linear_object(const Iterator & object,
3909  const Point<spacedim> &trial_point)
3910  {
3911  // let's look at this for simplicity for a quad (structdim==2) in a
3912  // space with spacedim>2 (notate trial_point by y): all points on the
3913  // surface are given by
3914  // x(\xi) = sum_i v_i phi_x(\xi)
3915  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3916  // reference coordinates of the quad. so what we are trying to do is
3917  // find a point x on the surface that is closest to the point y. there
3918  // are different ways to solve this problem, but in the end it's a
3919  // nonlinear problem and we have to find reference coordinates \xi so
3920  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3921  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3922  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3923  // have to use a Newton method to find the answer. This leads to the
3924  // following formulation of Newton steps:
3925  //
3926  // Given \xi_k, find \delta\xi_k so that
3927  // H_k \delta\xi_k = - F_k
3928  // where H_k is an approximation to the second derivatives of J at
3929  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3930  // number of times until the right hand side is small enough. As a
3931  // stopping criterion, we terminate if ||\delta\xi||<eps.
3932  //
3933  // As for the Hessian, the best choice would be
3934  // H_k = J''(\xi_k)
3935  // but we'll opt for the simpler Gauss-Newton form
3936  // H_k = A^T A
3937  // i.e.
3938  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3939  // \partial_n phi_i *
3940  // \partial_m phi_j
3941  // we start at xi=(0.5, 0.5).
3942  Point<structdim> xi;
3943  for (unsigned int d = 0; d < structdim; ++d)
3944  xi[d] = 0.5;
3945 
3946  Point<spacedim> x_k;
3947  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3948  x_k += object->vertex(i) *
3949  GeometryInfo<structdim>::d_linear_shape_function(xi, i);
3950 
3951  do
3952  {
3953  Tensor<1, structdim> F_k;
3954  for (const unsigned int i :
3955  GeometryInfo<structdim>::vertex_indices())
3956  F_k +=
3957  (x_k - trial_point) * object->vertex(i) *
3958  GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
3959  i);
3960 
3961  Tensor<2, structdim> H_k;
3962  for (const unsigned int i :
3963  GeometryInfo<structdim>::vertex_indices())
3964  for (const unsigned int j :
3965  GeometryInfo<structdim>::vertex_indices())
3966  {
3967  Tensor<2, structdim> tmp = outer_product(
3968  GeometryInfo<structdim>::d_linear_shape_function_gradient(
3969  xi, i),
3970  GeometryInfo<structdim>::d_linear_shape_function_gradient(
3971  xi, j));
3972  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3973  }
3974 
3975  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3976  xi += delta_xi;
3977 
3978  x_k = Point<spacedim>();
3979  for (const unsigned int i :
3980  GeometryInfo<structdim>::vertex_indices())
3981  x_k += object->vertex(i) *
3982  GeometryInfo<structdim>::d_linear_shape_function(xi, i);
3983 
3984  if (delta_xi.norm() < 1e-7)
3985  break;
3986  }
3987  while (true);
3988 
3989  return x_k;
3990  }
3991  } // namespace ProjectToObject
3992  } // namespace internal
3993 
3994 
3995 
3996  namespace internal
3997  {
3998  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3999  // inside the project_to_object function below.
4000  template <int structdim>
4001  inline bool
4002  weights_are_ok(
4003  const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v)
4004  {
4005  // clang has trouble figuring out structdim here, so define it
4006  // again:
4007  static const std::size_t n_vertices_per_cell =
4008  Tensor<1, GeometryInfo<structdim>::vertices_per_cell>::
4009  n_independent_components;
4010  std::array<double, n_vertices_per_cell> copied_weights;
4011  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4012  {
4013  copied_weights[i] = v[i];
4014  if (v[i] < 0.0 || v[i] > 1.0)
4015  return false;
4016  }
4017 
4018  // check the sum: try to avoid some roundoff errors by summing in order
4019  std::sort(copied_weights.begin(), copied_weights.end());
4020  const double sum =
4021  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4022  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
4023  }
4024  } // namespace internal
4025 
4026  template <typename Iterator>
4027  Point<Iterator::AccessorType::space_dimension>
4028  project_to_object(
4029  const Iterator & object,
4030  const Point<Iterator::AccessorType::space_dimension> &trial_point)
4031  {
4032  const int spacedim = Iterator::AccessorType::space_dimension;
4033  const int structdim = Iterator::AccessorType::structure_dimension;
4034 
4035  Point<spacedim> projected_point = trial_point;
4036 
4037  if (structdim >= spacedim)
4038  return projected_point;
4039  else if (structdim == 1 || structdim == 2)
4040  {
4041  using namespace internal::ProjectToObject;
4042  // Try to use the special flat algorithm for quads (this is better
4043  // than the general algorithm in 3D). This does not take into account
4044  // whether projected_point is outside the quad, but we optimize along
4045  // lines below anyway:
4046  const int dim = Iterator::AccessorType::dimension;
4047  const Manifold<dim, spacedim> &manifold = object->get_manifold();
4048  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
4049  &manifold) != nullptr)
4050  {
4051  projected_point =
4052  project_to_d_linear_object<Iterator, spacedim, structdim>(
4053  object, trial_point);
4054  }
4055  else
4056  {
4057  // We want to find a point on the convex hull (defined by the
4058  // vertices of the object and the manifold description) that is
4059  // relatively close to the trial point. This has a few issues:
4060  //
4061  // 1. For a general convex hull we are not guaranteed that a unique
4062  // minimum exists.
4063  // 2. The independent variables in the optimization process are the
4064  // weights given to Manifold::get_new_point, which must sum to 1,
4065  // so we cannot use standard finite differences to approximate a
4066  // gradient.
4067  //
4068  // There is not much we can do about 1., but for 2. we can derive
4069  // finite difference stencils that work on a structdim-dimensional
4070  // simplex and rewrite the optimization problem to use those
4071  // instead. Consider the structdim 2 case and let
4072  //
4073  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
4074  // c2, c3})
4075  //
4076  // where {c0, c1, c2, c3} are the weights for the four vertices on
4077  // the quadrilateral. We seek to minimize the Euclidean distance
4078  // between F(...) and trial_point. We can solve for c3 in terms of
4079  // the other weights and get, for one coordinate direction
4080  //
4081  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
4082  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
4083  //
4084  // where we substitute back in for c3 after taking the
4085  // derivative. We can compute a stencil for the cross derivative
4086  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
4087  // (and gradient_entry computes the sum over the independent
4088  // variables). Below, we somewhat arbitrarily pick the last
4089  // component as the dependent one.
4090  //
4091  // Since we can now calculate derivatives of the objective
4092  // function we can use gradient descent to minimize it.
4093  //
4094  // Of course, this is much simpler in the structdim = 1 case (we
4095  // could rewrite the projection as a 1D optimization problem), but
4096  // to reduce the potential for bugs we use the same code in both
4097  // cases.
4098  const double step_size = object->diameter() / 64.0;
4099 
4100  constexpr unsigned int n_vertices_per_cell =
4101  GeometryInfo<structdim>::vertices_per_cell;
4102 
4103  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
4104  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4105  ++vertex_n)
4106  vertices[vertex_n] = object->vertex(vertex_n);
4107 
4108  auto get_point_from_weights =
4109  [&](const Tensor<1, n_vertices_per_cell> &weights)
4110  -> Point<spacedim> {
4111  return object->get_manifold().get_new_point(
4112  make_array_view(vertices.begin(), vertices.end()),
4113  make_array_view(weights.begin_raw(), weights.end_raw()));
4114  };
4115 
4116  // pick the initial weights as (normalized) inverse distances from
4117  // the trial point:
4118  Tensor<1, n_vertices_per_cell> guess_weights;
4119  double guess_weights_sum = 0.0;
4120  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4121  ++vertex_n)
4122  {
4123  const double distance =
4124  vertices[vertex_n].distance(trial_point);
4125  if (distance == 0.0)
4126  {
4127  guess_weights = 0.0;
4128  guess_weights[vertex_n] = 1.0;
4129  guess_weights_sum = 1.0;
4130  break;
4131  }
4132  else
4133  {
4134  guess_weights[vertex_n] = 1.0 / distance;
4135  guess_weights_sum += guess_weights[vertex_n];
4136  }
4137  }
4138  guess_weights /= guess_weights_sum;
4139  Assert(internal::weights_are_ok<structdim>(guess_weights),
4140  ExcInternalError());
4141 
4142  // The optimization algorithm consists of two parts:
4143  //
4144  // 1. An outer loop where we apply the gradient descent algorithm.
4145  // 2. An inner loop where we do a line search to find the optimal
4146  // length of the step one should take in the gradient direction.
4147  //
4148  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4149  {
4150  const unsigned int dependent_direction =
4151  n_vertices_per_cell - 1;
4152  Tensor<1, n_vertices_per_cell> current_gradient;
4153  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4154  ++row_n)
4155  {
4156  if (row_n != dependent_direction)
4157  {
4158  current_gradient[row_n] =
4159  gradient_entry<spacedim, structdim>(
4160  row_n,
4161  dependent_direction,
4162  trial_point,
4163  guess_weights,
4164  step_size,
4165  get_point_from_weights);
4166 
4167  current_gradient[dependent_direction] -=
4168  current_gradient[row_n];
4169  }
4170  }
4171 
4172  // We need to travel in the -gradient direction, as noted
4173  // above, but we may not want to take a full step in that
4174  // direction; instead, guess that we will go -0.5*gradient and
4175  // do quasi-Newton iteration to pick the best multiplier. The
4176  // goal is to find a scalar alpha such that
4177  //
4178  // F(x - alpha g)
4179  //
4180  // is minimized, where g is the gradient and F is the
4181  // objective function. To find the optimal value we find roots
4182  // of the derivative of the objective function with respect to
4183  // alpha by Newton iteration, where we approximate the first
4184  // and second derivatives of F(x - alpha g) with centered
4185  // finite differences.
4186  double gradient_weight = -0.5;
4187  auto gradient_weight_objective_function =
4188  [&](const double gradient_weight_guess) -> double {
4189  return (trial_point -
4190  get_point_from_weights(guess_weights +
4191  gradient_weight_guess *
4192  current_gradient))
4193  .norm_square();
4194  };
4195 
4196  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4197  {
4198  const double update_numerator = centered_first_difference(
4199  gradient_weight,
4200  step_size,
4201  gradient_weight_objective_function);
4202  const double update_denominator =
4203  centered_second_difference(
4204  gradient_weight,
4205  step_size,
4206  gradient_weight_objective_function);
4207 
4208  // avoid division by zero. Note that we limit the gradient
4209  // weight below
4210  if (std::abs(update_denominator) == 0.0)
4211  break;
4212  gradient_weight =
4213  gradient_weight - update_numerator / update_denominator;
4214 
4215  // Put a fairly lenient bound on the largest possible
4216  // gradient (things tend to be locally flat, so the gradient
4217  // itself is usually small)
4218  if (std::abs(gradient_weight) > 10)
4219  {
4220  gradient_weight = -10.0;
4221  break;
4222  }
4223  }
4224 
4225  // It only makes sense to take convex combinations with weights
4226  // between zero and one. If the update takes us outside of this
4227  // region then rescale the update to stay within the region and
4228  // try again
4229  Tensor<1, n_vertices_per_cell> tentative_weights =
4230  guess_weights + gradient_weight * current_gradient;
4231 
4232  double new_gradient_weight = gradient_weight;
4233  for (unsigned int iteration_count = 0; iteration_count < 40;
4234  ++iteration_count)
4235  {
4236  if (internal::weights_are_ok<structdim>(tentative_weights))
4237  break;
4238 
4239  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4240  {
4241  if (tentative_weights[i] < 0.0)
4242  {
4243  tentative_weights -=
4244  (tentative_weights[i] / current_gradient[i]) *
4245  current_gradient;
4246  }
4247  if (tentative_weights[i] < 0.0 ||
4248  1.0 < tentative_weights[i])
4249  {
4250  new_gradient_weight /= 2.0;
4251  tentative_weights =
4252  guess_weights +
4253  new_gradient_weight * current_gradient;
4254  }
4255  }
4256  }
4257 
4258  // the update might still send us outside the valid region, so
4259  // check again and quit if the update is still not valid
4260  if (!internal::weights_are_ok<structdim>(tentative_weights))
4261  break;
4262 
4263  // if we cannot get closer by traveling in the gradient
4264  // direction then quit
4265  if (get_point_from_weights(tentative_weights)
4266  .distance(trial_point) <
4267  get_point_from_weights(guess_weights).distance(trial_point))
4268  guess_weights = tentative_weights;
4269  else
4270  break;
4271  Assert(internal::weights_are_ok<structdim>(guess_weights),
4272  ExcInternalError());
4273  }
4274  Assert(internal::weights_are_ok<structdim>(guess_weights),
4275  ExcInternalError());
4276  projected_point = get_point_from_weights(guess_weights);
4277  }
4278 
4279  // if structdim == 2 and the optimal point is not on the interior then
4280  // we may be able to get a more accurate result by projecting onto the
4281  // lines.
4282  if (structdim == 2)
4283  {
4284  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4285  line_projections;
4286  for (unsigned int line_n = 0;
4287  line_n < GeometryInfo<structdim>::lines_per_cell;
4288  ++line_n)
4289  {
4290  line_projections[line_n] =
4291  project_to_object(object->line(line_n), trial_point);
4292  }
4293  std::sort(line_projections.begin(),
4294  line_projections.end(),
4295  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4296  return a.distance(trial_point) <
4297  b.distance(trial_point);
4298  });
4299  if (line_projections[0].distance(trial_point) <
4300  projected_point.distance(trial_point))
4301  projected_point = line_projections[0];
4302  }
4303  }
4304  else
4305  {
4306  Assert(false, ExcNotImplemented());
4307  return projected_point;
4308  }
4309 
4310  return projected_point;
4311  }
4312 
4313 
4314 
4315  namespace internal
4316  {
4317  template <typename DataType,
4318  typename MeshType,
4319  typename MeshCellIteratorType>
4320  inline void
4321  exchange_cell_data(
4322  const MeshType &mesh,
4323  const std::function<
4324  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4325  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4326  & unpack,
4327  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4328  const std::function<void(
4329  const std::function<void(const MeshCellIteratorType &,
4330  const types::subdomain_id)> &)> &process_cells,
4331  const std::function<std::set<types::subdomain_id>(
4332  const parallel::TriangulationBase<MeshType::dimension,
4333  MeshType::space_dimension> &)>
4334  &compute_ghost_owners)
4335  {
4336 # ifndef DEAL_II_WITH_MPI
4337  (void)mesh;
4338  (void)pack;
4339  (void)unpack;
4340  (void)cell_filter;
4341  (void)process_cells;
4342  (void)compute_ghost_owners;
4343  Assert(false, ExcNeedsMPI());
4344 # else
4345  constexpr int dim = MeshType::dimension;
4346  constexpr int spacedim = MeshType::space_dimension;
4347  auto tria =
4348  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4349  &mesh.get_triangulation());
4350  Assert(
4351  tria != nullptr,
4352  ExcMessage(
4353  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4354 
4355  if (const auto tria = dynamic_cast<
4356  const parallel::shared::Triangulation<dim, spacedim> *>(
4357  &mesh.get_triangulation()))
4358  {
4359  Assert(
4360  tria->with_artificial_cells(),
4361  ExcMessage(
4362  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4363  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4364  "operate on a single layer of ghost cells. However, you have "
4365  "given a Triangulation object of type "
4366  "parallel::shared::Triangulation without artificial cells "
4367  "resulting in arbitrary numbers of ghost layers."));
4368  }
4369 
4370  // build list of cells to request for each neighbor
4371  std::set<::types::subdomain_id> ghost_owners =
4372  compute_ghost_owners(*tria);
4373  std::map<::types::subdomain_id,
4374  std::vector<typename CellId::binary_type>>
4375  neighbor_cell_list;
4376 
4377  for (const auto ghost_owner : ghost_owners)
4378  neighbor_cell_list[ghost_owner] = {};
4379 
4380  process_cells([&](const auto &cell, const auto key) -> void {
4381  if (cell_filter(cell))
4382  {
4383  constexpr int spacedim = MeshType::space_dimension;
4384  neighbor_cell_list[key].emplace_back(
4385  cell->id().template to_binary<spacedim>());
4386  }
4387  });
4388 
4389  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4390  ExcInternalError());
4391 
4392 
4393  // Before sending & receiving, make sure we protect this section with
4394  // a mutex:
4395  static Utilities::MPI::CollectiveMutex mutex;
4396  Utilities::MPI::CollectiveMutex::ScopedLock lock(
4397  mutex, tria->get_communicator());
4398 
4399  const int mpi_tag =
4400  Utilities::MPI::internal::Tags::exchange_cell_data_request;
4401  const int mpi_tag_reply =
4402  Utilities::MPI::internal::Tags::exchange_cell_data_reply;
4403 
4404  // send our requests
4405  std::vector<MPI_Request> requests(ghost_owners.size());
4406  {
4407  unsigned int idx = 0;
4408  for (const auto &it : neighbor_cell_list)
4409  {
4410  // send the data about the relevant cells
4411  const int ierr = MPI_Isend(it.second.data(),
4412  it.second.size() * sizeof(it.second[0]),
4413  MPI_BYTE,
4414  it.first,
4415  mpi_tag,
4416  tria->get_communicator(),
4417  &requests[idx]);
4418  AssertThrowMPI(ierr);
4419  ++idx;
4420  }
4421  }
4422 
4423  // receive requests and reply with the results
4424  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4425  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4426 
4427  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4428  {
4429  MPI_Status status;
4430  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4431  mpi_tag,
4432  tria->get_communicator(),
4433  &status);
4434  AssertThrowMPI(ierr);
4435 
4436  int len;
4437  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4438  AssertThrowMPI(ierr);
4439  Assert(len % sizeof(typename CellId::binary_type) == 0,
4440  ExcInternalError());
4441 
4442  const unsigned int n_cells =
4443  len / sizeof(typename CellId::binary_type);
4444  std::vector<typename CellId::binary_type> cells_with_requests(
4445  n_cells);
4446  std::vector<DataType> data_to_send;
4447  data_to_send.reserve(n_cells);
4448  std::vector<bool> cell_carries_data(n_cells, false);
4449 
4450  ierr = MPI_Recv(cells_with_requests.data(),
4451  len,
4452  MPI_BYTE,
4453  status.MPI_SOURCE,
4454  status.MPI_TAG,
4455  tria->get_communicator(),
4456  &status);
4457  AssertThrowMPI(ierr);
4458 
4459  // store data for each cell
4460  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4461  {
4462  const auto cell =
4463  tria->create_cell_iterator(CellId(cells_with_requests[c]));
4464 
4465  MeshCellIteratorType mesh_it(tria,
4466  cell->level(),
4467  cell->index(),
4468  &mesh);
4469 
4470  const std_cxx17::optional<DataType> data = pack(mesh_it);
4471  if (data)
4472  {
4473  data_to_send.emplace_back(*data);
4474  cell_carries_data[c] = true;
4475  }
4476  }
4477 
4478  // collect data for sending the reply in a buffer
4479 
4480  // (a) make room for storing the local offsets in case we receive
4481  // other data
4482  sendbuffers[idx].resize(sizeof(std::size_t));
4483 
4484  // (b) append the actual data and store how much memory it
4485  // corresponds to, which we then insert into the leading position of
4486  // the sendbuffer
4487  std::size_t size_of_send =
4488  Utilities::pack(data_to_send,
4489  sendbuffers[idx],
4490  /*enable_compression*/ false);
4491  std::memcpy(sendbuffers[idx].data(),
4492  &size_of_send,
4493  sizeof(std::size_t));
4494 
4495  // (c) append information of certain cells that got left out in case
4496  // we need it
4497  if (data_to_send.size() < n_cells)
4498  Utilities::pack(cell_carries_data,
4499  sendbuffers[idx],
4500  /*enable_compression*/ false);
4501 
4502  // send data
4503  ierr = MPI_Isend(sendbuffers[idx].data(),
4504  sendbuffers[idx].size(),
4505  MPI_BYTE,
4506  status.MPI_SOURCE,
4507  mpi_tag_reply,
4508  tria->get_communicator(),
4509  &reply_requests[idx]);
4510  AssertThrowMPI(ierr);
4511  }
4512 
4513  // finally receive the replies
4514  std::vector<char> receive;
4515  for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
4516  {
4517  MPI_Status status;
4518  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4519  mpi_tag_reply,
4520  tria->get_communicator(),
4521  &status);
4522  AssertThrowMPI(ierr);
4523 
4524  int len;
4525  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4526  AssertThrowMPI(ierr);
4527 
4528  receive.resize(len);
4529 
4530  ierr = MPI_Recv(receive.data(),
4531  len,
4532  MPI_BYTE,
4533  status.MPI_SOURCE,
4534  status.MPI_TAG,
4535  tria->get_communicator(),
4536  &status);
4537  AssertThrowMPI(ierr);
4538 
4539  // (a) first determine the length of the data section in the
4540  // received buffer
4541  auto data_iterator = receive.begin();
4542  std::size_t size_of_received_data =
4543  Utilities::unpack<std::size_t>(data_iterator,
4544  data_iterator + sizeof(std::size_t));
4545  data_iterator += sizeof(std::size_t);
4546 
4547  // (b) unpack the data section in the indicated region
4548  auto received_data = Utilities::unpack<std::vector<DataType>>(
4549  data_iterator,
4550  data_iterator + size_of_received_data,
4551  /*enable_compression*/ false);
4552  data_iterator += size_of_received_data;
4553 
4554  // (c) check if the received data contained fewer entries than the
4555  // number of cells we identified in the beginning, in which case we
4556  // need to extract the boolean vector with the relevant information
4557  const std::vector<typename CellId::binary_type> &this_cell_list =
4558  neighbor_cell_list[status.MPI_SOURCE];
4559  AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
4560  std::vector<bool> cells_with_data;
4561  if (received_data.size() < this_cell_list.size())
4562  {
4563  cells_with_data = Utilities::unpack<std::vector<bool>>(
4564  data_iterator, receive.end(), /*enable_compression*/ false);
4565  AssertDimension(cells_with_data.size(), this_cell_list.size());
4566  }
4567 
4568  // (d) go through the received data and call the user-provided
4569  // unpack function
4570  auto received_data_iterator = received_data.begin();
4571  for (unsigned int c = 0; c < this_cell_list.size(); ++c)
4572  if (cells_with_data.empty() || cells_with_data[c])
4573  {
4574  const typename Triangulation<dim, spacedim>::cell_iterator
4575  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
4576 
4577  MeshCellIteratorType cell(tria,
4578  tria_cell->level(),
4579  tria_cell->index(),
4580  &mesh);
4581 
4582  unpack(cell, *received_data_iterator);
4583  ++received_data_iterator;
4584  }
4585  }
4586 
4587  // make sure that all communication is finished
4588  // when we leave this function.
4589  if (requests.size() > 0)
4590  {
4591  const int ierr =
4592  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4593  AssertThrowMPI(ierr);
4594  }
4595  if (reply_requests.size() > 0)
4596  {
4597  const int ierr = MPI_Waitall(reply_requests.size(),
4598  reply_requests.data(),
4599  MPI_STATUSES_IGNORE);
4600  AssertThrowMPI(ierr);
4601  }
4602 
4603 
4604 # endif // DEAL_II_WITH_MPI
4605  }
4606 
4607  } // namespace internal
4608 
4609  template <typename DataType, typename MeshType>
4610  inline void
4611  exchange_cell_data_to_ghosts(
4612  const MeshType & mesh,
4613  const std::function<std_cxx17::optional<DataType>(
4614  const typename MeshType::active_cell_iterator &)> &pack,
4615  const std::function<void(const typename MeshType::active_cell_iterator &,
4616  const DataType &)> & unpack,
4617  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4618  &cell_filter)
4619  {
4620 # ifndef DEAL_II_WITH_MPI
4621  (void)mesh;
4622  (void)pack;
4623  (void)unpack;
4624  (void)cell_filter;
4625  Assert(false, ExcNeedsMPI());
4626 # else
4627  internal::exchange_cell_data<DataType,
4628  MeshType,
4629  typename MeshType::active_cell_iterator>(
4630  mesh,
4631  pack,
4632  unpack,
4633  cell_filter,
4634  [&](const auto &process) {
4635  for (const auto &cell : mesh.active_cell_iterators())
4636  if (cell->is_ghost())
4637  process(cell, cell->subdomain_id());
4638  },
4639  [](const auto &tria) { return tria.ghost_owners(); });
4640 # endif
4641  }
4642 
4643 
4644 
4645  template <typename DataType, typename MeshType>
4646  inline void
4647  exchange_cell_data_to_level_ghosts(
4648  const MeshType & mesh,
4649  const std::function<std_cxx17::optional<DataType>(
4650  const typename MeshType::level_cell_iterator &)> &pack,
4651  const std::function<void(const typename MeshType::level_cell_iterator &,
4652  const DataType &)> & unpack,
4653  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4654  &cell_filter)
4655  {
4656 # ifndef DEAL_II_WITH_MPI
4657  (void)mesh;
4658  (void)pack;
4659  (void)unpack;
4660  (void)cell_filter;
4661  Assert(false, ExcNeedsMPI());
4662 # else
4663  internal::exchange_cell_data<DataType,
4664  MeshType,
4665  typename MeshType::level_cell_iterator>(
4666  mesh,
4667  pack,
4668  unpack,
4669  cell_filter,
4670  [&](const auto &process) {
4671  for (const auto &cell : mesh.cell_iterators())
4672  if (cell->level_subdomain_id() !=
4673  ::numbers::artificial_subdomain_id &&
4674  !cell->is_locally_owned_on_level())
4675  process(cell, cell->level_subdomain_id());
4676  },
4677  [](const auto &tria) { return tria.level_ghost_owners(); });
4678 # endif
4679  }
4680 } // namespace GridTools
4681 
4682 #endif // DOXYGEN
4683 
4684 DEAL_II_NAMESPACE_CLOSE
4685 
4686 #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:699
bounding_box.h
bounding_box.h
CellId
Definition: cell_id.h:71
CellId::binary_type
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
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 >> &cells) const
Definition: grid_tools.cc:6642
GridTools::MarchingCubeAlgorithm::value_type
typename VectorType::value_type value_type
Definition: grid_tools.h:3414
GridTools::MarchingCubeAlgorithm::n_subdivisions
const unsigned int n_subdivisions
Definition: grid_tools.h:3495
GridTools::MarchingCubeAlgorithm::MarchingCubeAlgorithm
MarchingCubeAlgorithm(const Mapping< dim, dim > &mapping, const FiniteElement< dim, dim > &fe, const unsigned int n_subdivisions=1, const double tolerance=1e-10)
Definition: grid_tools.cc:6579
GridTools::MarchingCubeAlgorithm::create_quadrature_rule
static Quadrature< dim > create_quadrature_rule(const unsigned int n_subdivisions)
Definition: grid_tools.cc:6596
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 >> &cells) const
Definition: grid_tools.cc:6626
GridTools::MarchingCubeAlgorithm::process_sub_cell
void process_sub_cell(const std::vector< value_type > &ls_values, const std::vector< Point< 2 >> &points, const std::vector< unsigned int > mask, const double iso_level, std::vector< Point< 2 >> &vertices, std::vector< CellData< 1 >> &cells) const
Definition: grid_tools.cc:6807
Point
Definition: point.h:111
Point::distance
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
SymmetricTensor::invert
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:3419
SymmetricTensor::outer_product
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition: symmetric_tensor.h:3470
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
Triangulation::begin
cell_iterator begin(const unsigned int level=0) const
Vector
Definition: vector.h:109
internal::ActiveCellIterator::type
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:106
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition: config.h:164
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
center
Point< 3 > center
Definition: data_out_base.cc:187
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:192
dof_handler.h
fe_values.h
geometry_info.h
level
unsigned int level
Definition: grid_out.cc:4606
vertex_indices
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1293
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
DeclException0
#define DeclException0(Exception0)
Definition: exceptions.h:464
ExcMeshNotOrientable
static ::ExceptionBase & ExcMeshNotOrientable()
GridTools::ExcScalingFactorNotPositive
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNeedsMPI
static ::ExceptionBase & ExcNeedsMPI()
GridTools::ExcVertexNotUsed
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
GridTools::ExcTriangulationHasBeenRefined
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
DeclExceptionMsg
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
GridTools::ExcInvalidNumberOfPartitions
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
DeclException1
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
GridTools::ExcNonExistentSubdomain
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
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:6398
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:4925
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:5017
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:4950
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:6521
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::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:6239
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::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:6335
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:5045
dof_handler.h
la_parallel_vector.h
la_vector.h
mapping_q1.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:5930
GridTools::get_face_connectivity_of_cells
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3732
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:289
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:3085
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2084
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:4072
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:2227
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:4883
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:1434
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:6194
GridTools::transform
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)
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:936
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:3003
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:436
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:741
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:6175
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:1105
GridTools::get_locally_owned_vertices
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4363
GridTools::cell_measure
double cell_measure(const std::vector< Point< dim >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< dim >::vertices_per_cell])
GridTools::project_to_object
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
GridTools::regularize_corner_cells
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:5227
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:5198
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:4390
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2050
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:737
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:2187
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:877
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:5572
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:5165
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:855
GridTools::partition_multigrid_levels
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4178
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:2738
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:869
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:3232
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:4347
GridTools::consistently_order_cells
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1959
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:5543
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:1697
GridTools::rotate
void rotate(const double angle, Triangulation< dim > &triangulation)
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:774
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:3831
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:313
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:3330
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:139
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:1021
GridTools::delete_unused_vertices
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:636
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:379
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:1210
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:1387
GridTools::get_vertex_connectivity_of_cells
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3766
GridTools::get_active_child_cells
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
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:4277
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:2259
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:1480
GridTools::diameter
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:83
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:3795
GridTools::find_cells_adjacent_to_vertex
std::vector< typename MeshType< dim, spacedim >::active_cell_iterator > find_cells_adjacent_to_vertex(const MeshType< dim, spacedim > &container, const unsigned int vertex_index)
Definition: grid_tools.cc:2610
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:3379
GridTools::compute_bounding_box
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:395
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:4417
GridTools::get_longest_direction
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5133
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:598
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:535
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:2413
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)
Definition: grid_tools.cc:5740
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:58
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::VectorRelations::angle
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
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::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:1483
Utilities::unpack
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1648
hp
Definition: hp.h:118
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13734
numbers::artificial_subdomain_id
const types::subdomain_id artificial_subdomain_id
Definition: types.h:298
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
numbers::flat_manifold_id
const types::manifold_id flat_manifold_id
Definition: types.h:269
types::manifold_id
unsigned int manifold_id
Definition: types.h:141
types::global_dof_index
unsigned int global_dof_index
Definition: types.h:76
types::subdomain_id
unsigned int subdomain_id
Definition: types.h:43
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5572
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:146
shared_tria.h
sparsity_tools.h
GeometryInfo::d_linear_shape_function
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
GeometryInfo::d_linear_shape_function_gradient
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
GridTools::CellDataTransferBuffer::load
void load(Archive &ar, const unsigned int)
Definition: grid_tools.h:3354
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:3328
GridTools::CellDataTransferBuffer::cell_ids
std::vector< CellId > cell_ids
Definition: grid_tools.h:3311
GridTools::PeriodicFacePair::orientation
std::bitset< 3 > orientation
Definition: grid_tools.h:2738
GridTools::PeriodicFacePair::matrix
FullMatrix< double > matrix
Definition: grid_tools.h:2752
GridTools::PeriodicFacePair::face_idx
unsigned int face_idx[2]
Definition: grid_tools.h:2731
GridTools::internal::DistributedComputePointLocationsInternal::recv_components
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1183
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:1158
always_return
Definition: template_constraints.h:358
tria_accessor.h
vertex_to_coinciding_vertex_group
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
Definition: tria_description.cc:731
coinciding_vertex_groups
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
Definition: tria_description.cc:730
tria
const ::Triangulation< dim, spacedim > & tria
Definition: tria_description.cc:717
tria_iterator.h
trilinos_vector.h
Generated by doxygen 1.9.1