Reference documentation for deal.II version GIT 1e9e64a91d 2022-09-28 19:20:02+00:00
\(\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 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_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 
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 
328 
346  virtual Point<spacedim>
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>
410  const ArrayView<const Point<spacedim>> &surrounding_points,
411  const Point<spacedim> & candidate) const;
412 
426  virtual Point<spacedim>
428  const typename Triangulation<dim, spacedim>::line_iterator &line) const;
429 
447  virtual Point<spacedim>
449  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const;
450 
469  virtual Point<spacedim>
471  const typename Triangulation<dim, spacedim>::hex_iterator &hex) const;
472 
473 
482  const typename Triangulation<dim, spacedim>::face_iterator &face) const;
483 
484 
493  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
494 
538  virtual Tensor<1, spacedim>
540  const Point<spacedim> &x2) const;
541 
597  virtual Tensor<1, spacedim>
599  const typename Triangulation<dim, spacedim>::face_iterator &face,
600  const Point<spacedim> & p) const;
601 
616  virtual void
618  const typename Triangulation<dim, spacedim>::face_iterator &face,
619  FaceVertexNormals &face_vertex_normals) const;
620 
622 };
623 
624 
634 template <int dim, int spacedim = dim>
635 class FlatManifold : public Manifold<dim, spacedim>
636 {
637 public:
666  const double tolerance = 1e-10);
667 
671  virtual std::unique_ptr<Manifold<dim, spacedim>>
672  clone() const override;
673 
695  virtual Point<spacedim>
696  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
697  const ArrayView<const double> & weights) const override;
698 
699 
710  virtual void
711  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
712  const Table<2, double> & weights,
713  ArrayView<Point<spacedim>> new_points) const override;
714 
722  virtual Point<spacedim>
724  const Point<spacedim> &candidate) const override;
725 
747  virtual Tensor<1, spacedim>
749  const Point<spacedim> &x2) const override;
750 
758  virtual Tensor<1, spacedim>
760  const typename Triangulation<dim, spacedim>::face_iterator &face,
761  const Point<spacedim> &p) const override;
762 
771  virtual void
773  const typename Triangulation<dim, spacedim>::face_iterator &face,
774  typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
775  const override;
776 
780  const Tensor<1, spacedim> &
782 
783 private:
798 
800  int,
802  double,
803  << "The component number " << arg1 << " of the point [ "
804  << arg2 << " ] is not in the interval [ 0, " << arg3
805  << "), bailing out.");
806 
811  const double tolerance;
812 };
813 
814 
901 template <int dim, int spacedim = dim, int chartdim = dim>
902 class ChartManifold : public Manifold<dim, spacedim>
903 {
904 public:
905  // explicitly check for sensible template arguments
906  static_assert(dim <= spacedim,
907  "The dimension <dim> of a ChartManifold must be less than or "
908  "equal to the space dimension <spacedim> in which it lives.");
909 
924  ChartManifold(const Tensor<1, chartdim> &periodicity = Tensor<1, chartdim>());
925 
930  virtual ~ChartManifold() override = default;
931 
936  virtual Point<spacedim>
938  const Point<spacedim> &p2,
939  const double w) const override;
940 
945  virtual Point<spacedim>
946  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
947  const ArrayView<const double> & weights) const override;
948 
970  virtual void
971  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
972  const Table<2, double> & weights,
973  ArrayView<Point<spacedim>> new_points) const override;
980  virtual Point<chartdim>
981  pull_back(const Point<spacedim> &space_point) const = 0;
982 
989  virtual Point<spacedim>
990  push_forward(const Point<chartdim> &chart_point) const = 0;
991 
1009  push_forward_gradient(const Point<chartdim> &chart_point) const;
1010 
1066  virtual Tensor<1, spacedim>
1068  const Point<spacedim> &x2) const override;
1069 
1073  const Tensor<1, chartdim> &
1074  get_periodicity() const;
1075 
1076 private:
1089 };
1090 
1091 
1092 /* -------------- declaration of explicit specializations ------------- */
1093 
1094 #ifndef DOXYGEN
1095 
1096 template <>
1097 Point<1>
1099  const Triangulation<1, 1>::face_iterator &) const;
1100 
1101 template <>
1102 Point<2>
1104  const Triangulation<1, 2>::face_iterator &) const;
1105 
1106 
1107 template <>
1108 Point<3>
1110  const Triangulation<1, 3>::face_iterator &) const;
1111 
1112 
1113 template <>
1114 Point<1>
1116  const Triangulation<1, 1>::quad_iterator &) const;
1117 
1118 template <>
1119 Point<2>
1121  const Triangulation<1, 2>::quad_iterator &) const;
1122 
1123 
1124 template <>
1125 Point<3>
1127  const Triangulation<1, 3>::quad_iterator &) const;
1128 
1129 
1130 template <>
1131 Point<3>
1133  const Triangulation<3, 3>::hex_iterator &) const;
1134 
1135 /*---Templated functions---*/
1136 
1137 namespace Manifolds
1138 {
1139  template <typename MeshIteratorType>
1140  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
1141  n_default_points_per_cell<MeshIteratorType>()>,
1142  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
1143  get_default_points_and_weights(const MeshIteratorType &iterator,
1144  const bool with_interpolation)
1145  {
1146  const int dim = MeshIteratorType::AccessorType::structure_dimension;
1147  const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1148  constexpr std::size_t points_per_cell =
1149  n_default_points_per_cell<MeshIteratorType>();
1150 
1151  std::pair<std::array<Point<spacedim>, points_per_cell>,
1152  std::array<double, points_per_cell>>
1153  points_weights;
1154 
1155 
1156  // note that the exact weights are chosen such as to minimize the
1157  // distortion of the four new quads from the optimal shape; their
1158  // derivation and values is copied over from the
1159  // interpolation function in the mapping
1160  switch (dim)
1161  {
1162  case 1:
1163  Assert(points_weights.first.size() == 2, ExcInternalError());
1164  Assert(points_weights.second.size() == 2, ExcInternalError());
1165  points_weights.first[0] = iterator->vertex(0);
1166  points_weights.second[0] = .5;
1167  points_weights.first[1] = iterator->vertex(1);
1168  points_weights.second[1] = .5;
1169  break;
1170  case 2:
1171  Assert(points_weights.first.size() == 8, ExcInternalError());
1172  Assert(points_weights.second.size() == 8, ExcInternalError());
1173 
1174  for (unsigned int i = 0; i < 4; ++i)
1175  {
1176  points_weights.first[i] = iterator->vertex(i);
1177  points_weights.first[4 + i] =
1178  (iterator->line(i)->has_children() ?
1179  iterator->line(i)->child(0)->vertex(1) :
1180  iterator->line(i)->get_manifold().get_new_point_on_line(
1181  iterator->line(i)));
1182  }
1183 
1184  if (with_interpolation)
1185  {
1186  std::fill(points_weights.second.begin(),
1187  points_weights.second.begin() + 4,
1188  -0.25);
1189  std::fill(points_weights.second.begin() + 4,
1190  points_weights.second.end(),
1191  0.5);
1192  }
1193  else
1194  std::fill(points_weights.second.begin(),
1195  points_weights.second.end(),
1196  1.0 / 8.0);
1197  break;
1198  case 3:
1199  {
1201  static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
1202  const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
1205  Assert(points_weights.first.size() == np, ExcInternalError());
1206  Assert(points_weights.second.size() == np, ExcInternalError());
1207  auto *sp3 = reinterpret_cast<
1208  std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
1209  *>(&points_weights.first);
1210 
1211  unsigned int j = 0;
1212 
1213  // note that the exact weights are chosen such as to minimize the
1214  // distortion of the eight new hexes from the optimal shape through
1215  // transfinite interpolation from the faces and vertices, see
1216  // TransfiniteInterpolationManifold for a deeper explanation of the
1217  // mechanisms
1218  if (with_interpolation)
1219  {
1220  for (unsigned int i = 0;
1221  i < GeometryInfo<dim>::vertices_per_cell;
1222  ++i, ++j)
1223  {
1224  (*sp3)[j] = hex->vertex(i);
1225  points_weights.second[j] = 1.0 / 8.0;
1226  }
1227  for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
1228  ++i, ++j)
1229  {
1230  (*sp3)[j] =
1231  (hex->line(i)->has_children() ?
1232  hex->line(i)->child(0)->vertex(1) :
1233  hex->line(i)->get_manifold().get_new_point_on_line(
1234  hex->line(i)));
1235  points_weights.second[j] = -1.0 / 4.0;
1236  }
1237  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1238  ++i, ++j)
1239  {
1240  (*sp3)[j] =
1241  (hex->quad(i)->has_children() ?
1242  hex->quad(i)->isotropic_child(0)->vertex(3) :
1243  hex->quad(i)->get_manifold().get_new_point_on_quad(
1244  hex->quad(i)));
1245  points_weights.second[j] = 1.0 / 2.0;
1246  }
1247  }
1248  else
1249  // Overwrite the weights with 1/np if we don't want to use
1250  // interpolation.
1251  std::fill(points_weights.second.begin(),
1252  points_weights.second.end(),
1253  1.0 / np);
1254  }
1255  break;
1256  default:
1257  Assert(false, ExcInternalError());
1258  break;
1259  }
1260  return points_weights;
1261  }
1262 } // namespace Manifolds
1263 
1264 #endif // DOXYGEN
1265 
1267 
1268 #endif
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:1036
const FlatManifold< chartdim, chartdim > sub_manifold
Definition: manifold.h:1088
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:1012
virtual Point< chartdim > pull_back(const Point< spacedim > &space_point) const =0
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const =0
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1099
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition: manifold.cc:1086
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1131
virtual ~ChartManifold() override=default
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:1057
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:1021
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
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 > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
const Tensor< 1, spacedim > & get_periodicity() const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:797
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
const double tolerance
Definition: manifold.h:811
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:125
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:65
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const
Definition: manifold.cc:40
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:339
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:368
virtual ~Manifold() override=default
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:307
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
Definition: manifold.cc:488
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:275
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Definition: manifold.cc:52
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:382
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Definition: manifold.cc:514
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:402
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:354
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1503
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1479
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)
constexpr std::size_t n_default_points_per_cell()
Definition: manifold.h:54
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)