Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
grid_tools.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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 
25 # include <deal.II/boost_adaptors/bounding_box.h>
26 
27 # include <deal.II/dofs/dof_handler.h>
28 
29 # include <deal.II/fe/mapping.h>
30 # include <deal.II/fe/mapping_q1.h>
31 
32 # include <deal.II/grid/manifold.h>
33 # include <deal.II/grid/tria.h>
34 # include <deal.II/grid/tria_accessor.h>
35 # include <deal.II/grid/tria_iterator.h>
36 
37 # include <deal.II/hp/dof_handler.h>
38 
39 # include <deal.II/lac/sparsity_tools.h>
40 
41 # include <deal.II/numerics/rtree.h>
42 
43 # include <boost/archive/binary_iarchive.hpp>
44 # include <boost/archive/binary_oarchive.hpp>
45 # include <boost/geometry/index/detail/serialization.hpp>
46 # include <boost/geometry/index/rtree.hpp>
47 # include <boost/optional.hpp>
48 # include <boost/serialization/array.hpp>
49 # include <boost/serialization/vector.hpp>
50 
51 # ifdef DEAL_II_WITH_ZLIB
52 # include <boost/iostreams/device/back_inserter.hpp>
53 # include <boost/iostreams/filter/gzip.hpp>
54 # include <boost/iostreams/filtering_stream.hpp>
55 # include <boost/iostreams/stream.hpp>
56 # endif
57 
58 
59 # include <bitset>
60 # include <list>
61 # include <set>
62 
63 DEAL_II_NAMESPACE_OPEN
64 
65 namespace parallel
66 {
67  namespace distributed
68  {
69  template <int, int>
70  class Triangulation;
71  }
72 } // namespace parallel
73 
74 namespace hp
75 {
76  template <int, int>
77  class MappingCollection;
78 }
79 
80 class SparsityPattern;
81 
82 namespace internal
83 {
84  template <int dim, int spacedim, class MeshType>
85  class ActiveCellIterator
86  {
87  public:
88 # ifndef _MSC_VER
89  using type = typename MeshType::active_cell_iterator;
90 # else
92 # endif
93  };
94 
95 # ifdef _MSC_VER
96  template <int dim, int spacedim>
97  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
98  {
99  public:
100  using type = TriaActiveIterator<
102  };
103 
104  template <int dim, int spacedim>
105  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim>>
106  {
107  public:
108  using type = TriaActiveIterator<
110  };
111 # endif
112 } // namespace internal
113 
122 namespace GridTools
123 {
124  template <int dim, int spacedim>
125  class Cache;
126 
131 
138  template <int dim, int spacedim>
139  double
141 
168  template <int dim, int spacedim>
169  double
171  const Mapping<dim, spacedim> & mapping =
173 
184  template <int dim, int spacedim>
185  double
187  const Mapping<dim, spacedim> & mapping =
189 
200  template <int dim, int spacedim>
201  double
203  const Mapping<dim, spacedim> & mapping =
205 
215  template <int dim>
216  double
217  cell_measure(
218  const std::vector<Point<dim>> &all_vertices,
219  const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]);
220 
226  template <int dim, typename T>
227  double
228  cell_measure(const T &, ...);
229 
243  template <int dim, int spacedim>
246 
266  template <typename Iterator>
269  const Iterator & object,
271 
280  template <int dim, int spacedim>
281  std::
282  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
284 
290 
307  template <int dim, int spacedim>
308  void
309  delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
310  std::vector<CellData<dim>> & cells,
311  SubCellData & subcelldata);
312 
329  template <int dim, int spacedim>
330  void
331  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
332  std::vector<CellData<dim>> & cells,
333  SubCellData & subcelldata,
334  std::vector<unsigned int> & considered_vertices,
335  const double tol = 1e-12);
336 
342 
377  template <int dim, typename Transformation, int spacedim>
378  void
379  transform(const Transformation & transformation,
380  Triangulation<dim, spacedim> &triangulation);
381 
387  template <int dim, int spacedim>
388  void
389  shift(const Tensor<1, spacedim> & shift_vector,
390  Triangulation<dim, spacedim> &triangulation);
391 
392 
400  void
401  rotate(const double angle, Triangulation<2> &triangulation);
402 
415  template <int dim>
416  void
417  rotate(const double angle,
418  const unsigned int axis,
419  Triangulation<dim, 3> &triangulation);
420 
478  template <int dim>
479  void
480  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
481  Triangulation<dim> & tria,
482  const Function<dim, double> *coefficient = nullptr,
483  const bool solve_for_absolute_positions = false);
484 
490  template <int dim, int spacedim>
491  std::map<unsigned int, Point<spacedim>>
493 
501  template <int dim, int spacedim>
502  void
503  scale(const double scaling_factor,
504  Triangulation<dim, spacedim> &triangulation);
505 
516  template <int dim, int spacedim>
517  void
518  distort_random(const double factor,
519  Triangulation<dim, spacedim> &triangulation,
520  const bool keep_boundary = true);
521 
557  template <int dim, int spacedim>
558  void
560  const bool isotropic = false,
561  const unsigned int max_iterations = 100);
562 
589  template <int dim, int spacedim>
590  void
592  const double max_ratio = 1.6180339887,
593  const unsigned int max_iterations = 5);
594 
686  template <int dim, int spacedim>
687  void
689  const double limit_angle_fraction = .75);
690 
696 
748  template <int dim, int spacedim>
749 # ifndef DOXYGEN
750  std::tuple<
751  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
752  std::vector<std::vector<Point<dim>>>,
753  std::vector<std::vector<unsigned int>>>
754 # else
755  return_type
756 # endif
758  const Cache<dim, spacedim> & cache,
759  const std::vector<Point<spacedim>> &points,
761  &cell_hint =
763 
792  template <int dim, int spacedim>
793 # ifndef DOXYGEN
794  std::tuple<
795  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
796  std::vector<std::vector<Point<dim>>>,
797  std::vector<std::vector<unsigned int>>,
798  std::vector<unsigned int>>
799 # else
800  return_type
801 # endif
803  const Cache<dim, spacedim> & cache,
804  const std::vector<Point<spacedim>> &points,
806  &cell_hint =
808 
873  template <int dim, int spacedim>
874 # ifndef DOXYGEN
875  std::tuple<
876  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
877  std::vector<std::vector<Point<dim>>>,
878  std::vector<std::vector<unsigned int>>,
879  std::vector<std::vector<Point<spacedim>>>,
880  std::vector<std::vector<unsigned int>>>
881 # else
882  return_type
883 # endif
885  const GridTools::Cache<dim, spacedim> & cache,
886  const std::vector<Point<spacedim>> & local_points,
887  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes);
888 
927  template <int dim, int spacedim>
928  std::map<unsigned int, Point<spacedim>>
930  const Mapping<dim, spacedim> & mapping =
932 
944  template <int spacedim>
945  unsigned int
946  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
947  const Point<spacedim> & p);
948 
975  template <int dim, template <int, int> class MeshType, int spacedim>
976  unsigned int
977  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
978  const Point<spacedim> & p,
979  const std::vector<bool> & marked_vertices = {});
980 
1006  template <int dim, template <int, int> class MeshType, int spacedim>
1007  unsigned int
1009  const MeshType<dim, spacedim> &mesh,
1010  const Point<spacedim> & p,
1011  const std::vector<bool> & marked_vertices = {});
1012 
1013 
1038  template <int dim, template <int, int> class MeshType, int spacedim>
1039 # ifndef _MSC_VER
1040  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1041 # else
1042  std::vector<
1043  typename ::internal::
1044  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1045 # endif
1046  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1047  const unsigned int vertex_index);
1048 
1049 
1074  template <int dim, template <int, int> class MeshType, int spacedim>
1075 # ifndef _MSC_VER
1076  typename MeshType<dim, spacedim>::active_cell_iterator
1077 # else
1078  typename ::internal::ActiveCellIterator<dim,
1079  spacedim,
1080  MeshType<dim, spacedim>>::type
1081 # endif
1082  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1083  const Point<spacedim> & p,
1084  const std::vector<bool> &marked_vertices = {});
1085 
1173  template <int dim, template <int, int> class MeshType, int spacedim>
1174 # ifndef _MSC_VER
1175  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1176 # else
1177  std::pair<typename ::internal::
1178  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1179  Point<dim>>
1180 # endif
1182  const MeshType<dim, spacedim> &mesh,
1183  const Point<spacedim> & p,
1184  const std::vector<bool> &marked_vertices = {});
1185 
1197  template <int dim, template <int, int> class MeshType, int spacedim>
1198 # ifndef _MSC_VER
1199  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1200 # else
1201  std::pair<typename ::internal::
1202  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1203  Point<dim>>
1204 # endif
1206  const Mapping<dim, spacedim> & mapping,
1207  const MeshType<dim, spacedim> &mesh,
1208  const Point<spacedim> & p,
1209  const std::vector<
1210  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1211  & vertex_to_cell_map,
1212  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1213  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1214  typename MeshType<dim, spacedim>::active_cell_iterator(),
1215  const std::vector<bool> & marked_vertices = {},
1216  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1217  RTree<std::pair<Point<spacedim>, unsigned int>>{});
1218 
1240  template <int dim, int spacedim>
1241  std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator,
1242  Point<dim>>
1244  const hp::MappingCollection<dim, spacedim> &mapping,
1245  const hp::DoFHandler<dim, spacedim> & mesh,
1246  const Point<spacedim> & p);
1247 
1254  template <int dim, int spacedim>
1255  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1256  Point<dim>>
1258  const Cache<dim, spacedim> &cache,
1259  const Point<spacedim> & p,
1262  const std::vector<bool> &marked_vertices = {});
1263 
1281  template <int dim, template <int, int> class MeshType, int spacedim>
1282 # ifndef _MSC_VER
1283  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1284  Point<dim>>>
1285 # else
1286  std::vector<std::pair<
1287  typename ::internal::
1288  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1289  Point<dim>>>
1290 # endif
1292  const Mapping<dim, spacedim> & mapping,
1293  const MeshType<dim, spacedim> &mesh,
1294  const Point<spacedim> & p,
1295  const double tolerance = 1e-12,
1296  const std::vector<bool> & marked_vertices = {});
1297 
1319  template <class MeshType>
1320  std::vector<typename MeshType::active_cell_iterator>
1321  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1322 
1340  template <class MeshType>
1341  void
1343  const typename MeshType::active_cell_iterator & cell,
1344  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1345 
1398  template <class MeshType>
1399  std::vector<typename MeshType::active_cell_iterator>
1401  const MeshType &mesh,
1402  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1403  &predicate);
1404 
1405 
1413  template <class MeshType>
1414  std::vector<typename MeshType::cell_iterator>
1416  const MeshType &mesh,
1417  const std::function<bool(const typename MeshType::cell_iterator &)>
1418  & predicate,
1419  const unsigned int level);
1420 
1421 
1437  template <class MeshType>
1438  std::vector<typename MeshType::active_cell_iterator>
1439  compute_ghost_cell_halo_layer(const MeshType &mesh);
1440 
1492  template <class MeshType>
1493  std::vector<typename MeshType::active_cell_iterator>
1495  const MeshType &mesh,
1496  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1497  & predicate,
1498  const double layer_thickness);
1499 
1525  template <class MeshType>
1526  std::vector<typename MeshType::active_cell_iterator>
1527  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1528  const double layer_thickness);
1529 
1545  template <class MeshType>
1546  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1548  const MeshType &mesh,
1549  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1550  &predicate);
1551 
1604  template <class MeshType>
1605  std::vector<BoundingBox<MeshType::space_dimension>>
1607  const MeshType &mesh,
1608  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1609  & predicate,
1610  const unsigned int refinement_level = 0,
1611  const bool allow_merge = false,
1612  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1613 
1643  template <int spacedim>
1644 # ifndef DOXYGEN
1645  std::tuple<std::vector<std::vector<unsigned int>>,
1646  std::map<unsigned int, unsigned int>,
1647  std::map<unsigned int, std::vector<unsigned int>>>
1648 # else
1649  return_type
1650 # endif
1652  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1653  const std::vector<Point<spacedim>> & points);
1654 
1655 
1690  template <int spacedim>
1691 # ifndef DOXYGEN
1692  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1693  std::map<unsigned int, unsigned int>,
1694  std::map<unsigned int, std::vector<unsigned int>>>
1695 # else
1696  return_type
1697 # endif
1699  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1700  const std::vector<Point<spacedim>> & points);
1701 
1702 
1711  template <int dim, int spacedim>
1712  std::vector<
1713  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1714  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1715 
1730  template <int dim, int spacedim>
1731  std::vector<std::vector<Tensor<1, spacedim>>>
1733  const Triangulation<dim, spacedim> &mesh,
1734  const std::vector<
1736  &vertex_to_cells);
1737 
1738 
1745  template <int dim, int spacedim>
1746  unsigned int
1749  const Point<spacedim> &position);
1750 
1762  template <int dim, int spacedim>
1763  std::map<unsigned int, types::global_vertex_index>
1766 
1780  template <int dim, int spacedim>
1781  std::pair<unsigned int, double>
1784 
1790 
1799  template <int dim, int spacedim>
1800  void
1802  const Triangulation<dim, spacedim> &triangulation,
1803  DynamicSparsityPattern & connectivity);
1804 
1813  template <int dim, int spacedim>
1814  void
1816  const Triangulation<dim, spacedim> &triangulation,
1817  DynamicSparsityPattern & connectivity);
1818 
1827  template <int dim, int spacedim>
1828  void
1830  const Triangulation<dim, spacedim> &triangulation,
1831  const unsigned int level,
1832  DynamicSparsityPattern & connectivity);
1833 
1854  template <int dim, int spacedim>
1855  void
1856  partition_triangulation(const unsigned int n_partitions,
1857  Triangulation<dim, spacedim> & triangulation,
1858  const SparsityTools::Partitioner partitioner =
1860 
1871  template <int dim, int spacedim>
1872  void
1873  partition_triangulation(const unsigned int n_partitions,
1874  const std::vector<unsigned int> &cell_weights,
1875  Triangulation<dim, spacedim> & triangulation,
1876  const SparsityTools::Partitioner partitioner =
1878 
1924  template <int dim, int spacedim>
1925  void
1926  partition_triangulation(const unsigned int n_partitions,
1927  const SparsityPattern & cell_connection_graph,
1928  Triangulation<dim, spacedim> &triangulation,
1929  const SparsityTools::Partitioner partitioner =
1931 
1942  template <int dim, int spacedim>
1943  void
1944  partition_triangulation(const unsigned int n_partitions,
1945  const std::vector<unsigned int> &cell_weights,
1946  const SparsityPattern & cell_connection_graph,
1947  Triangulation<dim, spacedim> &triangulation,
1948  const SparsityTools::Partitioner partitioner =
1950 
1958  template <int dim, int spacedim>
1959  void
1960  partition_triangulation_zorder(const unsigned int n_partitions,
1961  Triangulation<dim, spacedim> &triangulation);
1962 
1974  template <int dim, int spacedim>
1975  void
1977 
1988  template <int dim, int spacedim>
1989  void
1990  get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
1991  std::vector<types::subdomain_id> & subdomain);
1992 
2007  template <int dim, int spacedim>
2008  unsigned int
2010  const Triangulation<dim, spacedim> &triangulation,
2011  const types::subdomain_id subdomain);
2012 
2013 
2043  template <int dim, int spacedim>
2044  std::vector<bool>
2046 
2052 
2086  template <typename MeshType>
2087  std::list<std::pair<typename MeshType::cell_iterator,
2088  typename MeshType::cell_iterator>>
2089  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2090 
2100  template <int dim, int spacedim>
2101  bool
2103  const Triangulation<dim, spacedim> &mesh_2);
2104 
2114  template <typename MeshType>
2115  bool
2116  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2117 
2123 
2139  template <int dim, int spacedim>
2143  & distorted_cells,
2144  Triangulation<dim, spacedim> &triangulation);
2145 
2146 
2147 
2156 
2157 
2197  template <class MeshType>
2198  std::vector<typename MeshType::active_cell_iterator>
2199  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2200 
2201 
2225  template <class Container>
2226  std::vector<typename Container::cell_iterator>
2228  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2229 
2298  template <class Container>
2299  void
2301  const std::vector<typename Container::active_cell_iterator> &patch,
2303  &local_triangulation,
2304  std::map<
2305  typename Triangulation<Container::dimension,
2306  Container::space_dimension>::active_cell_iterator,
2307  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2308 
2345  template <class DoFHandlerType>
2346  std::map<types::global_dof_index,
2347  std::vector<typename DoFHandlerType::active_cell_iterator>>
2348  get_dof_to_support_patch_map(DoFHandlerType &dof_handler);
2349 
2350 
2357 
2363  template <typename CellIterator>
2364  struct PeriodicFacePair
2365  {
2369  CellIterator cell[2];
2370 
2375  unsigned int face_idx[2];
2376 
2382  std::bitset<3> orientation;
2383 
2397  };
2398 
2399 
2465  template <typename FaceIterator>
2466  bool orthogonal_equality(
2467  std::bitset<3> & orientation,
2468  const FaceIterator & face1,
2469  const FaceIterator & face2,
2470  const int direction,
2473  const FullMatrix<double> &matrix = FullMatrix<double>());
2474 
2475 
2479  template <typename FaceIterator>
2480  bool
2482  const FaceIterator & face1,
2483  const FaceIterator & face2,
2484  const int direction,
2487  const FullMatrix<double> &matrix = FullMatrix<double>());
2488 
2489 
2548  template <typename MeshType>
2549  void
2551  const MeshType & mesh,
2552  const types::boundary_id b_id1,
2553  const types::boundary_id b_id2,
2554  const int direction,
2556  & matched_pairs,
2557  const Tensor<1, MeshType::space_dimension> &offset =
2559  const FullMatrix<double> &matrix = FullMatrix<double>());
2560 
2561 
2584  template <typename MeshType>
2585  void
2587  const MeshType & mesh,
2588  const types::boundary_id b_id,
2589  const int direction,
2591  & matched_pairs,
2592  const ::Tensor<1, MeshType::space_dimension> &offset =
2594  const FullMatrix<double> &matrix = FullMatrix<double>());
2595 
2601 
2624  template <int dim, int spacedim>
2625  void
2627  const bool reset_boundary_ids = false);
2628 
2652  template <int dim, int spacedim>
2653  void
2655  const std::vector<types::boundary_id> &src_boundary_ids,
2656  const std::vector<types::manifold_id> &dst_manifold_ids,
2658  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2659 
2691  template <int dim, int spacedim>
2692  void
2694  const bool compute_face_ids = false);
2695 
2722  template <int dim, int spacedim>
2723  void
2726  const std::function<types::manifold_id(
2727  const std::set<types::manifold_id> &)> &disambiguation_function =
2728  [](const std::set<types::manifold_id> &manifold_ids) {
2729  if (manifold_ids.size() == 1)
2730  return *manifold_ids.begin();
2731  else
2733  },
2734  bool overwrite_only_flat_manifold_ids = true);
2819  template <typename DataType, typename MeshType>
2820  void
2822  const MeshType & mesh,
2823  const std::function<boost::optional<DataType>(
2824  const typename MeshType::active_cell_iterator &)> &pack,
2825  const std::function<void(const typename MeshType::active_cell_iterator &,
2826  const DataType &)> & unpack);
2827 
2828  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2829  * boxes @p local_bboxes.
2830  *
2831  * This function is meant to exchange bounding boxes describing the locally
2832  * owned cells in a distributed triangulation obtained with the function
2833  * GridTools::compute_mesh_predicate_bounding_box .
2834  *
2835  * The output vector's size is the number of processes of the MPI
2836  * communicator:
2837  * its i-th entry contains the vector @p local_bboxes of the i-th process.
2838  */
2839  template <int spacedim>
2840  std::vector<std::vector<BoundingBox<spacedim>>>
2841  exchange_local_bounding_boxes(
2842  const std::vector<BoundingBox<spacedim>> &local_bboxes,
2843  MPI_Comm mpi_communicator);
2844 
2879  template <int spacedim>
2880  RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
2882  const std::vector<BoundingBox<spacedim>> &local_description,
2883  MPI_Comm mpi_communicator);
2884 
2897  template <int dim, typename T>
2899  {
2903  std::vector<CellId> cell_ids;
2904 
2908  std::vector<T> data;
2909 
2917  template <class Archive>
2918  void
2919  save(Archive &ar, const unsigned int version) const;
2920 
2925  template <class Archive>
2926  void
2927  load(Archive &ar, const unsigned int version);
2928 
2929  BOOST_SERIALIZATION_SPLIT_MEMBER()
2930  };
2931 
2936 
2941  int,
2942  << "The number of partitions you gave is " << arg1
2943  << ", but must be greater than zero.");
2948  int,
2949  << "The subdomain id " << arg1
2950  << " has no cells associated with it.");
2955 
2960  double,
2961  << "The scaling factor must be positive, but it is " << arg1
2962  << ".");
2966  template <int N>
2968  Point<N>,
2969  << "The point <" << arg1
2970  << "> could not be found inside any of the "
2971  << "coarse grid cells.");
2975  template <int N>
2977  Point<N>,
2978  << "The point <" << arg1
2979  << "> could not be found inside any of the "
2980  << "subcells of a coarse grid cell.");
2981 
2986  unsigned int,
2987  << "The given vertex with index " << arg1
2988  << " is not used in the given triangulation.");
2989 
2990 
2993 } /*namespace GridTools*/
2994 
2995 
2996 
2997 /* ----------------- Template function --------------- */
2998 
2999 # ifndef DOXYGEN
3000 
3001 namespace GridTools
3002 {
3003  template <int dim, typename T>
3004  double
3005  cell_measure(const T &, ...)
3006  {
3007  Assert(false, ExcNotImplemented());
3008  return std::numeric_limits<double>::quiet_NaN();
3009  }
3010 
3011  template <int dim, typename Predicate, int spacedim>
3012  void
3013  transform(const Predicate & predicate,
3014  Triangulation<dim, spacedim> &triangulation)
3015  {
3016  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3017 
3018  // loop over all active cells, and
3019  // transform those vertices that
3020  // have not yet been touched. note
3021  // that we get to all vertices in
3022  // the triangulation by only
3023  // visiting the active cells.
3025  cell = triangulation.begin_active(),
3026  endc = triangulation.end();
3027  for (; cell != endc; ++cell)
3028  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
3029  if (treated_vertices[cell->vertex_index(v)] == false)
3030  {
3031  // transform this vertex
3032  cell->vertex(v) = predicate(cell->vertex(v));
3033  // and mark it as treated
3034  treated_vertices[cell->vertex_index(v)] = true;
3035  };
3036 
3037 
3038  // now fix any vertices on hanging nodes so that we don't create any holes
3039  if (dim == 2)
3040  {
3042  cell = triangulation.begin_active(),
3043  endc = triangulation.end();
3044  for (; cell != endc; ++cell)
3045  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
3046  ++face)
3047  if (cell->face(face)->has_children() &&
3048  !cell->face(face)->at_boundary())
3049  {
3050  // this line has children
3051  cell->face(face)->child(0)->vertex(1) =
3052  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3053  2;
3054  }
3055  }
3056  else if (dim == 3)
3057  {
3059  cell = triangulation.begin_active(),
3060  endc = triangulation.end();
3061  for (; cell != endc; ++cell)
3062  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
3063  ++face)
3064  if (cell->face(face)->has_children() &&
3065  !cell->face(face)->at_boundary())
3066  {
3067  // this face has hanging nodes
3068  cell->face(face)->child(0)->vertex(1) =
3069  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3070  2.0;
3071  cell->face(face)->child(0)->vertex(2) =
3072  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3073  2.0;
3074  cell->face(face)->child(1)->vertex(3) =
3075  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3076  2.0;
3077  cell->face(face)->child(2)->vertex(3) =
3078  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3079  2.0;
3080 
3081  // center of the face
3082  cell->face(face)->child(0)->vertex(3) =
3083  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3084  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3085  4.0;
3086  }
3087  }
3088 
3089  // Make sure FEValues notices that the mesh has changed
3090  triangulation.signals.mesh_movement();
3091  }
3092 
3093 
3094 
3095  template <class MeshType>
3096  std::vector<typename MeshType::active_cell_iterator>
3097  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3098  {
3099  std::vector<typename MeshType::active_cell_iterator> child_cells;
3100 
3101  if (cell->has_children())
3102  {
3103  for (unsigned int child = 0; child < cell->n_children(); ++child)
3104  if (cell->child(child)->has_children())
3105  {
3106  const std::vector<typename MeshType::active_cell_iterator>
3107  children = get_active_child_cells<MeshType>(cell->child(child));
3108  child_cells.insert(child_cells.end(),
3109  children.begin(),
3110  children.end());
3111  }
3112  else
3113  child_cells.push_back(cell->child(child));
3114  }
3115 
3116  return child_cells;
3117  }
3118 
3119 
3120 
3121  template <class MeshType>
3122  void
3124  const typename MeshType::active_cell_iterator & cell,
3125  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3126  {
3127  active_neighbors.clear();
3128  for (unsigned int n = 0;
3129  n < GeometryInfo<MeshType::dimension>::faces_per_cell;
3130  ++n)
3131  if (!cell->at_boundary(n))
3132  {
3133  if (MeshType::dimension == 1)
3134  {
3135  // check children of neighbor. note
3136  // that in 1d children of the neighbor
3137  // may be further refined. In 1d the
3138  // case is simple since we know what
3139  // children bound to the present cell
3140  typename MeshType::cell_iterator neighbor_child =
3141  cell->neighbor(n);
3142  if (!neighbor_child->active())
3143  {
3144  while (neighbor_child->has_children())
3145  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3146 
3147  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3148  ExcInternalError());
3149  }
3150  active_neighbors.push_back(neighbor_child);
3151  }
3152  else
3153  {
3154  if (cell->face(n)->has_children())
3155  // this neighbor has children. find
3156  // out which border to the present
3157  // cell
3158  for (unsigned int c = 0;
3159  c < cell->face(n)->number_of_children();
3160  ++c)
3161  active_neighbors.push_back(
3162  cell->neighbor_child_on_subface(n, c));
3163  else
3164  {
3165  // the neighbor must be active
3166  // himself
3167  Assert(cell->neighbor(n)->active(), ExcInternalError());
3168  active_neighbors.push_back(cell->neighbor(n));
3169  }
3170  }
3171  }
3172  }
3173 
3174 
3175 
3176  namespace internal
3177  {
3178  namespace ProjectToObject
3179  {
3192  struct CrossDerivative
3193  {
3194  const unsigned int direction_0;
3195  const unsigned int direction_1;
3196 
3197  CrossDerivative(const unsigned int d0, const unsigned int d1);
3198  };
3199 
3200  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3201  const unsigned int d1)
3202  : direction_0(d0)
3203  , direction_1(d1)
3204  {}
3205 
3206 
3207 
3212  template <typename F>
3213  inline auto
3214  centered_first_difference(const double center,
3215  const double step,
3216  const F &f) -> decltype(f(center) - f(center))
3217  {
3218  return (f(center + step) - f(center - step)) / (2.0 * step);
3219  }
3220 
3221 
3222 
3227  template <typename F>
3228  inline auto
3229  centered_second_difference(const double center,
3230  const double step,
3231  const F &f) -> decltype(f(center) - f(center))
3232  {
3233  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3234  (step * step);
3235  }
3236 
3237 
3238 
3248  template <int structdim, typename F>
3249  inline auto
3250  cross_stencil(
3251  const CrossDerivative cross_derivative,
3253  const double step,
3254  const F &f) -> decltype(f(center) - f(center))
3255  {
3257  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3258  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3259  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3260  1.0 / 3.0 * f(center - simplex_vector) +
3261  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3262  step;
3263  }
3264 
3265 
3266 
3273  template <int spacedim, int structdim, typename F>
3274  inline double
3275  gradient_entry(
3276  const unsigned int row_n,
3277  const unsigned int dependent_direction,
3278  const Point<spacedim> &p0,
3280  const double step,
3281  const F & f)
3282  {
3284  dependent_direction <
3286  ExcMessage("This function assumes that the last weight is a "
3287  "dependent variable (and hence we cannot take its "
3288  "derivative directly)."));
3289  Assert(row_n != dependent_direction,
3290  ExcMessage(
3291  "We cannot differentiate with respect to the variable "
3292  "that is assumed to be dependent."));
3293 
3294  const Point<spacedim> manifold_point = f(center);
3295  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3296  {row_n, dependent_direction}, center, step, f);
3297  double entry = 0.0;
3298  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3299  entry +=
3300  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3301  return entry;
3302  }
3303 
3309  template <typename Iterator, int spacedim, int structdim>
3311  project_to_d_linear_object(const Iterator & object,
3312  const Point<spacedim> &trial_point)
3313  {
3314  // let's look at this for simplicity for a quad (structdim==2) in a
3315  // space with spacedim>2 (notate trial_point by y): all points on the
3316  // surface are given by
3317  // x(\xi) = sum_i v_i phi_x(\xi)
3318  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3319  // reference coordinates of the quad. so what we are trying to do is
3320  // find a point x on the surface that is closest to the point y. there
3321  // are different ways to solve this problem, but in the end it's a
3322  // nonlinear problem and we have to find reference coordinates \xi so
3323  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3324  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3325  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3326  // have to use a Newton method to find the answer. This leads to the
3327  // following formulation of Newton steps:
3328  //
3329  // Given \xi_k, find \delta\xi_k so that
3330  // H_k \delta\xi_k = - F_k
3331  // where H_k is an approximation to the second derivatives of J at
3332  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3333  // number of times until the right hand side is small enough. As a
3334  // stopping criterion, we terminate if ||\delta\xi||<eps.
3335  //
3336  // As for the Hessian, the best choice would be
3337  // H_k = J''(\xi_k)
3338  // but we'll opt for the simpler Gauss-Newton form
3339  // H_k = A^T A
3340  // i.e.
3341  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3342  // \partial_n phi_i *
3343  // \partial_m phi_j
3344  // we start at xi=(0.5, 0.5).
3345  Point<structdim> xi;
3346  for (unsigned int d = 0; d < structdim; ++d)
3347  xi[d] = 0.5;
3348 
3349  Point<spacedim> x_k;
3350  for (unsigned int i = 0; i < GeometryInfo<structdim>::vertices_per_cell;
3351  ++i)
3352  x_k += object->vertex(i) *
3354 
3355  do
3356  {
3358  for (unsigned int i = 0;
3359  i < GeometryInfo<structdim>::vertices_per_cell;
3360  ++i)
3361  F_k +=
3362  (x_k - trial_point) * object->vertex(i) *
3364  i);
3365 
3367  for (unsigned int i = 0;
3368  i < GeometryInfo<structdim>::vertices_per_cell;
3369  ++i)
3370  for (unsigned int j = 0;
3371  j < GeometryInfo<structdim>::vertices_per_cell;
3372  ++j)
3373  {
3376  xi, i),
3378  xi, j));
3379  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3380  }
3381 
3382  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3383  xi += delta_xi;
3384 
3385  x_k = Point<spacedim>();
3386  for (unsigned int i = 0;
3387  i < GeometryInfo<structdim>::vertices_per_cell;
3388  ++i)
3389  x_k += object->vertex(i) *
3391 
3392  if (delta_xi.norm() < 1e-7)
3393  break;
3394  }
3395  while (true);
3396 
3397  return x_k;
3398  }
3399  } // namespace ProjectToObject
3400  } // namespace internal
3401 
3402 
3403 
3404  namespace internal
3405  {
3406  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3407  // inside the project_to_object function below.
3408  template <int structdim>
3409  inline bool
3410  weights_are_ok(
3412  {
3413  // clang has trouble figuring out structdim here, so define it
3414  // again:
3415  static const std::size_t n_vertices_per_cell =
3417  n_independent_components;
3418  std::array<double, n_vertices_per_cell> copied_weights;
3419  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3420  {
3421  copied_weights[i] = v[i];
3422  if (v[i] < 0.0 || v[i] > 1.0)
3423  return false;
3424  }
3425 
3426  // check the sum: try to avoid some roundoff errors by summing in order
3427  std::sort(copied_weights.begin(), copied_weights.end());
3428  const double sum =
3429  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3430  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3431  }
3432  } // namespace internal
3433 
3434  template <typename Iterator>
3437  const Iterator & object,
3439  {
3440  const int spacedim = Iterator::AccessorType::space_dimension;
3441  const int structdim = Iterator::AccessorType::structure_dimension;
3442 
3443  Point<spacedim> projected_point = trial_point;
3444 
3445  if (structdim >= spacedim)
3446  return projected_point;
3447  else if (structdim == 1 || structdim == 2)
3448  {
3449  using namespace internal::ProjectToObject;
3450  // Try to use the special flat algorithm for quads (this is better
3451  // than the general algorithm in 3D). This does not take into account
3452  // whether projected_point is outside the quad, but we optimize along
3453  // lines below anyway:
3454  const int dim = Iterator::AccessorType::dimension;
3455  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3456  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3457  &manifold) != nullptr)
3458  {
3459  projected_point =
3460  project_to_d_linear_object<Iterator, spacedim, structdim>(
3461  object, trial_point);
3462  }
3463  else
3464  {
3465  // We want to find a point on the convex hull (defined by the
3466  // vertices of the object and the manifold description) that is
3467  // relatively close to the trial point. This has a few issues:
3468  //
3469  // 1. For a general convex hull we are not guaranteed that a unique
3470  // minimum exists.
3471  // 2. The independent variables in the optimization process are the
3472  // weights given to Manifold::get_new_point, which must sum to 1,
3473  // so we cannot use standard finite differences to approximate a
3474  // gradient.
3475  //
3476  // There is not much we can do about 1., but for 2. we can derive
3477  // finite difference stencils that work on a structdim-dimensional
3478  // simplex and rewrite the optimization problem to use those
3479  // instead. Consider the structdim 2 case and let
3480  //
3481  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3482  // c2, c3})
3483  //
3484  // where {c0, c1, c2, c3} are the weights for the four vertices on
3485  // the quadrilateral. We seek to minimize the Euclidean distance
3486  // between F(...) and trial_point. We can solve for c3 in terms of
3487  // the other weights and get, for one coordinate direction
3488  //
3489  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3490  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3491  //
3492  // where we substitute back in for c3 after taking the
3493  // derivative. We can compute a stencil for the cross derivative
3494  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3495  // (and gradient_entry computes the sum over the independent
3496  // variables). Below, we somewhat arbitrarily pick the last
3497  // component as the dependent one.
3498  //
3499  // Since we can now calculate derivatives of the objective
3500  // function we can use gradient descent to minimize it.
3501  //
3502  // Of course, this is much simpler in the structdim = 1 case (we
3503  // could rewrite the projection as a 1D optimization problem), but
3504  // to reduce the potential for bugs we use the same code in both
3505  // cases.
3506  const double step_size = object->diameter() / 64.0;
3507 
3508  constexpr unsigned int n_vertices_per_cell =
3510 
3511  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3512  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3513  ++vertex_n)
3514  vertices[vertex_n] = object->vertex(vertex_n);
3515 
3516  auto get_point_from_weights =
3517  [&](const Tensor<1, n_vertices_per_cell> &weights)
3518  -> Point<spacedim> {
3519  return object->get_manifold().get_new_point(
3520  make_array_view(vertices.begin(), vertices.end()),
3521  make_array_view(weights.begin_raw(), weights.end_raw()));
3522  };
3523 
3524  // pick the initial weights as (normalized) inverse distances from
3525  // the trial point:
3526  Tensor<1, n_vertices_per_cell> guess_weights;
3527  double guess_weights_sum = 0.0;
3528  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3529  ++vertex_n)
3530  {
3531  const double distance =
3532  vertices[vertex_n].distance(trial_point);
3533  if (distance == 0.0)
3534  {
3535  guess_weights = 0.0;
3536  guess_weights[vertex_n] = 1.0;
3537  guess_weights_sum = 1.0;
3538  break;
3539  }
3540  else
3541  {
3542  guess_weights[vertex_n] = 1.0 / distance;
3543  guess_weights_sum += guess_weights[vertex_n];
3544  }
3545  }
3546  guess_weights /= guess_weights_sum;
3547  Assert(internal::weights_are_ok<structdim>(guess_weights),
3548  ExcInternalError());
3549 
3550  // The optimization algorithm consists of two parts:
3551  //
3552  // 1. An outer loop where we apply the gradient descent algorithm.
3553  // 2. An inner loop where we do a line search to find the optimal
3554  // length of the step one should take in the gradient direction.
3555  //
3556  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3557  {
3558  const unsigned int dependent_direction =
3559  n_vertices_per_cell - 1;
3560  Tensor<1, n_vertices_per_cell> current_gradient;
3561  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3562  ++row_n)
3563  {
3564  if (row_n != dependent_direction)
3565  {
3566  current_gradient[row_n] =
3567  gradient_entry<spacedim, structdim>(
3568  row_n,
3569  dependent_direction,
3570  trial_point,
3571  guess_weights,
3572  step_size,
3573  get_point_from_weights);
3574 
3575  current_gradient[dependent_direction] -=
3576  current_gradient[row_n];
3577  }
3578  }
3579 
3580  // We need to travel in the -gradient direction, as noted
3581  // above, but we may not want to take a full step in that
3582  // direction; instead, guess that we will go -0.5*gradient and
3583  // do quasi-Newton iteration to pick the best multiplier. The
3584  // goal is to find a scalar alpha such that
3585  //
3586  // F(x - alpha g)
3587  //
3588  // is minimized, where g is the gradient and F is the
3589  // objective function. To find the optimal value we find roots
3590  // of the derivative of the objective function with respect to
3591  // alpha by Newton iteration, where we approximate the first
3592  // and second derivatives of F(x - alpha g) with centered
3593  // finite differences.
3594  double gradient_weight = -0.5;
3595  auto gradient_weight_objective_function =
3596  [&](const double gradient_weight_guess) -> double {
3597  return (trial_point -
3598  get_point_from_weights(guess_weights +
3599  gradient_weight_guess *
3600  current_gradient))
3601  .norm_square();
3602  };
3603 
3604  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3605  {
3606  const double update_numerator = centered_first_difference(
3607  gradient_weight,
3608  step_size,
3609  gradient_weight_objective_function);
3610  const double update_denominator =
3611  centered_second_difference(
3612  gradient_weight,
3613  step_size,
3614  gradient_weight_objective_function);
3615 
3616  // avoid division by zero. Note that we limit the gradient
3617  // weight below
3618  if (std::abs(update_denominator) == 0.0)
3619  break;
3620  gradient_weight =
3621  gradient_weight - update_numerator / update_denominator;
3622 
3623  // Put a fairly lenient bound on the largest possible
3624  // gradient (things tend to be locally flat, so the gradient
3625  // itself is usually small)
3626  if (std::abs(gradient_weight) > 10)
3627  {
3628  gradient_weight = -10.0;
3629  break;
3630  }
3631  }
3632 
3633  // It only makes sense to take convex combinations with weights
3634  // between zero and one. If the update takes us outside of this
3635  // region then rescale the update to stay within the region and
3636  // try again
3637  Tensor<1, n_vertices_per_cell> tentative_weights =
3638  guess_weights + gradient_weight * current_gradient;
3639 
3640  double new_gradient_weight = gradient_weight;
3641  for (unsigned int iteration_count = 0; iteration_count < 40;
3642  ++iteration_count)
3643  {
3644  if (internal::weights_are_ok<structdim>(tentative_weights))
3645  break;
3646 
3647  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3648  {
3649  if (tentative_weights[i] < 0.0)
3650  {
3651  tentative_weights -=
3652  (tentative_weights[i] / current_gradient[i]) *
3653  current_gradient;
3654  }
3655  if (tentative_weights[i] < 0.0 ||
3656  1.0 < tentative_weights[i])
3657  {
3658  new_gradient_weight /= 2.0;
3659  tentative_weights =
3660  guess_weights +
3661  new_gradient_weight * current_gradient;
3662  }
3663  }
3664  }
3665 
3666  // the update might still send us outside the valid region, so
3667  // check again and quit if the update is still not valid
3668  if (!internal::weights_are_ok<structdim>(tentative_weights))
3669  break;
3670 
3671  // if we cannot get closer by traveling in the gradient
3672  // direction then quit
3673  if (get_point_from_weights(tentative_weights)
3674  .distance(trial_point) <
3675  get_point_from_weights(guess_weights).distance(trial_point))
3676  guess_weights = tentative_weights;
3677  else
3678  break;
3679  Assert(internal::weights_are_ok<structdim>(guess_weights),
3680  ExcInternalError());
3681  }
3682  Assert(internal::weights_are_ok<structdim>(guess_weights),
3683  ExcInternalError());
3684  projected_point = get_point_from_weights(guess_weights);
3685  }
3686 
3687  // if structdim == 2 and the optimal point is not on the interior then
3688  // we may be able to get a more accurate result by projecting onto the
3689  // lines.
3690  if (structdim == 2)
3691  {
3692  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
3693  line_projections;
3694  for (unsigned int line_n = 0;
3695  line_n < GeometryInfo<structdim>::lines_per_cell;
3696  ++line_n)
3697  {
3698  line_projections[line_n] =
3699  project_to_object(object->line(line_n), trial_point);
3700  }
3701  std::sort(line_projections.begin(),
3702  line_projections.end(),
3703  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
3704  return a.distance(trial_point) <
3705  b.distance(trial_point);
3706  });
3707  if (line_projections[0].distance(trial_point) <
3708  projected_point.distance(trial_point))
3709  projected_point = line_projections[0];
3710  }
3711  }
3712  else
3713  {
3714  Assert(false, ExcNotImplemented());
3715  return projected_point;
3716  }
3717 
3718  return projected_point;
3719  }
3720 
3721 
3722 
3723  template <int dim, typename T>
3724  template <class Archive>
3725  void
3727  const unsigned int /*version*/) const
3728  {
3729  Assert(cell_ids.size() == data.size(),
3730  ExcDimensionMismatch(cell_ids.size(), data.size()));
3731  // archive the cellids in an efficient binary format
3732  const std::size_t n_cells = cell_ids.size();
3733  ar & n_cells;
3734  for (const auto &id : cell_ids)
3735  {
3736  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3737  ar & binary_cell_id;
3738  }
3739 
3740  ar &data;
3741  }
3742 
3743 
3744 
3745  template <int dim, typename T>
3746  template <class Archive>
3747  void
3749  const unsigned int /*version*/)
3750  {
3751  std::size_t n_cells;
3752  ar & n_cells;
3753  cell_ids.clear();
3754  cell_ids.reserve(n_cells);
3755  for (unsigned int c = 0; c < n_cells; ++c)
3756  {
3757  CellId::binary_type value;
3758  ar & value;
3759  cell_ids.emplace_back(value);
3760  }
3761  ar &data;
3762  }
3763 
3764 
3765 
3766  template <typename DataType, typename MeshType>
3767  void
3769  const MeshType & mesh,
3770  const std::function<boost::optional<DataType>(
3771  const typename MeshType::active_cell_iterator &)> &pack,
3772  const std::function<void(const typename MeshType::active_cell_iterator &,
3773  const DataType &)> & unpack)
3774  {
3775 # ifndef DEAL_II_WITH_MPI
3776  (void)mesh;
3777  (void)pack;
3778  (void)unpack;
3779  Assert(false,
3780  ExcMessage(
3781  "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
3782 # else
3783  constexpr int dim = MeshType::dimension;
3784  constexpr int spacedim = MeshType::space_dimension;
3785  auto tria = static_cast<const parallel::Triangulation<dim, spacedim> *>(
3786  &mesh.get_triangulation());
3787  Assert(
3788  tria != nullptr,
3789  ExcMessage(
3790  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3791 
3792  // map neighbor_id -> data_buffer where we accumulate the data to send
3793  using DestinationToBufferMap =
3794  std::map<::types::subdomain_id,
3796  DestinationToBufferMap destination_to_data_buffer_map;
3797 
3798  std::map<unsigned int, std::set<::types::subdomain_id>>
3799  vertices_with_ghost_neighbors =
3800  tria->compute_vertices_with_ghost_neighbors();
3801 
3802  for (const auto &cell : tria->active_cell_iterators())
3803  if (cell->is_locally_owned())
3804  {
3805  std::set<::types::subdomain_id> send_to;
3806  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
3807  ++v)
3808  {
3809  const std::map<unsigned int,
3810  std::set<::types::subdomain_id>>::
3811  const_iterator neighbor_subdomains_of_vertex =
3812  vertices_with_ghost_neighbors.find(cell->vertex_index(v));
3813 
3814  if (neighbor_subdomains_of_vertex ==
3815  vertices_with_ghost_neighbors.end())
3816  continue;
3817 
3818  Assert(neighbor_subdomains_of_vertex->second.size() != 0,
3819  ExcInternalError());
3820 
3821  send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
3822  neighbor_subdomains_of_vertex->second.end());
3823  }
3824 
3825  if (send_to.size() > 0)
3826  {
3827  // this cell's data needs to be sent to someone
3828  typename MeshType::active_cell_iterator mesh_it(tria,
3829  cell->level(),
3830  cell->index(),
3831  &mesh);
3832 
3833  const boost::optional<DataType> data = pack(mesh_it);
3834 
3835  if (data)
3836  {
3837  const CellId cellid = cell->id();
3838 
3839  for (const auto subdomain : send_to)
3840  {
3841  // find the data buffer for proc "subdomain" if it exists
3842  // or create an empty one otherwise
3843  typename DestinationToBufferMap::iterator p =
3844  destination_to_data_buffer_map
3845  .insert(std::make_pair(
3846  subdomain, CellDataTransferBuffer<dim, DataType>()))
3847  .first;
3848 
3849  p->second.cell_ids.emplace_back(cellid);
3850  p->second.data.emplace_back(data.get());
3851  }
3852  }
3853  }
3854  }
3855 
3856 
3857  // 2. send our messages
3858  std::set<::types::subdomain_id> ghost_owners = tria->ghost_owners();
3859  const unsigned int n_ghost_owners = ghost_owners.size();
3860  std::vector<std::vector<char>> sendbuffers(n_ghost_owners);
3861  std::vector<MPI_Request> requests(n_ghost_owners);
3862 
3863  unsigned int idx = 0;
3864  for (auto it = ghost_owners.begin(); it != ghost_owners.end(); ++it, ++idx)
3865  {
3866  CellDataTransferBuffer<dim, DataType> &data =
3867  destination_to_data_buffer_map[*it];
3868 
3869  // pack all the data into the buffer for this recipient and send it.
3870  // keep data around till we can make sure that the packet has been
3871  // received
3872  sendbuffers[idx] = Utilities::pack(data);
3873  const int ierr = MPI_Isend(sendbuffers[idx].data(),
3874  sendbuffers[idx].size(),
3875  MPI_BYTE,
3876  *it,
3877  786,
3878  tria->get_communicator(),
3879  &requests[idx]);
3880  AssertThrowMPI(ierr);
3881  }
3882 
3883  // 3. receive messages
3884  std::vector<char> receive;
3885  for (unsigned int idx = 0; idx < n_ghost_owners; ++idx)
3886  {
3887  MPI_Status status;
3888  int len;
3889  int ierr =
3890  MPI_Probe(MPI_ANY_SOURCE, 786, tria->get_communicator(), &status);
3891  AssertThrowMPI(ierr);
3892  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3893  AssertThrowMPI(ierr);
3894 
3895  receive.resize(len);
3896 
3897  char *ptr = receive.data();
3898  ierr = MPI_Recv(ptr,
3899  len,
3900  MPI_BYTE,
3901  status.MPI_SOURCE,
3902  status.MPI_TAG,
3903  tria->get_communicator(),
3904  &status);
3905  AssertThrowMPI(ierr);
3906 
3907  auto cellinfo =
3908  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(receive);
3909 
3910  DataType *data = cellinfo.data.data();
3911  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
3912  {
3914  tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
3915 
3916  const typename MeshType::active_cell_iterator cell(
3917  tria, tria_cell->level(), tria_cell->index(), &mesh);
3918 
3919  unpack(cell, *data);
3920  }
3921  }
3922 
3923  // make sure that all communication is finished
3924  // when we leave this function.
3925  if (requests.size())
3926  {
3927  const int ierr =
3928  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
3929  AssertThrowMPI(ierr);
3930  }
3931 # endif // DEAL_II_WITH_MPI
3932  }
3933 } // namespace GridTools
3934 
3935 # endif
3936 
3937 DEAL_II_NAMESPACE_CLOSE
3938 
3939 /*---------------------------- grid_tools.h ---------------------------*/
3940 /* end of #ifndef dealii_grid_tools_h */
3941 #endif
3942 /*---------------------------- grid_tools.h ---------------------------*/
void remove_hanging_nodes(Triangulation< dim, spacedim > &tria, const bool isotropic=false, const unsigned int max_iterations=100)
Definition: grid_tools.cc:3937
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:3715
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)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:211
unsigned int n_vertices() const
const types::manifold_id flat_manifold_id
Definition: types.h:246
static const unsigned int invalid_unsigned_int
Definition: types.h:173
unsigned int manifold_id
Definition: types.h:123
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:1055
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3690
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:2054
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:75
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true)
Definition: grid_tools.cc:1088
static ::ExceptionBase & ExcPointNotFoundInCoarseGrid(Point< N > arg1)
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: grid_tools.cc:5309
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:3905
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:420
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position)
Definition: grid_tools.cc:1824
void load(Archive &ar, const unsigned int version)
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:4007
bool orthogonal_equality(std::bitset< 3 > &orientation, const FaceIterator &face1, const FaceIterator &face2, const int direction, const Tensor< 1, FaceIterator::AccessorType::space_dimension > &offset=Tensor< 1, FaceIterator::AccessorType::space_dimension >(), const FullMatrix< double > &matrix=FullMatrix< double >())
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12055
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2152
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
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)
Definition: grid_tools.cc:5075
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11883
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:2203
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1318
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, MPI_Comm mpi_communicator)
Definition: grid_tools.cc:5467
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:73
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)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3156
std::bitset< 3 > orientation
Definition: grid_tools.h:2382
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3009
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3128
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2906
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:2679
__global__ void scale(Number *val, const Number *V_val, const size_type N)
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:612
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
MeshType< dim, spacedim >::active_cell_iterator find_active_cell_around_point(const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={})
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2608
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)
cell_iterator end() const
Definition: tria.cc:11949
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:501
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
Definition: tria.h:81
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2643
#define Assert(cond, exc)
Definition: exceptions.h:1407
Signals signals
Definition: tria.h:2394
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:2903
Abstract base class for mapping classes.
Definition: dof_tools.h:57
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
#define DeclException0(Exception0)
Definition: exceptions.h:473
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:1907
void save(Archive &ar, const unsigned int version) const
void remove_anisotropy(Triangulation< dim, spacedim > &tria, const double max_ratio=1.6180339887, const unsigned int max_iterations=5)
Definition: grid_tools.cc:3974
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)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:3789
void rotate(const double angle, Triangulation< 2 > &triangulation)
Definition: grid_tools.cc:887
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:3651
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1174
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: hp.h:117
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_halo_layer(const MeshType &mesh)
Definition: cell_id.h:63
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:1597
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
unsigned int global_dof_index
Definition: types.h:89
__global__ void set(Number *val, const Number s, const size_type N)
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:4364
double cell_measure(const std::vector< Point< dim >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< dim >::vertices_per_cell])
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1695
std::map< types::global_dof_index, std::vector< typename DoFHandlerType::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandlerType &dof_handler)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
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:1464
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2560
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id1, const types::boundary_id b_id2, const 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 >())
FullMatrix< double > matrix
Definition: grid_tools.h:2396
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:3817
static ::ExceptionBase & ExcNotImplemented()
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:4328
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3072
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< boost::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack)
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
unsigned int boundary_id
Definition: types.h:111
static ::ExceptionBase & ExcPointNotFound(Point< N > arg1)
std::vector< typename MeshType::active_cell_iterator > compute_active_cell_layer_within_distance(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const double layer_thickness)
std::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=1e-12, const std::vector< bool > &marked_vertices={})
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:3053
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:5328
std::vector< typename MeshType::active_cell_iterator > get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
static ::ExceptionBase & ExcInternalError()
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:713