Reference documentation for deal.II version 9.5.0
Searching...
No Matches
mapping_c1.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2001 - 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
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
18
22
23#include <cmath>
24#include <memory>
25
27
28
29
30template <int dim, int spacedim>
32 : MappingQ<dim, spacedim>(3)
33{
34 Assert(dim > 1, ExcImpossibleInDim(dim));
35}
36
37
38
39template <>
40void
42 std::vector<Point<1>> &) const
43{
44 const unsigned int dim = 1;
45 (void)dim;
46 Assert(dim > 1, ExcImpossibleInDim(dim));
47}
48
49
50
51template <>
52void
55 std::vector<Point<2>> & a) const
56{
57 const unsigned int dim = 2;
58 const std::array<double, 2> interior_gl_points{
59 {0.5 - 0.5 * std::sqrt(1.0 / 5.0), 0.5 + 0.5 * std::sqrt(1.0 / 5.0)}};
60
61 // loop over each of the lines, and if it is at the boundary, then first get
62 // the boundary description and second compute the points on it. if not at
63 // the boundary, get the respective points from another function
64 for (unsigned int line_no = 0; line_no < GeometryInfo<dim>::lines_per_cell;
65 ++line_no)
66 {
67 const Triangulation<dim>::line_iterator line = cell->line(line_no);
68
69 if (line->at_boundary())
70 {
71 // first get the normal vectors at the two vertices of this line
72 // from the boundary description
73 const Manifold<dim> &manifold = line->get_manifold();
74
75 Manifold<dim>::FaceVertexNormals face_vertex_normals;
76 manifold.get_normals_at_vertices(line, face_vertex_normals);
77
78 // then transform them into interpolation points for a cubic
79 // polynomial
80 //
81 // for this, note that if we describe the boundary curve as a
82 // polynomial in tangential coordinate @p{t=0..1} (along the line)
83 // and @p{s} in normal direction, then the cubic mapping is such
84 // that @p{s = a*t**3 + b*t**2 + c*t + d}, and we want to determine
85 // the interpolation points at @p{t=0.276} and @p{t=0.724}
86 // (Gauss-Lobatto points). Since at @p{t=0,1} we want a vertex which
87 // is actually at the boundary, we know that @p{d=0} and @p{a=-b-c},
88 // which gives @p{s(0.276)} and @p{s(0.726)} in terms of @p{b,c}. As
89 // side-conditions, we want that the derivatives at @p{t=0} and
90 // @p{t=1}, i.e. at the vertices match those returned by the
91 // boundary.
92 //
93 // The task is then first to determine the coefficients from the
94 // tangentials. for that, first rotate the tangents of @p{s(t)} into
95 // the global coordinate system. they are @p{A (1,c)} and @p{A
96 // (1,-b-2c)} with @p{A} the rotation matrix, since the tangentials
97 // in the coordinate system relative to the line are @p{(1,c)} and
98 // @p{(1,-b-2c)} at the two vertices, respectively. We then have to
99 // make sure by matching @p{b,c} that these tangentials are
100 // orthogonal to the normals returned by the boundary object
101 const Tensor<1, 2> coordinate_vector =
102 line->vertex(1) - line->vertex(0);
103 const double h = std::sqrt(coordinate_vector * coordinate_vector);
104 Tensor<1, 2> coordinate_axis = coordinate_vector;
105 coordinate_axis /= h;
106
107 const double alpha =
108 std::atan2(coordinate_axis[1], coordinate_axis[0]);
109 const double c = -((face_vertex_normals[0][1] * std::sin(alpha) +
110 face_vertex_normals[0][0] * std::cos(alpha)) /
111 (face_vertex_normals[0][1] * std::cos(alpha) -
112 face_vertex_normals[0][0] * std::sin(alpha)));
113 const double b = ((face_vertex_normals[1][1] * std::sin(alpha) +
114 face_vertex_normals[1][0] * std::cos(alpha)) /
115 (face_vertex_normals[1][1] * std::cos(alpha) -
116 face_vertex_normals[1][0] * std::sin(alpha))) -
117 2 * c;
118
119 const double t1 = interior_gl_points[0];
120 const double t2 = interior_gl_points[1];
121 const double s_t1 = (((-b - c) * t1 + b) * t1 + c) * t1;
122 const double s_t2 = (((-b - c) * t2 + b) * t2 + c) * t2;
123
124 // next evaluate the so determined cubic polynomial at the points
125 // 1/3 and 2/3, first in unit coordinates
126 const Point<2> new_unit_points[2] = {Point<2>(t1, s_t1),
127 Point<2>(t2, s_t2)};
128 // then transform these points to real coordinates by rotating,
129 // scaling and shifting
130 for (const auto &new_unit_point : new_unit_points)
131 {
132 Point<2> real_point(std::cos(alpha) * new_unit_point[0] -
133 std::sin(alpha) * new_unit_point[1],
134 std::sin(alpha) * new_unit_point[0] +
135 std::cos(alpha) * new_unit_point[1]);
136 real_point *= h;
137 real_point += line->vertex(0);
138 a.push_back(real_point);
139 }
140 }
141 else
142 // not at boundary, so just use scaled Gauss-Lobatto points (i.e., use
143 // plain straight lines).
144 {
145 // Note that the zeroth Gauss-Lobatto point is a boundary point, so
146 // we push back mapped versions of the first and second.
147 a.push_back((1.0 - interior_gl_points[0]) * line->vertex(0) +
148 (interior_gl_points[0] * line->vertex(1)));
149 a.push_back((1.0 - interior_gl_points[1]) * line->vertex(0) +
150 (interior_gl_points[1] * line->vertex(1)));
151 }
152 }
153}
154
155
156
157template <int dim, int spacedim>
158void
160 const typename Triangulation<dim>::cell_iterator &,
161 std::vector<Point<dim>> &) const
162{
163 Assert(false, ExcNotImplemented());
164}
165
166
167
168template <>
169void
171 std::vector<Point<1>> &) const
172{
173 const unsigned int dim = 1;
174 (void)dim;
175 Assert(dim > 2, ExcImpossibleInDim(dim));
176}
177
178
179
180template <>
181void
183 std::vector<Point<2>> &) const
184{
185 const unsigned int dim = 2;
186 (void)dim;
187 Assert(dim > 2, ExcImpossibleInDim(dim));
188}
189
190
191
192template <int dim, int spacedim>
193void
195 const typename Triangulation<dim>::cell_iterator &,
196 std::vector<Point<dim>> &) const
197{
198 Assert(false, ExcNotImplemented());
199}
200
201
202
203template <int dim, int spacedim>
204std::unique_ptr<Mapping<dim, spacedim>>
206{
207 return std::make_unique<MappingC1<dim, spacedim>>();
208}
209
210
211
212// explicit instantiations
213#include "mapping_c1.inst"
214
215
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition manifold.h:307
virtual void add_line_support_points(const typename Triangulation< dim >::cell_iterator &cell, std::vector< Point< dim > > &a) const override
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
virtual void add_quad_support_points(const typename Triangulation< dim >::cell_iterator &cell, std::vector< Point< dim > > &a) const override
Definition point.h:112
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1457
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)