Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\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 # include <deal.II/base/std_cxx17/optional.h>
25 
26 # include <deal.II/boost_adaptors/bounding_box.h>
27 
28 # include <deal.II/dofs/dof_handler.h>
29 
30 # include <deal.II/fe/mapping.h>
31 # include <deal.II/fe/mapping_q1.h>
32 
33 # include <deal.II/grid/manifold.h>
34 # include <deal.II/grid/tria.h>
35 # include <deal.II/grid/tria_accessor.h>
36 # include <deal.II/grid/tria_iterator.h>
37 
38 # include <deal.II/hp/dof_handler.h>
39 
40 # include <deal.II/lac/sparsity_tools.h>
41 
42 # include <deal.II/numerics/rtree.h>
43 
44 # include <boost/archive/binary_iarchive.hpp>
45 # include <boost/archive/binary_oarchive.hpp>
46 # include <boost/geometry/index/detail/serialization.hpp>
47 # include <boost/geometry/index/rtree.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 // Forward declarations
66 # ifndef DOXYGEN
67 namespace parallel
68 {
69  namespace distributed
70  {
71  template <int, int>
72  class Triangulation;
73  }
74 } // namespace parallel
75 
76 namespace hp
77 {
78  template <int, int>
79  class MappingCollection;
80 }
81 
82 class SparsityPattern;
83 # endif
84 
85 namespace internal
86 {
87  template <int dim, int spacedim, class MeshType>
88  class ActiveCellIterator
89  {
90  public:
91 # ifndef _MSC_VER
92  using type = typename MeshType::active_cell_iterator;
93 # else
95 # endif
96  };
97 
98 # ifdef _MSC_VER
99  template <int dim, int spacedim>
100  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
101  {
102  public:
103  using type = TriaActiveIterator<
105  };
106 
107  template <int dim, int spacedim>
108  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim>>
109  {
110  public:
111  using type = TriaActiveIterator<
113  };
114 # endif
115 } // namespace internal
116 
125 namespace GridTools
126 {
127  template <int dim, int spacedim>
128  class Cache;
129 
134 
141  template <int dim, int spacedim>
142  double
144 
171  template <int dim, int spacedim>
172  double
174  const Mapping<dim, spacedim> & mapping =
176 
187  template <int dim, int spacedim>
188  double
190  const Mapping<dim, spacedim> & mapping =
192 
203  template <int dim, int spacedim>
204  double
206  const Mapping<dim, spacedim> & mapping =
208 
218  template <int dim>
219  double
220  cell_measure(
221  const std::vector<Point<dim>> &all_vertices,
222  const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]);
223 
229  template <int dim, typename T>
230  double
231  cell_measure(const T &, ...);
232 
261  template <int dim>
264  const Mapping<dim> & mapping,
265  const Quadrature<dim> & quadrature);
266 
276  template <int dim>
277  double
279  const Mapping<dim> & mapping,
280  const Quadrature<dim> & quadrature);
281 
295  template <int dim, int spacedim>
298 
318  template <typename Iterator>
321  const Iterator & object,
323 
336  template <int dim, int spacedim>
337  std::
338  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
340 
346 
363  template <int dim, int spacedim>
364  void
365  delete_unused_vertices(std::vector<Point<spacedim>> &vertices,
366  std::vector<CellData<dim>> & cells,
367  SubCellData & subcelldata);
368 
385  template <int dim, int spacedim>
386  void
387  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
388  std::vector<CellData<dim>> & cells,
389  SubCellData & subcelldata,
390  std::vector<unsigned int> & considered_vertices,
391  const double tol = 1e-12);
392 
398 
433  template <int dim, typename Transformation, int spacedim>
434  void
435  transform(const Transformation & transformation,
436  Triangulation<dim, spacedim> &triangulation);
437 
443  template <int dim, int spacedim>
444  void
445  shift(const Tensor<1, spacedim> & shift_vector,
446  Triangulation<dim, spacedim> &triangulation);
447 
448 
456  void
457  rotate(const double angle, Triangulation<2> &triangulation);
458 
471  template <int dim>
472  void
473  rotate(const double angle,
474  const unsigned int axis,
475  Triangulation<dim, 3> &triangulation);
476 
534  template <int dim>
535  void
536  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
537  Triangulation<dim> & tria,
538  const Function<dim, double> *coefficient = nullptr,
539  const bool solve_for_absolute_positions = false);
540 
546  template <int dim, int spacedim>
547  std::map<unsigned int, Point<spacedim>>
549 
557  template <int dim, int spacedim>
558  void
559  scale(const double scaling_factor,
560  Triangulation<dim, spacedim> &triangulation);
561 
572  template <int dim, int spacedim>
573  void
574  distort_random(const double factor,
575  Triangulation<dim, spacedim> &triangulation,
576  const bool keep_boundary = true);
577 
613  template <int dim, int spacedim>
614  void
616  const bool isotropic = false,
617  const unsigned int max_iterations = 100);
618 
645  template <int dim, int spacedim>
646  void
648  const double max_ratio = 1.6180339887,
649  const unsigned int max_iterations = 5);
650 
742  template <int dim, int spacedim>
743  void
745  const double limit_angle_fraction = .75);
746 
752 
804  template <int dim, int spacedim>
805 # ifndef DOXYGEN
806  std::tuple<
807  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
808  std::vector<std::vector<Point<dim>>>,
809  std::vector<std::vector<unsigned int>>>
810 # else
811  return_type
812 # endif
814  const Cache<dim, spacedim> & cache,
815  const std::vector<Point<spacedim>> &points,
817  &cell_hint =
819 
848  template <int dim, int spacedim>
849 # ifndef DOXYGEN
850  std::tuple<
851  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
852  std::vector<std::vector<Point<dim>>>,
853  std::vector<std::vector<unsigned int>>,
854  std::vector<unsigned int>>
855 # else
856  return_type
857 # endif
859  const Cache<dim, spacedim> & cache,
860  const std::vector<Point<spacedim>> &points,
862  &cell_hint =
864 
929  template <int dim, int spacedim>
930 # ifndef DOXYGEN
931  std::tuple<
932  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
933  std::vector<std::vector<Point<dim>>>,
934  std::vector<std::vector<unsigned int>>,
935  std::vector<std::vector<Point<spacedim>>>,
936  std::vector<std::vector<unsigned int>>>
937 # else
938  return_type
939 # endif
941  const GridTools::Cache<dim, spacedim> & cache,
942  const std::vector<Point<spacedim>> & local_points,
943  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes);
944 
983  template <int dim, int spacedim>
984  std::map<unsigned int, Point<spacedim>>
986  const Mapping<dim, spacedim> & mapping =
988 
1000  template <int spacedim>
1001  unsigned int
1002  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1003  const Point<spacedim> & p);
1004 
1031  template <int dim, template <int, int> class MeshType, int spacedim>
1032  unsigned int
1033  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1034  const Point<spacedim> & p,
1035  const std::vector<bool> & marked_vertices = {});
1036 
1062  template <int dim, template <int, int> class MeshType, int spacedim>
1063  unsigned int
1065  const MeshType<dim, spacedim> &mesh,
1066  const Point<spacedim> & p,
1067  const std::vector<bool> & marked_vertices = {});
1068 
1069 
1094  template <int dim, template <int, int> class MeshType, int spacedim>
1095 # ifndef _MSC_VER
1096  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1097 # else
1098  std::vector<
1099  typename ::internal::
1100  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1101 # endif
1102  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1103  const unsigned int vertex_index);
1104 
1105 
1130  template <int dim, template <int, int> class MeshType, int spacedim>
1131 # ifndef _MSC_VER
1132  typename MeshType<dim, spacedim>::active_cell_iterator
1133 # else
1134  typename ::internal::ActiveCellIterator<dim,
1135  spacedim,
1136  MeshType<dim, spacedim>>::type
1137 # endif
1138  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1139  const Point<spacedim> & p,
1140  const std::vector<bool> &marked_vertices = {});
1141 
1229  template <int dim, template <int, int> class MeshType, int spacedim>
1230 # ifndef _MSC_VER
1231  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1232 # else
1233  std::pair<typename ::internal::
1234  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1235  Point<dim>>
1236 # endif
1238  const MeshType<dim, spacedim> &mesh,
1239  const Point<spacedim> & p,
1240  const std::vector<bool> &marked_vertices = {});
1241 
1253  template <int dim, template <int, int> class MeshType, int spacedim>
1254 # ifndef _MSC_VER
1255  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1256 # else
1257  std::pair<typename ::internal::
1258  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1259  Point<dim>>
1260 # endif
1262  const Mapping<dim, spacedim> & mapping,
1263  const MeshType<dim, spacedim> &mesh,
1264  const Point<spacedim> & p,
1265  const std::vector<
1266  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1267  & vertex_to_cell_map,
1268  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1269  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1270  typename MeshType<dim, spacedim>::active_cell_iterator(),
1271  const std::vector<bool> & marked_vertices = {},
1272  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1273  RTree<std::pair<Point<spacedim>, unsigned int>>{});
1274 
1296  template <int dim, int spacedim>
1297  std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator,
1298  Point<dim>>
1300  const hp::MappingCollection<dim, spacedim> &mapping,
1301  const hp::DoFHandler<dim, spacedim> & mesh,
1302  const Point<spacedim> & p);
1303 
1310  template <int dim, int spacedim>
1311  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1312  Point<dim>>
1314  const Cache<dim, spacedim> &cache,
1315  const Point<spacedim> & p,
1318  const std::vector<bool> &marked_vertices = {});
1319 
1337  template <int dim, template <int, int> class MeshType, int spacedim>
1338 # ifndef _MSC_VER
1339  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1340  Point<dim>>>
1341 # else
1342  std::vector<std::pair<
1343  typename ::internal::
1344  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1345  Point<dim>>>
1346 # endif
1348  const Mapping<dim, spacedim> & mapping,
1349  const MeshType<dim, spacedim> &mesh,
1350  const Point<spacedim> & p,
1351  const double tolerance = 1e-12,
1352  const std::vector<bool> & marked_vertices = {});
1353 
1375  template <class MeshType>
1376  std::vector<typename MeshType::active_cell_iterator>
1377  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1378 
1396  template <class MeshType>
1397  void
1399  const typename MeshType::active_cell_iterator & cell,
1400  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1401 
1454  template <class MeshType>
1455  std::vector<typename MeshType::active_cell_iterator>
1457  const MeshType &mesh,
1458  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1459  &predicate);
1460 
1461 
1469  template <class MeshType>
1470  std::vector<typename MeshType::cell_iterator>
1472  const MeshType &mesh,
1473  const std::function<bool(const typename MeshType::cell_iterator &)>
1474  & predicate,
1475  const unsigned int level);
1476 
1477 
1493  template <class MeshType>
1494  std::vector<typename MeshType::active_cell_iterator>
1495  compute_ghost_cell_halo_layer(const MeshType &mesh);
1496 
1548  template <class MeshType>
1549  std::vector<typename MeshType::active_cell_iterator>
1551  const MeshType &mesh,
1552  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1553  & predicate,
1554  const double layer_thickness);
1555 
1581  template <class MeshType>
1582  std::vector<typename MeshType::active_cell_iterator>
1583  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1584  const double layer_thickness);
1585 
1601  template <class MeshType>
1602  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1604  const MeshType &mesh,
1605  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1606  &predicate);
1607 
1660  template <class MeshType>
1661  std::vector<BoundingBox<MeshType::space_dimension>>
1663  const MeshType &mesh,
1664  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1665  & predicate,
1666  const unsigned int refinement_level = 0,
1667  const bool allow_merge = false,
1668  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1669 
1699  template <int spacedim>
1700 # ifndef DOXYGEN
1701  std::tuple<std::vector<std::vector<unsigned int>>,
1702  std::map<unsigned int, unsigned int>,
1703  std::map<unsigned int, std::vector<unsigned int>>>
1704 # else
1705  return_type
1706 # endif
1708  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1709  const std::vector<Point<spacedim>> & points);
1710 
1711 
1746  template <int spacedim>
1747 # ifndef DOXYGEN
1748  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1749  std::map<unsigned int, unsigned int>,
1750  std::map<unsigned int, std::vector<unsigned int>>>
1751 # else
1752  return_type
1753 # endif
1755  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1756  const std::vector<Point<spacedim>> & points);
1757 
1758 
1767  template <int dim, int spacedim>
1768  std::vector<
1769  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1770  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1771 
1786  template <int dim, int spacedim>
1787  std::vector<std::vector<Tensor<1, spacedim>>>
1789  const Triangulation<dim, spacedim> &mesh,
1790  const std::vector<
1792  &vertex_to_cells);
1793 
1794 
1801  template <int dim, int spacedim>
1802  unsigned int
1805  const Point<spacedim> &position);
1806 
1818  template <int dim, int spacedim>
1819  std::map<unsigned int, types::global_vertex_index>
1822 
1836  template <int dim, int spacedim>
1837  std::pair<unsigned int, double>
1840 
1846 
1855  template <int dim, int spacedim>
1856  void
1858  const Triangulation<dim, spacedim> &triangulation,
1859  DynamicSparsityPattern & connectivity);
1860 
1869  template <int dim, int spacedim>
1870  void
1872  const Triangulation<dim, spacedim> &triangulation,
1873  DynamicSparsityPattern & connectivity);
1874 
1883  template <int dim, int spacedim>
1884  void
1886  const Triangulation<dim, spacedim> &triangulation,
1887  const unsigned int level,
1888  DynamicSparsityPattern & connectivity);
1889 
1910  template <int dim, int spacedim>
1911  void
1912  partition_triangulation(const unsigned int n_partitions,
1913  Triangulation<dim, spacedim> & triangulation,
1914  const SparsityTools::Partitioner partitioner =
1916 
1927  template <int dim, int spacedim>
1928  void
1929  partition_triangulation(const unsigned int n_partitions,
1930  const std::vector<unsigned int> &cell_weights,
1931  Triangulation<dim, spacedim> & triangulation,
1932  const SparsityTools::Partitioner partitioner =
1934 
1980  template <int dim, int spacedim>
1981  void
1982  partition_triangulation(const unsigned int n_partitions,
1983  const SparsityPattern & cell_connection_graph,
1984  Triangulation<dim, spacedim> &triangulation,
1985  const SparsityTools::Partitioner partitioner =
1987 
1998  template <int dim, int spacedim>
1999  void
2000  partition_triangulation(const unsigned int n_partitions,
2001  const std::vector<unsigned int> &cell_weights,
2002  const SparsityPattern & cell_connection_graph,
2003  Triangulation<dim, spacedim> &triangulation,
2004  const SparsityTools::Partitioner partitioner =
2006 
2014  template <int dim, int spacedim>
2015  void
2016  partition_triangulation_zorder(const unsigned int n_partitions,
2017  Triangulation<dim, spacedim> &triangulation);
2018 
2030  template <int dim, int spacedim>
2031  void
2033 
2044  template <int dim, int spacedim>
2045  void
2046  get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
2047  std::vector<types::subdomain_id> & subdomain);
2048 
2063  template <int dim, int spacedim>
2064  unsigned int
2066  const Triangulation<dim, spacedim> &triangulation,
2067  const types::subdomain_id subdomain);
2068 
2069 
2099  template <int dim, int spacedim>
2100  std::vector<bool>
2102 
2108 
2142  template <typename MeshType>
2143  std::list<std::pair<typename MeshType::cell_iterator,
2144  typename MeshType::cell_iterator>>
2145  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2146 
2156  template <int dim, int spacedim>
2157  bool
2159  const Triangulation<dim, spacedim> &mesh_2);
2160 
2170  template <typename MeshType>
2171  bool
2172  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2173 
2179 
2195  template <int dim, int spacedim>
2199  & distorted_cells,
2200  Triangulation<dim, spacedim> &triangulation);
2201 
2202 
2203 
2212 
2213 
2253  template <class MeshType>
2254  std::vector<typename MeshType::active_cell_iterator>
2255  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2256 
2257 
2281  template <class Container>
2282  std::vector<typename Container::cell_iterator>
2284  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2285 
2354  template <class Container>
2355  void
2357  const std::vector<typename Container::active_cell_iterator> &patch,
2359  &local_triangulation,
2360  std::map<
2361  typename Triangulation<Container::dimension,
2362  Container::space_dimension>::active_cell_iterator,
2363  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2364 
2401  template <class DoFHandlerType>
2402  std::map<types::global_dof_index,
2403  std::vector<typename DoFHandlerType::active_cell_iterator>>
2404  get_dof_to_support_patch_map(DoFHandlerType &dof_handler);
2405 
2406 
2413 
2419  template <typename CellIterator>
2421  {
2425  CellIterator cell[2];
2426 
2431  unsigned int face_idx[2];
2432 
2438  std::bitset<3> orientation;
2439 
2453  };
2454 
2455 
2521  template <typename FaceIterator>
2522  bool orthogonal_equality(
2523  std::bitset<3> & orientation,
2524  const FaceIterator & face1,
2525  const FaceIterator & face2,
2526  const int direction,
2529  const FullMatrix<double> &matrix = FullMatrix<double>());
2530 
2531 
2535  template <typename FaceIterator>
2536  bool
2538  const FaceIterator & face1,
2539  const FaceIterator & face2,
2540  const int direction,
2543  const FullMatrix<double> &matrix = FullMatrix<double>());
2544 
2545 
2604  template <typename MeshType>
2605  void
2607  const MeshType & mesh,
2608  const types::boundary_id b_id1,
2609  const types::boundary_id b_id2,
2610  const int direction,
2612  & matched_pairs,
2613  const Tensor<1, MeshType::space_dimension> &offset =
2615  const FullMatrix<double> &matrix = FullMatrix<double>());
2616 
2617 
2642  template <typename MeshType>
2643  void
2645  const MeshType & mesh,
2646  const types::boundary_id b_id,
2647  const int direction,
2649  & matched_pairs,
2650  const ::Tensor<1, MeshType::space_dimension> &offset =
2652  const FullMatrix<double> &matrix = FullMatrix<double>());
2653 
2659 
2682  template <int dim, int spacedim>
2683  void
2685  const bool reset_boundary_ids = false);
2686 
2710  template <int dim, int spacedim>
2711  void
2713  const std::vector<types::boundary_id> &src_boundary_ids,
2714  const std::vector<types::manifold_id> &dst_manifold_ids,
2716  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2717 
2749  template <int dim, int spacedim>
2750  void
2752  const bool compute_face_ids = false);
2753 
2780  template <int dim, int spacedim>
2781  void
2784  const std::function<types::manifold_id(
2785  const std::set<types::manifold_id> &)> &disambiguation_function =
2786  [](const std::set<types::manifold_id> &manifold_ids) {
2787  if (manifold_ids.size() == 1)
2788  return *manifold_ids.begin();
2789  else
2791  },
2792  bool overwrite_only_flat_manifold_ids = true);
2877  template <typename DataType, typename MeshType>
2878  void
2880  const MeshType & mesh,
2881  const std::function<std_cxx17::optional<DataType>(
2882  const typename MeshType::active_cell_iterator &)> &pack,
2883  const std::function<void(const typename MeshType::active_cell_iterator &,
2884  const DataType &)> & unpack);
2885 
2886  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2887  * boxes @p local_bboxes.
2888  *
2889  * This function is meant to exchange bounding boxes describing the locally
2890  * owned cells in a distributed triangulation obtained with the function
2891  * GridTools::compute_mesh_predicate_bounding_box .
2892  *
2893  * The output vector's size is the number of processes of the MPI
2894  * communicator:
2895  * its i-th entry contains the vector @p local_bboxes of the i-th process.
2896  */
2897  template <int spacedim>
2898  std::vector<std::vector<BoundingBox<spacedim>>>
2899  exchange_local_bounding_boxes(
2900  const std::vector<BoundingBox<spacedim>> &local_bboxes,
2901  MPI_Comm mpi_communicator);
2902 
2937  template <int spacedim>
2938  RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
2940  const std::vector<BoundingBox<spacedim>> &local_description,
2941  MPI_Comm mpi_communicator);
2942 
2955  template <int dim, typename T>
2957  {
2961  std::vector<CellId> cell_ids;
2962 
2966  std::vector<T> data;
2967 
2975  template <class Archive>
2976  void
2977  save(Archive &ar, const unsigned int version) const;
2978 
2983  template <class Archive>
2984  void
2985  load(Archive &ar, const unsigned int version);
2986 
2987  BOOST_SERIALIZATION_SPLIT_MEMBER()
2988  };
2989 
2994 
2999  int,
3000  << "The number of partitions you gave is " << arg1
3001  << ", but must be greater than zero.");
3006  int,
3007  << "The subdomain id " << arg1
3008  << " has no cells associated with it.");
3013 
3018  double,
3019  << "The scaling factor must be positive, but it is " << arg1
3020  << ".");
3024  template <int N>
3026  Point<N>,
3027  << "The point <" << arg1
3028  << "> could not be found inside any of the "
3029  << "coarse grid cells.");
3033  template <int N>
3035  Point<N>,
3036  << "The point <" << arg1
3037  << "> could not be found inside any of the "
3038  << "subcells of a coarse grid cell.");
3039 
3044  unsigned int,
3045  << "The given vertex with index " << arg1
3046  << " is not used in the given triangulation.");
3047 
3048 
3051 } /*namespace GridTools*/
3052 
3053 
3054 
3055 /* ----------------- Template function --------------- */
3056 
3057 # ifndef DOXYGEN
3058 
3059 namespace GridTools
3060 {
3061  // declare specializations
3062  template <>
3063  double
3064  cell_measure<1>(const std::vector<Point<1>> &,
3065  const unsigned int (&)[GeometryInfo<1>::vertices_per_cell]);
3066 
3067  template <>
3068  double
3069  cell_measure<2>(const std::vector<Point<2>> &,
3070  const unsigned int (&)[GeometryInfo<2>::vertices_per_cell]);
3071 
3072  template <>
3073  double
3074  cell_measure<3>(const std::vector<Point<3>> &,
3075  const unsigned int (&)[GeometryInfo<3>::vertices_per_cell]);
3076 
3077 
3078 
3079  template <int dim, typename T>
3080  double
3081  cell_measure(const T &, ...)
3082  {
3083  Assert(false, ExcNotImplemented());
3084  return std::numeric_limits<double>::quiet_NaN();
3085  }
3086 
3087 
3088 
3089  template <int dim, typename Predicate, int spacedim>
3090  void
3091  transform(const Predicate & predicate,
3092  Triangulation<dim, spacedim> &triangulation)
3093  {
3094  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3095 
3096  // loop over all active cells, and
3097  // transform those vertices that
3098  // have not yet been touched. note
3099  // that we get to all vertices in
3100  // the triangulation by only
3101  // visiting the active cells.
3103  cell = triangulation.begin_active(),
3104  endc = triangulation.end();
3105  for (; cell != endc; ++cell)
3106  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
3107  if (treated_vertices[cell->vertex_index(v)] == false)
3108  {
3109  // transform this vertex
3110  cell->vertex(v) = predicate(cell->vertex(v));
3111  // and mark it as treated
3112  treated_vertices[cell->vertex_index(v)] = true;
3113  };
3114 
3115 
3116  // now fix any vertices on hanging nodes so that we don't create any holes
3117  if (dim == 2)
3118  {
3120  cell = triangulation.begin_active(),
3121  endc = triangulation.end();
3122  for (; cell != endc; ++cell)
3123  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
3124  ++face)
3125  if (cell->face(face)->has_children() &&
3126  !cell->face(face)->at_boundary())
3127  {
3128  // this line has children
3129  cell->face(face)->child(0)->vertex(1) =
3130  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3131  2;
3132  }
3133  }
3134  else if (dim == 3)
3135  {
3137  cell = triangulation.begin_active(),
3138  endc = triangulation.end();
3139  for (; cell != endc; ++cell)
3140  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
3141  ++face)
3142  if (cell->face(face)->has_children() &&
3143  !cell->face(face)->at_boundary())
3144  {
3145  // this face has hanging nodes
3146  cell->face(face)->child(0)->vertex(1) =
3147  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3148  2.0;
3149  cell->face(face)->child(0)->vertex(2) =
3150  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3151  2.0;
3152  cell->face(face)->child(1)->vertex(3) =
3153  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3154  2.0;
3155  cell->face(face)->child(2)->vertex(3) =
3156  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3157  2.0;
3158 
3159  // center of the face
3160  cell->face(face)->child(0)->vertex(3) =
3161  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3162  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3163  4.0;
3164  }
3165  }
3166 
3167  // Make sure FEValues notices that the mesh has changed
3168  triangulation.signals.mesh_movement();
3169  }
3170 
3171 
3172 
3173  template <class MeshType>
3174  std::vector<typename MeshType::active_cell_iterator>
3175  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3176  {
3177  std::vector<typename MeshType::active_cell_iterator> child_cells;
3178 
3179  if (cell->has_children())
3180  {
3181  for (unsigned int child = 0; child < cell->n_children(); ++child)
3182  if (cell->child(child)->has_children())
3183  {
3184  const std::vector<typename MeshType::active_cell_iterator>
3185  children = get_active_child_cells<MeshType>(cell->child(child));
3186  child_cells.insert(child_cells.end(),
3187  children.begin(),
3188  children.end());
3189  }
3190  else
3191  child_cells.push_back(cell->child(child));
3192  }
3193 
3194  return child_cells;
3195  }
3196 
3197 
3198 
3199  template <class MeshType>
3200  void
3202  const typename MeshType::active_cell_iterator & cell,
3203  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3204  {
3205  active_neighbors.clear();
3206  for (unsigned int n = 0;
3207  n < GeometryInfo<MeshType::dimension>::faces_per_cell;
3208  ++n)
3209  if (!cell->at_boundary(n))
3210  {
3211  if (MeshType::dimension == 1)
3212  {
3213  // check children of neighbor. note
3214  // that in 1d children of the neighbor
3215  // may be further refined. In 1d the
3216  // case is simple since we know what
3217  // children bound to the present cell
3218  typename MeshType::cell_iterator neighbor_child =
3219  cell->neighbor(n);
3220  if (!neighbor_child->active())
3221  {
3222  while (neighbor_child->has_children())
3223  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3224 
3225  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3226  ExcInternalError());
3227  }
3228  active_neighbors.push_back(neighbor_child);
3229  }
3230  else
3231  {
3232  if (cell->face(n)->has_children())
3233  // this neighbor has children. find
3234  // out which border to the present
3235  // cell
3236  for (unsigned int c = 0;
3237  c < cell->face(n)->number_of_children();
3238  ++c)
3239  active_neighbors.push_back(
3240  cell->neighbor_child_on_subface(n, c));
3241  else
3242  {
3243  // the neighbor must be active
3244  // himself
3245  Assert(cell->neighbor(n)->active(), ExcInternalError());
3246  active_neighbors.push_back(cell->neighbor(n));
3247  }
3248  }
3249  }
3250  }
3251 
3252 
3253 
3254  namespace internal
3255  {
3256  namespace ProjectToObject
3257  {
3270  struct CrossDerivative
3271  {
3272  const unsigned int direction_0;
3273  const unsigned int direction_1;
3274 
3275  CrossDerivative(const unsigned int d0, const unsigned int d1);
3276  };
3277 
3278  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3279  const unsigned int d1)
3280  : direction_0(d0)
3281  , direction_1(d1)
3282  {}
3283 
3284 
3285 
3290  template <typename F>
3291  inline auto
3292  centered_first_difference(const double center,
3293  const double step,
3294  const F &f) -> decltype(f(center) - f(center))
3295  {
3296  return (f(center + step) - f(center - step)) / (2.0 * step);
3297  }
3298 
3299 
3300 
3305  template <typename F>
3306  inline auto
3307  centered_second_difference(const double center,
3308  const double step,
3309  const F &f) -> decltype(f(center) - f(center))
3310  {
3311  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3312  (step * step);
3313  }
3314 
3315 
3316 
3326  template <int structdim, typename F>
3327  inline auto
3328  cross_stencil(
3329  const CrossDerivative cross_derivative,
3331  const double step,
3332  const F &f) -> decltype(f(center) - f(center))
3333  {
3335  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3336  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3337  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3338  1.0 / 3.0 * f(center - simplex_vector) +
3339  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3340  step;
3341  }
3342 
3343 
3344 
3351  template <int spacedim, int structdim, typename F>
3352  inline double
3353  gradient_entry(
3354  const unsigned int row_n,
3355  const unsigned int dependent_direction,
3356  const Point<spacedim> &p0,
3358  const double step,
3359  const F & f)
3360  {
3362  dependent_direction <
3364  ExcMessage("This function assumes that the last weight is a "
3365  "dependent variable (and hence we cannot take its "
3366  "derivative directly)."));
3367  Assert(row_n != dependent_direction,
3368  ExcMessage(
3369  "We cannot differentiate with respect to the variable "
3370  "that is assumed to be dependent."));
3371 
3372  const Point<spacedim> manifold_point = f(center);
3373  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3374  {row_n, dependent_direction}, center, step, f);
3375  double entry = 0.0;
3376  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3377  entry +=
3378  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3379  return entry;
3380  }
3381 
3387  template <typename Iterator, int spacedim, int structdim>
3389  project_to_d_linear_object(const Iterator & object,
3390  const Point<spacedim> &trial_point)
3391  {
3392  // let's look at this for simplicity for a quad (structdim==2) in a
3393  // space with spacedim>2 (notate trial_point by y): all points on the
3394  // surface are given by
3395  // x(\xi) = sum_i v_i phi_x(\xi)
3396  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3397  // reference coordinates of the quad. so what we are trying to do is
3398  // find a point x on the surface that is closest to the point y. there
3399  // are different ways to solve this problem, but in the end it's a
3400  // nonlinear problem and we have to find reference coordinates \xi so
3401  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3402  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3403  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3404  // have to use a Newton method to find the answer. This leads to the
3405  // following formulation of Newton steps:
3406  //
3407  // Given \xi_k, find \delta\xi_k so that
3408  // H_k \delta\xi_k = - F_k
3409  // where H_k is an approximation to the second derivatives of J at
3410  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3411  // number of times until the right hand side is small enough. As a
3412  // stopping criterion, we terminate if ||\delta\xi||<eps.
3413  //
3414  // As for the Hessian, the best choice would be
3415  // H_k = J''(\xi_k)
3416  // but we'll opt for the simpler Gauss-Newton form
3417  // H_k = A^T A
3418  // i.e.
3419  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3420  // \partial_n phi_i *
3421  // \partial_m phi_j
3422  // we start at xi=(0.5, 0.5).
3423  Point<structdim> xi;
3424  for (unsigned int d = 0; d < structdim; ++d)
3425  xi[d] = 0.5;
3426 
3427  Point<spacedim> x_k;
3428  for (unsigned int i = 0; i < GeometryInfo<structdim>::vertices_per_cell;
3429  ++i)
3430  x_k += object->vertex(i) *
3432 
3433  do
3434  {
3436  for (unsigned int i = 0;
3437  i < GeometryInfo<structdim>::vertices_per_cell;
3438  ++i)
3439  F_k +=
3440  (x_k - trial_point) * object->vertex(i) *
3442  i);
3443 
3445  for (unsigned int i = 0;
3446  i < GeometryInfo<structdim>::vertices_per_cell;
3447  ++i)
3448  for (unsigned int j = 0;
3449  j < GeometryInfo<structdim>::vertices_per_cell;
3450  ++j)
3451  {
3454  xi, i),
3456  xi, j));
3457  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3458  }
3459 
3460  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3461  xi += delta_xi;
3462 
3463  x_k = Point<spacedim>();
3464  for (unsigned int i = 0;
3465  i < GeometryInfo<structdim>::vertices_per_cell;
3466  ++i)
3467  x_k += object->vertex(i) *
3469 
3470  if (delta_xi.norm() < 1e-7)
3471  break;
3472  }
3473  while (true);
3474 
3475  return x_k;
3476  }
3477  } // namespace ProjectToObject
3478  } // namespace internal
3479 
3480 
3481 
3482  namespace internal
3483  {
3484  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3485  // inside the project_to_object function below.
3486  template <int structdim>
3487  inline bool
3488  weights_are_ok(
3490  {
3491  // clang has trouble figuring out structdim here, so define it
3492  // again:
3493  static const std::size_t n_vertices_per_cell =
3495  n_independent_components;
3496  std::array<double, n_vertices_per_cell> copied_weights;
3497  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3498  {
3499  copied_weights[i] = v[i];
3500  if (v[i] < 0.0 || v[i] > 1.0)
3501  return false;
3502  }
3503 
3504  // check the sum: try to avoid some roundoff errors by summing in order
3505  std::sort(copied_weights.begin(), copied_weights.end());
3506  const double sum =
3507  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3508  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3509  }
3510  } // namespace internal
3511 
3512  template <typename Iterator>
3515  const Iterator & object,
3517  {
3518  const int spacedim = Iterator::AccessorType::space_dimension;
3519  const int structdim = Iterator::AccessorType::structure_dimension;
3520 
3521  Point<spacedim> projected_point = trial_point;
3522 
3523  if (structdim >= spacedim)
3524  return projected_point;
3525  else if (structdim == 1 || structdim == 2)
3526  {
3527  using namespace internal::ProjectToObject;
3528  // Try to use the special flat algorithm for quads (this is better
3529  // than the general algorithm in 3D). This does not take into account
3530  // whether projected_point is outside the quad, but we optimize along
3531  // lines below anyway:
3532  const int dim = Iterator::AccessorType::dimension;
3533  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3534  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3535  &manifold) != nullptr)
3536  {
3537  projected_point =
3538  project_to_d_linear_object<Iterator, spacedim, structdim>(
3539  object, trial_point);
3540  }
3541  else
3542  {
3543  // We want to find a point on the convex hull (defined by the
3544  // vertices of the object and the manifold description) that is
3545  // relatively close to the trial point. This has a few issues:
3546  //
3547  // 1. For a general convex hull we are not guaranteed that a unique
3548  // minimum exists.
3549  // 2. The independent variables in the optimization process are the
3550  // weights given to Manifold::get_new_point, which must sum to 1,
3551  // so we cannot use standard finite differences to approximate a
3552  // gradient.
3553  //
3554  // There is not much we can do about 1., but for 2. we can derive
3555  // finite difference stencils that work on a structdim-dimensional
3556  // simplex and rewrite the optimization problem to use those
3557  // instead. Consider the structdim 2 case and let
3558  //
3559  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3560  // c2, c3})
3561  //
3562  // where {c0, c1, c2, c3} are the weights for the four vertices on
3563  // the quadrilateral. We seek to minimize the Euclidean distance
3564  // between F(...) and trial_point. We can solve for c3 in terms of
3565  // the other weights and get, for one coordinate direction
3566  //
3567  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3568  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3569  //
3570  // where we substitute back in for c3 after taking the
3571  // derivative. We can compute a stencil for the cross derivative
3572  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3573  // (and gradient_entry computes the sum over the independent
3574  // variables). Below, we somewhat arbitrarily pick the last
3575  // component as the dependent one.
3576  //
3577  // Since we can now calculate derivatives of the objective
3578  // function we can use gradient descent to minimize it.
3579  //
3580  // Of course, this is much simpler in the structdim = 1 case (we
3581  // could rewrite the projection as a 1D optimization problem), but
3582  // to reduce the potential for bugs we use the same code in both
3583  // cases.
3584  const double step_size = object->diameter() / 64.0;
3585 
3586  constexpr unsigned int n_vertices_per_cell =
3588 
3589  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3590  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3591  ++vertex_n)
3592  vertices[vertex_n] = object->vertex(vertex_n);
3593 
3594  auto get_point_from_weights =
3595  [&](const Tensor<1, n_vertices_per_cell> &weights)
3596  -> Point<spacedim> {
3597  return object->get_manifold().get_new_point(
3598  make_array_view(vertices.begin(), vertices.end()),
3599  make_array_view(weights.begin_raw(), weights.end_raw()));
3600  };
3601 
3602  // pick the initial weights as (normalized) inverse distances from
3603  // the trial point:
3604  Tensor<1, n_vertices_per_cell> guess_weights;
3605  double guess_weights_sum = 0.0;
3606  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3607  ++vertex_n)
3608  {
3609  const double distance =
3610  vertices[vertex_n].distance(trial_point);
3611  if (distance == 0.0)
3612  {
3613  guess_weights = 0.0;
3614  guess_weights[vertex_n] = 1.0;
3615  guess_weights_sum = 1.0;
3616  break;
3617  }
3618  else
3619  {
3620  guess_weights[vertex_n] = 1.0 / distance;
3621  guess_weights_sum += guess_weights[vertex_n];
3622  }
3623  }
3624  guess_weights /= guess_weights_sum;
3625  Assert(internal::weights_are_ok<structdim>(guess_weights),
3626  ExcInternalError());
3627 
3628  // The optimization algorithm consists of two parts:
3629  //
3630  // 1. An outer loop where we apply the gradient descent algorithm.
3631  // 2. An inner loop where we do a line search to find the optimal
3632  // length of the step one should take in the gradient direction.
3633  //
3634  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3635  {
3636  const unsigned int dependent_direction =
3637  n_vertices_per_cell - 1;
3638  Tensor<1, n_vertices_per_cell> current_gradient;
3639  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3640  ++row_n)
3641  {
3642  if (row_n != dependent_direction)
3643  {
3644  current_gradient[row_n] =
3645  gradient_entry<spacedim, structdim>(
3646  row_n,
3647  dependent_direction,
3648  trial_point,
3649  guess_weights,
3650  step_size,
3651  get_point_from_weights);
3652 
3653  current_gradient[dependent_direction] -=
3654  current_gradient[row_n];
3655  }
3656  }
3657 
3658  // We need to travel in the -gradient direction, as noted
3659  // above, but we may not want to take a full step in that
3660  // direction; instead, guess that we will go -0.5*gradient and
3661  // do quasi-Newton iteration to pick the best multiplier. The
3662  // goal is to find a scalar alpha such that
3663  //
3664  // F(x - alpha g)
3665  //
3666  // is minimized, where g is the gradient and F is the
3667  // objective function. To find the optimal value we find roots
3668  // of the derivative of the objective function with respect to
3669  // alpha by Newton iteration, where we approximate the first
3670  // and second derivatives of F(x - alpha g) with centered
3671  // finite differences.
3672  double gradient_weight = -0.5;
3673  auto gradient_weight_objective_function =
3674  [&](const double gradient_weight_guess) -> double {
3675  return (trial_point -
3676  get_point_from_weights(guess_weights +
3677  gradient_weight_guess *
3678  current_gradient))
3679  .norm_square();
3680  };
3681 
3682  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3683  {
3684  const double update_numerator = centered_first_difference(
3685  gradient_weight,
3686  step_size,
3687  gradient_weight_objective_function);
3688  const double update_denominator =
3689  centered_second_difference(
3690  gradient_weight,
3691  step_size,
3692  gradient_weight_objective_function);
3693 
3694  // avoid division by zero. Note that we limit the gradient
3695  // weight below
3696  if (std::abs(update_denominator) == 0.0)
3697  break;
3698  gradient_weight =
3699  gradient_weight - update_numerator / update_denominator;
3700 
3701  // Put a fairly lenient bound on the largest possible
3702  // gradient (things tend to be locally flat, so the gradient
3703  // itself is usually small)
3704  if (std::abs(gradient_weight) > 10)
3705  {
3706  gradient_weight = -10.0;
3707  break;
3708  }
3709  }
3710 
3711  // It only makes sense to take convex combinations with weights
3712  // between zero and one. If the update takes us outside of this
3713  // region then rescale the update to stay within the region and
3714  // try again
3715  Tensor<1, n_vertices_per_cell> tentative_weights =
3716  guess_weights + gradient_weight * current_gradient;
3717 
3718  double new_gradient_weight = gradient_weight;
3719  for (unsigned int iteration_count = 0; iteration_count < 40;
3720  ++iteration_count)
3721  {
3722  if (internal::weights_are_ok<structdim>(tentative_weights))
3723  break;
3724 
3725  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3726  {
3727  if (tentative_weights[i] < 0.0)
3728  {
3729  tentative_weights -=
3730  (tentative_weights[i] / current_gradient[i]) *
3731  current_gradient;
3732  }
3733  if (tentative_weights[i] < 0.0 ||
3734  1.0 < tentative_weights[i])
3735  {
3736  new_gradient_weight /= 2.0;
3737  tentative_weights =
3738  guess_weights +
3739  new_gradient_weight * current_gradient;
3740  }
3741  }
3742  }
3743 
3744  // the update might still send us outside the valid region, so
3745  // check again and quit if the update is still not valid
3746  if (!internal::weights_are_ok<structdim>(tentative_weights))
3747  break;
3748 
3749  // if we cannot get closer by traveling in the gradient
3750  // direction then quit
3751  if (get_point_from_weights(tentative_weights)
3752  .distance(trial_point) <
3753  get_point_from_weights(guess_weights).distance(trial_point))
3754  guess_weights = tentative_weights;
3755  else
3756  break;
3757  Assert(internal::weights_are_ok<structdim>(guess_weights),
3758  ExcInternalError());
3759  }
3760  Assert(internal::weights_are_ok<structdim>(guess_weights),
3761  ExcInternalError());
3762  projected_point = get_point_from_weights(guess_weights);
3763  }
3764 
3765  // if structdim == 2 and the optimal point is not on the interior then
3766  // we may be able to get a more accurate result by projecting onto the
3767  // lines.
3768  if (structdim == 2)
3769  {
3770  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
3771  line_projections;
3772  for (unsigned int line_n = 0;
3773  line_n < GeometryInfo<structdim>::lines_per_cell;
3774  ++line_n)
3775  {
3776  line_projections[line_n] =
3777  project_to_object(object->line(line_n), trial_point);
3778  }
3779  std::sort(line_projections.begin(),
3780  line_projections.end(),
3781  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
3782  return a.distance(trial_point) <
3783  b.distance(trial_point);
3784  });
3785  if (line_projections[0].distance(trial_point) <
3786  projected_point.distance(trial_point))
3787  projected_point = line_projections[0];
3788  }
3789  }
3790  else
3791  {
3792  Assert(false, ExcNotImplemented());
3793  return projected_point;
3794  }
3795 
3796  return projected_point;
3797  }
3798 
3799 
3800 
3801  template <int dim, typename T>
3802  template <class Archive>
3803  void
3805  const unsigned int /*version*/) const
3806  {
3807  Assert(cell_ids.size() == data.size(),
3808  ExcDimensionMismatch(cell_ids.size(), data.size()));
3809  // archive the cellids in an efficient binary format
3810  const std::size_t n_cells = cell_ids.size();
3811  ar & n_cells;
3812  for (const auto &id : cell_ids)
3813  {
3814  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3815  ar & binary_cell_id;
3816  }
3817 
3818  ar &data;
3819  }
3820 
3821 
3822 
3823  template <int dim, typename T>
3824  template <class Archive>
3825  void
3827  const unsigned int /*version*/)
3828  {
3829  std::size_t n_cells;
3830  ar & n_cells;
3831  cell_ids.clear();
3832  cell_ids.reserve(n_cells);
3833  for (unsigned int c = 0; c < n_cells; ++c)
3834  {
3835  CellId::binary_type value;
3836  ar & value;
3837  cell_ids.emplace_back(value);
3838  }
3839  ar &data;
3840  }
3841 
3842 
3843 
3844  template <typename DataType, typename MeshType>
3845  void
3847  const MeshType & mesh,
3848  const std::function<std_cxx17::optional<DataType>(
3849  const typename MeshType::active_cell_iterator &)> &pack,
3850  const std::function<void(const typename MeshType::active_cell_iterator &,
3851  const DataType &)> & unpack)
3852  {
3853 # ifndef DEAL_II_WITH_MPI
3854  (void)mesh;
3855  (void)pack;
3856  (void)unpack;
3857  Assert(false,
3858  ExcMessage(
3859  "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
3860 # else
3861  constexpr int dim = MeshType::dimension;
3862  constexpr int spacedim = MeshType::space_dimension;
3863  auto tria = static_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3864  &mesh.get_triangulation());
3865  Assert(
3866  tria != nullptr,
3867  ExcMessage(
3868  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3869 
3870  // map neighbor_id -> data_buffer where we accumulate the data to send
3871  using DestinationToBufferMap =
3872  std::map<::types::subdomain_id,
3874  DestinationToBufferMap destination_to_data_buffer_map;
3875 
3876  std::map<unsigned int, std::set<::types::subdomain_id>>
3877  vertices_with_ghost_neighbors =
3878  tria->compute_vertices_with_ghost_neighbors();
3879 
3880  for (const auto &cell : tria->active_cell_iterators())
3881  if (cell->is_locally_owned())
3882  {
3883  std::set<::types::subdomain_id> send_to;
3884  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
3885  ++v)
3886  {
3887  const std::map<unsigned int,
3888  std::set<::types::subdomain_id>>::
3889  const_iterator neighbor_subdomains_of_vertex =
3890  vertices_with_ghost_neighbors.find(cell->vertex_index(v));
3891 
3892  if (neighbor_subdomains_of_vertex ==
3893  vertices_with_ghost_neighbors.end())
3894  continue;
3895 
3896  Assert(neighbor_subdomains_of_vertex->second.size() != 0,
3897  ExcInternalError());
3898 
3899  send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
3900  neighbor_subdomains_of_vertex->second.end());
3901  }
3902 
3903  if (send_to.size() > 0)
3904  {
3905  // this cell's data needs to be sent to someone
3906  typename MeshType::active_cell_iterator mesh_it(tria,
3907  cell->level(),
3908  cell->index(),
3909  &mesh);
3910 
3911  const std_cxx17::optional<DataType> data = pack(mesh_it);
3912 
3913  if (data)
3914  {
3915  const CellId cellid = cell->id();
3916 
3917  for (const auto subdomain : send_to)
3918  {
3919  // find the data buffer for proc "subdomain" if it exists
3920  // or create an empty one otherwise
3921  typename DestinationToBufferMap::iterator p =
3922  destination_to_data_buffer_map
3923  .insert(std::make_pair(
3924  subdomain, CellDataTransferBuffer<dim, DataType>()))
3925  .first;
3926 
3927  p->second.cell_ids.emplace_back(cellid);
3928  p->second.data.emplace_back(*data);
3929  }
3930  }
3931  }
3932  }
3933 
3934 
3935  // 2. send our messages
3936  std::set<::types::subdomain_id> ghost_owners = tria->ghost_owners();
3937  const unsigned int n_ghost_owners = ghost_owners.size();
3938  std::vector<std::vector<char>> sendbuffers(n_ghost_owners);
3939  std::vector<MPI_Request> requests(n_ghost_owners);
3940 
3941  unsigned int idx = 0;
3942  for (auto it = ghost_owners.begin(); it != ghost_owners.end(); ++it, ++idx)
3943  {
3944  CellDataTransferBuffer<dim, DataType> &data =
3945  destination_to_data_buffer_map[*it];
3946 
3947  // pack all the data into the buffer for this recipient and send it.
3948  // keep data around till we can make sure that the packet has been
3949  // received
3950  sendbuffers[idx] = Utilities::pack(data, /*enable_compression*/ false);
3951  const int ierr = MPI_Isend(sendbuffers[idx].data(),
3952  sendbuffers[idx].size(),
3953  MPI_BYTE,
3954  *it,
3955  786,
3956  tria->get_communicator(),
3957  &requests[idx]);
3958  AssertThrowMPI(ierr);
3959  }
3960 
3961  // 3. receive messages
3962  std::vector<char> receive;
3963  for (unsigned int idx = 0; idx < n_ghost_owners; ++idx)
3964  {
3965  MPI_Status status;
3966  int len;
3967  int ierr =
3968  MPI_Probe(MPI_ANY_SOURCE, 786, tria->get_communicator(), &status);
3969  AssertThrowMPI(ierr);
3970  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3971  AssertThrowMPI(ierr);
3972 
3973  receive.resize(len);
3974 
3975  char *ptr = receive.data();
3976  ierr = MPI_Recv(ptr,
3977  len,
3978  MPI_BYTE,
3979  status.MPI_SOURCE,
3980  status.MPI_TAG,
3981  tria->get_communicator(),
3982  &status);
3983  AssertThrowMPI(ierr);
3984 
3985  auto cellinfo =
3986  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
3987  receive, /*enable_compression*/ false);
3988 
3989  DataType *data = cellinfo.data.data();
3990  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
3991  {
3993  tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
3994 
3995  const typename MeshType::active_cell_iterator cell(
3996  tria, tria_cell->level(), tria_cell->index(), &mesh);
3997 
3998  unpack(cell, *data);
3999  }
4000  }
4001 
4002  // make sure that all communication is finished
4003  // when we leave this function.
4004  if (requests.size())
4005  {
4006  const int ierr =
4007  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4008  AssertThrowMPI(ierr);
4009  }
4010 # endif // DEAL_II_WITH_MPI
4011  }
4012 } // namespace GridTools
4013 
4014 # endif
4015 
4016 DEAL_II_NAMESPACE_CLOSE
4017 
4018 /*---------------------------- grid_tools.h ---------------------------*/
4019 /* end of #ifndef dealii_grid_tools_h */
4020 #endif
4021 /*---------------------------- 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:4055
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:3833
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:267
static const unsigned int invalid_unsigned_int
Definition: types.h:187
unsigned int manifold_id
Definition: types.h:137
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:1181
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3808
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:2182
Definition: tria.h:136
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:76
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true)
Definition: grid_tools.cc:1214
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:5427
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:4023
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:494
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position)
Definition: grid_tools.cc:1952
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:4125
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:134
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12111
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2280
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:5193
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11939
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:2331
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:1416
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:5585
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:78
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:3274
std::bitset< 3 > orientation
Definition: grid_tools.h:2438
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3127
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3246
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3024
double compute_maximum_aspect_ratio(const Triangulation< dim > &triangulation, const Mapping< dim > &mapping, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:479
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:2797
__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:738
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:2726
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:12005
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:636
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
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:2761
#define Assert(cond, exc)
Definition: exceptions.h:1407
Signals signals
Definition: tria.h:2411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:2961
Abstract base class for mapping classes.
Definition: mapping.h:302
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:2035
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:4092
DEAL_II_CONSTEXPR SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
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:3907
void rotate(const double angle, Triangulation< 2 > &triangulation)
Definition: grid_tools.cc:1013
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:3769
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1170
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:68
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:1725
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:4482
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:1592
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2678
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:2452
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:3935
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:4446
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3190
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)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
Vector< double > compute_aspect_ratio_of_cells(const Triangulation< dim > &triangulation, const Mapping< dim > &mapping, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:420
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
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)
unsigned int boundary_id
Definition: types.h:125
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:3171
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:5446
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:839