Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_tools.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 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 
116 
117 # ifdef _MSC_VER
118  template <int dim, int spacedim>
119  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim>>
120  {
121  public:
122  using type =
124  };
125 # endif
126 } // namespace internal
127 
136 namespace GridTools
137 {
138  template <int dim, int spacedim>
139  class Cache;
140 
145 
152  template <int dim, int spacedim>
153  double
155 
182  template <int dim, int spacedim>
183  double
185  const Mapping<dim, spacedim> & mapping =
186  (ReferenceCells::get_hypercube<dim>()
187  .template get_default_linear_mapping<dim, spacedim>()));
188 
199  template <int dim, int spacedim>
200  double
203  const Mapping<dim, spacedim> & mapping =
204  (ReferenceCells::get_hypercube<dim>()
205  .template get_default_linear_mapping<dim, spacedim>()));
206 
217  template <int dim, int spacedim>
218  double
220  const Triangulation<dim, spacedim> &triangulation,
221  const Mapping<dim, spacedim> & mapping =
222  (ReferenceCells::get_hypercube<dim>()
223  .template get_default_linear_mapping<dim, spacedim>()));
224 
236  template <int dim>
238  cell_measure(
239  const std::vector<Point<dim>> &all_vertices,
241 
258  template <int dim>
259  double
260  cell_measure(const std::vector<Point<dim>> & all_vertices,
262 
285  template <int dim, int spacedim>
286  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
288 
315  template <int dim>
316  Vector<double>
318  const Triangulation<dim> &triangulation,
319  const Quadrature<dim> & quadrature);
320 
328  template <int dim>
329  double
331  const Triangulation<dim> &triangulation,
332  const Quadrature<dim> & quadrature);
333 
347  template <int dim, int spacedim>
350 
368  template <typename Iterator>
371  const Iterator & object,
373 
386  template <int dim, int spacedim>
387  std::
388  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
390 
396 
413  template <int dim, int spacedim>
414  void
416  std::vector<CellData<dim>> & cells,
417  SubCellData & subcelldata);
418 
436  template <int dim, int spacedim>
437  void
438  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
439  std::vector<CellData<dim>> & cells,
440  SubCellData & subcelldata,
441  std::vector<unsigned int> & considered_vertices,
442  const double tol = 1e-12);
443 
462  template <int dim, int spacedim>
463  void
465  const std::vector<Point<spacedim>> &all_vertices,
466  std::vector<CellData<dim>> & cells);
467 
477  template <int dim>
478  void
479  consistently_order_cells(std::vector<CellData<dim>> &cells);
480 
486 
540  template <int dim, typename Transformation, int spacedim>
541  void
542  transform(const Transformation & transformation,
543  Triangulation<dim, spacedim> &triangulation);
544 
550  template <int dim, int spacedim>
551  void
552  shift(const Tensor<1, spacedim> & shift_vector,
553  Triangulation<dim, spacedim> &triangulation);
554 
555 
565  template <int dim>
566  void
567  rotate(const double angle, Triangulation<dim> &triangulation);
568 
581  template <int dim>
582  void
583  rotate(const double angle,
584  const unsigned int axis,
585  Triangulation<dim, 3> &triangulation);
586 
644  template <int dim>
645  void
646  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
647  Triangulation<dim> & tria,
648  const Function<dim, double> *coefficient = nullptr,
649  const bool solve_for_absolute_positions = false);
650 
656  template <int dim, int spacedim>
657  std::map<unsigned int, Point<spacedim>>
659 
667  template <int dim, int spacedim>
668  void
669  scale(const double scaling_factor,
670  Triangulation<dim, spacedim> &triangulation);
671 
686  template <int dim, int spacedim>
687  void
689  const double factor,
690  Triangulation<dim, spacedim> &triangulation,
691  const bool keep_boundary = true,
692  const unsigned int seed = boost::random::mt19937::default_seed);
693 
727  template <int dim, int spacedim>
728  void
730  const bool isotropic = false,
731  const unsigned int max_iterations = 100);
732 
757  template <int dim, int spacedim>
758  void
760  const double max_ratio = 1.6180339887,
761  const unsigned int max_iterations = 5);
762 
852  template <int dim, int spacedim>
853  void
855  const double limit_angle_fraction = .75);
856 
862 
917  template <int dim, int spacedim>
918 # ifndef DOXYGEN
919  std::tuple<
920  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
921  std::vector<std::vector<Point<dim>>>,
922  std::vector<std::vector<unsigned int>>>
923 # else
924  return_type
925 # endif
927  const Cache<dim, spacedim> & cache,
928  const std::vector<Point<spacedim>> &points,
930  &cell_hint =
932 
966  template <int dim, int spacedim>
967 # ifndef DOXYGEN
968  std::tuple<
969  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
970  std::vector<std::vector<Point<dim>>>,
971  std::vector<std::vector<unsigned int>>,
972  std::vector<unsigned int>>
973 # else
974  return_type
975 # endif
977  const Cache<dim, spacedim> & cache,
978  const std::vector<Point<spacedim>> &points,
980  &cell_hint =
982 
1052  template <int dim, int spacedim>
1053 # ifndef DOXYGEN
1054  std::tuple<
1055  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1056  std::vector<std::vector<Point<dim>>>,
1057  std::vector<std::vector<unsigned int>>,
1058  std::vector<std::vector<Point<spacedim>>>,
1059  std::vector<std::vector<unsigned int>>>
1060 # else
1061  return_type
1062 # endif
1064  const GridTools::Cache<dim, spacedim> & cache,
1065  const std::vector<Point<spacedim>> & local_points,
1066  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1067  const double tolerance = 1e-10);
1068 
1069  namespace internal
1070  {
1085  template <int dim, int spacedim>
1087  {
1094  std::vector<std::tuple<std::pair<int, int>,
1095  unsigned int,
1096  unsigned int,
1097  Point<dim>,
1099  unsigned int>>
1101 
1105  std::vector<unsigned int> send_ranks;
1106 
1112  std::vector<unsigned int> send_ptrs;
1113 
1124  std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1126 
1130  std::vector<unsigned int> recv_ranks;
1131 
1137  std::vector<unsigned int> recv_ptrs;
1138  };
1139 
1149  template <int dim, int spacedim>
1152  const GridTools::Cache<dim, spacedim> & cache,
1153  const std::vector<Point<spacedim>> & points,
1154  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1155  const double tolerance,
1156  const bool perform_handshake,
1157  const bool enforce_unique_mapping = false);
1158 
1159  } // namespace internal
1160 
1197  template <int dim, int spacedim>
1198  std::map<unsigned int, Point<spacedim>>
1200  const Triangulation<dim, spacedim> &container,
1201  const Mapping<dim, spacedim> & mapping =
1202  (ReferenceCells::get_hypercube<dim>()
1203  .template get_default_linear_mapping<dim, spacedim>()));
1204 
1214  template <int spacedim>
1215  unsigned int
1216  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1217  const Point<spacedim> & p);
1218 
1242  template <int dim, template <int, int> class MeshType, int spacedim>
1243  unsigned int
1244  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1245  const Point<spacedim> & p,
1246  const std::vector<bool> & marked_vertices = {});
1247 
1271  template <int dim, template <int, int> class MeshType, int spacedim>
1272  unsigned int
1274  const MeshType<dim, spacedim> &mesh,
1275  const Point<spacedim> & p,
1276  const std::vector<bool> & marked_vertices = {});
1277 
1278 
1298  template <int dim, template <int, int> class MeshType, int spacedim>
1299 # ifndef _MSC_VER
1300  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1301 # else
1302  std::vector<
1303  typename ::internal::
1304  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1305 # endif
1306  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1307  const unsigned int vertex_index);
1308 
1371  template <int dim, template <int, int> class MeshType, int spacedim>
1372 # ifndef _MSC_VER
1373  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1374 # else
1375  std::pair<typename ::internal::
1376  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1377  Point<dim>>
1378 # endif
1380  const MeshType<dim, spacedim> &mesh,
1381  const Point<spacedim> & p,
1382  const std::vector<bool> &marked_vertices = {},
1383  const double tolerance = 1.e-10);
1384 
1392  template <int dim, template <int, int> class MeshType, int spacedim>
1393 # ifndef _MSC_VER
1394  typename MeshType<dim, spacedim>::active_cell_iterator
1395 # else
1396  typename ::internal::
1397  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1398 # endif
1399  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1400  const Point<spacedim> & p,
1401  const std::vector<bool> &marked_vertices = {},
1402  const double tolerance = 1.e-10);
1403 
1410  template <int dim, int spacedim>
1411  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1412  Point<dim>>
1414  const hp::MappingCollection<dim, spacedim> &mapping,
1415  const DoFHandler<dim, spacedim> & mesh,
1416  const Point<spacedim> & p,
1417  const double tolerance = 1.e-10);
1418 
1470  template <int dim, int spacedim>
1471  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1472  Point<dim>>
1474  const Cache<dim, spacedim> &cache,
1475  const Point<spacedim> & p,
1478  const std::vector<bool> &marked_vertices = {},
1479  const double tolerance = 1.e-10);
1480 
1494  template <int dim, template <int, int> class MeshType, int spacedim>
1495 # ifndef _MSC_VER
1496  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1497 # else
1498  std::pair<typename ::internal::
1499  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1500  Point<dim>>
1501 # endif
1503  const Mapping<dim, spacedim> & mapping,
1504  const MeshType<dim, spacedim> &mesh,
1505  const Point<spacedim> & p,
1506  const std::vector<
1507  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1509  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1510  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1511  typename MeshType<dim, spacedim>::active_cell_iterator(),
1512  const std::vector<bool> & marked_vertices = {},
1513  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1514  RTree<std::pair<Point<spacedim>, unsigned int>>{},
1515  const double tolerance = 1.e-10);
1516 
1538  template <int dim, template <int, int> class MeshType, int spacedim>
1539 # ifndef _MSC_VER
1540  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1541  Point<dim>>>
1542 # else
1543  std::vector<std::pair<
1544  typename ::internal::
1545  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1546  Point<dim>>>
1547 # endif
1549  const Mapping<dim, spacedim> & mapping,
1550  const MeshType<dim, spacedim> &mesh,
1551  const Point<spacedim> & p,
1552  const double tolerance,
1553  const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1554  Point<dim>> & first_cell);
1555 
1562  template <int dim, template <int, int> class MeshType, int spacedim>
1563 # ifndef _MSC_VER
1564  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1565  Point<dim>>>
1566 # else
1567  std::vector<std::pair<
1568  typename ::internal::
1569  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1570  Point<dim>>>
1571 # endif
1573  const Mapping<dim, spacedim> & mapping,
1574  const MeshType<dim, spacedim> &mesh,
1575  const Point<spacedim> & p,
1576  const double tolerance = 1e-10,
1577  const std::vector<bool> & marked_vertices = {});
1578 
1600  template <class MeshType>
1601  std::vector<typename MeshType::active_cell_iterator>
1602  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1603 
1628  template <class MeshType>
1629  void
1631  const typename MeshType::active_cell_iterator & cell,
1632  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1633 
1683  template <class MeshType>
1684  std::vector<typename MeshType::active_cell_iterator>
1686  const MeshType &mesh,
1687  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1688  &predicate);
1689 
1690 
1698  template <class MeshType>
1699  std::vector<typename MeshType::cell_iterator>
1701  const MeshType &mesh,
1702  const std::function<bool(const typename MeshType::cell_iterator &)>
1703  & predicate,
1704  const unsigned int level);
1705 
1706 
1719  template <class MeshType>
1720  std::vector<typename MeshType::active_cell_iterator>
1721  compute_ghost_cell_halo_layer(const MeshType &mesh);
1722 
1771  template <class MeshType>
1772  std::vector<typename MeshType::active_cell_iterator>
1774  const MeshType &mesh,
1775  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1776  & predicate,
1777  const double layer_thickness);
1778 
1801  template <class MeshType>
1802  std::vector<typename MeshType::active_cell_iterator>
1803  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1804  const double layer_thickness);
1805 
1821  template <class MeshType>
1822  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1824  const MeshType &mesh,
1825  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1826  &predicate);
1827 
1880  template <class MeshType>
1881  std::vector<BoundingBox<MeshType::space_dimension>>
1883  const MeshType &mesh,
1884  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1885  & predicate,
1886  const unsigned int refinement_level = 0,
1887  const bool allow_merge = false,
1888  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1889 
1917  template <int spacedim>
1918 # ifndef DOXYGEN
1919  std::tuple<std::vector<std::vector<unsigned int>>,
1920  std::map<unsigned int, unsigned int>,
1921  std::map<unsigned int, std::vector<unsigned int>>>
1922 # else
1923  return_type
1924 # endif
1926  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1927  const std::vector<Point<spacedim>> & points);
1928 
1929 
1964  template <int spacedim>
1965 # ifndef DOXYGEN
1966  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1967  std::map<unsigned int, unsigned int>,
1968  std::map<unsigned int, std::vector<unsigned int>>>
1969 # else
1970  return_type
1971 # endif
1973  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1974  const std::vector<Point<spacedim>> & points);
1975 
1976 
1985  template <int dim, int spacedim>
1986  std::vector<
1987  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1988  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1989 
2002  template <int dim, int spacedim>
2003  std::vector<std::vector<Tensor<1, spacedim>>>
2005  const Triangulation<dim, spacedim> &mesh,
2006  const std::vector<
2008  &vertex_to_cells);
2009 
2010 
2018  template <int dim, int spacedim>
2019  unsigned int
2022  const Point<spacedim> & position,
2023  const Mapping<dim, spacedim> & mapping =
2024  (ReferenceCells::get_hypercube<dim>()
2025  .template get_default_linear_mapping<dim, spacedim>()));
2026 
2038  template <int dim, int spacedim>
2039  std::map<unsigned int, types::global_vertex_index>
2042 
2054  template <int dim, int spacedim>
2055  std::pair<unsigned int, double>
2058 
2064 
2073  template <int dim, int spacedim>
2074  void
2076  const Triangulation<dim, spacedim> &triangulation,
2077  DynamicSparsityPattern & connectivity);
2078 
2087  template <int dim, int spacedim>
2088  void
2090  const Triangulation<dim, spacedim> &triangulation,
2091  DynamicSparsityPattern & connectivity);
2092 
2101  template <int dim, int spacedim>
2102  void
2104  const Triangulation<dim, spacedim> &triangulation,
2105  const unsigned int level,
2106  DynamicSparsityPattern & connectivity);
2107 
2128  template <int dim, int spacedim>
2129  void
2130  partition_triangulation(const unsigned int n_partitions,
2131  Triangulation<dim, spacedim> & triangulation,
2132  const SparsityTools::Partitioner partitioner =
2134 
2145  template <int dim, int spacedim>
2146  void
2147  partition_triangulation(const unsigned int n_partitions,
2148  const std::vector<unsigned int> &cell_weights,
2149  Triangulation<dim, spacedim> & triangulation,
2150  const SparsityTools::Partitioner partitioner =
2152 
2198  template <int dim, int spacedim>
2199  void
2200  partition_triangulation(const unsigned int n_partitions,
2201  const SparsityPattern & cell_connection_graph,
2202  Triangulation<dim, spacedim> &triangulation,
2203  const SparsityTools::Partitioner partitioner =
2205 
2216  template <int dim, int spacedim>
2217  void
2218  partition_triangulation(const unsigned int n_partitions,
2219  const std::vector<unsigned int> &cell_weights,
2220  const SparsityPattern & cell_connection_graph,
2221  Triangulation<dim, spacedim> &triangulation,
2222  const SparsityTools::Partitioner partitioner =
2224 
2239  template <int dim, int spacedim>
2240  void
2241  partition_triangulation_zorder(const unsigned int n_partitions,
2242  Triangulation<dim, spacedim> &triangulation,
2243  const bool group_siblings = true);
2244 
2256  template <int dim, int spacedim>
2257  void
2259 
2267  template <int dim, int spacedim>
2268  std::vector<types::subdomain_id>
2270  const std::vector<CellId> & cell_ids);
2271 
2282  template <int dim, int spacedim>
2283  void
2285  std::vector<types::subdomain_id> & subdomain);
2286 
2301  template <int dim, int spacedim>
2302  unsigned int
2304  const Triangulation<dim, spacedim> &triangulation,
2305  const types::subdomain_id subdomain);
2306 
2336  template <int dim, int spacedim>
2337  std::vector<bool>
2339 
2345 
2379  template <typename MeshType>
2380  std::list<std::pair<typename MeshType::cell_iterator,
2381  typename MeshType::cell_iterator>>
2382  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2383 
2393  template <int dim, int spacedim>
2394  bool
2396  const Triangulation<dim, spacedim> &mesh_2);
2397 
2407  template <typename MeshType>
2408  bool
2409  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2410 
2416 
2432  template <int dim, int spacedim>
2436  & distorted_cells,
2437  Triangulation<dim, spacedim> &triangulation);
2438 
2439 
2440 
2449 
2450 
2488  template <class MeshType>
2489  std::vector<typename MeshType::active_cell_iterator>
2490  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2491 
2492 
2514  template <class Container>
2515  std::vector<typename Container::cell_iterator>
2517  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2518 
2585  template <class Container>
2586  void
2588  const std::vector<typename Container::active_cell_iterator> &patch,
2590  &local_triangulation,
2591  std::map<
2592  typename Triangulation<Container::dimension,
2593  Container::space_dimension>::active_cell_iterator,
2594  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2595 
2627  template <int dim, int spacedim>
2628  std::map<
2630  std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2632 
2633 
2640 
2646  template <typename CellIterator>
2648  {
2652  CellIterator cell[2];
2653 
2658  unsigned int face_idx[2];
2659 
2665  std::bitset<3> orientation;
2666 
2680  };
2681 
2682 
2746  template <typename FaceIterator>
2747  bool orthogonal_equality(
2748  std::bitset<3> & orientation,
2749  const FaceIterator & face1,
2750  const FaceIterator & face2,
2751  const int direction,
2755 
2756 
2760  template <typename FaceIterator>
2761  bool
2763  const FaceIterator & face1,
2764  const FaceIterator & face2,
2765  const int direction,
2769 
2770 
2827  template <typename MeshType>
2828  void
2830  const MeshType & mesh,
2831  const types::boundary_id b_id1,
2832  const types::boundary_id b_id2,
2833  const int direction,
2835  & matched_pairs,
2836  const Tensor<1, MeshType::space_dimension> &offset =
2839 
2840 
2863  template <typename MeshType>
2864  void
2866  const MeshType & mesh,
2867  const types::boundary_id b_id,
2868  const int direction,
2870  & matched_pairs,
2871  const ::Tensor<1, MeshType::space_dimension> &offset =
2874 
2880 
2901  template <int dim, int spacedim>
2902  void
2904  const bool reset_boundary_ids = false);
2905 
2927  template <int dim, int spacedim>
2928  void
2930  const std::vector<types::boundary_id> &src_boundary_ids,
2931  const std::vector<types::manifold_id> &dst_manifold_ids,
2933  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2934 
2964  template <int dim, int spacedim>
2965  void
2967  const bool compute_face_ids = false);
2968 
2993  template <int dim, int spacedim>
2994  void
2997  const std::function<types::manifold_id(
2998  const std::set<types::manifold_id> &)> &disambiguation_function =
2999  [](const std::set<types::manifold_id> &manifold_ids) {
3000  if (manifold_ids.size() == 1)
3001  return *manifold_ids.begin();
3002  else
3004  },
3005  bool overwrite_only_flat_manifold_ids = true);
3092  template <typename DataType, typename MeshType>
3093  void
3095  const MeshType & mesh,
3096  const std::function<std_cxx17::optional<DataType>(
3097  const typename MeshType::active_cell_iterator &)> &pack,
3098  const std::function<void(const typename MeshType::active_cell_iterator &,
3099  const DataType &)> & unpack,
3100  const std::function<bool(const typename MeshType::active_cell_iterator &)>
3101  &cell_filter =
3103 
3114  template <typename DataType, typename MeshType>
3115  void
3117  const MeshType & mesh,
3118  const std::function<std_cxx17::optional<DataType>(
3119  const typename MeshType::level_cell_iterator &)> &pack,
3120  const std::function<void(const typename MeshType::level_cell_iterator &,
3121  const DataType &)> & unpack,
3122  const std::function<bool(const typename MeshType::level_cell_iterator &)> &
3124  true});
3125 
3126  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
3127  * boxes @p local_bboxes.
3128  *
3129  * This function is meant to exchange bounding boxes describing the locally
3130  * owned cells in a distributed triangulation obtained with the function
3131  * GridTools::compute_mesh_predicate_bounding_box .
3132  *
3133  * The output vector's size is the number of processes of the MPI
3134  * communicator:
3135  * its i-th entry contains the vector @p local_bboxes of the i-th process.
3136  */
3137  template <int spacedim>
3138  std::vector<std::vector<BoundingBox<spacedim>>>
3140  const std::vector<BoundingBox<spacedim>> &local_bboxes,
3141  const MPI_Comm & mpi_communicator);
3142 
3175  template <int spacedim>
3178  const std::vector<BoundingBox<spacedim>> &local_description,
3179  const MPI_Comm & mpi_communicator);
3180 
3198  template <int dim, int spacedim>
3199  void
3201  const Triangulation<dim, spacedim> & tria,
3202  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3203  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3204 
3211  template <int dim, int spacedim>
3212  std::map<unsigned int, std::set<::types::subdomain_id>>
3214  const Triangulation<dim, spacedim> &tria);
3215 
3228  template <int dim, typename T>
3230  {
3234  std::vector<CellId> cell_ids;
3235 
3239  std::vector<T> data;
3240 
3249  template <class Archive>
3250  void
3251  save(Archive &ar, const unsigned int version) const;
3252 
3259  template <class Archive>
3260  void
3261  load(Archive &ar, const unsigned int version);
3262 
3263 # ifdef DOXYGEN
3264 
3269  template <class Archive>
3270  void
3271  serialize(Archive &archive, const unsigned int version);
3272 # else
3273  // This macro defines the serialize() method that is compatible with
3274  // the templated save() and load() method that have been implemented.
3275  BOOST_SERIALIZATION_SPLIT_MEMBER()
3276 # endif
3277  };
3278 
3279 
3280 
3298  template <int dim, typename VectorType>
3300  {
3301  public:
3305  using value_type = typename VectorType::value_type;
3306 
3310  MarchingCubeAlgorithm(const Mapping<dim, dim> & mapping,
3311  const FiniteElement<dim, dim> &fe,
3312  const unsigned int n_subdivisions = 1,
3313  const double tolerance = 1e-10);
3314 
3319  void
3320  process(const DoFHandler<dim> & background_dof_handler,
3321  const VectorType & ls_vector,
3322  const double iso_level,
3323  std::vector<Point<dim>> & vertices,
3324  std::vector<CellData<dim - 1>> &cells) const;
3325 
3332  void
3333  process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
3334  const VectorType & ls_vector,
3335  const double iso_level,
3336  std::vector<Point<dim>> & vertices,
3337  std::vector<CellData<dim - 1>> &cells) const;
3338 
3339  private:
3344  static Quadrature<dim>
3345  create_quadrature_rule(const unsigned int n_subdivisions);
3346 
3350  void
3351  process_cell(std::vector<value_type> & ls_values,
3352  const std::vector<Point<dim>> & points,
3353  const double iso_level,
3354  std::vector<Point<dim>> & vertices,
3355  std::vector<CellData<dim - 1>> &cells) const;
3356 
3363  void
3364  process_sub_cell(const std::vector<value_type> & ls_values,
3365  const std::vector<Point<2>> & points,
3366  const std::vector<unsigned int> mask,
3367  const double iso_level,
3368  std::vector<Point<2>> & vertices,
3369  std::vector<CellData<1>> & cells) const;
3370 
3374  void
3375  process_sub_cell(const std::vector<value_type> & ls_values,
3376  const std::vector<Point<3>> & points,
3377  const std::vector<unsigned int> mask,
3378  const double iso_level,
3379  std::vector<Point<3>> & vertices,
3380  std::vector<CellData<2>> & cells) const;
3381 
3386  const unsigned int n_subdivisions;
3387 
3392  const double tolerance;
3393 
3399  };
3400 
3401 
3402 
3407 
3412  int,
3413  << "The number of partitions you gave is " << arg1
3414  << ", but must be greater than zero.");
3419  int,
3420  << "The subdomain id " << arg1
3421  << " has no cells associated with it.");
3426 
3431  double,
3432  << "The scaling factor must be positive, but it is " << arg1
3433  << ".");
3434 
3439  unsigned int,
3440  << "The given vertex with index " << arg1
3441  << " is not used in the given triangulation.");
3442 
3445 } /*namespace GridTools*/
3446 
3447 
3458  "The edges of the mesh are not consistently orientable.");
3459 
3460 
3461 
3462 /* ----------------- Template function --------------- */
3463 
3464 # ifndef DOXYGEN
3465 
3466 namespace GridTools
3467 {
3468  template <int dim>
3469  double
3470  cell_measure(
3471  const std::vector<Point<dim>> &all_vertices,
3472  const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
3473  {
3474  // We forward call to the ArrayView version:
3475  const ArrayView<const unsigned int> view(
3476  indices, GeometryInfo<dim>::vertices_per_cell);
3477  return cell_measure(all_vertices, view);
3478  }
3479 
3480 
3481 
3482  // This specialization is defined here so that the general template in the
3483  // source file doesn't need to have further 1D overloads for the internal
3484  // functions it calls.
3485  template <>
3489  {
3490  return {};
3491  }
3492 
3493 
3494 
3495  template <int dim, typename Predicate, int spacedim>
3496  void
3497  transform(const Predicate & predicate,
3498  Triangulation<dim, spacedim> &triangulation)
3499  {
3500  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3501 
3502  // loop over all active cells, and
3503  // transform those vertices that
3504  // have not yet been touched. note
3505  // that we get to all vertices in
3506  // the triangulation by only
3507  // visiting the active cells.
3509  cell = triangulation.begin_active(),
3510  endc = triangulation.end();
3511  for (; cell != endc; ++cell)
3512  for (const unsigned int v : cell->vertex_indices())
3513  if (treated_vertices[cell->vertex_index(v)] == false)
3514  {
3515  // transform this vertex
3516  cell->vertex(v) = predicate(cell->vertex(v));
3517  // and mark it as treated
3518  treated_vertices[cell->vertex_index(v)] = true;
3519  };
3520 
3521 
3522  // now fix any vertices on hanging nodes so that we don't create any holes
3523  if (dim == 2)
3524  {
3526  cell = triangulation.begin_active(),
3527  endc = triangulation.end();
3528  for (; cell != endc; ++cell)
3529  for (const unsigned int face : cell->face_indices())
3530  if (cell->face(face)->has_children() &&
3531  !cell->face(face)->at_boundary())
3532  {
3533  Assert(cell->reference_cell() ==
3534  ReferenceCells::get_hypercube<dim>(),
3535  ExcNotImplemented());
3536 
3537  // this line has children
3538  cell->face(face)->child(0)->vertex(1) =
3539  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3540  2;
3541  }
3542  }
3543  else if (dim == 3)
3544  {
3546  cell = triangulation.begin_active(),
3547  endc = triangulation.end();
3548  for (; cell != endc; ++cell)
3549  for (const unsigned int face : cell->face_indices())
3550  if (cell->face(face)->has_children() &&
3551  !cell->face(face)->at_boundary())
3552  {
3553  Assert(cell->reference_cell() ==
3554  ReferenceCells::get_hypercube<dim>(),
3555  ExcNotImplemented());
3556 
3557  // this face has hanging nodes
3558  cell->face(face)->child(0)->vertex(1) =
3559  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3560  2.0;
3561  cell->face(face)->child(0)->vertex(2) =
3562  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3563  2.0;
3564  cell->face(face)->child(1)->vertex(3) =
3565  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3566  2.0;
3567  cell->face(face)->child(2)->vertex(3) =
3568  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3569  2.0;
3570 
3571  // center of the face
3572  cell->face(face)->child(0)->vertex(3) =
3573  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3574  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3575  4.0;
3576  }
3577  }
3578 
3579  // Make sure FEValues notices that the mesh has changed
3580  triangulation.signals.mesh_movement();
3581  }
3582 
3583 
3584 
3585  template <class MeshType>
3586  std::vector<typename MeshType::active_cell_iterator>
3587  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3588  {
3589  std::vector<typename MeshType::active_cell_iterator> child_cells;
3590 
3591  if (cell->has_children())
3592  {
3593  for (unsigned int child = 0; child < cell->n_children(); ++child)
3594  if (cell->child(child)->has_children())
3595  {
3596  const std::vector<typename MeshType::active_cell_iterator>
3597  children = get_active_child_cells<MeshType>(cell->child(child));
3598  child_cells.insert(child_cells.end(),
3599  children.begin(),
3600  children.end());
3601  }
3602  else
3603  child_cells.push_back(cell->child(child));
3604  }
3605 
3606  return child_cells;
3607  }
3608 
3609 
3610 
3611  template <class MeshType>
3612  void
3614  const typename MeshType::active_cell_iterator & cell,
3615  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3616  {
3617  active_neighbors.clear();
3618  for (const unsigned int n : cell->face_indices())
3619  if (!cell->at_boundary(n))
3620  {
3621  if (MeshType::dimension == 1)
3622  {
3623  // check children of neighbor. note
3624  // that in 1d children of the neighbor
3625  // may be further refined. In 1d the
3626  // case is simple since we know what
3627  // children bound to the present cell
3628  typename MeshType::cell_iterator neighbor_child =
3629  cell->neighbor(n);
3630  if (!neighbor_child->is_active())
3631  {
3632  while (neighbor_child->has_children())
3633  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3634 
3635  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3636  ExcInternalError());
3637  }
3638  active_neighbors.push_back(neighbor_child);
3639  }
3640  else
3641  {
3642  if (cell->face(n)->has_children())
3643  // this neighbor has children. find
3644  // out which border to the present
3645  // cell
3646  for (unsigned int c = 0;
3647  c < cell->face(n)->n_active_descendants();
3648  ++c)
3649  active_neighbors.push_back(
3650  cell->neighbor_child_on_subface(n, c));
3651  else
3652  {
3653  // the neighbor must be active
3654  // himself
3655  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3656  active_neighbors.push_back(cell->neighbor(n));
3657  }
3658  }
3659  }
3660  }
3661 
3662 
3663 
3664  namespace internal
3665  {
3666  namespace ProjectToObject
3667  {
3680  struct CrossDerivative
3681  {
3682  const unsigned int direction_0;
3683  const unsigned int direction_1;
3684 
3685  CrossDerivative(const unsigned int d0, const unsigned int d1);
3686  };
3687 
3688  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3689  const unsigned int d1)
3690  : direction_0(d0)
3691  , direction_1(d1)
3692  {}
3693 
3694 
3695 
3700  template <typename F>
3701  inline auto
3702  centered_first_difference(const double center,
3703  const double step,
3704  const F &f) -> decltype(f(center) - f(center))
3705  {
3706  return (f(center + step) - f(center - step)) / (2.0 * step);
3707  }
3708 
3709 
3710 
3715  template <typename F>
3716  inline auto
3717  centered_second_difference(const double center,
3718  const double step,
3719  const F &f) -> decltype(f(center) - f(center))
3720  {
3721  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3722  (step * step);
3723  }
3724 
3725 
3726 
3736  template <int structdim, typename F>
3737  inline auto
3738  cross_stencil(
3739  const CrossDerivative cross_derivative,
3741  const double step,
3742  const F &f) -> decltype(f(center) - f(center))
3743  {
3745  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3746  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3747  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3748  1.0 / 3.0 * f(center - simplex_vector) +
3749  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3750  step;
3751  }
3752 
3753 
3754 
3761  template <int spacedim, int structdim, typename F>
3762  inline double
3763  gradient_entry(
3764  const unsigned int row_n,
3765  const unsigned int dependent_direction,
3766  const Point<spacedim> &p0,
3768  const double step,
3769  const F & f)
3770  {
3772  dependent_direction <
3774  ExcMessage("This function assumes that the last weight is a "
3775  "dependent variable (and hence we cannot take its "
3776  "derivative directly)."));
3777  Assert(row_n != dependent_direction,
3778  ExcMessage(
3779  "We cannot differentiate with respect to the variable "
3780  "that is assumed to be dependent."));
3781 
3782  const Point<spacedim> manifold_point = f(center);
3783  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3784  {row_n, dependent_direction}, center, step, f);
3785  double entry = 0.0;
3786  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3787  entry +=
3788  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3789  return entry;
3790  }
3791 
3797  template <typename Iterator, int spacedim, int structdim>
3799  project_to_d_linear_object(const Iterator & object,
3800  const Point<spacedim> &trial_point)
3801  {
3802  // let's look at this for simplicity for a quad (structdim==2) in a
3803  // space with spacedim>2 (notate trial_point by y): all points on the
3804  // surface are given by
3805  // x(\xi) = sum_i v_i phi_x(\xi)
3806  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3807  // reference coordinates of the quad. so what we are trying to do is
3808  // find a point x on the surface that is closest to the point y. there
3809  // are different ways to solve this problem, but in the end it's a
3810  // nonlinear problem and we have to find reference coordinates \xi so
3811  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3812  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3813  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3814  // have to use a Newton method to find the answer. This leads to the
3815  // following formulation of Newton steps:
3816  //
3817  // Given \xi_k, find \delta\xi_k so that
3818  // H_k \delta\xi_k = - F_k
3819  // where H_k is an approximation to the second derivatives of J at
3820  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3821  // number of times until the right hand side is small enough. As a
3822  // stopping criterion, we terminate if ||\delta\xi||<eps.
3823  //
3824  // As for the Hessian, the best choice would be
3825  // H_k = J''(\xi_k)
3826  // but we'll opt for the simpler Gauss-Newton form
3827  // H_k = A^T A
3828  // i.e.
3829  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3830  // \partial_n phi_i *
3831  // \partial_m phi_j
3832  // we start at xi=(0.5, 0.5).
3833  Point<structdim> xi;
3834  for (unsigned int d = 0; d < structdim; ++d)
3835  xi[d] = 0.5;
3836 
3837  Point<spacedim> x_k;
3838  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3839  x_k += object->vertex(i) *
3841 
3842  do
3843  {
3845  for (const unsigned int i :
3847  F_k +=
3848  (x_k - trial_point) * object->vertex(i) *
3850  i);
3851 
3853  for (const unsigned int i :
3855  for (const unsigned int j :
3857  {
3860  xi, i),
3862  xi, j));
3863  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3864  }
3865 
3866  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3867  xi += delta_xi;
3868 
3869  x_k = Point<spacedim>();
3870  for (const unsigned int i :
3872  x_k += object->vertex(i) *
3874 
3875  if (delta_xi.norm() < 1e-7)
3876  break;
3877  }
3878  while (true);
3879 
3880  return x_k;
3881  }
3882  } // namespace ProjectToObject
3883  } // namespace internal
3884 
3885 
3886 
3887  namespace internal
3888  {
3889  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3890  // inside the project_to_object function below.
3891  template <int structdim>
3892  inline bool
3893  weights_are_ok(
3895  {
3896  // clang has trouble figuring out structdim here, so define it
3897  // again:
3898  static const std::size_t n_vertices_per_cell =
3900  n_independent_components;
3901  std::array<double, n_vertices_per_cell> copied_weights;
3902  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3903  {
3904  copied_weights[i] = v[i];
3905  if (v[i] < 0.0 || v[i] > 1.0)
3906  return false;
3907  }
3908 
3909  // check the sum: try to avoid some roundoff errors by summing in order
3910  std::sort(copied_weights.begin(), copied_weights.end());
3911  const double sum =
3912  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3913  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3914  }
3915  } // namespace internal
3916 
3917  template <typename Iterator>
3920  const Iterator & object,
3922  {
3923  const int spacedim = Iterator::AccessorType::space_dimension;
3924  const int structdim = Iterator::AccessorType::structure_dimension;
3925 
3926  Point<spacedim> projected_point = trial_point;
3927 
3928  if (structdim >= spacedim)
3929  return projected_point;
3930  else if (structdim == 1 || structdim == 2)
3931  {
3932  using namespace internal::ProjectToObject;
3933  // Try to use the special flat algorithm for quads (this is better
3934  // than the general algorithm in 3D). This does not take into account
3935  // whether projected_point is outside the quad, but we optimize along
3936  // lines below anyway:
3937  const int dim = Iterator::AccessorType::dimension;
3938  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3939  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3940  &manifold) != nullptr)
3941  {
3942  projected_point =
3943  project_to_d_linear_object<Iterator, spacedim, structdim>(
3944  object, trial_point);
3945  }
3946  else
3947  {
3948  // We want to find a point on the convex hull (defined by the
3949  // vertices of the object and the manifold description) that is
3950  // relatively close to the trial point. This has a few issues:
3951  //
3952  // 1. For a general convex hull we are not guaranteed that a unique
3953  // minimum exists.
3954  // 2. The independent variables in the optimization process are the
3955  // weights given to Manifold::get_new_point, which must sum to 1,
3956  // so we cannot use standard finite differences to approximate a
3957  // gradient.
3958  //
3959  // There is not much we can do about 1., but for 2. we can derive
3960  // finite difference stencils that work on a structdim-dimensional
3961  // simplex and rewrite the optimization problem to use those
3962  // instead. Consider the structdim 2 case and let
3963  //
3964  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3965  // c2, c3})
3966  //
3967  // where {c0, c1, c2, c3} are the weights for the four vertices on
3968  // the quadrilateral. We seek to minimize the Euclidean distance
3969  // between F(...) and trial_point. We can solve for c3 in terms of
3970  // the other weights and get, for one coordinate direction
3971  //
3972  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3973  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3974  //
3975  // where we substitute back in for c3 after taking the
3976  // derivative. We can compute a stencil for the cross derivative
3977  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3978  // (and gradient_entry computes the sum over the independent
3979  // variables). Below, we somewhat arbitrarily pick the last
3980  // component as the dependent one.
3981  //
3982  // Since we can now calculate derivatives of the objective
3983  // function we can use gradient descent to minimize it.
3984  //
3985  // Of course, this is much simpler in the structdim = 1 case (we
3986  // could rewrite the projection as a 1D optimization problem), but
3987  // to reduce the potential for bugs we use the same code in both
3988  // cases.
3989  const double step_size = object->diameter() / 64.0;
3990 
3991  constexpr unsigned int n_vertices_per_cell =
3993 
3994  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3995  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3996  ++vertex_n)
3997  vertices[vertex_n] = object->vertex(vertex_n);
3998 
3999  auto get_point_from_weights =
4000  [&](const Tensor<1, n_vertices_per_cell> &weights)
4001  -> Point<spacedim> {
4002  return object->get_manifold().get_new_point(
4003  make_array_view(vertices.begin(), vertices.end()),
4004  make_array_view(weights.begin_raw(), weights.end_raw()));
4005  };
4006 
4007  // pick the initial weights as (normalized) inverse distances from
4008  // the trial point:
4009  Tensor<1, n_vertices_per_cell> guess_weights;
4010  double guess_weights_sum = 0.0;
4011  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4012  ++vertex_n)
4013  {
4014  const double distance =
4015  vertices[vertex_n].distance(trial_point);
4016  if (distance == 0.0)
4017  {
4018  guess_weights = 0.0;
4019  guess_weights[vertex_n] = 1.0;
4020  guess_weights_sum = 1.0;
4021  break;
4022  }
4023  else
4024  {
4025  guess_weights[vertex_n] = 1.0 / distance;
4026  guess_weights_sum += guess_weights[vertex_n];
4027  }
4028  }
4029  guess_weights /= guess_weights_sum;
4030  Assert(internal::weights_are_ok<structdim>(guess_weights),
4031  ExcInternalError());
4032 
4033  // The optimization algorithm consists of two parts:
4034  //
4035  // 1. An outer loop where we apply the gradient descent algorithm.
4036  // 2. An inner loop where we do a line search to find the optimal
4037  // length of the step one should take in the gradient direction.
4038  //
4039  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4040  {
4041  const unsigned int dependent_direction =
4042  n_vertices_per_cell - 1;
4043  Tensor<1, n_vertices_per_cell> current_gradient;
4044  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
4045  ++row_n)
4046  {
4047  if (row_n != dependent_direction)
4048  {
4049  current_gradient[row_n] =
4050  gradient_entry<spacedim, structdim>(
4051  row_n,
4052  dependent_direction,
4053  trial_point,
4054  guess_weights,
4055  step_size,
4056  get_point_from_weights);
4057 
4058  current_gradient[dependent_direction] -=
4059  current_gradient[row_n];
4060  }
4061  }
4062 
4063  // We need to travel in the -gradient direction, as noted
4064  // above, but we may not want to take a full step in that
4065  // direction; instead, guess that we will go -0.5*gradient and
4066  // do quasi-Newton iteration to pick the best multiplier. The
4067  // goal is to find a scalar alpha such that
4068  //
4069  // F(x - alpha g)
4070  //
4071  // is minimized, where g is the gradient and F is the
4072  // objective function. To find the optimal value we find roots
4073  // of the derivative of the objective function with respect to
4074  // alpha by Newton iteration, where we approximate the first
4075  // and second derivatives of F(x - alpha g) with centered
4076  // finite differences.
4077  double gradient_weight = -0.5;
4078  auto gradient_weight_objective_function =
4079  [&](const double gradient_weight_guess) -> double {
4080  return (trial_point -
4081  get_point_from_weights(guess_weights +
4082  gradient_weight_guess *
4083  current_gradient))
4084  .norm_square();
4085  };
4086 
4087  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4088  {
4089  const double update_numerator = centered_first_difference(
4090  gradient_weight,
4091  step_size,
4092  gradient_weight_objective_function);
4093  const double update_denominator =
4094  centered_second_difference(
4095  gradient_weight,
4096  step_size,
4097  gradient_weight_objective_function);
4098 
4099  // avoid division by zero. Note that we limit the gradient
4100  // weight below
4101  if (std::abs(update_denominator) == 0.0)
4102  break;
4103  gradient_weight =
4104  gradient_weight - update_numerator / update_denominator;
4105 
4106  // Put a fairly lenient bound on the largest possible
4107  // gradient (things tend to be locally flat, so the gradient
4108  // itself is usually small)
4109  if (std::abs(gradient_weight) > 10)
4110  {
4111  gradient_weight = -10.0;
4112  break;
4113  }
4114  }
4115 
4116  // It only makes sense to take convex combinations with weights
4117  // between zero and one. If the update takes us outside of this
4118  // region then rescale the update to stay within the region and
4119  // try again
4120  Tensor<1, n_vertices_per_cell> tentative_weights =
4121  guess_weights + gradient_weight * current_gradient;
4122 
4123  double new_gradient_weight = gradient_weight;
4124  for (unsigned int iteration_count = 0; iteration_count < 40;
4125  ++iteration_count)
4126  {
4127  if (internal::weights_are_ok<structdim>(tentative_weights))
4128  break;
4129 
4130  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
4131  {
4132  if (tentative_weights[i] < 0.0)
4133  {
4134  tentative_weights -=
4135  (tentative_weights[i] / current_gradient[i]) *
4136  current_gradient;
4137  }
4138  if (tentative_weights[i] < 0.0 ||
4139  1.0 < tentative_weights[i])
4140  {
4141  new_gradient_weight /= 2.0;
4142  tentative_weights =
4143  guess_weights +
4144  new_gradient_weight * current_gradient;
4145  }
4146  }
4147  }
4148 
4149  // the update might still send us outside the valid region, so
4150  // check again and quit if the update is still not valid
4151  if (!internal::weights_are_ok<structdim>(tentative_weights))
4152  break;
4153 
4154  // if we cannot get closer by traveling in the gradient
4155  // direction then quit
4156  if (get_point_from_weights(tentative_weights)
4157  .distance(trial_point) <
4158  get_point_from_weights(guess_weights).distance(trial_point))
4159  guess_weights = tentative_weights;
4160  else
4161  break;
4162  Assert(internal::weights_are_ok<structdim>(guess_weights),
4163  ExcInternalError());
4164  }
4165  Assert(internal::weights_are_ok<structdim>(guess_weights),
4166  ExcInternalError());
4167  projected_point = get_point_from_weights(guess_weights);
4168  }
4169 
4170  // if structdim == 2 and the optimal point is not on the interior then
4171  // we may be able to get a more accurate result by projecting onto the
4172  // lines.
4173  if (structdim == 2)
4174  {
4175  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
4176  line_projections;
4177  for (unsigned int line_n = 0;
4178  line_n < GeometryInfo<structdim>::lines_per_cell;
4179  ++line_n)
4180  {
4181  line_projections[line_n] =
4182  project_to_object(object->line(line_n), trial_point);
4183  }
4184  std::sort(line_projections.begin(),
4185  line_projections.end(),
4186  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
4187  return a.distance(trial_point) <
4188  b.distance(trial_point);
4189  });
4190  if (line_projections[0].distance(trial_point) <
4191  projected_point.distance(trial_point))
4192  projected_point = line_projections[0];
4193  }
4194  }
4195  else
4196  {
4197  Assert(false, ExcNotImplemented());
4198  return projected_point;
4199  }
4200 
4201  return projected_point;
4202  }
4203 
4204 
4205 
4206  template <int dim, typename T>
4207  template <class Archive>
4208  void
4210  const unsigned int /*version*/) const
4211  {
4212  Assert(cell_ids.size() == data.size(),
4213  ExcDimensionMismatch(cell_ids.size(), data.size()));
4214  // archive the cellids in an efficient binary format
4215  const std::size_t n_cells = cell_ids.size();
4216  ar & n_cells;
4217  for (const auto &id : cell_ids)
4218  {
4219  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
4220  ar & binary_cell_id;
4221  }
4222 
4223  ar &data;
4224  }
4225 
4226 
4227 
4228  template <int dim, typename T>
4229  template <class Archive>
4230  void
4232  const unsigned int /*version*/)
4233  {
4234  std::size_t n_cells;
4235  ar & n_cells;
4236  cell_ids.clear();
4237  cell_ids.reserve(n_cells);
4238  for (unsigned int c = 0; c < n_cells; ++c)
4239  {
4240  CellId::binary_type value;
4241  ar & value;
4242  cell_ids.emplace_back(value);
4243  }
4244  ar &data;
4245  }
4246 
4247 
4248  namespace internal
4249  {
4250  template <typename DataType,
4251  typename MeshType,
4252  typename MeshCellIteratorType>
4253  inline void
4254  exchange_cell_data(
4255  const MeshType &mesh,
4256  const std::function<
4257  std_cxx17::optional<DataType>(const MeshCellIteratorType &)> &pack,
4258  const std::function<void(const MeshCellIteratorType &, const DataType &)>
4259  & unpack,
4260  const std::function<bool(const MeshCellIteratorType &)> & cell_filter,
4261  const std::function<void(
4262  const std::function<void(const MeshCellIteratorType &,
4263  const types::subdomain_id)> &)> &process_cells,
4264  const std::function<std::set<types::subdomain_id>(
4265  const parallel::TriangulationBase<MeshType::dimension,
4266  MeshType::space_dimension> &)>
4267  &compute_ghost_owners)
4268  {
4269 # ifndef DEAL_II_WITH_MPI
4270  (void)mesh;
4271  (void)pack;
4272  (void)unpack;
4273  (void)cell_filter;
4274  (void)process_cells;
4275  (void)compute_ghost_owners;
4276  Assert(false, ExcNeedsMPI());
4277 # else
4278  constexpr int dim = MeshType::dimension;
4279  constexpr int spacedim = MeshType::space_dimension;
4280  auto tria =
4281  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
4282  &mesh.get_triangulation());
4283  Assert(
4284  tria != nullptr,
4285  ExcMessage(
4286  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4287 
4288  if (const auto tria = dynamic_cast<
4290  &mesh.get_triangulation()))
4291  {
4292  Assert(
4293  tria->with_artificial_cells(),
4294  ExcMessage(
4295  "The functions GridTools::exchange_cell_data_to_ghosts() and "
4296  "GridTools::exchange_cell_data_to_level_ghosts() can only "
4297  "operate on a single layer ghost cells. However, you have "
4298  "given a Triangulation object of type "
4299  "parallel::shared::Triangulation without artificial cells "
4300  "resulting in arbitrary numbers of ghost layers."));
4301  }
4302 
4303  // build list of cells to request for each neighbor
4304  std::set<::types::subdomain_id> ghost_owners =
4305  compute_ghost_owners(*tria);
4306  std::map<::types::subdomain_id,
4307  std::vector<typename CellId::binary_type>>
4308  neighbor_cell_list;
4309 
4310  for (const auto ghost_owner : ghost_owners)
4311  neighbor_cell_list[ghost_owner] = {};
4312 
4313  process_cells([&](const auto &cell, const auto key) {
4314  if (cell_filter(cell))
4315  neighbor_cell_list[key].emplace_back(
4316  cell->id().template to_binary<spacedim>());
4317  });
4318 
4319  Assert(ghost_owners.size() == neighbor_cell_list.size(),
4320  ExcInternalError());
4321 
4322 
4323  // Before sending & receiving, make sure we protect this section with
4324  // a mutex:
4325  static Utilities::MPI::CollectiveMutex mutex;
4327  mutex, tria->get_communicator());
4328 
4329  const int mpi_tag =
4331  const int mpi_tag_reply =
4333 
4334  // send our requests:
4335  std::vector<MPI_Request> requests(ghost_owners.size());
4336  {
4337  unsigned int idx = 0;
4338  for (const auto &it : neighbor_cell_list)
4339  {
4340  // send the data about the relevant cells
4341  const int ierr = MPI_Isend(it.second.data(),
4342  it.second.size() * sizeof(it.second[0]),
4343  MPI_BYTE,
4344  it.first,
4345  mpi_tag,
4346  tria->get_communicator(),
4347  &requests[idx]);
4348  AssertThrowMPI(ierr);
4349  ++idx;
4350  }
4351  }
4352 
4353  using DestinationToBufferMap =
4354  std::map<::types::subdomain_id,
4356  DestinationToBufferMap destination_to_data_buffer_map;
4357 
4358  // receive requests and reply with the ghost indices
4359  std::vector<std::vector<typename CellId::binary_type>> cell_data_to_send(
4360  ghost_owners.size());
4361  std::vector<std::vector<types::global_dof_index>>
4362  send_dof_numbers_and_indices(ghost_owners.size());
4363  std::vector<MPI_Request> reply_requests(ghost_owners.size());
4364  std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4365 
4366  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4367  {
4368  MPI_Status status;
4369  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4370  mpi_tag,
4371  tria->get_communicator(),
4372  &status);
4373  AssertThrowMPI(ierr);
4374 
4375  int len;
4376  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4377  AssertThrowMPI(ierr);
4378  Assert(len % sizeof(cell_data_to_send[idx][0]) == 0,
4379  ExcInternalError());
4380 
4381  const unsigned int n_cells =
4382  len / sizeof(typename CellId::binary_type);
4383  cell_data_to_send[idx].resize(n_cells);
4384 
4385  ierr = MPI_Recv(cell_data_to_send[idx].data(),
4386  len,
4387  MPI_BYTE,
4388  status.MPI_SOURCE,
4389  status.MPI_TAG,
4390  tria->get_communicator(),
4391  &status);
4392  AssertThrowMPI(ierr);
4393 
4394  // store data for each cell
4395  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
4396  {
4397  const auto cell =
4398  tria->create_cell_iterator(CellId(cell_data_to_send[idx][c]));
4399 
4400  MeshCellIteratorType mesh_it(tria,
4401  cell->level(),
4402  cell->index(),
4403  &mesh);
4404  const std_cxx17::optional<DataType> data = pack(mesh_it);
4405 
4406  if (data)
4407  {
4408  typename DestinationToBufferMap::iterator p =
4409  destination_to_data_buffer_map
4410  .insert(std::make_pair(
4411  idx,
4412  GridTools::CellDataTransferBuffer<dim, DataType>()))
4413  .first;
4414 
4415  p->second.cell_ids.emplace_back(cell->id());
4416  p->second.data.emplace_back(*data);
4417  }
4418  }
4419 
4420  // send reply
4421  GridTools::CellDataTransferBuffer<dim, DataType> &data =
4422  destination_to_data_buffer_map[idx];
4423 
4424  sendbuffers[idx] =
4425  Utilities::pack(data, /*enable_compression*/ false);
4426  ierr = MPI_Isend(sendbuffers[idx].data(),
4427  sendbuffers[idx].size(),
4428  MPI_BYTE,
4429  status.MPI_SOURCE,
4430  mpi_tag_reply,
4431  tria->get_communicator(),
4432  &reply_requests[idx]);
4433  AssertThrowMPI(ierr);
4434  }
4435 
4436  // finally receive the replies
4437  std::vector<char> receive;
4438  for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4439  {
4440  MPI_Status status;
4441  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4442  mpi_tag_reply,
4443  tria->get_communicator(),
4444  &status);
4445  AssertThrowMPI(ierr);
4446 
4447  int len;
4448  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4449  AssertThrowMPI(ierr);
4450 
4451  receive.resize(len);
4452 
4453  char *ptr = receive.data();
4454  ierr = MPI_Recv(ptr,
4455  len,
4456  MPI_BYTE,
4457  status.MPI_SOURCE,
4458  status.MPI_TAG,
4459  tria->get_communicator(),
4460  &status);
4461  AssertThrowMPI(ierr);
4462 
4463  auto cellinfo =
4464  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4465  receive, /*enable_compression*/ false);
4466 
4467  DataType *data = cellinfo.data.data();
4468  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4469  {
4471  tria_cell = tria->create_cell_iterator(cellinfo.cell_ids[c]);
4472 
4473  MeshCellIteratorType cell(tria,
4474  tria_cell->level(),
4475  tria_cell->index(),
4476  &mesh);
4477 
4478  unpack(cell, *data);
4479  }
4480  }
4481 
4482  // make sure that all communication is finished
4483  // when we leave this function.
4484  if (requests.size() > 0)
4485  {
4486  const int ierr =
4487  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4488  AssertThrowMPI(ierr);
4489  }
4490  if (reply_requests.size() > 0)
4491  {
4492  const int ierr = MPI_Waitall(reply_requests.size(),
4493  reply_requests.data(),
4494  MPI_STATUSES_IGNORE);
4495  AssertThrowMPI(ierr);
4496  }
4497 
4498 
4499 # endif // DEAL_II_WITH_MPI
4500  }
4501 
4502  } // namespace internal
4503 
4504  template <typename DataType, typename MeshType>
4505  inline void
4507  const MeshType & mesh,
4508  const std::function<std_cxx17::optional<DataType>(
4509  const typename MeshType::active_cell_iterator &)> &pack,
4510  const std::function<void(const typename MeshType::active_cell_iterator &,
4511  const DataType &)> & unpack,
4512  const std::function<bool(const typename MeshType::active_cell_iterator &)>
4513  &cell_filter)
4514  {
4515 # ifndef DEAL_II_WITH_MPI
4516  (void)mesh;
4517  (void)pack;
4518  (void)unpack;
4519  (void)cell_filter;
4520  Assert(false, ExcNeedsMPI());
4521 # else
4522  internal::exchange_cell_data<DataType,
4523  MeshType,
4524  typename MeshType::active_cell_iterator>(
4525  mesh,
4526  pack,
4527  unpack,
4528  cell_filter,
4529  [&](const auto &process) {
4530  for (const auto &cell : mesh.active_cell_iterators())
4531  if (cell->is_ghost())
4532  process(cell, cell->subdomain_id());
4533  },
4534  [](const auto &tria) { return tria.ghost_owners(); });
4535 # endif
4536  }
4537 
4538 
4539 
4540  template <typename DataType, typename MeshType>
4541  inline void
4543  const MeshType & mesh,
4544  const std::function<std_cxx17::optional<DataType>(
4545  const typename MeshType::level_cell_iterator &)> &pack,
4546  const std::function<void(const typename MeshType::level_cell_iterator &,
4547  const DataType &)> & unpack,
4548  const std::function<bool(const typename MeshType::level_cell_iterator &)>
4549  &cell_filter)
4550  {
4551 # ifndef DEAL_II_WITH_MPI
4552  (void)mesh;
4553  (void)pack;
4554  (void)unpack;
4555  (void)cell_filter;
4556  Assert(false, ExcNeedsMPI());
4557 # else
4558  internal::exchange_cell_data<DataType,
4559  MeshType,
4560  typename MeshType::level_cell_iterator>(
4561  mesh,
4562  pack,
4563  unpack,
4564  cell_filter,
4565  [&](const auto &process) {
4566  for (const auto &cell : mesh.cell_iterators())
4567  if (cell->level_subdomain_id() !=
4569  !cell->is_locally_owned_on_level())
4570  process(cell, cell->level_subdomain_id());
4571  },
4572  [](const auto &tria) { return tria.level_ghost_owners(); });
4573 # endif
4574  }
4575 } // namespace GridTools
4576 
4577 # endif
4578 
4580 
4581 /*---------------------------- grid_tools.h ---------------------------*/
4582 /* end of #ifndef dealii_grid_tools_h */
4583 #endif
4584 /*---------------------------- 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:5048
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:4833
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:6325
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:3386
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:2185
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:4808
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:1921
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:81
virtual MPI_Comm get_communicator() const
Definition: tria.cc:10124
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:312
Contents is actually a matrix.
#define DEAL_II_DEPRECATED_EARLY
Definition: config.h:165
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
Definition: grid_tools.cc:5016
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:2934
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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)
Definition: grid_tools.cc:2759
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:3211
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:5110
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2042
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:3261
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:11957
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:3310
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:2665
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:6165
grid_tools.h: exchange_cell_ghosts()
Definition: mpi_tags.h:66
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4094
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
Definition: rtree.h:145
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
const double angle
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:408
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:6229
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:626
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3697
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:3908
cell_iterator end() const
Definition: tria.cc:12048
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:534
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:4273
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:6709
std::vector< std::tuple< std::pair< int, int >, unsigned int, unsigned int, Point< dim >, Point< spacedim >, unsigned int > > send_components
Definition: grid_tools.h:1100
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::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:515
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:3726
typename VectorType::value_type value_type
Definition: grid_tools.h:3305
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:1465
Signals signals
Definition: tria.h:2295
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3234
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:4300
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1256
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:697
#define DeclException0(Exception0)
Definition: exceptions.h:470
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:3016
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:6616
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
unsigned int level
Definition: grid_out.cc:4590
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:5081
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:6388
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim >> &vertices)
Definition: grid_tools.cc:288
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:4900
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:4766
cell_iterator create_cell_iterator(const CellId &cell_id) const
Definition: tria.cc:12023
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:12587
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:2696
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:5454
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1746
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)
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:445
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
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:2217
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:137
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:2568
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:378
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:4216
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:3663
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:5779
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:2679
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:4928
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:2031
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:5425
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:4246
std::vector< std::tuple< unsigned int, unsigned int, unsigned int > > recv_components
Definition: grid_tools.h:1125
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:3989
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:2506
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 invert_all_negative_measure_cells(const std::vector< Point< spacedim >> &all_vertices, std::vector< CellData< dim >> &cells)
Definition: grid_tools.cc:863
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:2022
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:4230
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:731
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:6511