Reference documentation for deal.II version 9.2.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 - 2020 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 
30 # include <deal.II/fe/mapping.h>
31 # include <deal.II/fe/mapping_q1.h>
32 
33 # include <deal.II/grid/manifold.h>
34 # include <deal.II/grid/tria.h>
37 
38 # include <deal.II/hp/dof_handler.h>
39 
41 
42 # include <deal.II/numerics/rtree.h>
43 
45 # include <boost/archive/binary_iarchive.hpp>
46 # include <boost/archive/binary_oarchive.hpp>
47 # include <boost/geometry/index/rtree.hpp>
48 # include <boost/serialization/array.hpp>
49 # include <boost/serialization/vector.hpp>
50 
51 # ifdef DEAL_II_WITH_ZLIB
52 # include <boost/iostreams/device/back_inserter.hpp>
53 # include <boost/iostreams/filter/gzip.hpp>
54 # include <boost/iostreams/filtering_stream.hpp>
55 # include <boost/iostreams/stream.hpp>
56 # endif
58 
59 # include <bitset>
60 # include <list>
61 # include <set>
62 
64 
65 // Forward declarations
66 # ifndef DOXYGEN
67 namespace parallel
68 {
69  namespace distributed
70  {
71  template <int, int>
72  class Triangulation;
73  }
74 } // namespace parallel
75 
76 namespace hp
77 {
78  template <int, int>
79  class MappingCollection;
80 }
81 
82 class SparsityPattern;
83 # endif
84 
85 namespace internal
86 {
87  template <int dim, int spacedim, class MeshType>
89  {
90  public:
91 # ifndef _MSC_VER
92  using type = typename MeshType::active_cell_iterator;
93 # else
95 # endif
96  };
97 
98 # ifdef _MSC_VER
99  template <int dim, int spacedim>
100  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim>>
101  {
102  public:
103  using type = TriaActiveIterator<
105  };
106 
107  template <int dim, int spacedim>
108  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim>>
109  {
110  public:
111  using type = TriaActiveIterator<
113  };
114 # endif
115 } // namespace internal
116 
125 namespace GridTools
126 {
127  template <int dim, int spacedim>
128  class Cache;
129 
134 
141  template <int dim, int spacedim>
142  double
144 
171  template <int dim, int spacedim>
172  double
174  const Mapping<dim, spacedim> & mapping =
176 
187  template <int dim, int spacedim>
188  double
190  const Mapping<dim, spacedim> & mapping =
192 
203  template <int dim, int spacedim>
204  double
206  const Mapping<dim, spacedim> & mapping =
208 
218  template <int dim>
219  double
220  cell_measure(
221  const std::vector<Point<dim>> &all_vertices,
223 
229  template <int dim, typename T>
230  double
231  cell_measure(const T &, ...);
232 
261  template <int dim>
264  const Triangulation<dim> &triangulation,
265  const Quadrature<dim> & quadrature);
266 
276  template <int dim>
277  double
279  const Triangulation<dim> &triangulation,
280  const Quadrature<dim> & quadrature);
281 
295  template <int dim, int spacedim>
298 
318  template <typename Iterator>
321  const Iterator & object,
323 
336  template <int dim, int spacedim>
337  std::
338  tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
340 
346 
363  template <int dim, int spacedim>
364  void
366  std::vector<CellData<dim>> & cells,
367  SubCellData & subcelldata);
368 
386  template <int dim, int spacedim>
387  void
388  delete_duplicated_vertices(std::vector<Point<spacedim>> &all_vertices,
389  std::vector<CellData<dim>> & cells,
390  SubCellData & subcelldata,
391  std::vector<unsigned int> & considered_vertices,
392  const double tol = 1e-12);
393 
399 
453  template <int dim, typename Transformation, int spacedim>
454  void
455  transform(const Transformation & transformation,
456  Triangulation<dim, spacedim> &triangulation);
457 
463  template <int dim, int spacedim>
464  void
465  shift(const Tensor<1, spacedim> & shift_vector,
466  Triangulation<dim, spacedim> &triangulation);
467 
468 
478  template <int dim>
479  void
480  rotate(const double angle, Triangulation<dim> &triangulation);
481 
494  template <int dim>
495  void
496  rotate(const double angle,
497  const unsigned int axis,
498  Triangulation<dim, 3> &triangulation);
499 
557  template <int dim>
558  void
559  laplace_transform(const std::map<unsigned int, Point<dim>> &new_points,
560  Triangulation<dim> & tria,
561  const Function<dim, double> *coefficient = nullptr,
562  const bool solve_for_absolute_positions = false);
563 
569  template <int dim, int spacedim>
570  std::map<unsigned int, Point<spacedim>>
572 
580  template <int dim, int spacedim>
581  void
582  scale(const double scaling_factor,
583  Triangulation<dim, spacedim> &triangulation);
584 
595  template <int dim, int spacedim>
596  void
597  distort_random(const double factor,
598  Triangulation<dim, spacedim> &triangulation,
599  const bool keep_boundary = true);
600 
636  template <int dim, int spacedim>
637  void
639  const bool isotropic = false,
640  const unsigned int max_iterations = 100);
641 
668  template <int dim, int spacedim>
669  void
671  const double max_ratio = 1.6180339887,
672  const unsigned int max_iterations = 5);
673 
765  template <int dim, int spacedim>
766  void
768  const double limit_angle_fraction = .75);
769 
775 
827  template <int dim, int spacedim>
828 # ifndef DOXYGEN
829  std::tuple<
830  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
831  std::vector<std::vector<Point<dim>>>,
832  std::vector<std::vector<unsigned int>>>
833 # else
834  return_type
835 # endif
837  const Cache<dim, spacedim> & cache,
838  const std::vector<Point<spacedim>> &points,
840  &cell_hint =
842 
871  template <int dim, int spacedim>
872 # ifndef DOXYGEN
873  std::tuple<
874  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
875  std::vector<std::vector<Point<dim>>>,
876  std::vector<std::vector<unsigned int>>,
877  std::vector<unsigned int>>
878 # else
879  return_type
880 # endif
882  const Cache<dim, spacedim> & cache,
883  const std::vector<Point<spacedim>> &points,
885  &cell_hint =
887 
952  template <int dim, int spacedim>
953 # ifndef DOXYGEN
954  std::tuple<
955  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
956  std::vector<std::vector<Point<dim>>>,
957  std::vector<std::vector<unsigned int>>,
958  std::vector<std::vector<Point<spacedim>>>,
959  std::vector<std::vector<unsigned int>>>
960 # else
961  return_type
962 # endif
964  const GridTools::Cache<dim, spacedim> & cache,
965  const std::vector<Point<spacedim>> & local_points,
966  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes);
967 
1006  template <int dim, int spacedim>
1007  std::map<unsigned int, Point<spacedim>>
1009  const Mapping<dim, spacedim> & mapping =
1011 
1023  template <int spacedim>
1024  unsigned int
1025  find_closest_vertex(const std::map<unsigned int, Point<spacedim>> &vertices,
1026  const Point<spacedim> & p);
1027 
1054  template <int dim, template <int, int> class MeshType, int spacedim>
1055  unsigned int
1056  find_closest_vertex(const MeshType<dim, spacedim> &mesh,
1057  const Point<spacedim> & p,
1058  const std::vector<bool> & marked_vertices = {});
1059 
1085  template <int dim, template <int, int> class MeshType, int spacedim>
1086  unsigned int
1088  const MeshType<dim, spacedim> &mesh,
1089  const Point<spacedim> & p,
1090  const std::vector<bool> & marked_vertices = {});
1091 
1092 
1117  template <int dim, template <int, int> class MeshType, int spacedim>
1118 # ifndef _MSC_VER
1119  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1120 # else
1121  std::vector<
1122  typename ::internal::
1123  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1124 # endif
1125  find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &container,
1126  const unsigned int vertex_index);
1127 
1128 
1153  template <int dim, template <int, int> class MeshType, int spacedim>
1154 # ifndef _MSC_VER
1155  typename MeshType<dim, spacedim>::active_cell_iterator
1156 # else
1157  typename ::internal::ActiveCellIterator<dim,
1158  spacedim,
1159  MeshType<dim, spacedim>>::type
1160 # endif
1161  find_active_cell_around_point(const MeshType<dim, spacedim> &mesh,
1162  const Point<spacedim> & p,
1163  const std::vector<bool> &marked_vertices = {});
1164 
1252  template <int dim, template <int, int> class MeshType, int spacedim>
1253 # ifndef _MSC_VER
1254  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1255 # else
1256  std::pair<typename ::internal::
1257  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1258  Point<dim>>
1259 # endif
1261  const MeshType<dim, spacedim> &mesh,
1262  const Point<spacedim> & p,
1263  const std::vector<bool> &marked_vertices = {});
1264 
1276  template <int dim, template <int, int> class MeshType, int spacedim>
1277 # ifndef _MSC_VER
1278  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim>>
1279 # else
1280  std::pair<typename ::internal::
1281  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1282  Point<dim>>
1283 # endif
1285  const Mapping<dim, spacedim> & mapping,
1286  const MeshType<dim, spacedim> &mesh,
1287  const Point<spacedim> & p,
1288  const std::vector<
1289  std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
1291  const std::vector<std::vector<Tensor<1, spacedim>>> &vertex_to_cell_centers,
1292  const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1293  typename MeshType<dim, spacedim>::active_cell_iterator(),
1294  const std::vector<bool> & marked_vertices = {},
1295  const RTree<std::pair<Point<spacedim>, unsigned int>> &used_vertices_rtree =
1296  RTree<std::pair<Point<spacedim>, unsigned int>>{});
1297 
1319  template <int dim, int spacedim>
1320  std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator,
1321  Point<dim>>
1323  const hp::MappingCollection<dim, spacedim> &mapping,
1324  const hp::DoFHandler<dim, spacedim> & mesh,
1325  const Point<spacedim> & p);
1326 
1333  template <int dim, int spacedim>
1334  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1335  Point<dim>>
1337  const Cache<dim, spacedim> &cache,
1338  const Point<spacedim> & p,
1341  const std::vector<bool> &marked_vertices = {});
1342 
1360  template <int dim, template <int, int> class MeshType, int spacedim>
1361 # ifndef _MSC_VER
1362  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1363  Point<dim>>>
1364 # else
1365  std::vector<std::pair<
1366  typename ::internal::
1367  ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1368  Point<dim>>>
1369 # endif
1371  const Mapping<dim, spacedim> & mapping,
1372  const MeshType<dim, spacedim> &mesh,
1373  const Point<spacedim> & p,
1374  const double tolerance = 1e-10,
1375  const std::vector<bool> & marked_vertices = {});
1376 
1398  template <class MeshType>
1399  std::vector<typename MeshType::active_cell_iterator>
1400  get_active_child_cells(const typename MeshType::cell_iterator &cell);
1401 
1426  template <class MeshType>
1427  void
1429  const typename MeshType::active_cell_iterator & cell,
1430  std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1431 
1484  template <class MeshType>
1485  std::vector<typename MeshType::active_cell_iterator>
1487  const MeshType &mesh,
1488  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1489  &predicate);
1490 
1491 
1499  template <class MeshType>
1500  std::vector<typename MeshType::cell_iterator>
1502  const MeshType &mesh,
1503  const std::function<bool(const typename MeshType::cell_iterator &)>
1504  & predicate,
1505  const unsigned int level);
1506 
1507 
1523  template <class MeshType>
1524  std::vector<typename MeshType::active_cell_iterator>
1525  compute_ghost_cell_halo_layer(const MeshType &mesh);
1526 
1578  template <class MeshType>
1579  std::vector<typename MeshType::active_cell_iterator>
1581  const MeshType &mesh,
1582  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1583  & predicate,
1584  const double layer_thickness);
1585 
1611  template <class MeshType>
1612  std::vector<typename MeshType::active_cell_iterator>
1613  compute_ghost_cell_layer_within_distance(const MeshType &mesh,
1614  const double layer_thickness);
1615 
1631  template <class MeshType>
1632  std::pair<Point<MeshType::space_dimension>, Point<MeshType::space_dimension>>
1634  const MeshType &mesh,
1635  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1636  &predicate);
1637 
1690  template <class MeshType>
1691  std::vector<BoundingBox<MeshType::space_dimension>>
1693  const MeshType &mesh,
1694  const std::function<bool(const typename MeshType::active_cell_iterator &)>
1695  & predicate,
1696  const unsigned int refinement_level = 0,
1697  const bool allow_merge = false,
1698  const unsigned int max_boxes = numbers::invalid_unsigned_int);
1699 
1729  template <int spacedim>
1730 # ifndef DOXYGEN
1731  std::tuple<std::vector<std::vector<unsigned int>>,
1732  std::map<unsigned int, unsigned int>,
1733  std::map<unsigned int, std::vector<unsigned int>>>
1734 # else
1735  return_type
1736 # endif
1738  const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
1739  const std::vector<Point<spacedim>> & points);
1740 
1741 
1776  template <int spacedim>
1777 # ifndef DOXYGEN
1778  std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1779  std::map<unsigned int, unsigned int>,
1780  std::map<unsigned int, std::vector<unsigned int>>>
1781 # else
1782  return_type
1783 # endif
1785  const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
1786  const std::vector<Point<spacedim>> & points);
1787 
1788 
1797  template <int dim, int spacedim>
1798  std::vector<
1799  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1800  vertex_to_cell_map(const Triangulation<dim, spacedim> &triangulation);
1801 
1816  template <int dim, int spacedim>
1817  std::vector<std::vector<Tensor<1, spacedim>>>
1819  const Triangulation<dim, spacedim> &mesh,
1820  const std::vector<
1822  &vertex_to_cells);
1823 
1824 
1834  template <int dim, int spacedim>
1835  unsigned int
1838  const Point<spacedim> & position,
1839  const Mapping<dim, spacedim> & mapping =
1841 
1853  template <int dim, int spacedim>
1854  std::map<unsigned int, types::global_vertex_index>
1857 
1871  template <int dim, int spacedim>
1872  std::pair<unsigned int, double>
1875 
1881 
1890  template <int dim, int spacedim>
1891  void
1893  const Triangulation<dim, spacedim> &triangulation,
1894  DynamicSparsityPattern & connectivity);
1895 
1904  template <int dim, int spacedim>
1905  void
1907  const Triangulation<dim, spacedim> &triangulation,
1908  DynamicSparsityPattern & connectivity);
1909 
1918  template <int dim, int spacedim>
1919  void
1921  const Triangulation<dim, spacedim> &triangulation,
1922  const unsigned int level,
1923  DynamicSparsityPattern & connectivity);
1924 
1945  template <int dim, int spacedim>
1946  void
1947  partition_triangulation(const unsigned int n_partitions,
1948  Triangulation<dim, spacedim> & triangulation,
1949  const SparsityTools::Partitioner partitioner =
1951 
1962  template <int dim, int spacedim>
1963  void
1964  partition_triangulation(const unsigned int n_partitions,
1965  const std::vector<unsigned int> &cell_weights,
1966  Triangulation<dim, spacedim> & triangulation,
1967  const SparsityTools::Partitioner partitioner =
1969 
2015  template <int dim, int spacedim>
2016  void
2017  partition_triangulation(const unsigned int n_partitions,
2018  const SparsityPattern & cell_connection_graph,
2019  Triangulation<dim, spacedim> &triangulation,
2020  const SparsityTools::Partitioner partitioner =
2022 
2033  template <int dim, int spacedim>
2034  void
2035  partition_triangulation(const unsigned int n_partitions,
2036  const std::vector<unsigned int> &cell_weights,
2037  const SparsityPattern & cell_connection_graph,
2038  Triangulation<dim, spacedim> &triangulation,
2039  const SparsityTools::Partitioner partitioner =
2041 
2056  template <int dim, int spacedim>
2057  void
2058  partition_triangulation_zorder(const unsigned int n_partitions,
2059  Triangulation<dim, spacedim> &triangulation,
2060  const bool group_siblings = true);
2061 
2073  template <int dim, int spacedim>
2074  void
2076 
2087  template <int dim, int spacedim>
2088  void
2090  std::vector<types::subdomain_id> & subdomain);
2091 
2106  template <int dim, int spacedim>
2107  unsigned int
2109  const Triangulation<dim, spacedim> &triangulation,
2110  const types::subdomain_id subdomain);
2111 
2112 
2142  template <int dim, int spacedim>
2143  std::vector<bool>
2145 
2151 
2185  template <typename MeshType>
2186  std::list<std::pair<typename MeshType::cell_iterator,
2187  typename MeshType::cell_iterator>>
2188  get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2);
2189 
2199  template <int dim, int spacedim>
2200  bool
2202  const Triangulation<dim, spacedim> &mesh_2);
2203 
2213  template <typename MeshType>
2214  bool
2215  have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2);
2216 
2222 
2238  template <int dim, int spacedim>
2242  & distorted_cells,
2243  Triangulation<dim, spacedim> &triangulation);
2244 
2245 
2246 
2255 
2256 
2296  template <class MeshType>
2297  std::vector<typename MeshType::active_cell_iterator>
2298  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
2299 
2300 
2324  template <class Container>
2325  std::vector<typename Container::cell_iterator>
2327  const std::vector<typename Container::active_cell_iterator> &patch_cells);
2328 
2397  template <class Container>
2398  void
2400  const std::vector<typename Container::active_cell_iterator> &patch,
2402  &local_triangulation,
2403  std::map<
2404  typename Triangulation<Container::dimension,
2405  Container::space_dimension>::active_cell_iterator,
2406  typename Container::active_cell_iterator> &patch_to_global_tria_map);
2407 
2444  template <class DoFHandlerType>
2445  std::map<types::global_dof_index,
2446  std::vector<typename DoFHandlerType::active_cell_iterator>>
2447  get_dof_to_support_patch_map(DoFHandlerType &dof_handler);
2448 
2449 
2456 
2462  template <typename CellIterator>
2464  {
2468  CellIterator cell[2];
2469 
2474  unsigned int face_idx[2];
2475 
2481  std::bitset<3> orientation;
2482 
2496  };
2497 
2498 
2564  template <typename FaceIterator>
2565  bool orthogonal_equality(
2566  std::bitset<3> & orientation,
2567  const FaceIterator & face1,
2568  const FaceIterator & face2,
2569  const int direction,
2573 
2574 
2578  template <typename FaceIterator>
2579  bool
2581  const FaceIterator & face1,
2582  const FaceIterator & face2,
2583  const int direction,
2587 
2588 
2647  template <typename MeshType>
2648  void
2650  const MeshType & mesh,
2651  const types::boundary_id b_id1,
2652  const types::boundary_id b_id2,
2653  const int direction,
2655  & matched_pairs,
2656  const Tensor<1, MeshType::space_dimension> &offset =
2659 
2660 
2685  template <typename MeshType>
2686  void
2688  const MeshType & mesh,
2689  const types::boundary_id b_id,
2690  const int direction,
2692  & matched_pairs,
2693  const ::Tensor<1, MeshType::space_dimension> &offset =
2696 
2702 
2725  template <int dim, int spacedim>
2726  void
2728  const bool reset_boundary_ids = false);
2729 
2753  template <int dim, int spacedim>
2754  void
2756  const std::vector<types::boundary_id> &src_boundary_ids,
2757  const std::vector<types::manifold_id> &dst_manifold_ids,
2759  const std::vector<types::boundary_id> &reset_boundary_ids = {});
2760 
2792  template <int dim, int spacedim>
2793  void
2795  const bool compute_face_ids = false);
2796 
2823  template <int dim, int spacedim>
2824  void
2827  const std::function<types::manifold_id(
2828  const std::set<types::manifold_id> &)> &disambiguation_function =
2829  [](const std::set<types::manifold_id> &manifold_ids) {
2830  if (manifold_ids.size() == 1)
2831  return *manifold_ids.begin();
2832  else
2834  },
2835  bool overwrite_only_flat_manifold_ids = true);
2920  template <typename DataType, typename MeshType>
2921  void
2923  const MeshType & mesh,
2924  const std::function<std_cxx17::optional<DataType>(
2925  const typename MeshType::active_cell_iterator &)> &pack,
2926  const std::function<void(const typename MeshType::active_cell_iterator &,
2927  const DataType &)> & unpack);
2928 
2929  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
2930  * boxes @p local_bboxes.
2931  *
2932  * This function is meant to exchange bounding boxes describing the locally
2933  * owned cells in a distributed triangulation obtained with the function
2934  * GridTools::compute_mesh_predicate_bounding_box .
2935  *
2936  * The output vector's size is the number of processes of the MPI
2937  * communicator:
2938  * its i-th entry contains the vector @p local_bboxes of the i-th process.
2939  */
2940  template <int spacedim>
2941  std::vector<std::vector<BoundingBox<spacedim>>>
2943  const std::vector<BoundingBox<spacedim>> &local_bboxes,
2944  MPI_Comm mpi_communicator);
2945 
2980  template <int spacedim>
2983  const std::vector<BoundingBox<spacedim>> &local_description,
2984  MPI_Comm mpi_communicator);
2985 
3005  template <int dim, int spacedim>
3006  void
3008  const Triangulation<dim, spacedim> & tria,
3009  std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
3010  std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
3011 
3018  template <int dim, int spacedim>
3019  std::map<unsigned int, std::set<::types::subdomain_id>>
3021  const Triangulation<dim, spacedim> &tria);
3022 
3035  template <int dim, typename T>
3037  {
3041  std::vector<CellId> cell_ids;
3042 
3046  std::vector<T> data;
3047 
3055  template <class Archive>
3056  void
3057  save(Archive &ar, const unsigned int version) const;
3058 
3063  template <class Archive>
3064  void
3065  load(Archive &ar, const unsigned int version);
3066 
3067 # ifdef DOXYGEN
3068 
3072  template <class Archive>
3073  void
3074  serialize(Archive &archive, const unsigned int version);
3075 # else
3076  // This macro defines the serialize() method that is compatible with
3077  // the templated save() and load() method that have been implemented.
3078  BOOST_SERIALIZATION_SPLIT_MEMBER()
3079 # endif
3080  };
3081 
3086 
3091  int,
3092  << "The number of partitions you gave is " << arg1
3093  << ", but must be greater than zero.");
3098  int,
3099  << "The subdomain id " << arg1
3100  << " has no cells associated with it.");
3105 
3110  double,
3111  << "The scaling factor must be positive, but it is " << arg1
3112  << ".");
3116  template <int N>
3118  Point<N>,
3119  << "The point <" << arg1
3120  << "> could not be found inside any of the "
3121  << "coarse grid cells.");
3125  template <int N>
3127  Point<N>,
3128  << "The point <" << arg1
3129  << "> could not be found inside any of the "
3130  << "subcells of a coarse grid cell.");
3131 
3136  unsigned int,
3137  << "The given vertex with index " << arg1
3138  << " is not used in the given triangulation.");
3139 
3140 
3143 } /*namespace GridTools*/
3144 
3145 
3146 
3147 /* ----------------- Template function --------------- */
3148 
3149 # ifndef DOXYGEN
3150 
3151 namespace GridTools
3152 {
3153  // declare specializations
3154  template <>
3155  double
3156  cell_measure<1>(const std::vector<Point<1>> &,
3157  const unsigned int (&)[GeometryInfo<1>::vertices_per_cell]);
3158 
3159  template <>
3160  double
3161  cell_measure<2>(const std::vector<Point<2>> &,
3162  const unsigned int (&)[GeometryInfo<2>::vertices_per_cell]);
3163 
3164  template <>
3165  double
3166  cell_measure<3>(const std::vector<Point<3>> &,
3167  const unsigned int (&)[GeometryInfo<3>::vertices_per_cell]);
3168 
3169 
3170 
3171  template <int dim, typename T>
3172  double
3173  cell_measure(const T &, ...)
3174  {
3175  Assert(false, ExcNotImplemented());
3176  return std::numeric_limits<double>::quiet_NaN();
3177  }
3178 
3179 
3180 
3181  // This specialization is defined here so that the general template in the
3182  // source file doesn't need to have further 1D overloads for the internal
3183  // functions it calls.
3184  template <>
3188  {
3189  return {};
3190  }
3191 
3192 
3193 
3194  template <int dim, typename Predicate, int spacedim>
3195  void
3196  transform(const Predicate & predicate,
3197  Triangulation<dim, spacedim> &triangulation)
3198  {
3199  std::vector<bool> treated_vertices(triangulation.n_vertices(), false);
3200 
3201  // loop over all active cells, and
3202  // transform those vertices that
3203  // have not yet been touched. note
3204  // that we get to all vertices in
3205  // the triangulation by only
3206  // visiting the active cells.
3208  cell = triangulation.begin_active(),
3209  endc = triangulation.end();
3210  for (; cell != endc; ++cell)
3211  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
3212  if (treated_vertices[cell->vertex_index(v)] == false)
3213  {
3214  // transform this vertex
3215  cell->vertex(v) = predicate(cell->vertex(v));
3216  // and mark it as treated
3217  treated_vertices[cell->vertex_index(v)] = true;
3218  };
3219 
3220 
3221  // now fix any vertices on hanging nodes so that we don't create any holes
3222  if (dim == 2)
3223  {
3225  cell = triangulation.begin_active(),
3226  endc = triangulation.end();
3227  for (; cell != endc; ++cell)
3228  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3229  if (cell->face(face)->has_children() &&
3230  !cell->face(face)->at_boundary())
3231  {
3232  // this line has children
3233  cell->face(face)->child(0)->vertex(1) =
3234  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3235  2;
3236  }
3237  }
3238  else if (dim == 3)
3239  {
3241  cell = triangulation.begin_active(),
3242  endc = triangulation.end();
3243  for (; cell != endc; ++cell)
3244  for (const unsigned int face : GeometryInfo<dim>::face_indices())
3245  if (cell->face(face)->has_children() &&
3246  !cell->face(face)->at_boundary())
3247  {
3248  // this face has hanging nodes
3249  cell->face(face)->child(0)->vertex(1) =
3250  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3251  2.0;
3252  cell->face(face)->child(0)->vertex(2) =
3253  (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3254  2.0;
3255  cell->face(face)->child(1)->vertex(3) =
3256  (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3257  2.0;
3258  cell->face(face)->child(2)->vertex(3) =
3259  (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3260  2.0;
3261 
3262  // center of the face
3263  cell->face(face)->child(0)->vertex(3) =
3264  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3265  cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3266  4.0;
3267  }
3268  }
3269 
3270  // Make sure FEValues notices that the mesh has changed
3271  triangulation.signals.mesh_movement();
3272  }
3273 
3274 
3275 
3276  template <class MeshType>
3277  std::vector<typename MeshType::active_cell_iterator>
3278  get_active_child_cells(const typename MeshType::cell_iterator &cell)
3279  {
3280  std::vector<typename MeshType::active_cell_iterator> child_cells;
3281 
3282  if (cell->has_children())
3283  {
3284  for (unsigned int child = 0; child < cell->n_children(); ++child)
3285  if (cell->child(child)->has_children())
3286  {
3287  const std::vector<typename MeshType::active_cell_iterator>
3288  children = get_active_child_cells<MeshType>(cell->child(child));
3289  child_cells.insert(child_cells.end(),
3290  children.begin(),
3291  children.end());
3292  }
3293  else
3294  child_cells.push_back(cell->child(child));
3295  }
3296 
3297  return child_cells;
3298  }
3299 
3300 
3301 
3302  template <class MeshType>
3303  void
3305  const typename MeshType::active_cell_iterator & cell,
3306  std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3307  {
3308  active_neighbors.clear();
3309  for (const unsigned int n :
3311  if (!cell->at_boundary(n))
3312  {
3313  if (MeshType::dimension == 1)
3314  {
3315  // check children of neighbor. note
3316  // that in 1d children of the neighbor
3317  // may be further refined. In 1d the
3318  // case is simple since we know what
3319  // children bound to the present cell
3320  typename MeshType::cell_iterator neighbor_child =
3321  cell->neighbor(n);
3322  if (!neighbor_child->is_active())
3323  {
3324  while (neighbor_child->has_children())
3325  neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3326 
3327  Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3328  ExcInternalError());
3329  }
3330  active_neighbors.push_back(neighbor_child);
3331  }
3332  else
3333  {
3334  if (cell->face(n)->has_children())
3335  // this neighbor has children. find
3336  // out which border to the present
3337  // cell
3338  for (unsigned int c = 0;
3339  c < cell->face(n)->number_of_children();
3340  ++c)
3341  active_neighbors.push_back(
3342  cell->neighbor_child_on_subface(n, c));
3343  else
3344  {
3345  // the neighbor must be active
3346  // himself
3347  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
3348  active_neighbors.push_back(cell->neighbor(n));
3349  }
3350  }
3351  }
3352  }
3353 
3354 
3355 
3356  namespace internal
3357  {
3358  namespace ProjectToObject
3359  {
3372  struct CrossDerivative
3373  {
3374  const unsigned int direction_0;
3375  const unsigned int direction_1;
3376 
3377  CrossDerivative(const unsigned int d0, const unsigned int d1);
3378  };
3379 
3380  inline CrossDerivative::CrossDerivative(const unsigned int d0,
3381  const unsigned int d1)
3382  : direction_0(d0)
3383  , direction_1(d1)
3384  {}
3385 
3386 
3387 
3392  template <typename F>
3393  inline auto
3394  centered_first_difference(const double center,
3395  const double step,
3396  const F &f) -> decltype(f(center) - f(center))
3397  {
3398  return (f(center + step) - f(center - step)) / (2.0 * step);
3399  }
3400 
3401 
3402 
3407  template <typename F>
3408  inline auto
3409  centered_second_difference(const double center,
3410  const double step,
3411  const F &f) -> decltype(f(center) - f(center))
3412  {
3413  return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3414  (step * step);
3415  }
3416 
3417 
3418 
3428  template <int structdim, typename F>
3429  inline auto
3430  cross_stencil(
3431  const CrossDerivative cross_derivative,
3433  const double step,
3434  const F &f) -> decltype(f(center) - f(center))
3435  {
3437  simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3438  simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3439  return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3440  1.0 / 3.0 * f(center - simplex_vector) +
3441  16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3442  step;
3443  }
3444 
3445 
3446 
3453  template <int spacedim, int structdim, typename F>
3454  inline double
3455  gradient_entry(
3456  const unsigned int row_n,
3457  const unsigned int dependent_direction,
3458  const Point<spacedim> &p0,
3460  const double step,
3461  const F & f)
3462  {
3464  dependent_direction <
3466  ExcMessage("This function assumes that the last weight is a "
3467  "dependent variable (and hence we cannot take its "
3468  "derivative directly)."));
3469  Assert(row_n != dependent_direction,
3470  ExcMessage(
3471  "We cannot differentiate with respect to the variable "
3472  "that is assumed to be dependent."));
3473 
3474  const Point<spacedim> manifold_point = f(center);
3475  const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
3476  {row_n, dependent_direction}, center, step, f);
3477  double entry = 0.0;
3478  for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3479  entry +=
3480  -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3481  return entry;
3482  }
3483 
3489  template <typename Iterator, int spacedim, int structdim>
3491  project_to_d_linear_object(const Iterator & object,
3492  const Point<spacedim> &trial_point)
3493  {
3494  // let's look at this for simplicity for a quad (structdim==2) in a
3495  // space with spacedim>2 (notate trial_point by y): all points on the
3496  // surface are given by
3497  // x(\xi) = sum_i v_i phi_x(\xi)
3498  // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
3499  // reference coordinates of the quad. so what we are trying to do is
3500  // find a point x on the surface that is closest to the point y. there
3501  // are different ways to solve this problem, but in the end it's a
3502  // nonlinear problem and we have to find reference coordinates \xi so
3503  // that J(\xi) = 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function
3504  // that is structdim-linear in \xi, so J(\xi) is a polynomial of degree
3505  // 2*structdim that we'd like to minimize. unless structdim==1, we'll
3506  // have to use a Newton method to find the answer. This leads to the
3507  // following formulation of Newton steps:
3508  //
3509  // Given \xi_k, find \delta\xi_k so that
3510  // H_k \delta\xi_k = - F_k
3511  // where H_k is an approximation to the second derivatives of J at
3512  // \xi_k, and F_k is the first derivative of J. We'll iterate this a
3513  // number of times until the right hand side is small enough. As a
3514  // stopping criterion, we terminate if ||\delta\xi||<eps.
3515  //
3516  // As for the Hessian, the best choice would be
3517  // H_k = J''(\xi_k)
3518  // but we'll opt for the simpler Gauss-Newton form
3519  // H_k = A^T A
3520  // i.e.
3521  // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
3522  // \partial_n phi_i *
3523  // \partial_m phi_j
3524  // we start at xi=(0.5, 0.5).
3525  Point<structdim> xi;
3526  for (unsigned int d = 0; d < structdim; ++d)
3527  xi[d] = 0.5;
3528 
3529  Point<spacedim> x_k;
3530  for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
3531  x_k += object->vertex(i) *
3533 
3534  do
3535  {
3537  for (const unsigned int i :
3539  F_k +=
3540  (x_k - trial_point) * object->vertex(i) *
3542  i);
3543 
3545  for (const unsigned int i :
3547  for (const unsigned int j :
3549  {
3552  xi, i),
3554  xi, j));
3555  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
3556  }
3557 
3558  const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
3559  xi += delta_xi;
3560 
3561  x_k = Point<spacedim>();
3562  for (const unsigned int i :
3564  x_k += object->vertex(i) *
3566 
3567  if (delta_xi.norm() < 1e-7)
3568  break;
3569  }
3570  while (true);
3571 
3572  return x_k;
3573  }
3574  } // namespace ProjectToObject
3575  } // namespace internal
3576 
3577 
3578 
3579  namespace internal
3580  {
3581  // We hit an internal compiler error in ICC 15 if we define this as a lambda
3582  // inside the project_to_object function below.
3583  template <int structdim>
3584  inline bool
3585  weights_are_ok(
3587  {
3588  // clang has trouble figuring out structdim here, so define it
3589  // again:
3590  static const std::size_t n_vertices_per_cell =
3592  n_independent_components;
3593  std::array<double, n_vertices_per_cell> copied_weights;
3594  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3595  {
3596  copied_weights[i] = v[i];
3597  if (v[i] < 0.0 || v[i] > 1.0)
3598  return false;
3599  }
3600 
3601  // check the sum: try to avoid some roundoff errors by summing in order
3602  std::sort(copied_weights.begin(), copied_weights.end());
3603  const double sum =
3604  std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3605  return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
3606  }
3607  } // namespace internal
3608 
3609  template <typename Iterator>
3612  const Iterator & object,
3614  {
3615  const int spacedim = Iterator::AccessorType::space_dimension;
3616  const int structdim = Iterator::AccessorType::structure_dimension;
3617 
3618  Point<spacedim> projected_point = trial_point;
3619 
3620  if (structdim >= spacedim)
3621  return projected_point;
3622  else if (structdim == 1 || structdim == 2)
3623  {
3624  using namespace internal::ProjectToObject;
3625  // Try to use the special flat algorithm for quads (this is better
3626  // than the general algorithm in 3D). This does not take into account
3627  // whether projected_point is outside the quad, but we optimize along
3628  // lines below anyway:
3629  const int dim = Iterator::AccessorType::dimension;
3630  const Manifold<dim, spacedim> &manifold = object->get_manifold();
3631  if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
3632  &manifold) != nullptr)
3633  {
3634  projected_point =
3635  project_to_d_linear_object<Iterator, spacedim, structdim>(
3636  object, trial_point);
3637  }
3638  else
3639  {
3640  // We want to find a point on the convex hull (defined by the
3641  // vertices of the object and the manifold description) that is
3642  // relatively close to the trial point. This has a few issues:
3643  //
3644  // 1. For a general convex hull we are not guaranteed that a unique
3645  // minimum exists.
3646  // 2. The independent variables in the optimization process are the
3647  // weights given to Manifold::get_new_point, which must sum to 1,
3648  // so we cannot use standard finite differences to approximate a
3649  // gradient.
3650  //
3651  // There is not much we can do about 1., but for 2. we can derive
3652  // finite difference stencils that work on a structdim-dimensional
3653  // simplex and rewrite the optimization problem to use those
3654  // instead. Consider the structdim 2 case and let
3655  //
3656  // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
3657  // c2, c3})
3658  //
3659  // where {c0, c1, c2, c3} are the weights for the four vertices on
3660  // the quadrilateral. We seek to minimize the Euclidean distance
3661  // between F(...) and trial_point. We can solve for c3 in terms of
3662  // the other weights and get, for one coordinate direction
3663  //
3664  // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
3665  // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
3666  //
3667  // where we substitute back in for c3 after taking the
3668  // derivative. We can compute a stencil for the cross derivative
3669  // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
3670  // (and gradient_entry computes the sum over the independent
3671  // variables). Below, we somewhat arbitrarily pick the last
3672  // component as the dependent one.
3673  //
3674  // Since we can now calculate derivatives of the objective
3675  // function we can use gradient descent to minimize it.
3676  //
3677  // Of course, this is much simpler in the structdim = 1 case (we
3678  // could rewrite the projection as a 1D optimization problem), but
3679  // to reduce the potential for bugs we use the same code in both
3680  // cases.
3681  const double step_size = object->diameter() / 64.0;
3682 
3683  constexpr unsigned int n_vertices_per_cell =
3685 
3686  std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3687  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3688  ++vertex_n)
3689  vertices[vertex_n] = object->vertex(vertex_n);
3690 
3691  auto get_point_from_weights =
3692  [&](const Tensor<1, n_vertices_per_cell> &weights)
3693  -> Point<spacedim> {
3694  return object->get_manifold().get_new_point(
3695  make_array_view(vertices.begin(), vertices.end()),
3696  make_array_view(weights.begin_raw(), weights.end_raw()));
3697  };
3698 
3699  // pick the initial weights as (normalized) inverse distances from
3700  // the trial point:
3701  Tensor<1, n_vertices_per_cell> guess_weights;
3702  double guess_weights_sum = 0.0;
3703  for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3704  ++vertex_n)
3705  {
3706  const double distance =
3707  vertices[vertex_n].distance(trial_point);
3708  if (distance == 0.0)
3709  {
3710  guess_weights = 0.0;
3711  guess_weights[vertex_n] = 1.0;
3712  guess_weights_sum = 1.0;
3713  break;
3714  }
3715  else
3716  {
3717  guess_weights[vertex_n] = 1.0 / distance;
3718  guess_weights_sum += guess_weights[vertex_n];
3719  }
3720  }
3721  guess_weights /= guess_weights_sum;
3722  Assert(internal::weights_are_ok<structdim>(guess_weights),
3723  ExcInternalError());
3724 
3725  // The optimization algorithm consists of two parts:
3726  //
3727  // 1. An outer loop where we apply the gradient descent algorithm.
3728  // 2. An inner loop where we do a line search to find the optimal
3729  // length of the step one should take in the gradient direction.
3730  //
3731  for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3732  {
3733  const unsigned int dependent_direction =
3734  n_vertices_per_cell - 1;
3735  Tensor<1, n_vertices_per_cell> current_gradient;
3736  for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
3737  ++row_n)
3738  {
3739  if (row_n != dependent_direction)
3740  {
3741  current_gradient[row_n] =
3742  gradient_entry<spacedim, structdim>(
3743  row_n,
3744  dependent_direction,
3745  trial_point,
3746  guess_weights,
3747  step_size,
3748  get_point_from_weights);
3749 
3750  current_gradient[dependent_direction] -=
3751  current_gradient[row_n];
3752  }
3753  }
3754 
3755  // We need to travel in the -gradient direction, as noted
3756  // above, but we may not want to take a full step in that
3757  // direction; instead, guess that we will go -0.5*gradient and
3758  // do quasi-Newton iteration to pick the best multiplier. The
3759  // goal is to find a scalar alpha such that
3760  //
3761  // F(x - alpha g)
3762  //
3763  // is minimized, where g is the gradient and F is the
3764  // objective function. To find the optimal value we find roots
3765  // of the derivative of the objective function with respect to
3766  // alpha by Newton iteration, where we approximate the first
3767  // and second derivatives of F(x - alpha g) with centered
3768  // finite differences.
3769  double gradient_weight = -0.5;
3770  auto gradient_weight_objective_function =
3771  [&](const double gradient_weight_guess) -> double {
3772  return (trial_point -
3773  get_point_from_weights(guess_weights +
3774  gradient_weight_guess *
3775  current_gradient))
3776  .norm_square();
3777  };
3778 
3779  for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3780  {
3781  const double update_numerator = centered_first_difference(
3782  gradient_weight,
3783  step_size,
3784  gradient_weight_objective_function);
3785  const double update_denominator =
3786  centered_second_difference(
3787  gradient_weight,
3788  step_size,
3789  gradient_weight_objective_function);
3790 
3791  // avoid division by zero. Note that we limit the gradient
3792  // weight below
3793  if (std::abs(update_denominator) == 0.0)
3794  break;
3795  gradient_weight =
3796  gradient_weight - update_numerator / update_denominator;
3797 
3798  // Put a fairly lenient bound on the largest possible
3799  // gradient (things tend to be locally flat, so the gradient
3800  // itself is usually small)
3801  if (std::abs(gradient_weight) > 10)
3802  {
3803  gradient_weight = -10.0;
3804  break;
3805  }
3806  }
3807 
3808  // It only makes sense to take convex combinations with weights
3809  // between zero and one. If the update takes us outside of this
3810  // region then rescale the update to stay within the region and
3811  // try again
3812  Tensor<1, n_vertices_per_cell> tentative_weights =
3813  guess_weights + gradient_weight * current_gradient;
3814 
3815  double new_gradient_weight = gradient_weight;
3816  for (unsigned int iteration_count = 0; iteration_count < 40;
3817  ++iteration_count)
3818  {
3819  if (internal::weights_are_ok<structdim>(tentative_weights))
3820  break;
3821 
3822  for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
3823  {
3824  if (tentative_weights[i] < 0.0)
3825  {
3826  tentative_weights -=
3827  (tentative_weights[i] / current_gradient[i]) *
3828  current_gradient;
3829  }
3830  if (tentative_weights[i] < 0.0 ||
3831  1.0 < tentative_weights[i])
3832  {
3833  new_gradient_weight /= 2.0;
3834  tentative_weights =
3835  guess_weights +
3836  new_gradient_weight * current_gradient;
3837  }
3838  }
3839  }
3840 
3841  // the update might still send us outside the valid region, so
3842  // check again and quit if the update is still not valid
3843  if (!internal::weights_are_ok<structdim>(tentative_weights))
3844  break;
3845 
3846  // if we cannot get closer by traveling in the gradient
3847  // direction then quit
3848  if (get_point_from_weights(tentative_weights)
3849  .distance(trial_point) <
3850  get_point_from_weights(guess_weights).distance(trial_point))
3851  guess_weights = tentative_weights;
3852  else
3853  break;
3854  Assert(internal::weights_are_ok<structdim>(guess_weights),
3855  ExcInternalError());
3856  }
3857  Assert(internal::weights_are_ok<structdim>(guess_weights),
3858  ExcInternalError());
3859  projected_point = get_point_from_weights(guess_weights);
3860  }
3861 
3862  // if structdim == 2 and the optimal point is not on the interior then
3863  // we may be able to get a more accurate result by projecting onto the
3864  // lines.
3865  if (structdim == 2)
3866  {
3867  std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
3868  line_projections;
3869  for (unsigned int line_n = 0;
3870  line_n < GeometryInfo<structdim>::lines_per_cell;
3871  ++line_n)
3872  {
3873  line_projections[line_n] =
3874  project_to_object(object->line(line_n), trial_point);
3875  }
3876  std::sort(line_projections.begin(),
3877  line_projections.end(),
3878  [&](const Point<spacedim> &a, const Point<spacedim> &b) {
3879  return a.distance(trial_point) <
3880  b.distance(trial_point);
3881  });
3882  if (line_projections[0].distance(trial_point) <
3883  projected_point.distance(trial_point))
3884  projected_point = line_projections[0];
3885  }
3886  }
3887  else
3888  {
3889  Assert(false, ExcNotImplemented());
3890  return projected_point;
3891  }
3892 
3893  return projected_point;
3894  }
3895 
3896 
3897 
3898  template <int dim, typename T>
3899  template <class Archive>
3900  void
3902  const unsigned int /*version*/) const
3903  {
3904  Assert(cell_ids.size() == data.size(),
3905  ExcDimensionMismatch(cell_ids.size(), data.size()));
3906  // archive the cellids in an efficient binary format
3907  const std::size_t n_cells = cell_ids.size();
3908  ar & n_cells;
3909  for (const auto &id : cell_ids)
3910  {
3911  CellId::binary_type binary_cell_id = id.template to_binary<dim>();
3912  ar & binary_cell_id;
3913  }
3914 
3915  ar &data;
3916  }
3917 
3918 
3919 
3920  template <int dim, typename T>
3921  template <class Archive>
3922  void
3924  const unsigned int /*version*/)
3925  {
3926  std::size_t n_cells;
3927  ar & n_cells;
3928  cell_ids.clear();
3929  cell_ids.reserve(n_cells);
3930  for (unsigned int c = 0; c < n_cells; ++c)
3931  {
3933  ar & value;
3934  cell_ids.emplace_back(value);
3935  }
3936  ar &data;
3937  }
3938 
3939 
3940 
3941  template <typename DataType, typename MeshType>
3942  void
3944  const MeshType & mesh,
3945  const std::function<std_cxx17::optional<DataType>(
3946  const typename MeshType::active_cell_iterator &)> &pack,
3947  const std::function<void(const typename MeshType::active_cell_iterator &,
3948  const DataType &)> & unpack)
3949  {
3950 # ifndef DEAL_II_WITH_MPI
3951  (void)mesh;
3952  (void)pack;
3953  (void)unpack;
3954  Assert(false,
3955  ExcMessage(
3956  "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
3957 # else
3958  constexpr int dim = MeshType::dimension;
3959  constexpr int spacedim = MeshType::space_dimension;
3960  auto tria =
3961  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
3962  &mesh.get_triangulation());
3963  Assert(
3964  tria != nullptr,
3965  ExcMessage(
3966  "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3967 
3968  // map neighbor_id -> data_buffer where we accumulate the data to send
3969  using DestinationToBufferMap =
3970  std::map<::types::subdomain_id,
3972  DestinationToBufferMap destination_to_data_buffer_map;
3973 
3974  std::map<unsigned int, std::set<::types::subdomain_id>>
3977 
3978  for (const auto &cell : tria->active_cell_iterators())
3979  if (cell->is_locally_owned())
3980  {
3981  std::set<::types::subdomain_id> send_to;
3982  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
3983  {
3984  const std::map<unsigned int,
3985  std::set<::types::subdomain_id>>::
3986  const_iterator neighbor_subdomains_of_vertex =
3987  vertices_with_ghost_neighbors.find(cell->vertex_index(v));
3988 
3989  if (neighbor_subdomains_of_vertex ==
3990  vertices_with_ghost_neighbors.end())
3991  continue;
3992 
3993  Assert(neighbor_subdomains_of_vertex->second.size() != 0,
3994  ExcInternalError());
3995 
3996  send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
3997  neighbor_subdomains_of_vertex->second.end());
3998  }
3999 
4000  if (send_to.size() > 0)
4001  {
4002  // this cell's data needs to be sent to someone
4003  typename MeshType::active_cell_iterator mesh_it(tria,
4004  cell->level(),
4005  cell->index(),
4006  &mesh);
4007 
4008  const std_cxx17::optional<DataType> data = pack(mesh_it);
4009 
4010  if (data)
4011  {
4012  const CellId cellid = cell->id();
4013 
4014  for (const auto subdomain : send_to)
4015  {
4016  // find the data buffer for proc "subdomain" if it exists
4017  // or create an empty one otherwise
4018  typename DestinationToBufferMap::iterator p =
4019  destination_to_data_buffer_map
4020  .insert(std::make_pair(
4021  subdomain, CellDataTransferBuffer<dim, DataType>()))
4022  .first;
4023 
4024  p->second.cell_ids.emplace_back(cellid);
4025  p->second.data.emplace_back(*data);
4026  }
4027  }
4028  }
4029  }
4030 
4031  // Protect the following communication:
4032  static Utilities::MPI::CollectiveMutex mutex;
4034  tria->get_communicator());
4035 
4036  const int mpi_tag =
4038 
4039  // 2. send our messages
4040  std::set<::types::subdomain_id> ghost_owners = tria->ghost_owners();
4041  const unsigned int n_ghost_owners = ghost_owners.size();
4042  std::vector<std::vector<char>> sendbuffers(n_ghost_owners);
4043  std::vector<MPI_Request> requests(n_ghost_owners);
4044 
4045  unsigned int idx = 0;
4046  for (auto it = ghost_owners.begin(); it != ghost_owners.end(); ++it, ++idx)
4047  {
4048  CellDataTransferBuffer<dim, DataType> &data =
4049  destination_to_data_buffer_map[*it];
4050 
4051  // pack all the data into the buffer for this recipient and send it.
4052  // keep data around till we can make sure that the packet has been
4053  // received
4054  sendbuffers[idx] = Utilities::pack(data, /*enable_compression*/ false);
4055  const int ierr = MPI_Isend(sendbuffers[idx].data(),
4056  sendbuffers[idx].size(),
4057  MPI_BYTE,
4058  *it,
4059  mpi_tag,
4060  tria->get_communicator(),
4061  &requests[idx]);
4062  AssertThrowMPI(ierr);
4063  }
4064 
4065  // 3. receive messages
4066  std::vector<char> receive;
4067  for (unsigned int idx = 0; idx < n_ghost_owners; ++idx)
4068  {
4069  MPI_Status status;
4070  int ierr =
4071  MPI_Probe(MPI_ANY_SOURCE, mpi_tag, tria->get_communicator(), &status);
4072  AssertThrowMPI(ierr);
4073 
4074  int len;
4075  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4076  AssertThrowMPI(ierr);
4077 
4078  receive.resize(len);
4079 
4080  char *ptr = receive.data();
4081  ierr = MPI_Recv(ptr,
4082  len,
4083  MPI_BYTE,
4084  status.MPI_SOURCE,
4085  status.MPI_TAG,
4086  tria->get_communicator(),
4087  &status);
4088  AssertThrowMPI(ierr);
4089 
4090  auto cellinfo =
4091  Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(
4092  receive, /*enable_compression*/ false);
4093 
4094  DataType *data = cellinfo.data.data();
4095  for (unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++data)
4096  {
4098  tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
4099 
4100  const typename MeshType::active_cell_iterator cell(
4101  tria, tria_cell->level(), tria_cell->index(), &mesh);
4102 
4103  unpack(cell, *data);
4104  }
4105  }
4106 
4107  // make sure that all communication is finished
4108  // when we leave this function.
4109  if (requests.size())
4110  {
4111  const int ierr =
4112  MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4113  AssertThrowMPI(ierr);
4114  }
4115 # endif // DEAL_II_WITH_MPI
4116  }
4117 } // namespace GridTools
4118 
4119 # endif
4120 
4122 
4123 /*---------------------------- grid_tools.h ---------------------------*/
4124 /* end of #ifndef dealii_grid_tools_h */
4125 #endif
4126 /*---------------------------- 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:3816
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:3600
void laplace_transform(const std::map< unsigned int, Point< dim >> &new_points, Triangulation< dim > &tria, const Function< dim, double > *coefficient=nullptr, const bool solve_for_absolute_positions=false)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:211
unsigned int n_vertices() const
const types::manifold_id flat_manifold_id
Definition: types.h:273
static const unsigned int invalid_unsigned_int
Definition: types.h:191
std::map< unsigned int, Point< spacedim > > get_all_vertices_at_boundary(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:980
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
Definition: grid_tools.cc:3575
double diameter(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:76
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true)
Definition: grid_tools.cc:1013
static ::ExceptionBase & ExcPointNotFoundInCoarseGrid(Point< N > arg1)
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:192
std::map< unsigned int, Point< spacedim > > extract_used_vertices(const Triangulation< dim, spacedim > &container, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: grid_tools.cc:5174
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:3784
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void load(Archive &ar, const unsigned int version)
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:2022
std::vector< typename MeshType::active_cell_iterator > compute_ghost_cell_layer_within_distance(const MeshType &mesh, const double layer_thickness)
std::vector< unsigned int > vertex_indices
Definition: tria.cc:2244
void regularize_corner_cells(Triangulation< dim, spacedim > &tria, const double limit_angle_fraction=.75)
Definition: grid_tools.cc:3878
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12185
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2072
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
std::vector< typename MeshType::active_cell_iterator > get_active_child_cells(const typename MeshType::cell_iterator &cell)
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &local_points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes)
Definition: grid_tools.cc:4940
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:12013
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:2123
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, MPI_Comm mpi_communicator)
Definition: grid_tools.cc:5332
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:79
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)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3062
std::bitset< 3 > orientation
Definition: grid_tools.h:2481
std::pair< typename Triangulation< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Cache< dim, spacedim > &cache, const Point< spacedim > &p, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator(), const std::vector< bool > &marked_vertices={})
Definition: grid_tools.cc:5210
double cell_measure< 3 >(const std::vector< Point< 3 >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< 3 >::vertices_per_cell])
void partition_multigrid_levels(Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2921
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:3033
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
const double angle
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:372
void delete_unused_vertices(std::vector< Point< spacedim >> &vertices, std::vector< CellData< dim >> &cells, SubCellData &subcelldata)
Definition: grid_tools.cc:513
void get_vertex_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2524
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:2735
cell_iterator end() const
Definition: tria.cc:12079
std::tuple< std::vector< Point< spacedim > >, std::vector< CellData< dim > >, SubCellData > get_coarse_mesh_description(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:415
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
typename MeshType::active_cell_iterator type
Definition: grid_tools.h:92
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
static const char T
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:2553
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:1419
Signals signals
Definition: tria.h:2222
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< CellId > cell_ids
Definition: grid_tools.h:3041
Abstract base class for mapping classes.
Definition: mapping.h:302
void get_active_neighbors(const typename MeshType::active_cell_iterator &cell, std::vector< typename MeshType::active_cell_iterator > &active_neighbors)
bool have_same_coarse_mesh(const MeshType &mesh_1, const MeshType &mesh_2)
boost::geometry::index::rtree< LeafType, IndexType > RTree
Definition: rtree.h:134
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:607
#define DeclException0(Exception0)
Definition: exceptions.h:473
std::vector< BoundingBox< MeshType::space_dimension > > compute_mesh_predicate_bounding_box(const MeshType &mesh, const std::function< bool(const typename MeshType::active_cell_iterator &)> &predicate, const unsigned int refinement_level=0, const bool allow_merge=false, const unsigned int max_boxes=numbers::invalid_unsigned_int)
Definition: grid_tools.cc:1827
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, MPI_Comm mpi_communicator)
Definition: grid_tools.cc:5236
unsigned int level
Definition: grid_out.cc:4355
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:3849
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:5395
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
Definition: grid_tools.cc:3668
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:3533
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1209
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:12618
Definition: cell_id.h:69
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:1510
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
__global__ void set(Number *val, const Number s, const size_type N)
double cell_measure< 2 >(const std::vector< Point< 2 >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< 2 >::vertices_per_cell])
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:4229
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
Definition: table.h:1913
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1707
std::map< types::global_dof_index, std::vector< typename DoFHandlerType::active_cell_iterator > > get_dof_to_support_patch_map(DoFHandlerType &dof_handler)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
GridTools::exchange_cell_data_to_ghosts():
Definition: mpi_tags.h:60
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:409
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
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:1378
double cell_measure(const T &,...)
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
Definition: grid_tools.cc:258
void get_subdomain_association(const Triangulation< dim, spacedim > &triangulation, std::vector< types::subdomain_id > &subdomain)
Definition: grid_tools.cc:2948
void get_face_connectivity_of_cells(const Triangulation< dim, spacedim > &triangulation, DynamicSparsityPattern &connectivity)
Definition: grid_tools.cc:2476
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1353
Point< 3 > center
FullMatrix< double > matrix
Definition: grid_tools.h:2495
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:3696
void rotate(const double angle, const unsigned int axis, Triangulation< dim, 3 > &triangulation)
Definition: grid_tools.cc:826
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:4193
std::vector< bool > get_locally_owned_vertices(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2978
void partition_triangulation_zorder(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const bool group_siblings=true)
Definition: grid_tools.cc:2816
std::vector< typename Container::cell_iterator > get_cells_at_coarsest_common_level(const std::vector< typename Container::active_cell_iterator > &patch_cells)
double cell_measure< 1 >(const std::vector< Point< 1 >> &all_vertices, const unsigned int(&vertex_indices)[GeometryInfo< 1 >::vertices_per_cell])
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_of_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const Point< spacedim > &position, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: grid_tools.cc:1744
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:1316
std::list< std::pair< typename MeshType::cell_iterator, typename MeshType::cell_iterator > > get_finest_common_cells(const MeshType &mesh_1, const MeshType &mesh_2)
std::map< unsigned int, std::set<::types::subdomain_id > > * vertices_with_ghost_neighbors
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)
static const bool value
static ::ExceptionBase & ExcPointNotFound(Point< N > arg1)
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:817
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:2962
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:614
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:5517