Reference documentation for deal.II version Git cf38156a17 2020-04-01 20:09:03 -0400
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
386  template <int dim, int spacedim>
387  void
388  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
389  std::vector<CellData<dim>> & cells,
390  SubCellData & subcelldata,
391  std::vector<unsigned int> & considered_vertices,
392  const double tol = 1e-12);
393 
399 
434  template <int dim, typename Transformation, int spacedim>
435  void
436  transform(const Transformation & transformation,
437  Triangulation<dim, spacedim> &triangulation);
438 
444  template <int dim, int spacedim>
445  void
446  shift(const Tensor<1, spacedim> & shift_vector,
447  Triangulation<dim, spacedim> &triangulation);
448 
449 
457  void
458  rotate(const double angle, Triangulation<2> &triangulation);
459 
472  template <int dim>
473  void
474  rotate(const double angle,
475  const unsigned int axis,
476  Triangulation<dim, 3> &triangulation);
477 
535  template <int dim>
536  void
537  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
538  Triangulation<dim> & tria,
539  const Function<dim, double> *coefficient = nullptr,
540  const bool solve_for_absolute_positions = false);
541 
547  template <int dim, int spacedim>
548  std::map<unsigned int, Point<spacedim>>
550 
558  template <int dim, int spacedim>
559  void
560  scale(const double scaling_factor,
561  Triangulation<dim, spacedim> &triangulation);
562 
573  template <int dim, int spacedim>
574  void
575  distort_random(const double factor,
576  Triangulation<dim, spacedim> &triangulation,
577  const bool keep_boundary = true);
578 
614  template <int dim, int spacedim>
615  void
617  const bool isotropic = false,
618  const unsigned int max_iterations = 100);
619 
646  template <int dim, int spacedim>
647  void
649  const double max_ratio = 1.6180339887,
650  const unsigned int max_iterations = 5);
651 
743  template <int dim, int spacedim>
744  void
746  const double limit_angle_fraction = .75);
747 
753 
805  template <int dim, int spacedim>
806 # ifndef DOXYGEN
807  std::tuple<
808  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
809  std::vector<std::vector<Point<dim>>>,
810  std::vector<std::vector<unsigned int>>>
811 # else
812  return_type
813 # endif
815  const Cache<dim, spacedim> & cache,
816  const std::vector<Point<spacedim>> &points,
818  &cell_hint =
820 
849  template <int dim, int spacedim>
850 # ifndef DOXYGEN
851  std::tuple<
852  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
853  std::vector<std::vector<Point<dim>>>,
854  std::vector<std::vector<unsigned int>>,
855  std::vector<unsigned int>>
856 # else
857  return_type
858 # endif
860  const Cache<dim, spacedim> & cache,
861  const std::vector<Point<spacedim>> &points,
863  &cell_hint =
865 
930  template <int dim, int spacedim>
931 # ifndef DOXYGEN
932  std::tuple<
933  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
934  std::vector<std::vector<Point<dim>>>,
935  std::vector<std::vector<unsigned int>>,
936  std::vector<std::vector<Point<spacedim>>>,
937  std::vector<std::vector<unsigned int>>>
938 # else
939  return_type
940 # endif
942  const GridTools::Cache<dim, spacedim> & cache,
943  const std::vector<Point<spacedim>> & local_points,
944  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes);
945 
984  template <int dim, int spacedim>
985  std::map<unsigned int, Point<spacedim>>
987  const Mapping<dim, spacedim> & mapping =
989 
1001  template <int spacedim>
1002  unsigned int
1003  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1004  const Point<spacedim> & p);
1005 
1032  template <int dim, template <int, int> class MeshType, int spacedim>
1033  unsigned int
1034  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1035  const Point<spacedim> & p,
1036  const std::vector<bool> & marked_vertices = {});
1037 
1063  template <int dim, template <int, int> class MeshType, int spacedim>
1064  unsigned int
1066  const MeshType<dim, spacedim> &mesh,
1067  const Point<spacedim> & p,
1068  const std::vector<bool> & marked_vertices = {});
1069 
1070 
1095  template <int dim, template <int, int> class MeshType, int spacedim>
1096 # ifndef _MSC_VER
1097  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1098 # else
1099  std::vector<
1100  typename ::internal::
1101  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1102 # endif
1103  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1104  const unsigned int vertex_index);
1105 
1106 
1131  template <int dim, template <int, int> class MeshType, int spacedim>
1132 # ifndef _MSC_VER
1133  typename MeshType<dim, spacedim>::active_cell_iterator
1134 # else
1135  typename ::internal::ActiveCellIterator<dim,
1136  spacedim,
1137  MeshType<dim, spacedim>>::type
1138 # endif
1139  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1140  const Point<spacedim> & p,
1141  const std::vector<bool> &marked_vertices = {});
1142 
1230  template <int dim, template <int, int> class MeshType, int spacedim>
1231 # ifndef _MSC_VER
1232  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1233 # else
1234  std::pair<typename ::internal::
1235  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1236  Point<dim>>
1237 # endif
1239  const MeshType<dim, spacedim> &mesh,
1240  const Point<spacedim> & p,
1241  const std::vector<bool> &marked_vertices = {});
1242 
1254  template <int dim, template <int, int> class MeshType, int spacedim>
1255 # ifndef _MSC_VER
1256  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1257 # else
1258  std::pair<typename ::internal::
1259  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1260  Point<dim>>
1261 # endif
1263  const Mapping<dim, spacedim> & mapping,
1264  const MeshType<dim, spacedim> &mesh,
1265  const Point<spacedim> & p,
1266  const std::vector<
1267  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1268  & vertex_to_cell_map,
1269  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1270  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1271  typename MeshType<dim, spacedim>::active_cell_iterator(),
1272  const std::vector<bool> & marked_vertices = {},
1273  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1274  RTree<std::pair<Point<spacedim>, unsigned int>>{});
1275 
1297  template <int dim, int spacedim>
1298  std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator,
1299  Point<dim>>
1301  const hp::MappingCollection<dim, spacedim> &mapping,
1302  const hp::DoFHandler<dim, spacedim> & mesh,
1303  const Point<spacedim> & p);
1304 
1311  template <int dim, int spacedim>
1312  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1313  Point<dim>>
1315  const Cache<dim, spacedim> &cache,
1316  const Point<spacedim> & p,
1319  const std::vector<bool> &marked_vertices = {});
1320 
1338  template <int dim, template <int, int> class MeshType, int spacedim>
1339 # ifndef _MSC_VER
1340  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1341  Point<dim>>>
1342 # else
1343  std::vector<std::pair<
1344  typename ::internal::
1345  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1346  Point<dim>>>
1347 # endif
1349  const Mapping<dim, spacedim> & mapping,
1350  const MeshType<dim, spacedim> &mesh,
1351  const Point<spacedim> & p,
1352  const double tolerance = 1e-10,
1353  const std::vector<bool> & marked_vertices = {});
1354 
1376  template <class MeshType>
1377  std::vector<typename MeshType::active_cell_iterator>
1378  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1379 
1404  template <class MeshType>
1405  void
1407  const typename MeshType::active_cell_iterator & cell,
1408  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1409 
1462  template <class MeshType>
1463  std::vector<typename MeshType::active_cell_iterator>
1465  const MeshType &mesh,
1466  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1467  &predicate);
1468 
1469 
1477  template <class MeshType>
1478  std::vector<typename MeshType::cell_iterator>
1480  const MeshType &mesh,
1481  const std::function<bool(const typename MeshType::cell_iterator &)>
1482  & predicate,
1483  const unsigned int level);
1484 
1485 
1501  template <class MeshType>
1502  std::vector<typename MeshType::active_cell_iterator>
1503  compute_ghost_cell_halo_layer(const MeshType &mesh);
1504 
1556  template <class MeshType>
1557  std::vector<typename MeshType::active_cell_iterator>
1559  const MeshType &mesh,
1560  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1561  & predicate,
1562  const double layer_thickness);
1563 
1589  template <class MeshType>
1590  std::vector<typename MeshType::active_cell_iterator>
1591  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1592  const double layer_thickness);
1593 
1609  template <class MeshType>
1610  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1612  const MeshType &mesh,
1613  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1614  &predicate);
1615 
1668  template <class MeshType>
1669  std::vector<BoundingBox<MeshType::space_dimension>>
1671  const MeshType &mesh,
1672  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1673  & predicate,
1674  const unsigned int refinement_level = 0,
1675  const bool allow_merge = false,
1676  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1677 
1707  template <int spacedim>
1708 # ifndef DOXYGEN
1709  std::tuple<std::vector<std::vector<unsigned int>>,
1710  std::map<unsigned int, unsigned int>,
1711  std::map<unsigned int, std::vector<unsigned int>>>
1712 # else
1713  return_type
1714 # endif
1716  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1717  const std::vector<Point<spacedim>> & points);
1718 
1719 
1754  template <int spacedim>
1755 # ifndef DOXYGEN
1756  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1757  std::map<unsigned int, unsigned int>,
1758  std::map<unsigned int, std::vector<unsigned int>>>
1759 # else
1760  return_type
1761 # endif
1763  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1764  const std::vector<Point<spacedim>> & points);
1765 
1766 
1775  template <int dim, int spacedim>
1776  std::vector<
1777  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1778  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1779 
1794  template <int dim, int spacedim>
1795  std::vector<std::vector<Tensor<1, spacedim>>>
1797  const Triangulation<dim, spacedim> &mesh,
1798  const std::vector<
1800  &vertex_to_cells);
1801 
1802 
1809  template <int dim, int spacedim>
1810  unsigned int
1813  const Point<spacedim> &position);
1814 
1826  template <int dim, int spacedim>
1827  std::map<unsigned int, types::global_vertex_index>
1830 
1844  template <int dim, int spacedim>
1845  std::pair<unsigned int, double>
1848 
1854 
1863  template <int dim, int spacedim>
1864  void
1866  const Triangulation<dim, spacedim> &triangulation,
1867  DynamicSparsityPattern & connectivity);
1868 
1877  template <int dim, int spacedim>
1878  void
1880  const Triangulation<dim, spacedim> &triangulation,
1881  DynamicSparsityPattern & connectivity);
1882 
1891  template <int dim, int spacedim>
1892  void
1894  const Triangulation<dim, spacedim> &triangulation,
1895  const unsigned int level,
1896  DynamicSparsityPattern & connectivity);
1897 
1918  template <int dim, int spacedim>
1919  void
1920  partition_triangulation(const unsigned int n_partitions,
1921  Triangulation<dim, spacedim> & triangulation,
1922  const SparsityTools::Partitioner partitioner =
1924 
1935  template <int dim, int spacedim>
1936  void
1937  partition_triangulation(const unsigned int n_partitions,
1938  const std::vector<unsigned int> &cell_weights,
1939  Triangulation<dim, spacedim> & triangulation,
1940  const SparsityTools::Partitioner partitioner =
1942 
1988  template <int dim, int spacedim>
1989  void
1990  partition_triangulation(const unsigned int n_partitions,
1991  const SparsityPattern & cell_connection_graph,
1992  Triangulation<dim, spacedim> &triangulation,
1993  const SparsityTools::Partitioner partitioner =
1995 
2006  template <int dim, int spacedim>
2007  void
2008  partition_triangulation(const unsigned int n_partitions,
2009  const std::vector<unsigned int> &cell_weights,
2010  const SparsityPattern & cell_connection_graph,
2011  Triangulation<dim, spacedim> &triangulation,
2012  const SparsityTools::Partitioner partitioner =
2014 
2029  template <int dim, int spacedim>
2030  void
2031  partition_triangulation_zorder(const unsigned int n_partitions,
2032  Triangulation<dim, spacedim> &triangulation,
2033  const bool group_siblings = true);
2034 
2046  template <int dim, int spacedim>
2047  void
2049 
2060  template <int dim, int spacedim>
2061  void
2062  get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
2063  std::vector<types::subdomain_id> & subdomain);
2064 
2079  template <int dim, int spacedim>
2080  unsigned int
2082  const Triangulation<dim, spacedim> &triangulation,
2083  const types::subdomain_id subdomain);
2084 
2085 
2115  template <int dim, int spacedim>
2116  std::vector<bool>
2118 
2124 
2158  template <typename MeshType>
2159  std::list<std::pair<typename MeshType::cell_iterator,
2160  typename MeshType::cell_iterator>>
2161  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2162 
2172  template <int dim, int spacedim>
2173  bool
2175  const Triangulation<dim, spacedim> &mesh_2);
2176 
2186  template <typename MeshType>
2187  bool
2188  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2189 
2195 
2211  template <int dim, int spacedim>
2215  & distorted_cells,
2216  Triangulation<dim, spacedim> &triangulation);
2217 
2218 
2219 
2228 
2229 
2269  template <class MeshType>
2270  std::vector<typename MeshType::active_cell_iterator>
2271  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2272 
2273 
2297  template <class Container>
2298  std::vector<typename Container::cell_iterator>
2300  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2301 
2370  template <class Container>
2371  void
2373  const std::vector<typename Container::active_cell_iterator> &patch,
2375  &local_triangulation,
2376  std::map<
2377  typename Triangulation<Container::dimension,
2378  Container::space_dimension>::active_cell_iterator,
2379  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2380 
2417  template <class DoFHandlerType>
2418  std::map<types::global_dof_index,
2419  std::vector<typename DoFHandlerType::active_cell_iterator>>
2420  get_dof_to_support_patch_map(DoFHandlerType &dof_handler);
2421 
2422 
2429 
2435  template <typename CellIterator>
2437  {
2441  CellIterator cell[2];
2442 
2447  unsigned int face_idx[2];
2448 
2454  std::bitset<3> orientation;
2455 
2469  };
2470 
2471 
2537  template <typename FaceIterator>
2538  bool orthogonal_equality(
2539  std::bitset<3> & orientation,
2540  const FaceIterator & face1,
2541  const FaceIterator & face2,
2542  const int direction,
2545  const FullMatrix<double> &matrix = FullMatrix<double>());
2546 
2547 
2551  template <typename FaceIterator>
2552  bool
2554  const FaceIterator & face1,
2555  const FaceIterator & face2,
2556  const int direction,
2559  const FullMatrix<double> &matrix = FullMatrix<double>());
2560 
2561 
2620  template <typename MeshType>
2621  void
2623  const MeshType & mesh,
2624  const types::boundary_id b_id1,
2625  const types::boundary_id b_id2,
2626  const int direction,
2628  & matched_pairs,
2629  const Tensor<1, MeshType::space_dimension> &offset =
2631  const FullMatrix<double> &matrix = FullMatrix<double>());
2632 
2633 
2658  template <typename MeshType>
2659  void
2661  const MeshType & mesh,
2662  const types::boundary_id b_id,
2663  const int direction,
2665  & matched_pairs,
2666  const ::Tensor<1, MeshType::space_dimension> &offset =
2668  const FullMatrix<double> &matrix = FullMatrix<double>());
2669 
2675 
2698  template <int dim, int spacedim>
2699  void
2701  const bool reset_boundary_ids = false);
2702 
2726  template <int dim, int spacedim>
2727  void
2729  const std::vector<types::boundary_id> &src_boundary_ids,
2730  const std::vector<types::manifold_id> &dst_manifold_ids,
2732  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2733 
2765  template <int dim, int spacedim>
2766  void
2768  const bool compute_face_ids = false);
2769 
2796  template <int dim, int spacedim>
2797  void
2800  const std::function<types::manifold_id(
2801  const std::set<types::manifold_id> &)> &disambiguation_function =
2802  [](const std::set<types::manifold_id> &manifold_ids) {
2803  if (manifold_ids.size() == 1)
2804  return *manifold_ids.begin();
2805  else
2807  },
2808  bool overwrite_only_flat_manifold_ids = true);
2893  template <typename DataType, typename MeshType>
2894  void
2895  exchange_cell_data_to_ghosts(
2896  const MeshType & mesh,
2897  const std::function<std_cxx17::optional<DataType>(
2898  const typename MeshType::active_cell_iterator &)> &pack,
2899  const std::function<void(const typename MeshType::active_cell_iterator &,
2900  const DataType &)> & unpack);
2901 
2902  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2903  * boxes @p local_bboxes.
2904  *
2905  * This function is meant to exchange bounding boxes describing the locally
2906  * owned cells in a distributed triangulation obtained with the function
2907  * GridTools::compute_mesh_predicate_bounding_box .
2908  *
2909  * The output vector's size is the number of processes of the MPI
2910  * communicator:
2911  * its i-th entry contains the vector @p local_bboxes of the i-th process.
2912  */
2913  template <int spacedim>
2914  std::vector<std::vector<BoundingBox<spacedim>>>
2915  exchange_local_bounding_boxes(
2916  const std::vector<BoundingBox<spacedim>> &local_bboxes,
2917  MPI_Comm mpi_communicator);
2918 
2953  template <int spacedim>
2954  RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
2956  const std::vector<BoundingBox<spacedim>> &local_description,
2957  MPI_Comm mpi_communicator);
2958 
2978  template <int dim, int spacedim>
2979  void
2981  const Triangulation<dim, spacedim> & tria,
2982  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
2983  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
2984 
2997  template <int dim, typename T>
2999  {
3003  std::vector<CellId> cell_ids;
3004 
3008  std::vector<T> data;
3009 
3017  template <class Archive>
3018  void
3019  save(Archive &ar, const unsigned int version) const;
3020 
3025  template <class Archive>
3026  void
3027  load(Archive &ar, const unsigned int version);
3028 
3029  BOOST_SERIALIZATION_SPLIT_MEMBER()
3030  };
3031 
3036 
3041  int,
3042  << "The number of partitions you gave is " << arg1
3043  << ", but must be greater than zero.");
3048  int,
3049  << "The subdomain id " << arg1
3050  << " has no cells associated with it.");
3055 
3060  double,
3061  << "The scaling factor must be positive, but it is " << arg1
3062  << ".");
3066  template <int N>
3068  Point<N>,
3069  << "The point <" << arg1
3070  << "> could not be found inside any of the "
3071  << "coarse grid cells.");
3075  template <int N>
3077  Point<N>,
3078  << "The point <" << arg1
3079  << "> could not be found inside any of the "
3080  << "subcells of a coarse grid cell.");
3081 
3086  unsigned int,
3087  << "The given vertex with index " << arg1
3088  << " is not used in the given triangulation.");
3089 
3090 
3093 } /*namespace GridTools*/
3094 
3095 
3096 
3097 /* ----------------- Template function --------------- */
3098 
3099 # ifndef DOXYGEN
3100 
3101 namespace GridTools
3102 {
3103  // declare specializations
3104  template <>
3105  double
3106  cell_measure<1>(const std::vector<Point<1>> &,
3107  const unsigned int (&)[GeometryInfo<1>::vertices_per_cell]);
3108 
3109  template <>
3110  double
3111  cell_measure<2>(const std::vector<Point<2>> &,
3112  const unsigned int (&)[GeometryInfo<2>::vertices_per_cell]);
3113 
3114  template <>
3115  double
3116  cell_measure<3>(const std::vector<Point<3>> &,
3117  const unsigned int (&)[GeometryInfo<3>::vertices_per_cell]);
3118 
3119 
3120 
3121  template <int dim, typename T>
3122  double
3123  cell_measure(const T &, ...)
3124  {
3125  Assert(false, ExcNotImplemented());
3126  return std::numeric_limits<double>::quiet_NaN();
3127  }
3128 
3129 
3130 
3131  template <int dim, typename Predicate, int spacedim>
3132  void
3133  transform(const Predicate & predicate,
3134  Triangulation<dim, spacedim> &triangulation)
3135  {
3136  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3137 
3138  // loop over all active cells, and
3139  // transform those vertices that
3140  // have not yet been touched. note
3141  // that we get to all vertices in
3142  // the triangulation by only
3143  // visiting the active cells.
3145  cell = triangulation.begin_active(),
3146  endc = triangulation.end();
3147  for (; cell != endc; ++cell)
3148  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
3149  if (treated_vertices[cell->vertex_index(v)] == false)
3150  {
3151  // transform this vertex
3152  cell->vertex(v) = predicate(cell->vertex(v));
3153  // and mark it as treated
3154  treated_vertices[cell->vertex_index(v)] = true;
3155  };
3156 
3157 
3158  // now fix any vertices on hanging nodes so that we don't create any holes
3159  if (dim == 2)
3160  {
3162  cell = triangulation.begin_active(),
3163  endc = triangulation.end();
3164  for (; cell != endc; ++cell)
3165  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3166  if (cell->face(face)->has_children() &&
3167  !cell->face(face)->at_boundary())
3168  {
3169  // this line has children
3170  cell->face(face)->child(0)->vertex(1) =
3171  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3172  2;
3173  }
3174  }
3175  else if (dim == 3)
3176  {
3178  cell = triangulation.begin_active(),
3179  endc = triangulation.end();
3180  for (; cell != endc; ++cell)
3181  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3182  if (cell->face(face)->has_children() &&
3183  !cell->face(face)->at_boundary())
3184  {
3185  // this face has hanging nodes
3186  cell->face(face)->child(0)->vertex(1) =
3187  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3188  2.0;
3189  cell->face(face)->child(0)->vertex(2) =
3190  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3191  2.0;
3192  cell->face(face)->child(1)->vertex(3) =
3193  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3194  2.0;
3195  cell->face(face)->child(2)->vertex(3) =
3196  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3197  2.0;
3198 
3199  // center of the face
3200  cell->face(face)->child(0)->vertex(3) =
3201  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3202  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3203  4.0;
3204  }
3205  }
3206 
3207  // Make sure FEValues notices that the mesh has changed
3208  triangulation.signals.mesh_movement();
3209  }
3210 
3211 
3212 
3213  template <class MeshType>
3214  std::vector<typename MeshType::active_cell_iterator>
3215  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3216  {
3217  std::vector<typename MeshType::active_cell_iterator> child_cells;
3218 
3219  if (cell->has_children())
3220  {
3221  for (unsigned int child = 0; child < cell->n_children(); ++child)
3222  if (cell->child(child)->has_children())
3223  {
3224  const std::vector<typename MeshType::active_cell_iterator>
3225  children = get_active_child_cells<MeshType>(cell->child(child));
3226  child_cells.insert(child_cells.end(),
3227  children.begin(),
3228  children.end());
3229  }
3230  else
3231  child_cells.push_back(cell->child(child));
3232  }
3233 
3234  return child_cells;
3235  }
3236 
3237 
3238 
3239  template <class MeshType>
3240  void
3242  const typename MeshType::active_cell_iterator & cell,
3243  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3244  {
3245  active_neighbors.clear();
3246  for (const unsigned int n :
3248  if (!cell->at_boundary(n))
3249  {
3250  if (MeshType::dimension == 1)
3251  {
3252  // check children of neighbor. note
3253  // that in 1d children of the neighbor
3254  // may be further refined. In 1d the
3255  // case is simple since we know what
3256  // children bound to the present cell
3257  typename MeshType::cell_iterator neighbor_child =
3258  cell->neighbor(n);
3259  if (!neighbor_child->is_active())
3260  {
3261  while (neighbor_child->has_children())
3262  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3263 
3264  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3265  ExcInternalError());
3266  }
3267  active_neighbors.push_back(neighbor_child);
3268  }
3269  else
3270  {
3271  if (cell->face(n)->has_children())
3272  // this neighbor has children. find
3273  // out which border to the present
3274  // cell
3275  for (unsigned int c = 0;
3276  c < cell->face(n)->number_of_children();
3277  ++c)
3278  active_neighbors.push_back(
3279  cell->neighbor_child_on_subface(n, c));
3280  else
3281  {
3282  // the neighbor must be active
3283  // himself
3284  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3285  active_neighbors.push_back(cell->neighbor(n));
3286  }
3287  }
3288  }
3289  }
3290 
3291 
3292 
3293  namespace internal
3294  {
3295  namespace ProjectToObject
3296  {
3309  struct CrossDerivative
3310  {
3311  const unsigned int direction_0;
3312  const unsigned int direction_1;
3313 
3314  CrossDerivative(const unsigned int d0, const unsigned int d1);
3315  };
3316 
3317  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3318  const unsigned int d1)
3319  : direction_0(d0)
3320  , direction_1(d1)
3321  {}
3322 
3323 
3324 
3329  template <typename F>
3330  inline auto
3331  centered_first_difference(const double center,
3332  const double step,
3333  const F &f) -> decltype(f(center) - f(center))
3334  {
3335  return (f(center + step) - f(center - step)) / (2.0 * step);
3336  }
3337 
3338 
3339 
3344  template <typename F>
3345  inline auto
3346  centered_second_difference(const double center,
3347  const double step,
3348  const F &f) -> decltype(f(center) - f(center))
3349  {
3350  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3351  (step * step);
3352  }
3353 
3354 
3355 
3365  template <int structdim, typename F>
3366  inline auto
3367  cross_stencil(
3368  const CrossDerivative cross_derivative,
3370  const double step,
3371  const F &f) -> decltype(f(center) - f(center))
3372  {
3374  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3375  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3376  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3377  1.0 / 3.0 * f(center - simplex_vector) +
3378  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3379  step;
3380  }
3381 
3382 
3383 
3390  template <int spacedim, int structdim, typename F>
3391  inline double
3392  gradient_entry(
3393  const unsigned int row_n,
3394  const unsigned int dependent_direction,
3395  const Point<spacedim> &p0,
3397  const double step,
3398  const F & f)
3399  {
3401  dependent_direction <
3403  ExcMessage("This function assumes that the last weight is a "
3404  "dependent variable (and hence we cannot take its "
3405  "derivative directly)."));
3406  Assert(row_n != dependent_direction,
3407  ExcMessage(
3408  "We cannot differentiate with respect to the variable "
3409  "that is assumed to be dependent."));
3410 
3411  const Point<spacedim> manifold_point = f(center);
3412  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3413  {row_n, dependent_direction}, center, step, f);
3414  double entry = 0.0;
3415  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3416  entry +=
3417  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3418  return entry;
3419  }
3420 
3426  template <typename Iterator, int spacedim, int structdim>
3428  project_to_d_linear_object(const Iterator & object,
3429  const Point<spacedim> &trial_point)
3430  {
3431  // let's look at this for simplicity for a quad (structdim==2) in a
3432  // space with spacedim>2 (notate trial_point by y): all points on the
3433  // surface are given by
3434  // x(\xi) = sum_i v_i phi_x(\xi)
3435  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3436  // reference coordinates of the quad. so what we are trying to do is
3437  // find a point x on the surface that is closest to the point y. there
3438  // are different ways to solve this problem, but in the end it's a
3439  // nonlinear problem and we have to find reference coordinates \xi so
3440  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3441  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3442  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3443  // have to use a Newton method to find the answer. This leads to the
3444  // following formulation of Newton steps:
3445  //
3446  // Given \xi_k, find \delta\xi_k so that
3447  // H_k \delta\xi_k = - F_k
3448  // where H_k is an approximation to the second derivatives of J at
3449  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3450  // number of times until the right hand side is small enough. As a
3451  // stopping criterion, we terminate if ||\delta\xi||<eps.
3452  //
3453  // As for the Hessian, the best choice would be
3454  // H_k = J''(\xi_k)
3455  // but we'll opt for the simpler Gauss-Newton form
3456  // H_k = A^T A
3457  // i.e.
3458  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3459  // \partial_n phi_i *
3460  // \partial_m phi_j
3461  // we start at xi=(0.5, 0.5).
3462  Point<structdim> xi;
3463  for (unsigned int d = 0; d < structdim; ++d)
3464  xi[d] = 0.5;
3465 
3466  Point<spacedim> x_k;
3467  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3468  x_k += object->vertex(i) *
3470 
3471  do
3472  {
3474  for (const unsigned int i :
3476  F_k +=
3477  (x_k - trial_point) * object->vertex(i) *
3479  i);
3480 
3482  for (const unsigned int i :
3484  for (const unsigned int j :
3486  {
3489  xi, i),
3491  xi, j));
3492  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3493  }
3494 
3495  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3496  xi += delta_xi;
3497 
3498  x_k = Point<spacedim>();
3499  for (const unsigned int i :
3501  x_k += object->vertex(i) *
3503 
3504  if (delta_xi.norm() < 1e-7)
3505  break;
3506  }
3507  while (true);
3508 
3509  return x_k;
3510  }
3511  } // namespace ProjectToObject
3512  } // namespace internal
3513 
3514 
3515 
3516  namespace internal
3517  {
3518  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3519  // inside the project_to_object function below.
3520  template <int structdim>
3521  inline bool
3522  weights_are_ok(
3524  {
3525  // clang has trouble figuring out structdim here, so define it
3526  // again:
3527  static const std::size_t n_vertices_per_cell =
3529  n_independent_components;
3530  std::array<double, n_vertices_per_cell> copied_weights;
3531  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3532  {
3533  copied_weights[i] = v[i];
3534  if (v[i] < 0.0 || v[i] > 1.0)
3535  return false;
3536  }
3537 
3538  // check the sum: try to avoid some roundoff errors by summing in order
3539  std::sort(copied_weights.begin(), copied_weights.end());
3540  const double sum =
3541  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3542  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3543  }
3544  } // namespace internal
3545 
3546  template <typename Iterator>
3549  const Iterator & object,
3551  {
3552  const int spacedim = Iterator::AccessorType::space_dimension;
3553  const int structdim = Iterator::AccessorType::structure_dimension;
3554 
3555  Point<spacedim> projected_point = trial_point;
3556 
3557  if (structdim >= spacedim)
3558  return projected_point;
3559  else if (structdim == 1 || structdim == 2)
3560  {
3561  using namespace internal::ProjectToObject;
3562  // Try to use the special flat algorithm for quads (this is better
3563  // than the general algorithm in 3D). This does not take into account
3564  // whether projected_point is outside the quad, but we optimize along
3565  // lines below anyway:
3566  const int dim = Iterator::AccessorType::dimension;
3567  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3568  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3569  &manifold) != nullptr)
3570  {
3571  projected_point =
3572  project_to_d_linear_object<Iterator, spacedim, structdim>(
3573  object, trial_point);
3574  }
3575  else
3576  {
3577  // We want to find a point on the convex hull (defined by the
3578  // vertices of the object and the manifold description) that is
3579  // relatively close to the trial point. This has a few issues:
3580  //
3581  // 1. For a general convex hull we are not guaranteed that a unique
3582  // minimum exists.
3583  // 2. The independent variables in the optimization process are the
3584  // weights given to Manifold::get_new_point, which must sum to 1,
3585  // so we cannot use standard finite differences to approximate a
3586  // gradient.
3587  //
3588  // There is not much we can do about 1., but for 2. we can derive
3589  // finite difference stencils that work on a structdim-dimensional
3590  // simplex and rewrite the optimization problem to use those
3591  // instead. Consider the structdim 2 case and let
3592  //
3593  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3594  // c2, c3})
3595  //
3596  // where {c0, c1, c2, c3} are the weights for the four vertices on
3597  // the quadrilateral. We seek to minimize the Euclidean distance
3598  // between F(...) and trial_point. We can solve for c3 in terms of
3599  // the other weights and get, for one coordinate direction
3600  //
3601  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3602  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3603  //
3604  // where we substitute back in for c3 after taking the
3605  // derivative. We can compute a stencil for the cross derivative
3606  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3607  // (and gradient_entry computes the sum over the independent
3608  // variables). Below, we somewhat arbitrarily pick the last
3609  // component as the dependent one.
3610  //
3611  // Since we can now calculate derivatives of the objective
3612  // function we can use gradient descent to minimize it.
3613  //
3614  // Of course, this is much simpler in the structdim = 1 case (we
3615  // could rewrite the projection as a 1D optimization problem), but
3616  // to reduce the potential for bugs we use the same code in both
3617  // cases.
3618  const double step_size = object->diameter() / 64.0;
3619 
3620  constexpr unsigned int n_vertices_per_cell =
3622 
3623  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3624  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3625  ++vertex_n)
3626  vertices[vertex_n] = object->vertex(vertex_n);
3627 
3628  auto get_point_from_weights =
3629  [&](const Tensor<1, n_vertices_per_cell> &weights)
3630  -> Point<spacedim> {
3631  return object->get_manifold().get_new_point(
3632  make_array_view(vertices.begin(), vertices.end()),
3633  make_array_view(weights.begin_raw(), weights.end_raw()));
3634  };
3635 
3636  // pick the initial weights as (normalized) inverse distances from
3637  // the trial point:
3638  Tensor<1, n_vertices_per_cell> guess_weights;
3639  double guess_weights_sum = 0.0;
3640  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3641  ++vertex_n)
3642  {
3643  const double distance =
3644  vertices[vertex_n].distance(trial_point);
3645  if (distance == 0.0)
3646  {
3647  guess_weights = 0.0;
3648  guess_weights[vertex_n] = 1.0;
3649  guess_weights_sum = 1.0;
3650  break;
3651  }
3652  else
3653  {
3654  guess_weights[vertex_n] = 1.0 / distance;
3655  guess_weights_sum += guess_weights[vertex_n];
3656  }
3657  }
3658  guess_weights /= guess_weights_sum;
3659  Assert(internal::weights_are_ok<structdim>(guess_weights),
3660  ExcInternalError());
3661 
3662  // The optimization algorithm consists of two parts:
3663  //
3664  // 1. An outer loop where we apply the gradient descent algorithm.
3665  // 2. An inner loop where we do a line search to find the optimal
3666  // length of the step one should take in the gradient direction.
3667  //
3668  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3669  {
3670  const unsigned int dependent_direction =
3671  n_vertices_per_cell - 1;
3672  Tensor<1, n_vertices_per_cell> current_gradient;
3673  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3674  ++row_n)
3675  {
3676  if (row_n != dependent_direction)
3677  {
3678  current_gradient[row_n] =
3679  gradient_entry<spacedim, structdim>(
3680  row_n,
3681  dependent_direction,
3682  trial_point,
3683  guess_weights,
3684  step_size,
3685  get_point_from_weights);
3686 
3687  current_gradient[dependent_direction] -=
3688  current_gradient[row_n];
3689  }
3690  }
3691 
3692  // We need to travel in the -gradient direction, as noted
3693  // above, but we may not want to take a full step in that
3694  // direction; instead, guess that we will go -0.5*gradient and
3695  // do quasi-Newton iteration to pick the best multiplier. The
3696  // goal is to find a scalar alpha such that
3697  //
3698  // F(x - alpha g)
3699  //
3700  // is minimized, where g is the gradient and F is the
3701  // objective function. To find the optimal value we find roots
3702  // of the derivative of the objective function with respect to
3703  // alpha by Newton iteration, where we approximate the first
3704  // and second derivatives of F(x - alpha g) with centered
3705  // finite differences.
3706  double gradient_weight = -0.5;
3707  auto gradient_weight_objective_function =
3708  [&](const double gradient_weight_guess) -> double {
3709  return (trial_point -
3710  get_point_from_weights(guess_weights +
3711  gradient_weight_guess *
3712  current_gradient))
3713  .norm_square();
3714  };
3715 
3716  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3717  {
3718  const double update_numerator = centered_first_difference(
3719  gradient_weight,
3720  step_size,
3721  gradient_weight_objective_function);
3722  const double update_denominator =
3723  centered_second_difference(
3724  gradient_weight,
3725  step_size,
3726  gradient_weight_objective_function);
3727 
3728  // avoid division by zero. Note that we limit the gradient
3729  // weight below
3730  if (std::abs(update_denominator) == 0.0)
3731  break;
3732  gradient_weight =
3733  gradient_weight - update_numerator / update_denominator;
3734 
3735  // Put a fairly lenient bound on the largest possible
3736  // gradient (things tend to be locally flat, so the gradient
3737  // itself is usually small)
3738  if (std::abs(gradient_weight) > 10)
3739  {
3740  gradient_weight = -10.0;
3741  break;
3742  }
3743  }
3744 
3745  // It only makes sense to take convex combinations with weights
3746  // between zero and one. If the update takes us outside of this
3747  // region then rescale the update to stay within the region and
3748  // try again
3749  Tensor<1, n_vertices_per_cell> tentative_weights =
3750  guess_weights + gradient_weight * current_gradient;
3751 
3752  double new_gradient_weight = gradient_weight;
3753  for (unsigned int iteration_count = 0; iteration_count < 40;
3754  ++iteration_count)
3755  {
3756  if (internal::weights_are_ok<structdim>(tentative_weights))
3757  break;
3758 
3759  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3760  {
3761  if (tentative_weights[i] < 0.0)
3762  {
3763  tentative_weights -=
3764  (tentative_weights[i] / current_gradient[i]) *
3765  current_gradient;
3766  }
3767  if (tentative_weights[i] < 0.0 ||
3768  1.0 < tentative_weights[i])
3769  {
3770  new_gradient_weight /= 2.0;
3771  tentative_weights =
3772  guess_weights +
3773  new_gradient_weight * current_gradient;
3774  }
3775  }
3776  }
3777 
3778  // the update might still send us outside the valid region, so
3779  // check again and quit if the update is still not valid
3780  if (!internal::weights_are_ok<structdim>(tentative_weights))
3781  break;
3782 
3783  // if we cannot get closer by traveling in the gradient
3784  // direction then quit
3785  if (get_point_from_weights(tentative_weights)
3786  .distance(trial_point) <
3787  get_point_from_weights(guess_weights).distance(trial_point))
3788  guess_weights = tentative_weights;
3789  else
3790  break;
3791  Assert(internal::weights_are_ok<structdim>(guess_weights),
3792  ExcInternalError());
3793  }
3794  Assert(internal::weights_are_ok<structdim>(guess_weights),
3795  ExcInternalError());
3796  projected_point = get_point_from_weights(guess_weights);
3797  }
3798 
3799  // if structdim == 2 and the optimal point is not on the interior then
3800  // we may be able to get a more accurate result by projecting onto the
3801  // lines.
3802  if (structdim == 2)
3803  {
3804  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
3805  line_projections;
3806  for (unsigned int line_n = 0;
3807  line_n < GeometryInfo<structdim>::lines_per_cell;
3808  ++line_n)
3809  {
3810  line_projections[line_n] =
3811  project_to_object(object->line(line_n), trial_point);
3812  }
3813  std::sort(line_projections.begin(),
3814  line_projections.end(),
3815  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
3816  return a.distance(trial_point) <
3817  b.distance(trial_point);
3818  });
3819  if (line_projections[0].distance(trial_point) <
3820  projected_point.distance(trial_point))
3821  projected_point = line_projections[0];
3822  }
3823  }
3824  else
3825  {
3826  Assert(false, ExcNotImplemented());
3827  return projected_point;
3828  }
3829 
3830  return projected_point;
3831  }
3832 
3833 
3834 
3835  template <int dim, typename T>
3836  template <class Archive>
3837  void
3839  const unsigned int /*version*/) const
3840  {
3841  Assert(cell_ids.size() == data.size(),
3842  ExcDimensionMismatch(cell_ids.size(), data.size()));
3843  // archive the cellids in an efficient binary format
3844  const std::size_t n_cells = cell_ids.size();
3845  ar & n_cells;
3846  for (const auto &id : cell_ids)
3847  {
3848  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3849  ar & binary_cell_id;
3850  }
3851 
3852  ar &data;
3853  }
3854 
3855 
3856 
3857  template <int dim, typename T>
3858  template <class Archive>
3859  void
3861  const unsigned int /*version*/)
3862  {
3863  std::size_t n_cells;
3864  ar & n_cells;
3865  cell_ids.clear();
3866  cell_ids.reserve(n_cells);
3867  for (unsigned int c = 0; c < n_cells; ++c)
3868  {
3869  CellId::binary_type value;
3870  ar & value;
3871  cell_ids.emplace_back(value);
3872  }
3873  ar &data;
3874  }
3875 
3876 
3877 
3878  template <typename DataType, typename MeshType>
3879  void
3881  const MeshType & mesh,
3882  const std::function<std_cxx17::optional<DataType>(
3883  const typename MeshType::active_cell_iterator &)> &pack,
3884  const std::function<void(const typename MeshType::active_cell_iterator &,
3885  const DataType &)> & unpack)
3886  {
3887 # ifndef DEAL_II_WITH_MPI
3888  (void)mesh;
3889  (void)pack;
3890  (void)unpack;
3891  Assert(false,
3892  ExcMessage(
3893  "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
3894 # else
3895  constexpr int dim = MeshType::dimension;
3896  constexpr int spacedim = MeshType::space_dimension;
3897  auto tria = static_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3898  &mesh.get_triangulation());
3899  Assert(
3900  tria != nullptr,
3901  ExcMessage(
3902  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3903 
3904  // map neighbor_id -> data_buffer where we accumulate the data to send
3905  using DestinationToBufferMap =
3906  std::map<::types::subdomain_id,
3908  DestinationToBufferMap destination_to_data_buffer_map;
3909 
3910  std::map<unsigned int, std::set<::types::subdomain_id>>
3911  vertices_with_ghost_neighbors =
3912  tria->compute_vertices_with_ghost_neighbors();
3913 
3914  for (const auto &cell : tria->active_cell_iterators())
3915  if (cell->is_locally_owned())
3916  {
3917  std::set<::types::subdomain_id> send_to;
3918  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
3919  {
3920  const std::map<unsigned int,
3921  std::set<::types::subdomain_id>>::
3922  const_iterator neighbor_subdomains_of_vertex =
3923  vertices_with_ghost_neighbors.find(cell->vertex_index(v));
3924 
3925  if (neighbor_subdomains_of_vertex ==
3926  vertices_with_ghost_neighbors.end())
3927  continue;
3928 
3929  Assert(neighbor_subdomains_of_vertex->second.size() != 0,
3930  ExcInternalError());
3931 
3932  send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
3933  neighbor_subdomains_of_vertex->second.end());
3934  }
3935 
3936  if (send_to.size() > 0)
3937  {
3938  // this cell's data needs to be sent to someone
3939  typename MeshType::active_cell_iterator mesh_it(tria,
3940  cell->level(),
3941  cell->index(),
3942  &mesh);
3943 
3944  const std_cxx17::optional<DataType> data = pack(mesh_it);
3945 
3946  if (data)
3947  {
3948  const CellId cellid = cell->id();
3949 
3950  for (const auto subdomain : send_to)
3951  {
3952  // find the data buffer for proc "subdomain" if it exists
3953  // or create an empty one otherwise
3954  typename DestinationToBufferMap::iterator p =
3955  destination_to_data_buffer_map
3956  .insert(std::make_pair(
3957  subdomain, CellDataTransferBuffer<dim, DataType>()))
3958  .first;
3959 
3960  p->second.cell_ids.emplace_back(cellid);
3961  p->second.data.emplace_back(*data);
3962  }
3963  }
3964  }
3965  }
3966 
3967  // Protect the following communcation:
3968  static Utilities::MPI::CollectiveMutex mutex;
3970  tria->get_communicator());
3971 
3972  const int mpi_tag =
3974 
3975  // 2. send our messages
3976  std::set<::types::subdomain_id> ghost_owners = tria->ghost_owners();
3977  const unsigned int n_ghost_owners = ghost_owners.size();
3978  std::vector<std::vector<char>> sendbuffers(n_ghost_owners);
3979  std::vector<MPI_Request> requests(n_ghost_owners);
3980 
3981  unsigned int idx = 0;
3982  for (auto it = ghost_owners.begin(); it != ghost_owners.end(); ++it, ++idx)
3983  {
3984  CellDataTransferBuffer<dim, DataType> &data =
3985  destination_to_data_buffer_map[*it];
3986 
3987  // pack all the data into the buffer for this recipient and send it.
3988  // keep data around till we can make sure that the packet has been
3989  // received
3990  sendbuffers[idx] = Utilities::pack(data, /*enable_compression*/ false);
3991  const int ierr = MPI_Isend(sendbuffers[idx].data(),
3992  sendbuffers[idx].size(),
3993  MPI_BYTE,
3994  *it,
3995  mpi_tag,
3996  tria->get_communicator(),
3997  &requests[idx]);
3998  AssertThrowMPI(ierr);
3999  }
4000 
4001  // 3. receive messages
4002  std::vector<char> receive;
4003  for (unsigned int idx = 0; idx < n_ghost_owners; ++idx)
4004  {
4005  MPI_Status status;
4006  int ierr =
4007  MPI_Probe(MPI_ANY_SOURCE, mpi_tag, tria->get_communicator(), &status);
4008  AssertThrowMPI(ierr);
4009 
4010  int len;
4011  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4012  AssertThrowMPI(ierr);
4013 
4014  receive.resize(len);
4015 
4016  char *ptr = receive.data();
4017  ierr = MPI_Recv(ptr,
4018  len,
4019  MPI_BYTE,
4020  status.MPI_SOURCE,
4021  status.MPI_TAG,
4022  tria->get_communicator(),
4023  &status);
4024  AssertThrowMPI(ierr);
4025 
4026  auto cellinfo =
4027  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4028  receive, /*enable_compression*/ false);
4029 
4030  DataType *data = cellinfo.data.data();
4031  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4032  {
4034  tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
4035 
4036  const typename MeshType::active_cell_iterator cell(
4037  tria, tria_cell->level(), tria_cell->index(), &mesh);
4038 
4039  unpack(cell, *data);
4040  }
4041  }
4042 
4043  // make sure that all communication is finished
4044  // when we leave this function.
4045  if (requests.size())
4046  {
4047  const int ierr =
4048  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4049  AssertThrowMPI(ierr);
4050  }
4051 # endif // DEAL_II_WITH_MPI
4052  }
4053 } // namespace GridTools
4054 
4055 # endif
4056 
4057 DEAL_II_NAMESPACE_CLOSE
4058 
4059 /*---------------------------- grid_tools.h ---------------------------*/
4060 /* end of #ifndef dealii_grid_tools_h */
4061 #endif
4062 /*---------------------------- 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:4068
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:3852
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:270
static const unsigned int invalid_unsigned_int
Definition: types.h:190
unsigned int manifold_id
Definition: types.h:140
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:1235
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3827
return_type guess_point_owner(const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:2221
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:1268
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:5426
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:4036
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:501
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position)
Definition: grid_tools.cc:1992
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:4130
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:12183
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2319
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:5192
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12011
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:2370
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:1442
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:5584
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:3308
std::bitset< 3 > orientation
Definition: grid_tools.h:2454
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3168
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3280
double compute_maximum_aspect_ratio(const Triangulation< dim > &triangulation, const Mapping< dim > &mapping, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:485
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:2836
__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:740
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:2771
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:12077
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:643
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:44
#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:2800
#define Assert(cond, exc)
Definition: exceptions.h:1419
Signals signals
Definition: tria.h:2220
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3003
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:2074
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:4101
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 collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
Definition: grid_tools.cc:5647
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:3920
void rotate(const double angle, Triangulation< 2 > &triangulation)
Definition: grid_tools.cc:1073
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:3788
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1221
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: hp.h:117
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-10, const std::vector< bool > &marked_vertices={})
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:1765
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:77
__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:4481
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:1707
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)
GridTools::exchange_cell_data_to_ghosts():
Definition: mpi_tags.h:60
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:1633
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2723
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:2468
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:3948
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:4445
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3225
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:3063
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:419
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:128
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)
unsigned int count_cells_with_subdomain_association(const Triangulation< dim, spacedim > &triangulation, const types::subdomain_id subdomain)
Definition: grid_tools.cc:3209
unsigned int find_closest_vertex(const std::map< unsigned int, Point< spacedim >> &vertices, const Point< spacedim > &p)
Definition: grid_tools.cc:5445
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:841