deal.II version GIT relicensing-1929-g70a5450eef 2024-10-03 18:00:00+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
785 double
787
788private:
803
805 int,
807 double,
808 << "The component number " << arg1 << " of the point [ "
809 << arg2 << " ] is not in the interval [ 0, " << arg3
810 << "), bailing out.");
811
816 const double tolerance;
817};
818
819
906template <int dim, int spacedim = dim, int chartdim = dim>
907class ChartManifold : public Manifold<dim, spacedim>
908{
909public:
910 // explicitly check for sensible template arguments
911 static_assert(dim <= spacedim,
912 "The dimension <dim> of a ChartManifold must be less than or "
913 "equal to the space dimension <spacedim> in which it lives.");
914
930
935 virtual ~ChartManifold() override = default;
936
941 virtual Point<spacedim>
943 const Point<spacedim> &p2,
944 const double w) const override;
945
950 virtual Point<spacedim>
951 get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
952 const ArrayView<const double> &weights) const override;
953
975 virtual void
976 get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
977 const Table<2, double> &weights,
978 ArrayView<Point<spacedim>> new_points) const override;
985 virtual Point<chartdim>
986 pull_back(const Point<spacedim> &space_point) const = 0;
987
994 virtual Point<spacedim>
995 push_forward(const Point<chartdim> &chart_point) const = 0;
996
1014 push_forward_gradient(const Point<chartdim> &chart_point) const;
1015
1071 virtual Tensor<1, spacedim>
1073 const Point<spacedim> &x2) const override;
1074
1078 const Tensor<1, chartdim> &
1079 get_periodicity() const;
1080
1081private:
1094};
1095
1096
1097/* -------------- declaration of explicit specializations ------------- */
1098
1099#ifndef DOXYGEN
1100
1101template <>
1105
1106template <>
1110
1111
1112template <>
1116
1117
1118template <>
1122
1123template <>
1127
1128
1129template <>
1133
1134
1135template <>
1138 const Triangulation<3, 3>::hex_iterator &) const;
1139
1140/*---Templated functions---*/
1141
1142template <int dim, int spacedim>
1143inline double
1145{
1146 return tolerance;
1147}
1148
1149
1150
1151namespace Manifolds
1152{
1153 template <typename MeshIteratorType>
1154 std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
1155 n_default_points_per_cell<MeshIteratorType>()>,
1156 std::array<double, n_default_points_per_cell<MeshIteratorType>()>>
1157 get_default_points_and_weights(const MeshIteratorType &iterator,
1158 const bool with_interpolation)
1159 {
1160 const int dim = MeshIteratorType::AccessorType::structure_dimension;
1161 const int spacedim = MeshIteratorType::AccessorType::space_dimension;
1162 constexpr std::size_t points_per_cell =
1163 n_default_points_per_cell<MeshIteratorType>();
1164
1165 std::pair<std::array<Point<spacedim>, points_per_cell>,
1166 std::array<double, points_per_cell>>
1167 points_weights;
1168
1169
1170 // note that the exact weights are chosen such as to minimize the
1171 // distortion of the four new quads from the optimal shape; their
1172 // derivation and values is copied over from the
1173 // interpolation function in the mapping
1174 switch (dim)
1175 {
1176 case 1:
1177 Assert(points_weights.first.size() == 2, ExcInternalError());
1178 Assert(points_weights.second.size() == 2, ExcInternalError());
1179 points_weights.first[0] = iterator->vertex(0);
1180 points_weights.second[0] = .5;
1181 points_weights.first[1] = iterator->vertex(1);
1182 points_weights.second[1] = .5;
1183 break;
1184 case 2:
1185 Assert(points_weights.first.size() == 8, ExcInternalError());
1186 Assert(points_weights.second.size() == 8, ExcInternalError());
1187
1188 for (unsigned int i = 0; i < 4; ++i)
1189 {
1190 points_weights.first[i] = iterator->vertex(i);
1191 points_weights.first[4 + i] =
1192 (iterator->line(i)->has_children() ?
1193 iterator->line(i)->child(0)->vertex(1) :
1194 iterator->line(i)->get_manifold().get_new_point_on_line(
1195 iterator->line(i)));
1196 }
1197
1198 if (with_interpolation)
1199 {
1200 std::fill(points_weights.second.begin(),
1201 points_weights.second.begin() + 4,
1202 -0.25);
1203 std::fill(points_weights.second.begin() + 4,
1204 points_weights.second.end(),
1205 0.5);
1206 }
1207 else
1208 std::fill(points_weights.second.begin(),
1209 points_weights.second.end(),
1210 1.0 / 8.0);
1211 break;
1212 case 3:
1213 {
1215 static_cast<TriaIterator<TriaAccessor<3, 3, 3>>>(iterator);
1216 const unsigned int np = GeometryInfo<dim>::vertices_per_cell +
1219 Assert(points_weights.first.size() == np, ExcInternalError());
1220 Assert(points_weights.second.size() == np, ExcInternalError());
1221 auto *sp3 = reinterpret_cast<
1222 std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()>
1223 *>(&points_weights.first);
1224
1225 unsigned int j = 0;
1226
1227 // note that the exact weights are chosen such as to minimize the
1228 // distortion of the eight new hexes from the optimal shape through
1229 // transfinite interpolation from the faces and vertices, see
1230 // TransfiniteInterpolationManifold for a deeper explanation of the
1231 // mechanisms
1232 if (with_interpolation)
1233 {
1234 for (unsigned int i = 0;
1235 i < GeometryInfo<dim>::vertices_per_cell;
1236 ++i, ++j)
1237 {
1238 (*sp3)[j] = hex->vertex(i);
1239 points_weights.second[j] = 1.0 / 8.0;
1240 }
1241 for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
1242 ++i, ++j)
1243 {
1244 (*sp3)[j] =
1245 (hex->line(i)->has_children() ?
1246 hex->line(i)->child(0)->vertex(1) :
1247 hex->line(i)->get_manifold().get_new_point_on_line(
1248 hex->line(i)));
1249 points_weights.second[j] = -1.0 / 4.0;
1250 }
1251 for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell;
1252 ++i, ++j)
1253 {
1254 (*sp3)[j] =
1255 (hex->quad(i)->has_children() ?
1256 hex->quad(i)->isotropic_child(0)->vertex(3) :
1257 hex->quad(i)->get_manifold().get_new_point_on_quad(
1258 hex->quad(i)));
1259 points_weights.second[j] = 1.0 / 2.0;
1260 }
1261 }
1262 else
1263 // Overwrite the weights with 1/np if we don't want to use
1264 // interpolation.
1265 std::fill(points_weights.second.begin(),
1266 points_weights.second.end(),
1267 1.0 / np);
1268 }
1269 break;
1270 default:
1272 break;
1273 }
1274 return points_weights;
1275 }
1276} // namespace Manifolds
1277
1278#endif // DOXYGEN
1279
1281
1282#endif
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const =0
const FlatManifold< chartdim, chartdim > sub_manifold
Definition manifold.h:1093
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
double get_tolerance() const
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:802
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
const double tolerance
Definition manifold.h:816
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:557
typename IteratorSelector::hex_iterator hex_iterator
Definition tria.h:1691
typename IteratorSelector::quad_iterator quad_iterator
Definition tria.h:1667
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1643
#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