Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
manifold.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2018 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_tria_manifold_h
17 #define dealii_tria_manifold_h
18 
19 
20 /*---------------------------- manifold.h ---------------------------*/
21 
22 #include <deal.II/base/config.h>
23 
24 #include <deal.II/base/array_view.h>
25 #include <deal.II/base/derivative_form.h>
26 #include <deal.II/base/point.h>
27 #include <deal.II/base/quadrature_lib.h>
28 #include <deal.II/base/subscriptor.h>
29 #include <deal.II/base/thread_management.h>
30 
31 #include <deal.II/grid/tria.h>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 // forward declaration
36 template <int, typename>
37 class Table;
38 
43 namespace Manifolds
44 {
51  template <typename MeshIteratorType>
52  inline constexpr std::size_t
54  {
55  // Note that in C++11 a constexpr function can only have a return
56  // statement, so we cannot alias the structure dimension
58  vertices_per_cell +
60  lines_per_cell +
62  quads_per_cell +
64  hexes_per_cell -
65  1; // don't count the cell itself, just the bounding objects
66  }
67 
110  template <typename MeshIteratorType>
112  get_default_quadrature(const MeshIteratorType &iterator,
113  const bool with_interpolation = false);
114 
157  template <typename MeshIteratorType>
158  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
159  n_default_points_per_cell<MeshIteratorType>()>,
160  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
161  get_default_points_and_weights(const MeshIteratorType &iterator,
162  const bool with_interpolation = false);
163 } // namespace Manifolds
164 
328 template <int dim, int spacedim = dim>
329 class Manifold : public Subscriptor
330 {
331 public:
332  // explicitly check for sensible template arguments
333  static_assert(dim <= spacedim,
334  "The dimension <dim> of a Manifold must be less than or "
335  "equal to the space dimension <spacedim> in which it lives.");
336 
337 
349  using FaceVertexNormals =
350  std::array<Tensor<1, spacedim>, GeometryInfo<dim>::vertices_per_face>;
351 
352 
357  virtual ~Manifold() override = default;
358 
364  virtual std::unique_ptr<Manifold<dim, spacedim>>
365  clone() const = 0;
366 
370 
389  virtual Point<spacedim>
391  const Point<spacedim> &p2,
392  const double w) const;
393 
410  virtual Point<spacedim>
411  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
412  const ArrayView<const double> & weights) const;
413 
414 
435  virtual void
436  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
437  const Table<2, double> & weights,
438  ArrayView<Point<spacedim>> new_points) const;
439 
451  virtual Point<spacedim>
453  const ArrayView<const Point<spacedim>> &surrounding_points,
454  const Point<spacedim> & candidate) const;
455 
469  virtual Point<spacedim>
471  const typename Triangulation<dim, spacedim>::line_iterator &line) const;
472 
490  virtual Point<spacedim>
492  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const;
493 
512  virtual Point<spacedim>
514  const typename Triangulation<dim, spacedim>::hex_iterator &hex) const;
515 
516 
525  const typename Triangulation<dim, spacedim>::face_iterator &face) const;
526 
527 
536  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
537 
539 
543 
581  virtual Tensor<1, spacedim>
583  const Point<spacedim> &x2) const;
584 
586 
590 
637  virtual Tensor<1, spacedim>
639  const typename Triangulation<dim, spacedim>::face_iterator &face,
640  const Point<spacedim> & p) const;
641 
656  virtual void
658  const typename Triangulation<dim, spacedim>::face_iterator &face,
659  FaceVertexNormals &face_vertex_normals) const;
660 
662 };
663 
664 
676 template <int dim, int spacedim = dim>
677 class FlatManifold : public Manifold<dim, spacedim>
678 {
679 public:
708  const double tolerance = 1e-10);
709 
713  virtual std::unique_ptr<Manifold<dim, spacedim>>
714  clone() const override;
715 
737  virtual Point<spacedim>
738  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
739  const ArrayView<const double> & weights) const override;
740 
741 
752  virtual void
753  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
754  const Table<2, double> & weights,
755  ArrayView<Point<spacedim>> new_points) const override;
756 
764  virtual Point<spacedim>
765  project_to_manifold(const ArrayView<const Point<spacedim>> &points,
766  const Point<spacedim> &candidate) const override;
767 
789  virtual Tensor<1, spacedim>
791  const Point<spacedim> &x2) const override;
792 
800  virtual Tensor<1, spacedim>
802  const typename Triangulation<dim, spacedim>::face_iterator &face,
803  const Point<spacedim> &p) const override;
804 
813  virtual void
815  const typename Triangulation<dim, spacedim>::face_iterator &face,
816  typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
817  const override;
818 
822  const Tensor<1, spacedim> &
823  get_periodicity() const;
824 
825 private:
840 
842  int,
844  double,
845  << "The component number " << arg1 << " of the point [ "
846  << arg2 << " ] is not in the interval [ 0, " << arg3
847  << "), bailing out.");
848 
853  const double tolerance;
854 };
855 
856 
945 template <int dim, int spacedim = dim, int chartdim = dim>
946 class ChartManifold : public Manifold<dim, spacedim>
947 {
948 public:
949  // explicitly check for sensible template arguments
950  static_assert(dim <= spacedim,
951  "The dimension <dim> of a ChartManifold must be less than or "
952  "equal to the space dimension <spacedim> in which it lives.");
953 
968  ChartManifold(const Tensor<1, chartdim> &periodicity = Tensor<1, chartdim>());
969 
974  virtual ~ChartManifold() override = default;
975 
980  virtual Point<spacedim>
982  const Point<spacedim> &p2,
983  const double w) const override;
984 
989  virtual Point<spacedim>
990  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
991  const ArrayView<const double> & weights) const override;
992 
1014  virtual void
1015  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
1016  const Table<2, double> & weights,
1017  ArrayView<Point<spacedim>> new_points) const override;
1024  virtual Point<chartdim>
1025  pull_back(const Point<spacedim> &space_point) const = 0;
1026 
1033  virtual Point<spacedim>
1034  push_forward(const Point<chartdim> &chart_point) const = 0;
1035 
1053  push_forward_gradient(const Point<chartdim> &chart_point) const;
1054 
1110  virtual Tensor<1, spacedim>
1112  const Point<spacedim> &x2) const override;
1113 
1117  const Tensor<1, chartdim> &
1118  get_periodicity() const;
1119 
1120 private:
1133 };
1134 
1135 
1136 /* -------------- declaration of explicit specializations ------------- */
1137 
1138 #ifndef DOXYGEN
1139 
1140 template <>
1141 Point<1>
1143  const Triangulation<1, 1>::face_iterator &) const;
1144 
1145 template <>
1146 Point<2>
1148  const Triangulation<1, 2>::face_iterator &) const;
1149 
1150 
1151 template <>
1152 Point<3>
1154  const Triangulation<1, 3>::face_iterator &) const;
1155 
1156 
1157 template <>
1158 Point<1>
1160  const Triangulation<1, 1>::quad_iterator &) const;
1161 
1162 template <>
1163 Point<2>
1165  const Triangulation<1, 2>::quad_iterator &) const;
1166 
1167 
1168 template <>
1169 Point<3>
1171  const Triangulation<1, 3>::quad_iterator &) const;
1172 
1173 
1174 template <>
1175 Point<3>
1177  const Triangulation<3, 3>::hex_iterator &) const;
1178 
1179 /*---Templated functions---*/
1180 
1181 namespace Manifolds
1182 {
1183  template <typename MeshIteratorType>
1185  get_default_quadrature(const MeshIteratorType &iterator,
1186  const bool with_interpolation)
1187  {
1188  const auto points_and_weights =
1189  get_default_points_and_weights(iterator, with_interpolation);
1190  static const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1191  return Quadrature<spacedim>(
1192  std::vector<Point<spacedim>>(points_and_weights.first.begin(),
1193  points_and_weights.first.end()),
1194  std::vector<double>(points_and_weights.second.begin(),
1195  points_and_weights.second.end()));
1196  }
1197 
1198 
1199 
1200  template <typename MeshIteratorType>
1201  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
1202  n_default_points_per_cell<MeshIteratorType>()>,
1203  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
1204  get_default_points_and_weights(const MeshIteratorType &iterator,
1205  const bool with_interpolation)
1206  {
1207  const int dim = MeshIteratorType::AccessorType::structure_dimension;
1208  const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1209  constexpr std::size_t points_per_cell =
1210  n_default_points_per_cell<MeshIteratorType>();
1211 
1212  std::pair<std::array<Point<spacedim>, points_per_cell>,
1213  std::array<double, points_per_cell>>
1214  points_weights;
1215 
1216 
1217  // note that the exact weights are chosen such as to minimize the
1218  // distortion of the four new quads from the optimal shape; their
1219  // derivation and values is copied over from the
1220  // interpolation function in the mapping
1221  switch (dim)
1222  {
1223  case 1:
1224  Assert(points_weights.first.size() == 2, ExcInternalError());
1225  Assert(points_weights.second.size() == 2, ExcInternalError());
1226  points_weights.first[0] = iterator->vertex(0);
1227  points_weights.second[0] = .5;
1228  points_weights.first[1] = iterator->vertex(1);
1229  points_weights.second[1] = .5;
1230  break;
1231  case 2:
1232  Assert(points_weights.first.size() == 8, ExcInternalError());
1233  Assert(points_weights.second.size() == 8, ExcInternalError());
1234 
1235  for (unsigned int i = 0; i < 4; ++i)
1236  {
1237  points_weights.first[i] = iterator->vertex(i);
1238  points_weights.first[4 + i] =
1239  (iterator->line(i)->has_children() ?
1240  iterator->line(i)->child(0)->vertex(1) :
1241  iterator->line(i)->get_manifold().get_new_point_on_line(
1242  iterator->line(i)));
1243  }
1244 
1245  if (with_interpolation)
1246  {
1247  std::fill(points_weights.second.begin(),
1248  points_weights.second.begin() + 4,
1249  -0.25);
1250  std::fill(points_weights.second.begin() + 4,
1251  points_weights.second.end(),
1252  0.5);
1253  }
1254  else
1255  std::fill(points_weights.second.begin(),
1256  points_weights.second.end(),
1257  1.0 / 8.0);
1258  break;
1259  case 3:
1260  {
1262  static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
1263  const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
1266  Assert(points_weights.first.size() == np, ExcInternalError());
1267  Assert(points_weights.second.size() == np, ExcInternalError());
1268  auto *sp3 = reinterpret_cast<
1269  std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
1270  *>(&points_weights.first);
1271 
1272  unsigned int j = 0;
1273 
1274  // note that the exact weights are chosen such as to minimize the
1275  // distortion of the eight new hexes from the optimal shape through
1276  // transfinite interpolation from the faces and vertices, see
1277  // TransfiniteInterpolationManifold for a deeper explanation of the
1278  // mechanisms
1279  if (with_interpolation)
1280  {
1281  for (unsigned int i = 0;
1282  i < GeometryInfo<dim>::vertices_per_cell;
1283  ++i, ++j)
1284  {
1285  (*sp3)[j] = hex->vertex(i);
1286  points_weights.second[j] = 1.0 / 8.0;
1287  }
1288  for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
1289  ++i, ++j)
1290  {
1291  (*sp3)[j] =
1292  (hex->line(i)->has_children() ?
1293  hex->line(i)->child(0)->vertex(1) :
1294  hex->line(i)->get_manifold().get_new_point_on_line(
1295  hex->line(i)));
1296  points_weights.second[j] = -1.0 / 4.0;
1297  }
1298  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1299  ++i, ++j)
1300  {
1301  (*sp3)[j] =
1302  (hex->quad(i)->has_children() ?
1303  hex->quad(i)->isotropic_child(0)->vertex(3) :
1304  hex->quad(i)->get_manifold().get_new_point_on_quad(
1305  hex->quad(i)));
1306  points_weights.second[j] = 1.0 / 2.0;
1307  }
1308  }
1309  else
1310  // Overwrite the weights with 1/np if we don't want to use
1311  // interpolation.
1312  std::fill(points_weights.second.begin(),
1313  points_weights.second.end(),
1314  1.0 / np);
1315  }
1316  break;
1317  default:
1318  Assert(false, ExcInternalError());
1319  break;
1320  }
1321  return points_weights;
1322  }
1323 } // namespace Manifolds
1324 
1325 #endif // DOXYGEN
1326 
1327 DEAL_II_NAMESPACE_CLOSE
1328 
1329 #endif
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
Definition: manifold.cc:511
std::pair< std::array< Point< MeshIteratorType::AccessorType::space_dimension >, n_default_points_per_cell< MeshIteratorType >)>, std::array< double, n_default_points_per_cell< MeshIteratorType >)> > get_default_points_and_weights(const MeshIteratorType &iterator, const bool with_interpolation=false)
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
Definition: manifold.cc:853
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1592
const Tensor< 1, spacedim > & get_periodicity() const
Definition: manifold.cc:686
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
Definition: manifold.cc:444
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1640
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const =0
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const
Definition: manifold.cc:63
virtual ~ChartManifold() override=default
constexpr std::size_t n_default_points_per_cell()
Definition: manifold.h:53
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:310
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Definition: manifold.cc:469
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition: manifold.cc:1033
Quadrature< MeshIteratorType::AccessorType::space_dimension > get_default_quadrature(const MeshIteratorType &iterator, const bool with_interpolation=false)
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const
Definition: manifold.cc:38
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1078
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
Definition: manifold.cc:675
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Definition: manifold.cc:522
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1616
const double tolerance
Definition: manifold.h:853
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Definition: manifold.cc:50
const FlatManifold< chartdim, chartdim > sub_manifold
Definition: manifold.h:1132
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
Definition: manifold.cc:1004
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1046
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const
Definition: manifold.cc:117
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:695
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:983
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:324
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:532
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const =0
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:350
Definition: table.h:37
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:231
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
Definition: manifold.cc:593
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:959
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:338
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:295
virtual Point< chartdim > pull_back(const Point< spacedim > &space_point) const =0
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:968
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:839
static ::ExceptionBase & ExcInternalError()
virtual ~Manifold() override=default
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:358