Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +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>
113  get_default_quadrature(const MeshIteratorType &iterator,
114  const bool with_interpolation = false);
115 
158  template <typename MeshIteratorType>
159  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
160  n_default_points_per_cell<MeshIteratorType>()>,
161  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
162  get_default_points_and_weights(const MeshIteratorType &iterator,
163  const bool with_interpolation = false);
164 } // namespace Manifolds
165 
166 
167 
333 template <int dim, int spacedim = dim>
334 class Manifold : public Subscriptor
335 {
336 public:
337  // explicitly check for sensible template arguments
338  static_assert(dim <= spacedim,
339  "The dimension <dim> of a Manifold must be less than or "
340  "equal to the space dimension <spacedim> in which it lives.");
341 
342 
354  using FaceVertexNormals =
355  std::array<Tensor<1, spacedim>, GeometryInfo<dim>::vertices_per_face>;
356 
357 
362  virtual ~Manifold() override = default;
363 
369  virtual std::unique_ptr<Manifold<dim, spacedim>>
370  clone() const = 0;
371 
375 
394  virtual Point<spacedim>
395  get_intermediate_point(const Point<spacedim> &p1,
396  const Point<spacedim> &p2,
397  const double w) const;
398 
415  virtual Point<spacedim>
416  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
417  const ArrayView<const double> & weights) const;
418 
419 
440  virtual void
441  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
442  const Table<2, double> & weights,
443  ArrayView<Point<spacedim>> new_points) const;
444 
456  virtual Point<spacedim>
457  project_to_manifold(
458  const ArrayView<const Point<spacedim>> &surrounding_points,
459  const Point<spacedim> & candidate) const;
460 
474  virtual Point<spacedim>
475  get_new_point_on_line(
476  const typename Triangulation<dim, spacedim>::line_iterator &line) const;
477 
495  virtual Point<spacedim>
496  get_new_point_on_quad(
497  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const;
498 
517  virtual Point<spacedim>
518  get_new_point_on_hex(
519  const typename Triangulation<dim, spacedim>::hex_iterator &hex) const;
520 
521 
529  get_new_point_on_face(
530  const typename Triangulation<dim, spacedim>::face_iterator &face) const;
531 
532 
540  get_new_point_on_cell(
541  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
542 
544 
548 
586  virtual Tensor<1, spacedim>
587  get_tangent_vector(const Point<spacedim> &x1,
588  const Point<spacedim> &x2) const;
589 
591 
595 
642  virtual Tensor<1, spacedim>
643  normal_vector(
644  const typename Triangulation<dim, spacedim>::face_iterator &face,
645  const Point<spacedim> & p) const;
646 
661  virtual void
662  get_normals_at_vertices(
663  const typename Triangulation<dim, spacedim>::face_iterator &face,
664  FaceVertexNormals &face_vertex_normals) const;
665 
667 };
668 
669 
681 template <int dim, int spacedim = dim>
682 class FlatManifold : public Manifold<dim, spacedim>
683 {
684 public:
713  const double tolerance = 1e-10);
714 
718  virtual std::unique_ptr<Manifold<dim, spacedim>>
719  clone() const override;
720 
742  virtual Point<spacedim>
743  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
744  const ArrayView<const double> & weights) const override;
745 
746 
757  virtual void
758  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
759  const Table<2, double> & weights,
760  ArrayView<Point<spacedim>> new_points) const override;
761 
769  virtual Point<spacedim>
770  project_to_manifold(const ArrayView<const Point<spacedim>> &points,
771  const Point<spacedim> &candidate) const override;
772 
794  virtual Tensor<1, spacedim>
795  get_tangent_vector(const Point<spacedim> &x1,
796  const Point<spacedim> &x2) const override;
797 
805  virtual Tensor<1, spacedim>
806  normal_vector(
807  const typename Triangulation<dim, spacedim>::face_iterator &face,
808  const Point<spacedim> &p) const override;
809 
818  virtual void
819  get_normals_at_vertices(
820  const typename Triangulation<dim, spacedim>::face_iterator &face,
821  typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
822  const override;
823 
827  const Tensor<1, spacedim> &
828  get_periodicity() const;
829 
830 private:
845 
846  DeclException3(ExcPeriodicBox,
847  int,
849  double,
850  << "The component number " << arg1 << " of the point [ "
851  << arg2 << " ] is not in the interval [ 0, " << arg3
852  << "), bailing out.");
853 
858  const double tolerance;
859 };
860 
861 
950 template <int dim, int spacedim = dim, int chartdim = dim>
951 class ChartManifold : public Manifold<dim, spacedim>
952 {
953 public:
954  // explicitly check for sensible template arguments
955  static_assert(dim <= spacedim,
956  "The dimension <dim> of a ChartManifold must be less than or "
957  "equal to the space dimension <spacedim> in which it lives.");
958 
973  ChartManifold(const Tensor<1, chartdim> &periodicity = Tensor<1, chartdim>());
974 
979  virtual ~ChartManifold() override = default;
980 
985  virtual Point<spacedim>
986  get_intermediate_point(const Point<spacedim> &p1,
987  const Point<spacedim> &p2,
988  const double w) const override;
989 
994  virtual Point<spacedim>
995  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
996  const ArrayView<const double> & weights) const override;
997 
1019  virtual void
1020  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
1021  const Table<2, double> & weights,
1022  ArrayView<Point<spacedim>> new_points) const override;
1029  virtual Point<chartdim>
1030  pull_back(const Point<spacedim> &space_point) const = 0;
1031 
1038  virtual Point<spacedim>
1039  push_forward(const Point<chartdim> &chart_point) const = 0;
1040 
1058  push_forward_gradient(const Point<chartdim> &chart_point) const;
1059 
1115  virtual Tensor<1, spacedim>
1116  get_tangent_vector(const Point<spacedim> &x1,
1117  const Point<spacedim> &x2) const override;
1118 
1122  const Tensor<1, chartdim> &
1123  get_periodicity() const;
1124 
1125 private:
1138 };
1139 
1140 
1141 /* -------------- declaration of explicit specializations ------------- */
1142 
1143 #ifndef DOXYGEN
1144 
1145 template <>
1146 Point<1>
1148  const Triangulation<1, 1>::face_iterator &) const;
1149 
1150 template <>
1151 Point<2>
1153  const Triangulation<1, 2>::face_iterator &) const;
1154 
1155 
1156 template <>
1157 Point<3>
1159  const Triangulation<1, 3>::face_iterator &) const;
1160 
1161 
1162 template <>
1163 Point<1>
1165  const Triangulation<1, 1>::quad_iterator &) const;
1166 
1167 template <>
1168 Point<2>
1170  const Triangulation<1, 2>::quad_iterator &) const;
1171 
1172 
1173 template <>
1174 Point<3>
1176  const Triangulation<1, 3>::quad_iterator &) const;
1177 
1178 
1179 template <>
1180 Point<3>
1182  const Triangulation<3, 3>::hex_iterator &) const;
1183 
1184 /*---Templated functions---*/
1185 
1186 namespace Manifolds
1187 {
1188  template <typename MeshIteratorType>
1190  get_default_quadrature(const MeshIteratorType &iterator,
1191  const bool with_interpolation)
1192  {
1193  const auto points_and_weights =
1194  get_default_points_and_weights(iterator, with_interpolation);
1195  static const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1196  return Quadrature<spacedim>(
1197  std::vector<Point<spacedim>>(points_and_weights.first.begin(),
1198  points_and_weights.first.end()),
1199  std::vector<double>(points_and_weights.second.begin(),
1200  points_and_weights.second.end()));
1201  }
1202 
1203 
1204 
1205  template <typename MeshIteratorType>
1206  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
1207  n_default_points_per_cell<MeshIteratorType>()>,
1208  std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
1209  get_default_points_and_weights(const MeshIteratorType &iterator,
1210  const bool with_interpolation)
1211  {
1212  const int dim = MeshIteratorType::AccessorType::structure_dimension;
1213  const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1214  constexpr std::size_t points_per_cell =
1215  n_default_points_per_cell<MeshIteratorType>();
1216 
1217  std::pair<std::array<Point<spacedim>, points_per_cell>,
1218  std::array<double, points_per_cell>>
1219  points_weights;
1220 
1221 
1222  // note that the exact weights are chosen such as to minimize the
1223  // distortion of the four new quads from the optimal shape; their
1224  // derivation and values is copied over from the
1225  // interpolation function in the mapping
1226  switch (dim)
1227  {
1228  case 1:
1229  Assert(points_weights.first.size() == 2, ExcInternalError());
1230  Assert(points_weights.second.size() == 2, ExcInternalError());
1231  points_weights.first[0] = iterator->vertex(0);
1232  points_weights.second[0] = .5;
1233  points_weights.first[1] = iterator->vertex(1);
1234  points_weights.second[1] = .5;
1235  break;
1236  case 2:
1237  Assert(points_weights.first.size() == 8, ExcInternalError());
1238  Assert(points_weights.second.size() == 8, ExcInternalError());
1239 
1240  for (unsigned int i = 0; i < 4; ++i)
1241  {
1242  points_weights.first[i] = iterator->vertex(i);
1243  points_weights.first[4 + i] =
1244  (iterator->line(i)->has_children() ?
1245  iterator->line(i)->child(0)->vertex(1) :
1246  iterator->line(i)->get_manifold().get_new_point_on_line(
1247  iterator->line(i)));
1248  }
1249 
1250  if (with_interpolation)
1251  {
1252  std::fill(points_weights.second.begin(),
1253  points_weights.second.begin() + 4,
1254  -0.25);
1255  std::fill(points_weights.second.begin() + 4,
1256  points_weights.second.end(),
1257  0.5);
1258  }
1259  else
1260  std::fill(points_weights.second.begin(),
1261  points_weights.second.end(),
1262  1.0 / 8.0);
1263  break;
1264  case 3:
1265  {
1267  static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
1268  const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
1271  Assert(points_weights.first.size() == np, ExcInternalError());
1272  Assert(points_weights.second.size() == np, ExcInternalError());
1273  auto *sp3 = reinterpret_cast<
1274  std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
1275  *>(&points_weights.first);
1276 
1277  unsigned int j = 0;
1278 
1279  // note that the exact weights are chosen such as to minimize the
1280  // distortion of the eight new hexes from the optimal shape through
1281  // transfinite interpolation from the faces and vertices, see
1282  // TransfiniteInterpolationManifold for a deeper explanation of the
1283  // mechanisms
1284  if (with_interpolation)
1285  {
1286  for (unsigned int i = 0;
1287  i < GeometryInfo<dim>::vertices_per_cell;
1288  ++i, ++j)
1289  {
1290  (*sp3)[j] = hex->vertex(i);
1291  points_weights.second[j] = 1.0 / 8.0;
1292  }
1293  for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
1294  ++i, ++j)
1295  {
1296  (*sp3)[j] =
1297  (hex->line(i)->has_children() ?
1298  hex->line(i)->child(0)->vertex(1) :
1299  hex->line(i)->get_manifold().get_new_point_on_line(
1300  hex->line(i)));
1301  points_weights.second[j] = -1.0 / 4.0;
1302  }
1303  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1304  ++i, ++j)
1305  {
1306  (*sp3)[j] =
1307  (hex->quad(i)->has_children() ?
1308  hex->quad(i)->isotropic_child(0)->vertex(3) :
1309  hex->quad(i)->get_manifold().get_new_point_on_quad(
1310  hex->quad(i)));
1311  points_weights.second[j] = 1.0 / 2.0;
1312  }
1313  }
1314  else
1315  // Overwrite the weights with 1/np if we don't want to use
1316  // interpolation.
1317  std::fill(points_weights.second.begin(),
1318  points_weights.second.end(),
1319  1.0 / np);
1320  }
1321  break;
1322  default:
1323  Assert(false, ExcInternalError());
1324  break;
1325  }
1326  return points_weights;
1327  }
1328 } // namespace Manifolds
1329 
1330 #endif // DOXYGEN
1331 
1333 
1334 #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:1421
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:1469
constexpr std::size_t n_default_points_per_cell()
Definition: manifold.h:54
Quadrature< MeshIteratorType::AccessorType::space_dimension > get_default_quadrature(const MeshIteratorType &iterator, const bool with_interpolation=false)
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1445
const double tolerance
Definition: manifold.h:858
const FlatManifold< chartdim, chartdim > sub_manifold
Definition: manifold.h:1137
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:361
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:355
Definition: table.h:695
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:344
#define DEAL_II_DEPRECATED
Definition: config.h:151
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:844
static ::ExceptionBase & ExcInternalError()