1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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
17 #include <deal.II/fe/mapping_c1.h>
18
19 #include <deal.II/grid/manifold.h>
20 #include <deal.II/grid/tria_accessor.h>
21 #include <deal.II/grid/tria_iterator.h>
22
23 #include <cmath>
24 #include <memory>
25
26 DEAL_II_NAMESPACE_OPEN
27
28
29
30 template <int dim, int spacedim>
MappingC1Generic()31 MappingC1<dim, spacedim>::MappingC1Generic::MappingC1Generic()
32 : MappingQGeneric<dim, spacedim>(3)
33 {}
34
35
36
37 template <int dim, int spacedim>
MappingC1()38 MappingC1<dim, spacedim>::MappingC1()
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
54 template <>
55 void
add_line_support_points(const Triangulation<1>::cell_iterator &,std::vector<Point<1>> &) const56 MappingC1<1>::MappingC1Generic::add_line_support_points(
57 const Triangulation<1>::cell_iterator &,
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
67 template <>
68 void
add_line_support_points(const Triangulation<2>::cell_iterator & cell,std::vector<Point<2>> & a) const69 MappingC1<2>::MappingC1Generic::add_line_support_points(
70 const Triangulation<2>::cell_iterator &cell,
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
173 template <int dim, int spacedim>
174 void
add_line_support_points(const typename Triangulation<dim>::cell_iterator &,std::vector<Point<dim>> &) const175 MappingC1<dim, spacedim>::MappingC1Generic::add_line_support_points(
176 const typename Triangulation<dim>::cell_iterator &,
177 std::vector<Point<dim>> &) const
178 {
179 Assert(false, ExcNotImplemented());
180 }
181
182
183
184 template <>
185 void
add_quad_support_points(const Triangulation<1>::cell_iterator &,std::vector<Point<1>> &) const186 MappingC1<1>::MappingC1Generic::add_quad_support_points(
187 const Triangulation<1>::cell_iterator &,
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
197 template <>
198 void
add_quad_support_points(const Triangulation<2>::cell_iterator &,std::vector<Point<2>> &) const199 MappingC1<2>::MappingC1Generic::add_quad_support_points(
200 const Triangulation<2>::cell_iterator &,
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
210 template <int dim, int spacedim>
211 void
add_quad_support_points(const typename Triangulation<dim>::cell_iterator &,std::vector<Point<dim>> &) const212 MappingC1<dim, spacedim>::MappingC1Generic::add_quad_support_points(
213 const typename Triangulation<dim>::cell_iterator &,
214 std::vector<Point<dim>> &) const
215 {
216 Assert(false, ExcNotImplemented());
217 }
218
219
220
221 template <int dim, int spacedim>
222 std::unique_ptr<Mapping<dim, spacedim>>
clone() const223 MappingC1<dim, spacedim>::clone() const
224 {
225 return std::make_unique<MappingC1<dim, spacedim>>();
226 }
227
228
229
230 // explicit instantiations
231 #include "mapping_c1.inst"
232
233
234 DEAL_II_NAMESPACE_CLOSE
235