Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+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\}}\)
Loading...
Searching...
No Matches
manifold.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2014 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_tria_manifold_h
16#define dealii_tria_manifold_h
17
18
19/*---------------------------- manifold.h ---------------------------*/
20
21#include <deal.II/base/config.h>
22
25#include <deal.II/base/point.h>
28
29#include <deal.II/grid/tria.h>
30
32
33// forward declaration
34#ifndef DOXYGEN
35template <int, typename>
36class Table;
37#endif
38
43namespace 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>
111 std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
112 n_default_points_per_cell<MeshIteratorType>()>,
113 std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
114 get_default_points_and_weights(const MeshIteratorType &iterator,
115 const bool with_interpolation = false);
116} // namespace Manifolds
117
118
119
284template <int dim, int spacedim = dim>
285class Manifold : public Subscriptor
286{
287public:
288 // explicitly check for sensible template arguments
289 static_assert(dim <= spacedim,
290 "The dimension <dim> of a Manifold must be less than or "
291 "equal to the space dimension <spacedim> in which it lives.");
292
293
306 std::array<Tensor<1, spacedim>, GeometryInfo<dim>::vertices_per_face>;
307
308
313 virtual ~Manifold() override = default;
314
320 virtual std::unique_ptr<Manifold<dim, spacedim>>
321 clone() const = 0;
322
327
345 virtual Point<spacedim>
347 const Point<spacedim> &p2,
348 const double w) const;
349
366 virtual Point<spacedim>
367 get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
368 const ArrayView<const double> &weights) const;
369
370
391 virtual void
392 get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
393 const Table<2, double> &weights,
394 ArrayView<Point<spacedim>> new_points) const;
395
407 virtual Point<spacedim>
409 const ArrayView<const Point<spacedim>> &surrounding_points,
410 const Point<spacedim> &candidate) const;
411
425 virtual Point<spacedim>
427 const typename Triangulation<dim, spacedim>::line_iterator &line) const;
428
446 virtual Point<spacedim>
448 const typename Triangulation<dim, spacedim>::quad_iterator &quad) const;
449
468 virtual Point<spacedim>
470 const typename Triangulation<dim, spacedim>::hex_iterator &hex) const;
471
472
481 const typename Triangulation<dim, spacedim>::face_iterator &face) const;
482
483
492 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const;
493
537 virtual Tensor<1, spacedim>
539 const Point<spacedim> &x2) const;
540
596 virtual Tensor<1, spacedim>
599 const Point<spacedim> &p) const;
600
615 virtual void
618 FaceVertexNormals &face_vertex_normals) const;
619
621};
622
623
633template <int dim, int spacedim = dim>
634class FlatManifold : public Manifold<dim, spacedim>
635{
636public:
665 const double tolerance = 1e-10);
666
670 virtual std::unique_ptr<Manifold<dim, spacedim>>
671 clone() const override;
672
694 virtual Point<spacedim>
695 get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
696 const ArrayView<const double> &weights) const override;
697
698
709 virtual void
710 get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
711 const Table<2, double> &weights,
712 ArrayView<Point<spacedim>> new_points) const override;
713
721 virtual Point<spacedim>
723 const Point<spacedim> &candidate) const override;
724
746 virtual Tensor<1, spacedim>
748 const Point<spacedim> &x2) const override;
749
757 virtual Tensor<1, spacedim>
760 const Point<spacedim> &p) const override;
761
770 virtual void
773 typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
774 const override;
775
779 const Tensor<1, spacedim> &
781
782private:
797
799 int,
801 double,
802 << "The component number " << arg1 << " of the point [ "
803 << arg2 << " ] is not in the interval [ 0, " << arg3
804 << "), bailing out.");
805
810 const double tolerance;
811};
812
813
900template <int dim, int spacedim = dim, int chartdim = dim>
901class ChartManifold : public Manifold<dim, spacedim>
902{
903public:
904 // explicitly check for sensible template arguments
905 static_assert(dim <= spacedim,
906 "The dimension <dim> of a ChartManifold must be less than or "
907 "equal to the space dimension <spacedim> in which it lives.");
908
924
929 virtual ~ChartManifold() override = default;
930
935 virtual Point<spacedim>
937 const Point<spacedim> &p2,
938 const double w) const override;
939
944 virtual Point<spacedim>
945 get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
946 const ArrayView<const double> &weights) const override;
947
969 virtual void
970 get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
971 const Table<2, double> &weights,
972 ArrayView<Point<spacedim>> new_points) const override;
979 virtual Point<chartdim>
980 pull_back(const Point<spacedim> &space_point) const = 0;
981
988 virtual Point<spacedim>
989 push_forward(const Point<chartdim> &chart_point) const = 0;
990
1008 push_forward_gradient(const Point<chartdim> &chart_point) const;
1009
1065 virtual Tensor<1, spacedim>
1067 const Point<spacedim> &x2) const override;
1068
1072 const Tensor<1, chartdim> &
1073 get_periodicity() const;
1074
1075private:
1088};
1089
1090
1091/* -------------- declaration of explicit specializations ------------- */
1092
1093#ifndef DOXYGEN
1094
1095template <>
1099
1100template <>
1104
1105
1106template <>
1110
1111
1112template <>
1116
1117template <>
1121
1122
1123template <>
1127
1128
1129template <>
1132 const Triangulation<3, 3>::hex_iterator &) const;
1133
1134/*---Templated functions---*/
1135
1136namespace Manifolds
1137{
1138 template <typename MeshIteratorType>
1139 std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
1140 n_default_points_per_cell<MeshIteratorType>()>,
1141 std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
1142 get_default_points_and_weights(const MeshIteratorType &iterator,
1143 const bool with_interpolation)
1144 {
1145 const int dim = MeshIteratorType::AccessorType::structure_dimension;
1146 const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1147 constexpr std::size_t points_per_cell =
1148 n_default_points_per_cell<MeshIteratorType>();
1149
1150 std::pair<std::array<Point<spacedim>, points_per_cell>,
1151 std::array<double, points_per_cell>>
1152 points_weights;
1153
1154
1155 // note that the exact weights are chosen such as to minimize the
1156 // distortion of the four new quads from the optimal shape; their
1157 // derivation and values is copied over from the
1158 // interpolation function in the mapping
1159 switch (dim)
1160 {
1161 case 1:
1162 Assert(points_weights.first.size() == 2, ExcInternalError());
1163 Assert(points_weights.second.size() == 2, ExcInternalError());
1164 points_weights.first[0] = iterator->vertex(0);
1165 points_weights.second[0] = .5;
1166 points_weights.first[1] = iterator->vertex(1);
1167 points_weights.second[1] = .5;
1168 break;
1169 case 2:
1170 Assert(points_weights.first.size() == 8, ExcInternalError());
1171 Assert(points_weights.second.size() == 8, ExcInternalError());
1172
1173 for (unsigned int i = 0; i < 4; ++i)
1174 {
1175 points_weights.first[i] = iterator->vertex(i);
1176 points_weights.first[4 + i] =
1177 (iterator->line(i)->has_children() ?
1178 iterator->line(i)->child(0)->vertex(1) :
1179 iterator->line(i)->get_manifold().get_new_point_on_line(
1180 iterator->line(i)));
1181 }
1182
1183 if (with_interpolation)
1184 {
1185 std::fill(points_weights.second.begin(),
1186 points_weights.second.begin() + 4,
1187 -0.25);
1188 std::fill(points_weights.second.begin() + 4,
1189 points_weights.second.end(),
1190 0.5);
1191 }
1192 else
1193 std::fill(points_weights.second.begin(),
1194 points_weights.second.end(),
1195 1.0 / 8.0);
1196 break;
1197 case 3:
1198 {
1200 static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
1201 const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
1204 Assert(points_weights.first.size() == np, ExcInternalError());
1205 Assert(points_weights.second.size() == np, ExcInternalError());
1206 auto *sp3 = reinterpret_cast<
1207 std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
1208 *>(&points_weights.first);
1209
1210 unsigned int j = 0;
1211
1212 // note that the exact weights are chosen such as to minimize the
1213 // distortion of the eight new hexes from the optimal shape through
1214 // transfinite interpolation from the faces and vertices, see
1215 // TransfiniteInterpolationManifold for a deeper explanation of the
1216 // mechanisms
1217 if (with_interpolation)
1218 {
1219 for (unsigned int i = 0;
1220 i < GeometryInfo<dim>::vertices_per_cell;
1221 ++i, ++j)
1222 {
1223 (*sp3)[j] = hex->vertex(i);
1224 points_weights.second[j] = 1.0 / 8.0;
1225 }
1226 for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
1227 ++i, ++j)
1228 {
1229 (*sp3)[j] =
1230 (hex->line(i)->has_children() ?
1231 hex->line(i)->child(0)->vertex(1) :
1232 hex->line(i)->get_manifold().get_new_point_on_line(
1233 hex->line(i)));
1234 points_weights.second[j] = -1.0 / 4.0;
1235 }
1236 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1237 ++i, ++j)
1238 {
1239 (*sp3)[j] =
1240 (hex->quad(i)->has_children() ?
1241 hex->quad(i)->isotropic_child(0)->vertex(3) :
1242 hex->quad(i)->get_manifold().get_new_point_on_quad(
1243 hex->quad(i)));
1244 points_weights.second[j] = 1.0 / 2.0;
1245 }
1246 }
1247 else
1248 // Overwrite the weights with 1/np if we don't want to use
1249 // interpolation.
1250 std::fill(points_weights.second.begin(),
1251 points_weights.second.end(),
1252 1.0 / np);
1253 }
1254 break;
1255 default:
1257 break;
1258 }
1259 return points_weights;
1260 }
1261} // namespace Manifolds
1262
1263#endif // DOXYGEN
1264
1266
1267#endif
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const =0
const FlatManifold< chartdim, chartdim > sub_manifold
Definition manifold.h:1087
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:1048
virtual Point< chartdim > pull_back(const Point< spacedim > &space_point) const =0
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition manifold.cc:1087
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition manifold.cc:1074
const Tensor< 1, chartdim > & get_periodicity() const
Definition manifold.cc:1119
virtual ~ChartManifold() override=default
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition manifold.cc:1013
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
Definition manifold.cc:1027
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) 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 Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &points, const Point< spacedim > &candidate) const override
const Tensor< 1, spacedim > & get_periodicity() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
const Tensor< 1, spacedim > periodicity
Definition manifold.h:796
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
const double tolerance
Definition manifold.h:810
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 Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &surrounding_points, const Point< spacedim > &candidate) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
virtual ~Manifold() override=default
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition manifold.h:306
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const =0
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:562
typename IteratorSelector::hex_iterator hex_iterator
Definition tria.h:1705
typename IteratorSelector::quad_iterator quad_iterator
Definition tria.h:1681
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1657
#define DEAL_II_ASSERT_UNREACHABLE()
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:53