Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 1998 - 2023 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
36template <int, typename>
37class Table;
38#endif
39
44namespace 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
285template <int dim, int spacedim = dim>
286class Manifold : public Subscriptor
287{
288public:
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>
600 const Point<spacedim> & p) const;
601
616 virtual void
619 FaceVertexNormals &face_vertex_normals) const;
620
622};
623
624
634template <int dim, int spacedim = dim>
635class FlatManifold : public Manifold<dim, spacedim>
636{
637public:
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>
761 const Point<spacedim> &p) const override;
762
771 virtual void
774 typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
775 const override;
776
780 const Tensor<1, spacedim> &
782
783private:
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
901template <int dim, int spacedim = dim, int chartdim = dim>
902class ChartManifold : public Manifold<dim, spacedim>
903{
904public:
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
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
1076private:
1089};
1090
1091
1092/* -------------- declaration of explicit specializations ------------- */
1093
1094#ifndef DOXYGEN
1095
1096template <>
1100
1101template <>
1105
1106
1107template <>
1111
1112
1113template <>
1117
1118template <>
1122
1123
1124template <>
1128
1129
1130template <>
1133 const Triangulation<3, 3>::hex_iterator &) const;
1134
1135/*---Templated functions---*/
1136
1137namespace 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 > push_forward(const Point< chartdim > &chart_point) const =0
const FlatManifold< chartdim, chartdim > sub_manifold
Definition manifold.h:1088
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:1053
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:1095
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition manifold.cc:1082
const Tensor< 1, chartdim > & get_periodicity() const
Definition manifold.cc:1127
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:1017
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
Definition manifold.cc:1032
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:797
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
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 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:307
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:112
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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:556
typename IteratorSelector::hex_iterator hex_iterator
Definition tria.h:1505
typename IteratorSelector::quad_iterator quad_iterator
Definition tria.h:1481
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1457
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