Reference documentation for deal.II version Git 497f915867 2021-09-17 22:46:48 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2021 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 
25 
27 
29 
31 
32 # include <deal.II/fe/fe_values.h>
33 # include <deal.II/fe/mapping.h>
34 # include <deal.II/fe/mapping_q1.h>
35 
36 # include <deal.II/grid/manifold.h>
37 # include <deal.II/grid/tria.h>
40 
41 # include <deal.II/hp/dof_handler.h>
42 
44 # include <deal.II/lac/la_vector.h>
48 
49 # include <deal.II/numerics/rtree.h>
50 
52 # include <boost/archive/binary_iarchive.hpp>
53 # include <boost/archive/binary_oarchive.hpp>
54 # include <boost/geometry/index/rtree.hpp>
55 # include <boost/random/mersenne_twister.hpp>
56 # include <boost/serialization/array.hpp>
57 # include <boost/serialization/vector.hpp>
58 
59 # ifdef DEAL_II_WITH_ZLIB
60 # include <boost/iostreams/device/back_inserter.hpp>
61 # include <boost/iostreams/filter/gzip.hpp>
62 # include <boost/iostreams/filtering_stream.hpp>
63 # include <boost/iostreams/stream.hpp>
64 # endif
66 
67 # include <bitset>
68 # include <list>
69 # include <set>
70 
72 
73 // Forward declarations
74 # ifndef DOXYGEN
75 namespace parallel
76 {
77  namespace distributed
78  {
79  template <int, int>
80  class Triangulation;
81  }
82 } // namespace parallel
83 
84 namespace hp
85 {
86  template <int, int>
87  class MappingCollection;
88 }
89 
90 class SparsityPattern;
91 # endif
92 
93 namespace internal
94 {
95  template <int dim, int spacedim, class MeshType>
97  {
98  public:
99 # ifndef _MSC_VER
100  using type = typename MeshType::active_cell_iterator;
101 # else
103 # endif
104  };
105 
106 # ifdef _MSC_VER
107  template <int dim, int spacedim>
108  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
109  {
110  public:
111  using type =
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 =
175  (ReferenceCells::get_hypercube<dim>()
176  .template get_default_linear_mapping<dim, spacedim>()));
177 
188  template <int dim, int spacedim>
189  double
192  const Mapping<dim, spacedim> & mapping =
193  (ReferenceCells::get_hypercube<dim>()
194  .template get_default_linear_mapping<dim, spacedim>()));
195 
206  template <int dim, int spacedim>
207  double
209  const Triangulation<dim, spacedim> &triangulation,
210  const Mapping<dim, spacedim> & mapping =
211  (ReferenceCells::get_hypercube<dim>()
212  .template get_default_linear_mapping<dim, spacedim>()));
213 
225  template <int dim>
226  DEAL_II_DEPRECATED double
227  cell_measure(
228  const std::vector<Point<dim>> &all_vertices,
230 
247  template <int dim>
248  double
249  cell_measure(const std::vector<Point<dim>> & all_vertices,
251 
274  template <int dim, int spacedim>
275  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
277 
304  template <int dim>
305  Vector<double>
307  const Triangulation<dim> &triangulation,
308  const Quadrature<dim> & quadrature);
309 
317  template <int dim>
318  double
320  const Triangulation<dim> &triangulation,
321  const Quadrature<dim> & quadrature);
322 
336  template <int dim, int spacedim>
339 
357  template <typename Iterator>
360  const Iterator & object,
362 
375  template <int dim, int spacedim>
376  std::
377  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
379 
385 
402  template <int dim, int spacedim>
403  void
405  std::vector<CellData<dim>> & cells,
406  SubCellData & subcelldata);
407 
425  template <int dim, int spacedim>
426  void
427  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
428  std::vector<CellData<dim>> & cells,
429  SubCellData & subcelldata,
430  std::vector<unsigned int> & considered_vertices,
431  const double tol = 1e-12);
432 
451  template <int dim, int spacedim>
452  void
454  const std::vector<Point<spacedim>> &all_vertices,
455  std::vector<CellData<dim>> & cells);
456 
466  template <int dim>
467  void
468  consistently_order_cells(std::vector<CellData<dim>> &cells);
469 
475 
529  template <int dim, typename Transformation, int spacedim>
530  void
531  transform(const Transformation & transformation,
532  Triangulation<dim, spacedim> &triangulation);
533 
539  template <int dim, int spacedim>
540  void
541  shift(const Tensor<1, spacedim> & shift_vector,
542  Triangulation<dim, spacedim> &triangulation);
543 
544 
554  template <int dim>
555  void
556  rotate(const double angle, Triangulation<dim> &triangulation);
557 
570  template <int dim>
571  void
572  rotate(const double angle,
573  const unsigned int axis,
574  Triangulation<dim, 3> &triangulation);
575 
633  template <int dim>
634  void
635  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
636  Triangulation<dim> & tria,
637  const Function<dim, double> *coefficient = nullptr,
638  const bool solve_for_absolute_positions = false);
639 
645  template <int dim, int spacedim>
646  std::map<unsigned int, Point<spacedim>>
648 
656  template <int dim, int spacedim>
657  void
658  scale(const double scaling_factor,
659  Triangulation<dim, spacedim> &triangulation);
660 
675  template <int dim, int spacedim>
676  void
678  const double factor,
679  Triangulation<dim, spacedim> &triangulation,
680  const bool keep_boundary = true,
681  const unsigned int seed = boost::random::mt19937::default_seed);
682 
716  template <int dim, int spacedim>
717  void
719  const bool isotropic = false,
720  const unsigned int max_iterations = 100);
721 
746  template <int dim, int spacedim>
747  void
749  const double max_ratio = 1.6180339887,
750  const unsigned int max_iterations = 5);
751 
841  template <int dim, int spacedim>
842  void
844  const double limit_angle_fraction = .75);
845 
851 
906  template <int dim, int spacedim>
907 # ifndef DOXYGEN
908  std::tuple<
909  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
910  std::vector<std::vector<Point<dim>>>,
911  std::vector<std::vector<unsigned int>>>
912 # else
913  return_type
914 # endif
916  const Cache<dim, spacedim> & cache,
917  const std::vector<Point<spacedim>> &points,
919  &cell_hint =
921 
955  template <int dim, int spacedim>
956 # ifndef DOXYGEN
957  std::tuple<
958  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
959  std::vector<std::vector<Point<dim>>>,
960  std::vector<std::vector<unsigned int>>,
961  std::vector<unsigned int>>
962 # else
963  return_type
964 # endif
966  const Cache<dim, spacedim> & cache,
967  const std::vector<Point<spacedim>> &points,
969  &cell_hint =
971 
1041  template <int dim, int spacedim>
1042 # ifndef DOXYGEN
1043  std::tuple<
1044  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1045  std::vector<std::vector<Point<dim>>>,
1046  std::vector<std::vector<unsigned int>>,
1047  std::vector<std::vector<Point<spacedim>>>,
1048  std::vector<std::vector<unsigned int>>>
1049 # else
1050  return_type
1051 # endif
1053  const GridTools::Cache<dim, spacedim> & cache,
1054  const std::vector<Point<spacedim>> & local_points,
1055  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1056  const double tolerance = 1e-10);
1057 
1058  namespace internal
1059  {
1074  template <int dim, int spacedim>
1076  {
1083  std::vector<std::tuple<std::pair<int, int>,
1084  unsigned int,
1085  unsigned int,
1086  Point<dim>,
1088  unsigned int>>
1090 
1094  std::vector<unsigned int> send_ranks;
1095 
1101  std::vector<unsigned int> send_ptrs;
1102 
1113  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1115 
1119  std::vector<unsigned int> recv_ranks;
1120 
1126  std::vector<unsigned int> recv_ptrs;
1127  };
1128 
1138  template <int dim, int spacedim>
1141  const GridTools::Cache<dim, spacedim> & cache,
1142  const std::vector<Point<spacedim>> & points,
1143  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1144  const double tolerance,
1145  const bool perform_handshake,
1146  const bool enforce_unique_mapping = false);
1147 
1148  } // namespace internal
1149 
1186  template <int dim, int spacedim>
1187  std::map<unsigned int, Point<spacedim>>
1189  const Triangulation<dim, spacedim> &container,
1190  const Mapping<dim, spacedim> & mapping =
1191  (ReferenceCells::get_hypercube<dim>()
1192  .template get_default_linear_mapping<dim, spacedim>()));
1193 
1203  template <int spacedim>
1204  unsigned int
1205  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1206  const Point<spacedim> & p);
1207 
1231  template <int dim, template <int, int> class MeshType, int spacedim>
1232  unsigned int
1233  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1234  const Point<spacedim> & p,
1235  const std::vector<bool> & marked_vertices = {});
1236 
1260  template <int dim, template <int, int> class MeshType, int spacedim>
1261  unsigned int
1263  const MeshType<dim, spacedim> &mesh,
1264  const Point<spacedim> & p,
1265  const std::vector<bool> & marked_vertices = {});
1266 
1267 
1287  template <int dim, template <int, int> class MeshType, int spacedim>
1288 # ifndef _MSC_VER
1289  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1290 # else
1291  std::vector<
1292  typename ::internal::
1293  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1294 # endif
1295  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1296  const unsigned int vertex_index);
1297 
1360  template <int dim, template <int, int> class MeshType, int spacedim>
1361 # ifndef _MSC_VER
1362  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1363 # else
1364  std::pair<typename ::internal::
1365  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1366  Point<dim>>
1367 # endif
1369  const MeshType<dim, spacedim> &mesh,
1370  const Point<spacedim> & p,
1371  const std::vector<bool> &marked_vertices = {},
1372  const double tolerance = 1.e-10);
1373 
1381  template <int dim, template <int, int> class MeshType, int spacedim>
1382 # ifndef _MSC_VER
1383  typename MeshType<dim, spacedim>::active_cell_iterator
1384 # else
1385  typename ::internal::
1386  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1387 # endif
1388  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1389  const Point<spacedim> & p,
1390  const std::vector<bool> &marked_vertices = {},
1391  const double tolerance = 1.e-10);
1392 
1399  template <int dim, int spacedim>
1400  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1401  Point<dim>>
1403  const hp::MappingCollection<dim, spacedim> &mapping,
1404  const DoFHandler<dim, spacedim> & mesh,
1405  const Point<spacedim> & p,
1406  const double tolerance = 1.e-10);
1407 
1459  template <int dim, int spacedim>
1460  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1461  Point<dim>>
1463  const Cache<dim, spacedim> &cache,
1464  const Point<spacedim> & p,
1467  const std::vector<bool> &marked_vertices = {},
1468  const double tolerance = 1.e-10);
1469 
1483  template <int dim, template <int, int> class MeshType, int spacedim>
1484 # ifndef _MSC_VER
1485  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1486 # else
1487  std::pair<typename ::internal::
1488  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1489  Point<dim>>
1490 # endif
1492  const Mapping<dim, spacedim> & mapping,
1493  const MeshType<dim, spacedim> &mesh,
1494  const Point<spacedim> & p,
1495  const std::vector<
1496  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1498  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1499  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1500  typename MeshType<dim, spacedim>::active_cell_iterator(),
1501  const std::vector<bool> & marked_vertices = {},
1502  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1503  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1504  const double tolerance = 1.e-10,
1505  const RTree<
1506  std::pair<BoundingBox<spacedim>,
1508  *relevant_cell_bounding_boxes_rtree = nullptr);
1509 
1531  template <int dim, template <int, int> class MeshType, int spacedim>
1532 # ifndef _MSC_VER
1533  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1534  Point<dim>>>
1535 # else
1536  std::vector<std::pair<
1537  typename ::internal::
1538  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1539  Point<dim>>>
1540 # endif
1542  const Mapping<dim, spacedim> & mapping,
1543  const MeshType<dim, spacedim> &mesh,
1544  const Point<spacedim> & p,
1545  const double tolerance,
1546  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1547  Point<dim>> & first_cell);
1548 
1555  template <int dim, template <int, int> class MeshType, int spacedim>
1556 # ifndef _MSC_VER
1557  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1558  Point<dim>>>
1559 # else
1560  std::vector<std::pair<
1561  typename ::internal::
1562  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1563  Point<dim>>>
1564 # endif
1566  const Mapping<dim, spacedim> & mapping,
1567  const MeshType<dim, spacedim> &mesh,
1568  const Point<spacedim> & p,
1569  const double tolerance = 1e-10,
1570  const std::vector<bool> & marked_vertices = {});
1571 
1593  template <class MeshType>
1594  std::vector<typename MeshType::active_cell_iterator>
1595  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1596 
1621  template <class MeshType>
1622  void
1624  const typename MeshType::active_cell_iterator & cell,
1625  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1626 
1676  template <class MeshType>
1677  std::vector<typename MeshType::active_cell_iterator>
1679  const MeshType &mesh,
1680  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1681  &predicate);
1682 
1683 
1691  template <class MeshType>
1692  std::vector<typename MeshType::cell_iterator>
1694  const MeshType &mesh,
1695  const std::function<bool(const typename MeshType::cell_iterator &)>
1696  & predicate,
1697  const unsigned int level);
1698 
1699 
1712  template <class MeshType>
1713  std::vector<typename MeshType::active_cell_iterator>
1714  compute_ghost_cell_halo_layer(const MeshType &mesh);
1715 
1764  template <class MeshType>
1765  std::vector<typename MeshType::active_cell_iterator>
1767  const MeshType &mesh,
1768  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1769  & predicate,
1770  const double layer_thickness);
1771 
1794  template <class MeshType>
1795  std::vector<typename MeshType::active_cell_iterator>
1796  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1797  const double layer_thickness);
1798 
1814  template <class MeshType>
1815  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1817  const MeshType &mesh,
1818  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1819  &predicate);
1820 
1873  template <class MeshType>
1874  std::vector<BoundingBox<MeshType::space_dimension>>
1876  const MeshType &mesh,
1877  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1878  & predicate,
1879  const unsigned int refinement_level = 0,
1880  const bool allow_merge = false,
1881  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1882 
1910  template <int spacedim>
1911 # ifndef DOXYGEN
1912  std::tuple<std::vector<std::vector<unsigned int>>,
1913  std::map<unsigned int, unsigned int>,
1914  std::map<unsigned int, std::vector<unsigned int>>>
1915 # else
1916  return_type
1917 # endif
1919  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1920  const std::vector<Point<spacedim>> & points);
1921 
1922 
1957  template <int spacedim>
1958 # ifndef DOXYGEN
1959  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1960  std::map<unsigned int, unsigned int>,
1961  std::map<unsigned int, std::vector<unsigned int>>>
1962 # else
1963  return_type
1964 # endif
1966  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1967  const std::vector<Point<spacedim>> & points);
1968 
1969 
1978  template <int dim, int spacedim>
1979  std::vector<
1980  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1981  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1982 
1995  template <int dim, int spacedim>
1996  std::vector<std::vector<Tensor<1, spacedim>>>
1998  const Triangulation<dim, spacedim> &mesh,
1999  const std::vector<
2001  &vertex_to_cells);
2002 
2003 
2011  template <int dim, int spacedim>
2012  unsigned int
2015  const Point<spacedim> & position,
2016  const Mapping<dim, spacedim> & mapping =
2017  (ReferenceCells::get_hypercube<dim>()
2018  .template get_default_linear_mapping<dim, spacedim>()));
2019 
2031  template <int dim, int spacedim>
2032  std::map<unsigned int, types::global_vertex_index>
2035 
2047  template <int dim, int spacedim>
2048  std::pair<unsigned int, double>
2051 
2057 
2066  template <int dim, int spacedim>
2067  void
2069  const Triangulation<dim, spacedim> &triangulation,
2070  DynamicSparsityPattern & connectivity);
2071 
2080  template <int dim, int spacedim>
2081  void
2083  const Triangulation<dim, spacedim> &triangulation,
2084  DynamicSparsityPattern & connectivity);
2085 
2094  template <int dim, int spacedim>
2095  void
2097  const Triangulation<dim, spacedim> &triangulation,
2098  const unsigned int level,
2099  DynamicSparsityPattern & connectivity);
2100 
2121  template <int dim, int spacedim>
2122  void
2123  partition_triangulation(const unsigned int n_partitions,
2124  Triangulation<dim, spacedim> & triangulation,
2125  const SparsityTools::Partitioner partitioner =
2127 
2138  template <int dim, int spacedim>
2139  void
2140  partition_triangulation(const unsigned int n_partitions,
2141  const std::vector<unsigned int> &cell_weights,
2142  Triangulation<dim, spacedim> & triangulation,
2143  const SparsityTools::Partitioner partitioner =
2145 
2191  template <int dim, int spacedim>
2192  void
2193  partition_triangulation(const unsigned int n_partitions,
2194  const SparsityPattern & cell_connection_graph,
2195  Triangulation<dim, spacedim> &triangulation,
2196  const SparsityTools::Partitioner partitioner =
2198 
2209  template <int dim, int spacedim>
2210  void
2211  partition_triangulation(const unsigned int n_partitions,
2212  const std::vector<unsigned int> &cell_weights,
2213  const SparsityPattern & cell_connection_graph,
2214  Triangulation<dim, spacedim> &triangulation,
2215  const SparsityTools::Partitioner partitioner =
2217 
2232  template <int dim, int spacedim>
2233  void
2234  partition_triangulation_zorder(const unsigned int n_partitions,
2235  Triangulation<dim, spacedim> &triangulation,
2236  const bool group_siblings = true);
2237 
2249  template <int dim, int spacedim>
2250  void
2252 
2260  template <int dim, int spacedim>
2261  std::vector<types::subdomain_id>
2263  const std::vector<CellId> & cell_ids);
2264 
2275  template <int dim, int spacedim>
2276  void
2278  std::vector<types::subdomain_id> & subdomain);
2279 
2294  template <int dim, int spacedim>
2295  unsigned int
2297  const Triangulation<dim, spacedim> &triangulation,
2298  const types::subdomain_id subdomain);
2299 
2329  template <int dim, int spacedim>
2330  std::vector<bool>
2332 
2338 
2372  template <typename MeshType>
2373  std::list<std::pair<typename MeshType::cell_iterator,
2374  typename MeshType::cell_iterator>>
2375  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2376 
2386  template <int dim, int spacedim>
2387  bool
2389  const Triangulation<dim, spacedim> &mesh_2);
2390 
2400  template <typename MeshType>
2401  bool
2402  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2403 
2409 
2425  template <int dim, int spacedim>
2429  & distorted_cells,
2430  Triangulation<dim, spacedim> &triangulation);
2431 
2432 
2433 
2442 
2443 
2481  template <class MeshType>
2482  std::vector<typename MeshType::active_cell_iterator>
2483  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2484 
2485 
2507  template <class Container>
2508  std::vector<typename Container::cell_iterator>
2510  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2511 
2578  template <class Container>
2579  void
2581  const std::vector<typename Container::active_cell_iterator> &patch,
2583  &local_triangulation,
2584  std::map<
2585  typename Triangulation<Container::dimension,
2586  Container::space_dimension>::active_cell_iterator,
2587  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2588 
2620  template <int dim, int spacedim>
2621  std::map<
2623  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2625 
2626 
2633 
2639  template <typename CellIterator>
2641  {
2645  CellIterator cell[2];
2646 
2651  unsigned int face_idx[2];
2652 
2658  std::bitset<3> orientation;
2659 
2673  };
2674 
2675 
2739  template <typename FaceIterator>
2740  bool
2742  std::bitset<3> & orientation,
2743  const FaceIterator & face1,
2744  const FaceIterator & face2,
2745  const int direction,
2749 
2750 
2754  template <typename FaceIterator>
2755  bool
2757  const FaceIterator & face1,
2758  const FaceIterator & face2,
2759  const int direction,
2763 
2764 
2821  template <typename MeshType>
2822  void
2824  const MeshType & mesh,
2825  const types::boundary_id b_id1,
2826  const types::boundary_id b_id2,
2827  const int direction,
2829  & matched_pairs,
2830  const Tensor<1, MeshType::space_dimension> &offset =
2833 
2834 
2857  template <typename MeshType>
2858  void
2860  const MeshType & mesh,
2861  const types::boundary_id b_id,
2862  const int direction,
2864  & matched_pairs,
2865  const ::Tensor<1, MeshType::space_dimension> &offset =
2868 
2874 
2895  template <int dim, int spacedim>
2896  void
2898  const bool reset_boundary_ids = false);
2899 
2921  template <int dim, int spacedim>
2922  void
2924  const std::vector<types::boundary_id> &src_boundary_ids,
2925  const std::vector<types::manifold_id> &dst_manifold_ids,
2927  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2928 
2958  template <int dim, int spacedim>
2959  void
2961  const bool compute_face_ids = false);
2962 
2987  template <int dim, int spacedim>
2988  void
2991  const std::function<types::manifold_id(
2992  const std::set<types::manifold_id> &)> &disambiguation_function =
2993  [](const std::set<types::manifold_id> &manifold_ids) {
2994  if (manifold_ids.size() == 1)
2995  return *manifold_ids.begin();
2996  else
2998  },
2999  bool overwrite_only_flat_manifold_ids = true);
3086  template <typename DataType, typename MeshType>
3087  void
3089  const MeshType & mesh,
3090  const std::function<std_cxx17::optional<DataType>(
3091  const typename MeshType::active_cell_iterator &)> &pack,
3092  const std::function<void(const typename MeshType::active_cell_iterator &,
3093  const DataType &)> & unpack,
3094  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3095  &cell_filter =
3097 
3108  template <typename DataType, typename MeshType>
3109  void
3111  const MeshType & mesh,
3112  const std::function<std_cxx17::optional<DataType>(
3113  const typename MeshType::level_cell_iterator &)> &pack,
3114  const std::function<void(const typename MeshType::level_cell_iterator &,
3115  const DataType &)> & unpack,
3116  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3118  true});
3119 
3120  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3121  * boxes @p local_bboxes.
3122  *
3123  * This function is meant to exchange bounding boxes describing the locally
3124  * owned cells in a distributed triangulation obtained with the function
3125  * GridTools::compute_mesh_predicate_bounding_box .
3126  *
3127  * The output vector's size is the number of processes of the MPI
3128  * communicator:
3129  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3130  */
3131  template <int spacedim>
3132  std::vector<std::vector<BoundingBox<spacedim>>>
3134  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3135  const MPI_Comm & mpi_communicator);
3136 
3169  template <int spacedim>
3172  const std::vector<BoundingBox<spacedim>> &local_description,
3173  const MPI_Comm & mpi_communicator);
3174 
3192  template <int dim, int spacedim>
3193  void
3195  const Triangulation<dim, spacedim> & tria,
3196  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3197  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3198 
3205  template <int dim, int spacedim>
3206  std::map<unsigned int, std::set<::types::subdomain_id>>
3208  const Triangulation<dim, spacedim> &tria);
3209 
3222  template <int dim, typename T>
3224  {
3228  std::vector<CellId> cell_ids;
3229 
3233  std::vector<T> data;
3234 
3243  template <class Archive>
3244  void
3245  save(Archive &ar, const unsigned int version) const;
3246 
3253  template <class Archive>
3254  void
3255  load(Archive &ar, const unsigned int version);
3256 
3257 # ifdef DOXYGEN
3258 
3263  template <class Archive>
3264  void
3265  serialize(Archive &archive, const unsigned int version);
3266 # else
3267  // This macro defines the serialize() method that is compatible with
3268  // the templated save() and load() method that have been implemented.
3269  BOOST_SERIALIZATION_SPLIT_MEMBER()
3270 # endif
3271  };
3272 
3273 
3274 
3292  template <int dim, typename VectorType>
3294  {
3295  public:
3299  using value_type = typename VectorType::value_type;
3300 
3304  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3305  const FiniteElement<dim, dim> &fe,
3306  const unsigned int n_subdivisions = 1,
3307  const double tolerance = 1e-10);
3308 
3313  void
3314  process(const DoFHandler<dim> & background_dof_handler,
3315  const VectorType & ls_vector,
3316  const double iso_level,
3317  std::vector<Point<dim>> & vertices,
3318  std::vector<CellData<dim - 1>> &cells) const;
3319 
3326  void
3327  process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
3328  const VectorType & ls_vector,
3329  const double iso_level,
3330  std::vector<Point<dim>> & vertices,
3331  std::vector<CellData<dim - 1>> &cells) const;
3332 
3333  private:
3338  static Quadrature<dim>
3339  create_quadrature_rule(const unsigned int n_subdivisions);
3340 
3344  void
3345  process_cell(std::vector<value_type> & ls_values,
3346  const std::vector<Point<dim>> & points,
3347  const double iso_level,
3348  std::vector<Point<dim>> & vertices,
3349  std::vector<CellData<dim - 1>> &cells) const;
3350 
3357  void
3358  process_sub_cell(const std::vector<value_type> & ls_values,
3359  const std::vector<Point<2>> & points,
3360  const std::vector<unsigned int> mask,
3361  const double iso_level,
3362  std::vector<Point<2>> & vertices,
3363  std::vector<CellData<1>> & cells) const;
3364 
3368  void
3369  process_sub_cell(const std::vector<value_type> & ls_values,
3370  const std::vector<Point<3>> & points,
3371  const std::vector<unsigned int> mask,
3372  const double iso_level,
3373  std::vector<Point<3>> & vertices,
3374  std::vector<CellData<2>> & cells) const;
3375 
3380  const unsigned int n_subdivisions;
3381 
3386  const double tolerance;
3387 
3393  };
3394 
3395 
3396 
3401 
3406  int,
3407  << "The number of partitions you gave is " << arg1
3408  << ", but must be greater than zero.");
3413  int,
3414  << "The subdomain id " << arg1
3415  << " has no cells associated with it.");
3420 
3425  double,
3426  << "The scaling factor must be positive, but it is " << arg1
3427  << ".");
3428 
3433  unsigned int,
3434  << "The given vertex with index " << arg1
3435  << " is not used in the given triangulation.");
3436 
3439 } /*namespace GridTools*/
3440 
3441 
3452  "The edges of the mesh are not consistently orientable.");
3453 
3454 
3455 
3456 /* ----------------- Template function --------------- */
3457 
3458 # ifndef DOXYGEN
3459 
3460 namespace GridTools
3461 {
3462  template <int dim>
3463  double
3464  cell_measure(
3465  const std::vector<Point<dim>> &all_vertices,
3466  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3467  {
3468  // We forward call to the ArrayView version:
3469  const ArrayView<const unsigned int> view(
3470  indices, GeometryInfo<dim>::vertices_per_cell);
3471  return cell_measure(all_vertices, view);
3472  }
3473 
3474 
3475 
3476  // This specialization is defined here so that the general template in the
3477  // source file doesn't need to have further 1D overloads for the internal
3478  // functions it calls.
3479  template <>
3483  {
3484  return {};
3485  }
3486 
3487 
3488 
3489  template <int dim, typename Predicate, int spacedim>
3490  void
3491  transform(const Predicate & predicate,
3492  Triangulation<dim, spacedim> &triangulation)
3493  {
3494  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3495 
3496  // loop over all active cells, and
3497  // transform those vertices that
3498  // have not yet been touched. note
3499  // that we get to all vertices in
3500  // the triangulation by only
3501  // visiting the active cells.
3503  cell = triangulation.begin_active(),
3504  endc = triangulation.end();
3505  for (; cell != endc; ++cell)
3506  for (const unsigned int v : cell->vertex_indices())
3507  if (treated_vertices[cell->vertex_index(v)] == false)
3508  {
3509  // transform this vertex
3510  cell->vertex(v) = predicate(cell->vertex(v));
3511  // and mark it as treated
3512  treated_vertices[cell->vertex_index(v)] = true;
3513  };
3514 
3515 
3516  // now fix any vertices on hanging nodes so that we don't create any holes
3517  if (dim == 2)
3518  {
3520  cell = triangulation.begin_active(),
3521  endc = triangulation.end();
3522  for (; cell != endc; ++cell)
3523  for (const unsigned int face : cell->face_indices())
3524  if (cell->face(face)->has_children() &&
3525  !cell->face(face)->at_boundary())
3526  {
3527  Assert(cell->reference_cell() ==
3528  ReferenceCells::get_hypercube<dim>(),
3529  ExcNotImplemented());
3530 
3531  // this line has children
3532  cell->face(face)->child(0)->vertex(1) =
3533  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3534  2;
3535  }
3536  }
3537  else if (dim == 3)
3538  {
3540  cell = triangulation.begin_active(),
3541  endc = triangulation.end();
3542  for (; cell != endc; ++cell)
3543  for (const unsigned int face : cell->face_indices())
3544  if (cell->face(face)->has_children() &&
3545  !cell->face(face)->at_boundary())
3546  {
3547  Assert(cell->reference_cell() ==
3548  ReferenceCells::get_hypercube<dim>(),
3549  ExcNotImplemented());
3550 
3551  // this face has hanging nodes
3552  cell->face(face)->child(0)->vertex(1) =
3553  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3554  2.0;
3555  cell->face(face)->child(0)->vertex(2) =
3556  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3557  2.0;
3558  cell->face(face)->child(1)->vertex(3) =
3559  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3560  2.0;
3561  cell->face(face)->child(2)->vertex(3) =
3562  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3563  2.0;
3564 
3565  // center of the face
3566  cell->face(face)->child(0)->vertex(3) =
3567  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3568  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3569  4.0;
3570  }
3571  }
3572 
3573  // Make sure FEValues notices that the mesh has changed
3574  triangulation.signals.mesh_movement();
3575  }
3576 
3577 
3578 
3579  template <class MeshType>
3580  std::vector<typename MeshType::active_cell_iterator>
3581  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3582  {
3583  std::vector<typename MeshType::active_cell_iterator> child_cells;
3584 
3585  if (cell->has_children())
3586  {
3587  for (unsigned int child = 0; child < cell->n_children(); ++child)
3588  if (cell->child(child)->has_children())
3589  {
3590  const std::vector<typename MeshType::active_cell_iterator>
3591  children = get_active_child_cells<MeshType>(cell->child(child));
3592  child_cells.insert(child_cells.end(),
3593  children.begin(),
3594  children.end());
3595  }
3596  else
3597  child_cells.push_back(cell->child(child));
3598  }
3599 
3600  return child_cells;
3601  }
3602 
3603 
3604 
3605  template <class MeshType>
3606  void
3608  const typename MeshType::active_cell_iterator & cell,
3609  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3610  {
3611  active_neighbors.clear();
3612  for (const unsigned int n : cell->face_indices())
3613  if (!cell->at_boundary(n))
3614  {
3615  if (MeshType::dimension == 1)
3616  {
3617  // check children of neighbor. note
3618  // that in 1d children of the neighbor
3619  // may be further refined. In 1d the
3620  // case is simple since we know what
3621  // children bound to the present cell
3622  typename MeshType::cell_iterator neighbor_child =
3623  cell->neighbor(n);
3624  if (!neighbor_child->is_active())
3625  {
3626  while (neighbor_child->has_children())
3627  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3628 
3629  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3630  ExcInternalError());
3631  }
3632  active_neighbors.push_back(neighbor_child);
3633  }
3634  else
3635  {
3636  if (cell->face(n)->has_children())
3637  // this neighbor has children. find
3638  // out which border to the present
3639  // cell
3640  for (unsigned int c = 0;
3641  c < cell->face(n)->n_active_descendants();
3642  ++c)
3643  active_neighbors.push_back(
3644  cell->neighbor_child_on_subface(n, c));
3645  else
3646  {
3647  // the neighbor must be active
3648  // himself
3649  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3650  active_neighbors.push_back(cell->neighbor(n));
3651  }
3652  }
3653  }
3654  }
3655 
3656 
3657 
3658  namespace internal
3659  {
3660  namespace ProjectToObject
3661  {
3674  struct CrossDerivative
3675  {
3676  const unsigned int direction_0;
3677  const unsigned int direction_1;
3678 
3679  CrossDerivative(const unsigned int d0, const unsigned int d1);
3680  };
3681 
3682  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3683  const unsigned int d1)
3684  : direction_0(d0)
3685  , direction_1(d1)
3686  {}
3687 
3688 
3689 
3694  template <typename F>
3695  inline auto
3696  centered_first_difference(const double center,
3697  const double step,
3698  const F &f) -> decltype(f(center) - f(center))
3699  {
3700  return (f(center + step) - f(center - step)) / (2.0 * step);
3701  }
3702 
3703 
3704 
3709  template <typename F>
3710  inline auto
3711  centered_second_difference(const double center,
3712  const double step,
3713  const F &f) -> decltype(f(center) - f(center))
3714  {
3715  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3716  (step * step);
3717  }
3718 
3719 
3720 
3730  template <int structdim, typename F>
3731  inline auto
3732  cross_stencil(
3733  const CrossDerivative cross_derivative,
3735  const double step,
3736  const F &f) -> decltype(f(center) - f(center))
3737  {
3739  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3740  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3741  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3742  1.0 / 3.0 * f(center - simplex_vector) +
3743  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3744  step;
3745  }
3746 
3747 
3748 
3755  template <int spacedim, int structdim, typename F>
3756  inline double
3757  gradient_entry(
3758  const unsigned int row_n,
3759  const unsigned int dependent_direction,
3760  const Point<spacedim> &p0,
3762  const double step,
3763  const F & f)
3764  {
3766  dependent_direction <
3768  ExcMessage("This function assumes that the last weight is a "
3769  "dependent variable (and hence we cannot take its "
3770  "derivative directly)."));
3771  Assert(row_n != dependent_direction,
3772  ExcMessage(
3773  "We cannot differentiate with respect to the variable "
3774  "that is assumed to be dependent."));
3775 
3776  const Point<spacedim> manifold_point = f(center);
3777  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3778  {row_n, dependent_direction}, center, step, f);
3779  double entry = 0.0;
3780  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3781  entry +=
3782  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3783  return entry;
3784  }
3785 
3791  template <typename Iterator, int spacedim, int structdim>
3793  project_to_d_linear_object(const Iterator & object,
3794  const Point<spacedim> &trial_point)
3795  {
3796  // let's look at this for simplicity for a quad (structdim==2) in a
3797  // space with spacedim>2 (notate trial_point by y): all points on the
3798  // surface are given by
3799  // x(\xi) = sum_i v_i phi_x(\xi)
3800  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3801  // reference coordinates of the quad. so what we are trying to do is
3802  // find a point x on the surface that is closest to the point y. there
3803  // are different ways to solve this problem, but in the end it's a
3804  // nonlinear problem and we have to find reference coordinates \xi so
3805  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3806  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3807  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3808  // have to use a Newton method to find the answer. This leads to the
3809  // following formulation of Newton steps:
3810  //
3811  // Given \xi_k, find \delta\xi_k so that
3812  // H_k \delta\xi_k = - F_k
3813  // where H_k is an approximation to the second derivatives of J at
3814  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3815  // number of times until the right hand side is small enough. As a
3816  // stopping criterion, we terminate if ||\delta\xi||<eps.
3817  //
3818  // As for the Hessian, the best choice would be
3819  // H_k = J''(\xi_k)
3820  // but we'll opt for the simpler Gauss-Newton form
3821  // H_k = A^T A
3822  // i.e.
3823  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3824  // \partial_n phi_i *
3825  // \partial_m phi_j
3826  // we start at xi=(0.5, 0.5).
3827  Point<structdim> xi;
3828  for (unsigned int d = 0; d < structdim; ++d)
3829  xi[d] = 0.5;
3830 
3831  Point<spacedim> x_k;
3832  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3833  x_k += object->vertex(i) *
3835 
3836  do
3837  {
3839  for (const unsigned int i :
3841  F_k +=
3842  (x_k - trial_point) * object->vertex(i) *
3844  i);
3845 
3847  for (const unsigned int i :
3849  for (const unsigned int j :
3851  {
3854  xi, i),
3856  xi, j));
3857  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3858  }
3859 
3860  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3861  xi += delta_xi;
3862 
3863  x_k = Point<spacedim>();
3864  for (const unsigned int i :
3866  x_k += object->vertex(i) *
3868 
3869  if (delta_xi.norm() < 1e-7)
3870  break;
3871  }
3872  while (true);
3873 
3874  return x_k;
3875  }
3876  } // namespace ProjectToObject
3877  } // namespace internal
3878 
3879 
3880 
3881  namespace internal
3882  {
3883  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3884  // inside the project_to_object function below.
3885  template <int structdim>
3886  inline bool
3887  weights_are_ok(
3889  {
3890  // clang has trouble figuring out structdim here, so define it
3891  // again:
3892  static const std::size_t n_vertices_per_cell =
3894  n_independent_components;
3895  std::array<double, n_vertices_per_cell> copied_weights;
3896  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3897  {
3898  copied_weights[i] = v[i];
3899  if (v[i] < 0.0 || v[i] > 1.0)
3900  return false;
3901  }
3902 
3903  // check the sum: try to avoid some roundoff errors by summing in order
3904  std::sort(copied_weights.begin(), copied_weights.end());
3905  const double sum =
3906  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3907  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3908  }
3909  } // namespace internal
3910 
3911  template <typename Iterator>
3914  const Iterator & object,
3916  {
3917  const int spacedim = Iterator::AccessorType::space_dimension;
3918  const int structdim = Iterator::AccessorType::structure_dimension;
3919 
3920  Point<spacedim> projected_point = trial_point;
3921 
3922  if (structdim >= spacedim)
3923  return projected_point;
3924  else if (structdim == 1 || structdim == 2)
3925  {
3926  using namespace internal::ProjectToObject;
3927  // Try to use the special flat algorithm for quads (this is better
3928  // than the general algorithm in 3D). This does not take into account
3929  // whether projected_point is outside the quad, but we optimize along
3930  // lines below anyway:
3931  const int dim = Iterator::AccessorType::dimension;
3932  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3933  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3934  &manifold) != nullptr)
3935  {
3936  projected_point =
3937  project_to_d_linear_object<Iterator, spacedim, structdim>(
3938  object, trial_point);
3939  }
3940  else
3941  {
3942  // We want to find a point on the convex hull (defined by the
3943  // vertices of the object and the manifold description) that is
3944  // relatively close to the trial point. This has a few issues:
3945  //
3946  // 1. For a general convex hull we are not guaranteed that a unique
3947  // minimum exists.
3948  // 2. The independent variables in the optimization process are the
3949  // weights given to Manifold::get_new_point, which must sum to 1,
3950  // so we cannot use standard finite differences to approximate a
3951  // gradient.
3952  //
3953  // There is not much we can do about 1., but for 2. we can derive
3954  // finite difference stencils that work on a structdim-dimensional
3955  // simplex and rewrite the optimization problem to use those
3956  // instead. Consider the structdim 2 case and let
3957  //
3958  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3959  // c2, c3})
3960  //
3961  // where {c0, c1, c2, c3} are the weights for the four vertices on
3962  // the quadrilateral. We seek to minimize the Euclidean distance
3963  // between F(...) and trial_point. We can solve for c3 in terms of
3964  // the other weights and get, for one coordinate direction
3965  //
3966  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3967  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3968  //
3969  // where we substitute back in for c3 after taking the
3970  // derivative. We can compute a stencil for the cross derivative
3971  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3972  // (and gradient_entry computes the sum over the independent
3973  // variables). Below, we somewhat arbitrarily pick the last
3974  // component as the dependent one.
3975  //
3976  // Since we can now calculate derivatives of the objective
3977  // function we can use gradient descent to minimize it.
3978  //
3979  // Of course, this is much simpler in the structdim = 1 case (we
3980  // could rewrite the projection as a 1D optimization problem), but
3981  // to reduce the potential for bugs we use the same code in both
3982  // cases.
3983  const double step_size = object->diameter() / 64.0;
3984 
3985  constexpr unsigned int n_vertices_per_cell =
3987 
3988  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3989  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3990  ++vertex_n)
3991  vertices[vertex_n] = object->vertex(vertex_n);
3992 
3993  auto get_point_from_weights =
3994  [&](const Tensor<1, n_vertices_per_cell> &weights)
3995  -> Point<spacedim> {
3996  return object->get_manifold().get_new_point(
3997  make_array_view(vertices.begin(), vertices.end()),
3998  make_array_view(weights.begin_raw(), weights.end_raw()));
3999  };
4000 
4001  // pick the initial weights as (normalized) inverse distances from
4002  // the trial point:
4003  Tensor<1, n_vertices_per_cell> guess_weights;
4004  double guess_weights_sum = 0.0;
4005  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4006  ++vertex_n)
4007  {
4008  const double distance =
4009  vertices[vertex_n].distance(trial_point);
4010  if (distance == 0.0)
4011  {
4012  guess_weights = 0.0;
4013  guess_weights[vertex_n] = 1.0;
4014  guess_weights_sum = 1.0;
4015  break;
4016  }
4017  else
4018  {
4019  guess_weights[vertex_n] = 1.0 / distance;
4020  guess_weights_sum += guess_weights[vertex_n];
4021  }
4022  }
4023  guess_weights /= guess_weights_sum;
4024  Assert(internal::weights_are_ok<structdim>(guess_weights),
4025  ExcInternalError());
4026 
4027  // The optimization algorithm consists of two parts:
4028  //
4029  // 1. An outer loop where we apply the gradient descent algorithm.
4030  // 2. An inner loop where we do a line search to find the optimal
4031  // length of the step one should take in the gradient direction.
4032  //
4033  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4034  {
4035  const unsigned int dependent_direction =
4036  n_vertices_per_cell - 1;
4037  Tensor<1, n_vertices_per_cell> current_gradient;
4038  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4039  ++row_n)
4040  {
4041  if (row_n != dependent_direction)
4042  {
4043  current_gradient[row_n] =
4044  gradient_entry<spacedim, structdim>(
4045  row_n,
4046  dependent_direction,
4047  trial_point,
4048  guess_weights,
4049  step_size,
4050  get_point_from_weights);
4051 
4052  current_gradient[dependent_direction] -=
4053  current_gradient[row_n];
4054  }
4055  }
4056 
4057  // We need to travel in the -gradient direction, as noted
4058  // above, but we may not want to take a full step in that
4059  // direction; instead, guess that we will go -0.5*gradient and
4060  // do quasi-Newton iteration to pick the best multiplier. The
4061  // goal is to find a scalar alpha such that
4062  //
4063  // F(x - alpha g)
4064  //
4065  // is minimized, where g is the gradient and F is the
4066  // objective function. To find the optimal value we find roots
4067  // of the derivative of the objective function with respect to
4068  // alpha by Newton iteration, where we approximate the first
4069  // and second derivatives of F(x - alpha g) with centered
4070  // finite differences.
4071  double gradient_weight = -0.5;
4072  auto gradient_weight_objective_function =
4073  [&](const double gradient_weight_guess) -> double {
4074  return (trial_point -
4075  get_point_from_weights(guess_weights +
4076  gradient_weight_guess *
4077  current_gradient))
4078  .norm_square();
4079  };
4080 
4081  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4082  {
4083  const double update_numerator = centered_first_difference(
4084  gradient_weight,
4085  step_size,
4086  gradient_weight_objective_function);
4087  const double update_denominator =
4088  centered_second_difference(
4089  gradient_weight,
4090  step_size,
4091  gradient_weight_objective_function);
4092 
4093  // avoid division by zero. Note that we limit the gradient
4094  // weight below
4095  if (std::abs(update_denominator) == 0.0)
4096  break;
4097  gradient_weight =
4098  gradient_weight - update_numerator / update_denominator;
4099 
4100  // Put a fairly lenient bound on the largest possible
4101  // gradient (things tend to be locally flat, so the gradient
4102  // itself is usually small)
4103  if (std::abs(gradient_weight) > 10)
4104  {
4105  gradient_weight = -10.0;
4106  break;
4107  }
4108  }
4109 
4110  // It only makes sense to take convex combinations with weights
4111  // between zero and one. If the update takes us outside of this
4112  // region then rescale the update to stay within the region and
4113  // try again
4114  Tensor<1, n_vertices_per_cell> tentative_weights =
4115  guess_weights + gradient_weight * current_gradient;
4116 
4117  double new_gradient_weight = gradient_weight;
4118  for (unsigned int iteration_count = 0; iteration_count < 40;
4119  ++iteration_count)
4120  {
4121  if (internal::weights_are_ok<structdim>(tentative_weights))
4122  break;
4123 
4124  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4125  {
4126  if (tentative_weights[i] < 0.0)
4127  {
4128  tentative_weights -=
4129  (tentative_weights[i] / current_gradient[i]) *
4130  current_gradient;
4131  }
4132  if (tentative_weights[i] < 0.0 ||
4133  1.0 < tentative_weights[i])
4134  {
4135  new_gradient_weight /= 2.0;
4136  tentative_weights =
4137  guess_weights +
4138  new_gradient_weight * current_gradient;
4139  }
4140  }
4141  }
4142 
4143  // the update might still send us outside the valid region, so
4144  // check again and quit if the update is still not valid
4145  if (!internal::weights_are_ok<structdim>(tentative_weights))
4146  break;
4147 
4148  // if we cannot get closer by traveling in the gradient
4149  // direction then quit
4150  if (get_point_from_weights(tentative_weights)
4151  .distance(trial_point) <
4152  get_point_from_weights(guess_weights).distance(trial_point))
4153  guess_weights = tentative_weights;
4154  else
4155  break;
4156  Assert(internal::weights_are_ok<structdim>(guess_weights),
4157  ExcInternalError());
4158  }
4159  Assert(internal::weights_are_ok<structdim>(guess_weights),
4160  ExcInternalError());
4161  projected_point = get_point_from_weights(guess_weights);
4162  }
4163 
4164  // if structdim == 2 and the optimal point is not on the interior then
4165  // we may be able to get a more accurate result by projecting onto the
4166  // lines.
4167  if (structdim == 2)
4168  {
4169  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4170  line_projections;
4171  for (unsigned int line_n = 0;
4172  line_n < GeometryInfo<structdim>::lines_per_cell;
4173  ++line_n)
4174  {
4175  line_projections[line_n] =
4176  project_to_object(object->line(line_n), trial_point);
4177  }
4178  std::sort(line_projections.begin(),
4179  line_projections.end(),
4180  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4181  return a.distance(trial_point) <
4182  b.distance(trial_point);
4183  });
4184  if (line_projections[0].distance(trial_point) <
4185  projected_point.distance(trial_point))
4186  projected_point = line_projections[0];
4187  }
4188  }
4189  else
4190  {
4191  Assert(false, ExcNotImplemented());
4192  return projected_point;
4193  }
4194 
4195  return projected_point;
4196  }
4197 
4198 
4199 
4200  template <int dim, typename T>
4201  template <class Archive>
4202  void
4204  const unsigned int /*version*/) const
4205  {
4206  Assert(cell_ids.size() == data.size(),
4207  ExcDimensionMismatch(cell_ids.size(), data.size()));
4208  // archive the cellids in an efficient binary format
4209  const std::size_t n_cells = cell_ids.size();
4210  ar & n_cells;
4211  for (const auto &id : cell_ids)
4212  {
4213  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
4214  ar & binary_cell_id;
4215  }
4216 
4217  ar &data;
4218  }
4219 
4220 
4221 
4222  template <int dim, typename T>
4223  template <class Archive>
4224  void
4226  const unsigned int /*version*/)
4227  {
4228  std::size_t n_cells;
4229  ar & n_cells;
4230  cell_ids.clear();
4231  cell_ids.reserve(n_cells);
4232  for (unsigned int c = 0; c < n_cells; ++c)
4233  {
4234  CellId::binary_type value;
4235  ar & value;
4236  cell_ids.emplace_back(value);
4237  }
4238  ar &data;
4239  }
4240 
4241 
4242  namespace internal
4243  {
4244  template <typename DataType,
4245  typename MeshType,
4246  typename MeshCellIteratorType>
4247  inline void
4248  exchange_cell_data(
4249  const MeshType &mesh,
4250  const std::function<
4251  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4252  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4253  & unpack,
4254  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4255  const std::function<void(
4256  const std::function<void(const MeshCellIteratorType &,
4257  const types::subdomain_id)> &)> &process_cells,
4258  const std::function<std::set<types::subdomain_id>(
4259  const parallel::TriangulationBase<MeshType::dimension,
4260  MeshType::space_dimension> &)>
4261  &compute_ghost_owners)
4262  {
4263 # ifndef DEAL_II_WITH_MPI
4264  (void)mesh;
4265  (void)pack;
4266  (void)unpack;
4267  (void)cell_filter;
4268  (void)process_cells;
4269  (void)compute_ghost_owners;
4270  Assert(false, ExcNeedsMPI());
4271 # else
4272  constexpr int dim = MeshType::dimension;
4273  constexpr int spacedim = MeshType::space_dimension;
4274  auto tria =
4275  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4276  &mesh.get_triangulation());
4277  Assert(
4278  tria != nullptr,
4279  ExcMessage(
4280  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4281 
4282  if (const auto tria = dynamic_cast<
4284  &mesh.get_triangulation()))
4285  {
4286  Assert(
4287  tria->with_artificial_cells(),
4288  ExcMessage(
4289  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4290  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4291  "operate on a single layer ghost cells. However, you have "
4292  "given a Triangulation object of type "
4293  "parallel::shared::Triangulation without artificial cells "
4294  "resulting in arbitrary numbers of ghost layers."));
4295  }
4296 
4297  // build list of cells to request for each neighbor
4298  std::set<::types::subdomain_id> ghost_owners =
4299  compute_ghost_owners(*tria);
4300  std::map<::types::subdomain_id,
4301  std::vector<typename CellId::binary_type>>
4302  neighbor_cell_list;
4303 
4304  for (const auto ghost_owner : ghost_owners)
4305  neighbor_cell_list[ghost_owner] = {};
4306 
4307  process_cells([&](const auto &cell, const auto key) {
4308  if (cell_filter(cell))
4309  neighbor_cell_list[key].emplace_back(
4310  cell->id().template to_binary<spacedim>());
4311  });
4312 
4313  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4314  ExcInternalError());
4315 
4316 
4317  // Before sending & receiving, make sure we protect this section with
4318  // a mutex:
4319  static Utilities::MPI::CollectiveMutex mutex;
4321  mutex, tria->get_communicator());
4322 
4323  const int mpi_tag =
4325  const int mpi_tag_reply =
4327 
4328  // send our requests:
4329  std::vector<MPI_Request> requests(ghost_owners.size());
4330  {
4331  unsigned int idx = 0;
4332  for (const auto &it : neighbor_cell_list)
4333  {
4334  // send the data about the relevant cells
4335  const int ierr = MPI_Isend(it.second.data(),
4336  it.second.size() * sizeof(it.second[0]),
4337  MPI_BYTE,
4338  it.first,
4339  mpi_tag,
4340  tria->get_communicator(),
4341  &requests[idx]);
4342  AssertThrowMPI(ierr);
4343  ++idx;
4344  }
4345  }
4346 
4347  using DestinationToBufferMap =
4348  std::map<::types::subdomain_id,
4350  DestinationToBufferMap destination_to_data_buffer_map;
4351 
4352  // receive requests and reply with the ghost indices
4353  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4354  ghost_owners.size());
4355  std::vector<std::vector<types::global_dof_index>>
4356  send_dof_numbers_and_indices(ghost_owners.size());
4357  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4358  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4359 
4360  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4361  {
4362  MPI_Status status;
4363  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4364  mpi_tag,
4365  tria->get_communicator(),
4366  &status);
4367  AssertThrowMPI(ierr);
4368 
4369  int len;
4370  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4371  AssertThrowMPI(ierr);
4372  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4373  ExcInternalError());
4374 
4375  const unsigned int n_cells =
4376  len / sizeof(typename CellId::binary_type);
4377  cell_data_to_send[idx].resize(n_cells);
4378 
4379  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4380  len,
4381  MPI_BYTE,
4382  status.MPI_SOURCE,
4383  status.MPI_TAG,
4384  tria->get_communicator(),
4385  &status);
4386  AssertThrowMPI(ierr);
4387 
4388  // store data for each cell
4389  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4390  {
4391  const auto cell =
4392  tria->create_cell_iterator(CellId(cell_data_to_send[idx][c]));
4393 
4394  MeshCellIteratorType mesh_it(tria,
4395  cell->level(),
4396  cell->index(),
4397  &mesh);
4398  const std_cxx17::optional<DataType> data = pack(mesh_it);
4399 
4400  if (data)
4401  {
4402  typename DestinationToBufferMap::iterator p =
4403  destination_to_data_buffer_map
4404  .insert(std::make_pair(
4405  idx,
4406  GridTools::CellDataTransferBuffer<dim, DataType>()))
4407  .first;
4408 
4409  p->second.cell_ids.emplace_back(cell->id());
4410  p->second.data.emplace_back(*data);
4411  }
4412  }
4413 
4414  // send reply
4415  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4416  destination_to_data_buffer_map[idx];
4417 
4418  sendbuffers[idx] =
4419  Utilities::pack(data, /*enable_compression*/ false);
4420  ierr = MPI_Isend(sendbuffers[idx].data(),
4421  sendbuffers[idx].size(),
4422  MPI_BYTE,
4423  status.MPI_SOURCE,
4424  mpi_tag_reply,
4425  tria->get_communicator(),
4426  &reply_requests[idx]);
4427  AssertThrowMPI(ierr);
4428  }
4429 
4430  // finally receive the replies
4431  std::vector<char> receive;
4432  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4433  {
4434  MPI_Status status;
4435  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4436  mpi_tag_reply,
4437  tria->get_communicator(),
4438  &status);
4439  AssertThrowMPI(ierr);
4440 
4441  int len;
4442  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4443  AssertThrowMPI(ierr);
4444 
4445  receive.resize(len);
4446 
4447  char *ptr = receive.data();
4448  ierr = MPI_Recv(ptr,
4449  len,
4450  MPI_BYTE,
4451  status.MPI_SOURCE,
4452  status.MPI_TAG,
4453  tria->get_communicator(),
4454  &status);
4455  AssertThrowMPI(ierr);
4456 
4457  auto cellinfo =
4458  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4459  receive, /*enable_compression*/ false);
4460 
4461  DataType *data = cellinfo.data.data();
4462  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4463  {
4465  tria_cell = tria->create_cell_iterator(cellinfo.cell_ids[c]);
4466 
4467  MeshCellIteratorType cell(tria,
4468  tria_cell->level(),
4469  tria_cell->index(),
4470  &mesh);
4471 
4472  unpack(cell, *data);
4473  }
4474  }
4475 
4476  // make sure that all communication is finished
4477  // when we leave this function.
4478  if (requests.size() > 0)
4479  {
4480  const int ierr =
4481  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4482  AssertThrowMPI(ierr);
4483  }
4484  if (reply_requests.size() > 0)
4485  {
4486  const int ierr = MPI_Waitall(reply_requests.size(),
4487  reply_requests.data(),
4488  MPI_STATUSES_IGNORE);
4489  AssertThrowMPI(ierr);
4490  }
4491 
4492 
4493 # endif // DEAL_II_WITH_MPI
4494  }
4495 
4496  } // namespace internal
4497 
4498  template <typename DataType, typename MeshType>
4499  inline void
4501  const MeshType & mesh,
4502  const std::function<std_cxx17::optional<DataType>(
4503  const typename MeshType::active_cell_iterator &)> &pack,
4504  const std::function<void(const typename MeshType::active_cell_iterator &,
4505  const DataType &)> & unpack,
4506  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4507  &cell_filter)
4508  {
4509 # ifndef DEAL_II_WITH_MPI
4510  (void)mesh;
4511  (void)pack;
4512  (void)unpack;
4513  (void)cell_filter;
4514  Assert(false, ExcNeedsMPI());
4515 # else
4516  internal::exchange_cell_data<DataType,
4517  MeshType,
4518  typename MeshType::active_cell_iterator>(
4519  mesh,
4520  pack,
4521  unpack,
4522  cell_filter,
4523  [&](const auto &process) {
4524  for (const auto &cell : mesh.active_cell_iterators())
4525  if (cell->is_ghost())
4526  process(cell, cell->subdomain_id());
4527  },
4528  [](const auto &tria) { return tria.ghost_owners(); });
4529 # endif
4530  }
4531 
4532 
4533 
4534  template <typename DataType, typename MeshType>
4535  inline void
4537  const MeshType & mesh,
4538  const std::function<std_cxx17::optional<DataType>(
4539  const typename MeshType::level_cell_iterator &)> &pack,
4540  const std::function<void(const typename MeshType::level_cell_iterator &,
4541  const DataType &)> & unpack,
4542  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4543  &cell_filter)
4544  {
4545 # ifndef DEAL_II_WITH_MPI
4546  (void)mesh;
4547  (void)pack;
4548  (void)unpack;
4549  (void)cell_filter;
4550  Assert(false, ExcNeedsMPI());
4551 # else
4552  internal::exchange_cell_data<DataType,
4553  MeshType,
4554  typename MeshType::level_cell_iterator>(
4555  mesh,
4556  pack,
4557  unpack,
4558  cell_filter,
4559  [&](const auto &process) {
4560  for (const auto &cell : mesh.cell_iterators())
4561  if (cell->level_subdomain_id() !=
4563  !cell->is_locally_owned_on_level())
4564  process(cell, cell->level_subdomain_id());
4565  },
4566  [](const auto &tria) { return tria.level_ghost_owners(); });
4567 # endif
4568  }
4569 } // namespace GridTools
4570 
4571 # endif
4572 
4574 
4575 /*---------------------------- grid_tools.h ---------------------------*/
4576 /* end of #ifndef dealii_grid_tools_h */
4577 #endif
4578 /*---------------------------- 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:5110
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:4895
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:6273
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:213
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
unsigned int n_vertices() const
const types::manifold_id flat_manifold_id
Definition: types.h:264
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static ::ExceptionBase & ExcNeedsMPI()
const unsigned int n_subdivisions
Definition: grid_tools.h:3380
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2184
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:4870
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1920
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:80
virtual MPI_Comm get_communicator() const
Definition: tria.cc:11233
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:310
Contents is actually a matrix.
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5078
unsigned int find_closest_vertex_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:2960
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void load(Archive &ar, const unsigned int version)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
return_type guess_point_owner(const RTree< std::pair< BoundingBox< spacedim >, unsigned int >> &covering_rtree, const std::vector< Point< spacedim >> &points)
Definition: grid_tools.cc:3237
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:5172
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2041
const ::Triangulation< dim, spacedim > & tria
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3287
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:69
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:13067
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:3336
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
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)
std::pair< Point< MeshType::space_dimension >, Point< MeshType::space_dimension > > compute_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate)
std::bitset< 3 > orientation
Definition: grid_tools.h:2658
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:6113
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4120
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:150
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
const double angle
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:414
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm &mpi_communicator)
Definition: grid_tools.cc:6177
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:624
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3723
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)
void partition_triangulation(const unsigned int n_partitions, const std::vector< unsigned int > &cell_weights, const SparsityPattern &cell_connection_graph, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)
Definition: grid_tools.cc:3934
cell_iterator end() const
Definition: tria.cc:13158
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:532
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:4335
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
void process_sub_cell(const std::array< unsigned int, n_configurations > &cut_line_table, const ndarray< unsigned int, n_configurations, n_cols > &new_line_table, const ndarray< unsigned int, n_lines, 2 > &line_to_vertex_table, const std::vector< value_type > &ls_values, const std::vector< Point< dim >> &points, const std::vector< unsigned int > &mask, const double iso_level, const double tolerance, std::vector< Point< dim >> &vertices, std::vector< CellData< dim - 1 >> &cells)
Definition: grid_tools.cc:6657
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1089
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:100
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< std::set< typename MeshType< dim, spacedim >::active_cell_iterator >> &vertex_to_cell_map, const std::vector< std::vector< Tensor< 1, spacedim >>> &vertex_to_cell_centers, const typename MeshType< dim, spacedim >::active_cell_iterator &cell_hint=typename MeshType< dim, spacedim >::active_cell_iterator(), const std::vector< bool > &marked_vertices={}, const RTree< std::pair< Point< spacedim >, unsigned int >> &used_vertices_rtree=RTree< std::pair< Point< spacedim >, unsigned int >>{}, const double tolerance=1.e-10, const RTree< std::pair< BoundingBox< spacedim >, typename Triangulation< dim, spacedim >::active_cell_iterator >> *relevant_cell_bounding_boxes_rtree=nullptr)
Definition: grid_tools.cc:2758
std::map< types::global_dof_index, std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandler< dim, spacedim > &dof_handler)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
T sum(const T &t, const MPI_Comm &mpi_communicator)
void get_vertex_connectivity_of_cells_on_level(const Triangulation< dim, spacedim > &triangulation, const unsigned int level, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3752
typename VectorType::value_type value_type
Definition: grid_tools.h:3299
bool orthogonal_equality(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 >())
#define Assert(cond, exc)
Definition: exceptions.h:1461
Signals signals
Definition: tria.h:2289
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3228
Abstract base class for mapping classes.
Definition: mapping.h:303
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:4362
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1254
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:698
#define DeclException0(Exception0)
Definition: exceptions.h:464
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:3042
void process(const DoFHandler< dim > &background_dof_handler, const VectorType &ls_vector, const double iso_level, std::vector< Point< dim >> &vertices, std::vector< CellData< dim - 1 >> &cells) const
Definition: grid_tools.cc:6564
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
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:5143
Point< 3 > vertices[4]
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:6336
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:286
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:4962
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:4828
cell_iterator create_cell_iterator(const CellId &cell_id) const
Definition: tria.cc:13133
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1218
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13697
Definition: cell_id.h:70
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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:2695
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:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
__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:5516
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1776
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
double cell_measure(const std::vector< Point< dim >> &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:452
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition: grid_tools.cc:2216
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:136
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:2567
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:376
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:4278
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3689
DistributedComputePointLocationsInternal< dim, spacedim > distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes, const double tolerance, const bool perform_handshake, const bool enforce_unique_mapping=false)
Definition: grid_tools.cc:5874
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1321
Point< 3 > center
FullMatrix< double > matrix
Definition: grid_tools.h:2672
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:4990
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:2030
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:5487
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4308
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1114
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:4015
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 > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
numbers::NumberTraits< Number >::real_type norm() const
unsigned int find_closest_vertex(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={})
Definition: grid_tools.cc:2505
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
#define DEAL_II_DEPRECATED
Definition: config.h:160
void invert_all_negative_measure_cells(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:861
static ::ExceptionBase & ExcMeshNotOrientable()
void collect_periodic_faces(const MeshType &mesh, const types::boundary_id b_id, 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 >())
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2021
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:4292
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:729
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6459