Reference documentation for deal.II version GIT b065060f03 2023-05-30 23:50:02+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\}}\)
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
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>
22 
23 #include <cmath>
24 #include <memory>
25 
27 
28 
29 
30 template <int dim, int spacedim>
32  : MappingQ<dim, spacedim>(3)
33 {
34  Assert(dim > 1, ExcImpossibleInDim(dim));
35 }
36 
37 
38 
39 template <>
40 void
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 
51 template <>
52 void
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 
157 template <int dim, int spacedim>
158 void
160  const typename Triangulation<dim>::cell_iterator &,
161  std::vector<Point<dim>> &) const
162 {
163  Assert(false, ExcNotImplemented());
164 }
165 
166 
167 
168 template <>
169 void
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 
180 template <>
181 void
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 
192 template <int dim, int spacedim>
193 void
195  const typename Triangulation<dim>::cell_iterator &,
196  std::vector<Point<dim>> &) const
197 {
198  Assert(false, ExcNotImplemented());
199 }
200 
201 
202 
203 template <int dim, int spacedim>
204 std::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
Definition: manifold.cc:338
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:307
virtual std::unique_ptr< Mapping< dim, spacedim > > clone() const override
Definition: mapping_c1.cc:205
virtual void add_quad_support_points(const typename Triangulation< dim >::cell_iterator &cell, std::vector< Point< dim >> &a) const override
Definition: mapping_c1.cc:194
virtual void add_line_support_points(const typename Triangulation< dim >::cell_iterator &cell, std::vector< Point< dim >> &a) const override
Definition: mapping_c1.cc:159
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1458
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)