Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
manifold.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 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_tria_manifold_h
17 #define dealii_tria_manifold_h
18 
19 
20 /*---------------------------- manifold.h ---------------------------*/
21 
22 #include <deal.II/base/config.h>
23 
26 #include <deal.II/base/point.h>
29 
30 #include <deal.II/grid/tria.h>
31 
33 
34 // forward declaration
35 #ifndef DOXYGEN
36 template <int, typename>
37 class Table;
38 #endif
39 
44 namespace Manifolds
45 {
52  template <typename MeshIteratorType>
53  inline constexpr std::size_t
55  {
56  // Note that in C++11 a constexpr function can only have a return
57  // statement, so we cannot alias the structure dimension
59  vertices_per_cell +
61  lines_per_cell +
63  quads_per_cell +
65  hexes_per_cell -
66  1; // don't count the cell itself, just the bounding objects
67  }
68 
111  template <typename MeshIteratorType>
112  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
113  n_default_points_per_cell<MeshIteratorType>()>,
114  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
115  get_default_points_and_weights(const MeshIteratorType &iterator,
116  const bool with_interpolation = false);
117 } // namespace Manifolds
118 
119 
120 
285 template <int dim, int spacedim = dim>
286 class Manifold : public Subscriptor
287 {
288 public:
289  // explicitly check for sensible template arguments
290  static_assert(dim <= spacedim,
291  "The dimension <dim> of a Manifold must be less than or "
292  "equal to the space dimension <spacedim> in which it lives.");
293 
294 
306  using FaceVertexNormals =
307  std::array<Tensor<1, spacedim>, GeometryInfo<dim>::vertices_per_face>;
308 
309 
314  virtual ~Manifold() override = default;
315 
321  virtual std::unique_ptr<Manifold<dim, spacedim>>
322  clone() const = 0;
323 
327 
346  virtual Point<spacedim>
347  get_intermediate_point(const Point<spacedim> &p1,
348  const Point<spacedim> &p2,
349  const double w) const;
350 
367  virtual Point<spacedim>
368  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
369  const ArrayView<const double> & weights) const;
370 
371 
392  virtual void
393  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
394  const Table<2, double> & weights,
395  ArrayView<Point<spacedim>> new_points) const;
396 
408  virtual Point<spacedim>
409  project_to_manifold(
410  const ArrayView<const Point<spacedim>> &surrounding_points,
411  const Point<spacedim> & candidate) const;
412 
426  virtual Point<spacedim>
427  get_new_point_on_line(
428  const typename Triangulation<dim, spacedim>::line_iterator &line) const;
429 
447  virtual Point<spacedim>
448  get_new_point_on_quad(
449  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const;
450 
469  virtual Point<spacedim>
470  get_new_point_on_hex(
471  const typename Triangulation<dim, spacedim>::hex_iterator &hex) const;
472 
473 
481  get_new_point_on_face(
482  const typename Triangulation<dim, spacedim>::face_iterator &face) const;
483 
484 
492  get_new_point_on_cell(
493  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
494 
496 
500 
538  virtual Tensor<1, spacedim>
539  get_tangent_vector(const Point<spacedim> &x1,
540  const Point<spacedim> &x2) const;
541 
543 
547 
594  virtual Tensor<1, spacedim>
595  normal_vector(
596  const typename Triangulation<dim, spacedim>::face_iterator &face,
597  const Point<spacedim> & p) const;
598 
613  virtual void
614  get_normals_at_vertices(
615  const typename Triangulation<dim, spacedim>::face_iterator &face,
616  FaceVertexNormals &face_vertex_normals) const;
617 
619 };
620 
621 
631 template <int dim, int spacedim = dim>
632 class FlatManifold : public Manifold<dim, spacedim>
633 {
634 public:
663  const double tolerance = 1e-10);
664 
668  virtual std::unique_ptr<Manifold<dim, spacedim>>
669  clone() const override;
670 
692  virtual Point<spacedim>
693  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
694  const ArrayView<const double> & weights) const override;
695 
696 
707  virtual void
708  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
709  const Table<2, double> & weights,
710  ArrayView<Point<spacedim>> new_points) const override;
711 
719  virtual Point<spacedim>
720  project_to_manifold(const ArrayView<const Point<spacedim>> &points,
721  const Point<spacedim> &candidate) const override;
722 
744  virtual Tensor<1, spacedim>
745  get_tangent_vector(const Point<spacedim> &x1,
746  const Point<spacedim> &x2) const override;
747 
755  virtual Tensor<1, spacedim>
756  normal_vector(
757  const typename Triangulation<dim, spacedim>::face_iterator &face,
758  const Point<spacedim> &p) const override;
759 
768  virtual void
769  get_normals_at_vertices(
770  const typename Triangulation<dim, spacedim>::face_iterator &face,
771  typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
772  const override;
773 
777  const Tensor<1, spacedim> &
778  get_periodicity() const;
779 
780 private:
795 
796  DeclException3(ExcPeriodicBox,
797  int,
799  double,
800  << "The component number " << arg1 << " of the point [ "
801  << arg2 << " ] is not in the interval [ 0, " << arg3
802  << "), bailing out.");
803 
808  const double tolerance;
809 };
810 
811 
898 template <int dim, int spacedim = dim, int chartdim = dim>
899 class ChartManifold : public Manifold<dim, spacedim>
900 {
901 public:
902  // explicitly check for sensible template arguments
903  static_assert(dim <= spacedim,
904  "The dimension <dim> of a ChartManifold must be less than or "
905  "equal to the space dimension <spacedim> in which it lives.");
906 
921  ChartManifold(const Tensor<1, chartdim> &periodicity = Tensor<1, chartdim>());
922 
927  virtual ~ChartManifold() override = default;
928 
933  virtual Point<spacedim>
934  get_intermediate_point(const Point<spacedim> &p1,
935  const Point<spacedim> &p2,
936  const double w) const override;
937 
942  virtual Point<spacedim>
943  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
944  const ArrayView<const double> & weights) const override;
945 
967  virtual void
968  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
969  const Table<2, double> & weights,
970  ArrayView<Point<spacedim>> new_points) const override;
977  virtual Point<chartdim>
978  pull_back(const Point<spacedim> &space_point) const = 0;
979 
986  virtual Point<spacedim>
987  push_forward(const Point<chartdim> &chart_point) const = 0;
988 
1006  push_forward_gradient(const Point<chartdim> &chart_point) const;
1007 
1063  virtual Tensor<1, spacedim>
1064  get_tangent_vector(const Point<spacedim> &x1,
1065  const Point<spacedim> &x2) const override;
1066 
1070  const Tensor<1, chartdim> &
1071  get_periodicity() const;
1072 
1073 private:
1086 };
1087 
1088 
1089 /* -------------- declaration of explicit specializations ------------- */
1090 
1091 #ifndef DOXYGEN
1092 
1093 template <>
1094 Point<1>
1096  const Triangulation<1, 1>::face_iterator &) const;
1097 
1098 template <>
1099 Point<2>
1101  const Triangulation<1, 2>::face_iterator &) const;
1102 
1103 
1104 template <>
1105 Point<3>
1107  const Triangulation<1, 3>::face_iterator &) const;
1108 
1109 
1110 template <>
1111 Point<1>
1113  const Triangulation<1, 1>::quad_iterator &) const;
1114 
1115 template <>
1116 Point<2>
1118  const Triangulation<1, 2>::quad_iterator &) const;
1119 
1120 
1121 template <>
1122 Point<3>
1124  const Triangulation<1, 3>::quad_iterator &) const;
1125 
1126 
1127 template <>
1128 Point<3>
1130  const Triangulation<3, 3>::hex_iterator &) const;
1131 
1132 /*---Templated functions---*/
1133 
1134 namespace Manifolds
1135 {
1136  template <typename MeshIteratorType>
1137  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
1138  n_default_points_per_cell<MeshIteratorType>()>,
1139  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
1140  get_default_points_and_weights(const MeshIteratorType &iterator,
1141  const bool with_interpolation)
1142  {
1143  const int dim = MeshIteratorType::AccessorType::structure_dimension;
1144  const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1145  constexpr std::size_t points_per_cell =
1146  n_default_points_per_cell<MeshIteratorType>();
1147 
1148  std::pair<std::array<Point<spacedim>, points_per_cell>,
1149  std::array<double, points_per_cell>>
1150  points_weights;
1151 
1152 
1153  // note that the exact weights are chosen such as to minimize the
1154  // distortion of the four new quads from the optimal shape; their
1155  // derivation and values is copied over from the
1156  // interpolation function in the mapping
1157  switch (dim)
1158  {
1159  case 1:
1160  Assert(points_weights.first.size() == 2, ExcInternalError());
1161  Assert(points_weights.second.size() == 2, ExcInternalError());
1162  points_weights.first[0] = iterator->vertex(0);
1163  points_weights.second[0] = .5;
1164  points_weights.first[1] = iterator->vertex(1);
1165  points_weights.second[1] = .5;
1166  break;
1167  case 2:
1168  Assert(points_weights.first.size() == 8, ExcInternalError());
1169  Assert(points_weights.second.size() == 8, ExcInternalError());
1170 
1171  for (unsigned int i = 0; i < 4; ++i)
1172  {
1173  points_weights.first[i] = iterator->vertex(i);
1174  points_weights.first[4 + i] =
1175  (iterator->line(i)->has_children() ?
1176  iterator->line(i)->child(0)->vertex(1) :
1177  iterator->line(i)->get_manifold().get_new_point_on_line(
1178  iterator->line(i)));
1179  }
1180 
1181  if (with_interpolation)
1182  {
1183  std::fill(points_weights.second.begin(),
1184  points_weights.second.begin() + 4,
1185  -0.25);
1186  std::fill(points_weights.second.begin() + 4,
1187  points_weights.second.end(),
1188  0.5);
1189  }
1190  else
1191  std::fill(points_weights.second.begin(),
1192  points_weights.second.end(),
1193  1.0 / 8.0);
1194  break;
1195  case 3:
1196  {
1198  static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
1199  const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
1202  Assert(points_weights.first.size() == np, ExcInternalError());
1203  Assert(points_weights.second.size() == np, ExcInternalError());
1204  auto *sp3 = reinterpret_cast<
1205  std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
1206  *>(&points_weights.first);
1207 
1208  unsigned int j = 0;
1209 
1210  // note that the exact weights are chosen such as to minimize the
1211  // distortion of the eight new hexes from the optimal shape through
1212  // transfinite interpolation from the faces and vertices, see
1213  // TransfiniteInterpolationManifold for a deeper explanation of the
1214  // mechanisms
1215  if (with_interpolation)
1216  {
1217  for (unsigned int i = 0;
1218  i < GeometryInfo<dim>::vertices_per_cell;
1219  ++i, ++j)
1220  {
1221  (*sp3)[j] = hex->vertex(i);
1222  points_weights.second[j] = 1.0 / 8.0;
1223  }
1224  for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
1225  ++i, ++j)
1226  {
1227  (*sp3)[j] =
1228  (hex->line(i)->has_children() ?
1229  hex->line(i)->child(0)->vertex(1) :
1230  hex->line(i)->get_manifold().get_new_point_on_line(
1231  hex->line(i)));
1232  points_weights.second[j] = -1.0 / 4.0;
1233  }
1234  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1235  ++i, ++j)
1236  {
1237  (*sp3)[j] =
1238  (hex->quad(i)->has_children() ?
1239  hex->quad(i)->isotropic_child(0)->vertex(3) :
1240  hex->quad(i)->get_manifold().get_new_point_on_quad(
1241  hex->quad(i)));
1242  points_weights.second[j] = 1.0 / 2.0;
1243  }
1244  }
1245  else
1246  // Overwrite the weights with 1/np if we don't want to use
1247  // interpolation.
1248  std::fill(points_weights.second.begin(),
1249  points_weights.second.end(),
1250  1.0 / np);
1251  }
1252  break;
1253  default:
1254  Assert(false, ExcInternalError());
1255  break;
1256  }
1257  return points_weights;
1258  }
1259 } // namespace Manifolds
1260 
1261 #endif // DOXYGEN
1262 
1264 
1265 #endif
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)
Point< dim > push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
Definition: utilities.cc:828
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1418
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
Definition: manifold.cc:450
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1466
constexpr std::size_t n_default_points_per_cell()
Definition: manifold.h:54
#define Assert(cond, exc)
Definition: exceptions.h:1403
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1442
const double tolerance
Definition: manifold.h:808
const FlatManifold< chartdim, chartdim > sub_manifold
Definition: manifold.h:1085
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:330
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:307
Definition: table.h:687
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:561
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:344
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:794
static ::ExceptionBase & ExcInternalError()