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