Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
mapping_c1.cc
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 #include <deal.II/base/std_cxx14/memory.h>
17 
18 #include <deal.II/fe/mapping_c1.h>
19 
20 #include <deal.II/grid/manifold.h>
21 #include <deal.II/grid/tria_accessor.h>
22 #include <deal.II/grid/tria_iterator.h>
23 
24 #include <cmath>
25 
26 DEAL_II_NAMESPACE_OPEN
27 
28 
29 
30 template <int dim, int spacedim>
32  : MappingQGeneric<dim, spacedim>(3)
33 {}
34 
35 
36 
37 template <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 
54 template <>
55 void
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
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
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
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
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
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>>
224 {
225  return std_cxx14::make_unique<MappingC1<dim, spacedim>>();
226 }
227 
228 
229 
230 // explicit instantiations
231 #include "mapping_c1.inst"
232 
233 
234 DEAL_II_NAMESPACE_CLOSE
std::shared_ptr< const MappingQGeneric< dim, spacedim > > qp_mapping
Definition: mapping_q.h:368
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_c1.cc:223
Definition: point.h:110
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
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
#define Assert(cond, exc)
Definition: exceptions.h:1407
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
Definition: tensor.h:422
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:351
static ::ExceptionBase & ExcNotImplemented()
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:301
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
Definition: tria.cc:10398