Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
polynomials_raviart_thomas.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2004 - 2022 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
20
21#include <iomanip>
22#include <iostream>
23#include <memory>
24
25
27
28
29namespace
30{
31 // Create nodal Raviart-Thomas polynomials as the tensor product of Lagrange
32 // polynomials on Gauss-Lobatto points of the given degrees in the normal and
33 // tangential directions, respectively (we could also choose Lagrange
34 // polynomials on Gauss points but those are slightly more expensive to handle
35 // in classes).
36 std::vector<std::vector<Polynomials::Polynomial<double>>>
37 create_rt_polynomials(const unsigned int dim,
38 const unsigned int normal_degree,
39 const unsigned int tangential_degree)
40 {
41 std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
42 if (normal_degree > 0)
44 QGaussLobatto<1>(normal_degree + 1).get_points());
45 else
47 QMidpoint<1>().get_points());
48 if (tangential_degree > 0)
49 for (unsigned int d = 1; d < dim; ++d)
51 QGaussLobatto<1>(tangential_degree + 1).get_points());
52 else
53 for (unsigned int d = 1; d < dim; ++d)
55 QMidpoint<1>().get_points());
56
57 return pols;
58 }
59} // namespace
60
61
62
63template <int dim>
65 const unsigned int normal_degree,
66 const unsigned int tangential_degree)
67 : TensorPolynomialsBase<dim>(std::min(normal_degree, tangential_degree),
68 n_polynomials(normal_degree, tangential_degree))
69 , normal_degree(normal_degree)
70 , tangential_degree(tangential_degree)
71 , polynomial_space(
72 create_rt_polynomials(dim, normal_degree, tangential_degree))
73{
74 // create renumbering of the unknowns from the lexicographic order to the
75 // actual order required by the finite element class with unknowns on
76 // faces placed first
77 const unsigned int n_pols = polynomial_space.n();
80
83
84 // since we only store an anisotropic polynomial for the first component,
85 // we set up a second numbering to switch out the actual coordinate
86 // direction
87 renumber_aniso[0].resize(n_pols);
88 for (unsigned int i = 0; i < n_pols; ++i)
89 renumber_aniso[0][i] = i;
90 if (dim == 2)
91 {
92 // switch x and y component (i and j loops)
93 renumber_aniso[1].resize(n_pols);
94 for (unsigned int j = 0; j < normal_degree + 1; ++j)
95 for (unsigned int i = 0; i < tangential_degree + 1; ++i)
96 renumber_aniso[1][j * (tangential_degree + 1) + i] =
97 j + i * (normal_degree + 1);
98 }
99 if (dim == 3)
100 {
101 // switch x, y, and z component (i, j, k) -> (j, k, i)
102 renumber_aniso[1].resize(n_pols);
103 for (unsigned int k = 0; k < tangential_degree + 1; ++k)
104 for (unsigned int j = 0; j < normal_degree + 1; ++j)
105 for (unsigned int i = 0; i < tangential_degree + 1; ++i)
106 renumber_aniso[1][(k * (normal_degree + 1) + j) *
107 (tangential_degree + 1) +
108 i] =
109 j + (normal_degree + 1) * (k + i * (tangential_degree + 1));
110
111 // switch x, y, and z component (i, j, k) -> (k, i, j)
112 renumber_aniso[2].resize(n_pols);
113 for (unsigned int k = 0; k < normal_degree + 1; ++k)
114 for (unsigned int j = 0; j < tangential_degree + 1; ++j)
115 for (unsigned int i = 0; i < tangential_degree + 1; ++i)
116 renumber_aniso[2][(k * (tangential_degree + 1) + j) *
117 (tangential_degree + 1) +
118 i] =
119 k + (normal_degree + 1) * (i + j * (tangential_degree + 1));
120 }
121}
122
123
124
125template <int dim>
127 : PolynomialsRaviartThomas(k + 1, k)
128{}
129
130
131
132template <int dim>
133void
135 const Point<dim> & unit_point,
136 std::vector<Tensor<1, dim>> &values,
137 std::vector<Tensor<2, dim>> &grads,
138 std::vector<Tensor<3, dim>> &grad_grads,
139 std::vector<Tensor<4, dim>> &third_derivatives,
140 std::vector<Tensor<5, dim>> &fourth_derivatives) const
141{
142 Assert(values.size() == this->n() || values.size() == 0,
143 ExcDimensionMismatch(values.size(), this->n()));
144 Assert(grads.size() == this->n() || grads.size() == 0,
145 ExcDimensionMismatch(grads.size(), this->n()));
146 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
147 ExcDimensionMismatch(grad_grads.size(), this->n()));
148 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
149 ExcDimensionMismatch(third_derivatives.size(), this->n()));
150 Assert(fourth_derivatives.size() == this->n() ||
151 fourth_derivatives.size() == 0,
152 ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
153
154 std::vector<double> p_values;
155 std::vector<Tensor<1, dim>> p_grads;
156 std::vector<Tensor<2, dim>> p_grad_grads;
157 std::vector<Tensor<3, dim>> p_third_derivatives;
158 std::vector<Tensor<4, dim>> p_fourth_derivatives;
159
160 const unsigned int n_sub = polynomial_space.n();
161 p_values.resize((values.size() == 0) ? 0 : n_sub);
162 p_grads.resize((grads.size() == 0) ? 0 : n_sub);
163 p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
164 p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
165 p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
166
167 for (unsigned int d = 0; d < dim; ++d)
168 {
169 // First we copy the point. The polynomial space for component d
170 // consists of polynomials of degree k in x_d and degree k+1 in the
171 // other variables. in order to simplify this, we use the same
172 // AnisotropicPolynomial space and simply rotate the coordinates
173 // through all directions.
174 Point<dim> p;
175 for (unsigned int c = 0; c < dim; ++c)
176 p(c) = unit_point((c + d) % dim);
177
178 polynomial_space.evaluate(p,
179 p_values,
180 p_grads,
181 p_grad_grads,
182 p_third_derivatives,
183 p_fourth_derivatives);
184
185 for (unsigned int i = 0; i < p_values.size(); ++i)
186 values[lexicographic_to_hierarchic[i + d * n_sub]][d] =
187 p_values[renumber_aniso[d][i]];
188
189 for (unsigned int i = 0; i < p_grads.size(); ++i)
190 for (unsigned int d1 = 0; d1 < dim; ++d1)
191 grads[lexicographic_to_hierarchic[i + d * n_sub]][d][(d1 + d) % dim] =
192 p_grads[renumber_aniso[d][i]][d1];
193
194 for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
195 for (unsigned int d1 = 0; d1 < dim; ++d1)
196 for (unsigned int d2 = 0; d2 < dim; ++d2)
197 grad_grads[lexicographic_to_hierarchic[i + d * n_sub]][d]
198 [(d1 + d) % dim][(d2 + d) % dim] =
199 p_grad_grads[renumber_aniso[d][i]][d1][d2];
200
201 for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
202 for (unsigned int d1 = 0; d1 < dim; ++d1)
203 for (unsigned int d2 = 0; d2 < dim; ++d2)
204 for (unsigned int d3 = 0; d3 < dim; ++d3)
205 third_derivatives[lexicographic_to_hierarchic[i + d * n_sub]][d]
206 [(d1 + d) % dim][(d2 + d) % dim]
207 [(d3 + d) % dim] =
208 p_third_derivatives[renumber_aniso[d][i]][d1]
209 [d2][d3];
210
211 for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
212 for (unsigned int d1 = 0; d1 < dim; ++d1)
213 for (unsigned int d2 = 0; d2 < dim; ++d2)
214 for (unsigned int d3 = 0; d3 < dim; ++d3)
215 for (unsigned int d4 = 0; d4 < dim; ++d4)
216 fourth_derivatives[lexicographic_to_hierarchic[i + d * n_sub]]
217 [d][(d1 + d) % dim][(d2 + d) % dim]
218 [(d3 + d) % dim][(d4 + d) % dim] =
219 p_fourth_derivatives[renumber_aniso[d][i]]
220 [d1][d2][d3][d4];
221 }
222}
223
224
225
226template <int dim>
227std::string
229{
230 return "RaviartThomas";
231}
232
233
234
235template <int dim>
236unsigned int
238{
239 return n_polynomials(degree + 1, degree);
240}
241
242
243
244template <int dim>
245unsigned int
247 const unsigned int normal_degree,
248 const unsigned int tangential_degree)
249{
250 return dim * (normal_degree + 1) *
251 Utilities::pow(tangential_degree + 1, dim - 1);
252}
253
254
255
256template <int dim>
257std::vector<unsigned int>
259 const unsigned int normal_degree,
260 const unsigned int tangential_degree)
261{
262 const unsigned int n_dofs_face =
263 Utilities::pow(tangential_degree + 1, dim - 1);
264 std::vector<unsigned int> lexicographic_numbering;
265 // component 1
266 for (unsigned int j = 0; j < n_dofs_face; ++j)
267 {
268 lexicographic_numbering.push_back(j);
269 if (normal_degree > 1)
270 for (unsigned int i = n_dofs_face * 2 * dim;
271 i < n_dofs_face * 2 * dim + normal_degree - 1;
272 ++i)
273 lexicographic_numbering.push_back(i + j * (normal_degree - 1));
274 lexicographic_numbering.push_back(n_dofs_face + j);
275 }
276
277 // component 2
278 unsigned int layers = (dim == 3) ? tangential_degree + 1 : 1;
279 for (unsigned int k = 0; k < layers; ++k)
280 {
281 unsigned int k_add = k * (tangential_degree + 1);
282 for (unsigned int j = n_dofs_face * 2;
283 j < n_dofs_face * 2 + tangential_degree + 1;
284 ++j)
285 lexicographic_numbering.push_back(j + k_add);
286
287 if (normal_degree > 1)
288 for (unsigned int i = n_dofs_face * (2 * dim + (normal_degree - 1));
289 i < n_dofs_face * (2 * dim + (normal_degree - 1)) +
290 (normal_degree - 1) * (tangential_degree + 1);
291 ++i)
292 {
293 lexicographic_numbering.push_back(i + k_add * tangential_degree);
294 }
295 for (unsigned int j = n_dofs_face * 3;
296 j < n_dofs_face * 3 + tangential_degree + 1;
297 ++j)
298 lexicographic_numbering.push_back(j + k_add);
299 }
300
301 // component 3
302 if (dim == 3)
303 {
304 for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; ++i)
305 lexicographic_numbering.push_back(i);
306 if (normal_degree > 1)
307 for (unsigned int i =
308 6 * n_dofs_face + n_dofs_face * 2 * (normal_degree - 1);
309 i < 6 * n_dofs_face + n_dofs_face * 3 * (normal_degree - 1);
310 ++i)
311 lexicographic_numbering.push_back(i);
312 for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; ++i)
313 lexicographic_numbering.push_back(i);
314 }
315
316 return lexicographic_numbering;
317}
318
319
320
321template <int dim>
322std::unique_ptr<TensorPolynomialsBase<dim>>
324{
325 return std::make_unique<PolynomialsRaviartThomas<dim>>(*this);
326}
327
328
329
330template <int dim>
331std::vector<Point<dim>>
333{
334 Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
335 const Quadrature<1> tangential(
336 tangential_degree == 0 ?
337 static_cast<Quadrature<1>>(QMidpoint<1>()) :
338 static_cast<Quadrature<1>>(QGaussLobatto<1>(tangential_degree + 1)));
339 const Quadrature<1> normal(
340 normal_degree == 0 ?
341 static_cast<Quadrature<1>>(QMidpoint<1>()) :
342 static_cast<Quadrature<1>>(QGaussLobatto<1>(normal_degree + 1)));
343 const QAnisotropic<dim> quad =
344 (dim == 1 ? QAnisotropic<dim>(normal) :
345 (dim == 2 ? QAnisotropic<dim>(normal, tangential) :
346 QAnisotropic<dim>(normal, tangential, tangential)));
347
348 const unsigned int n_sub = polynomial_space.n();
349 std::vector<Point<dim>> points(dim * n_sub);
350 points.resize(n_polynomials(normal_degree, tangential_degree));
351 for (unsigned int d = 0; d < dim; ++d)
352 for (unsigned int i = 0; i < n_sub; ++i)
353 for (unsigned int c = 0; c < dim; ++c)
354 points[lexicographic_to_hierarchic[i + d * n_sub]][(c + d) % dim] =
355 quad.point(renumber_aniso[d][i])[c];
356 return points;
357}
358
359
360
361template class PolynomialsRaviartThomas<1>;
362template class PolynomialsRaviartThomas<2>;
363template class PolynomialsRaviartThomas<3>;
364
365
Definition point.h:112
static std::vector< unsigned int > get_lexicographic_numbering(const unsigned int normal_degree, const unsigned int tangential_degree)
std::array< std::vector< unsigned int >, dim > renumber_aniso
void evaluate(const Point< dim > &unit_point, std::vector< Tensor< 1, dim > > &values, std::vector< Tensor< 2, dim > > &grads, std::vector< Tensor< 3, dim > > &grad_grads, std::vector< Tensor< 4, dim > > &third_derivatives, std::vector< Tensor< 5, dim > > &fourth_derivatives) const override
std::string name() const override
PolynomialsRaviartThomas(const unsigned int degree_normal, const unsigned int degree_tangential)
static unsigned int n_polynomials(const unsigned int normal_degree, const unsigned int tangential_degree)
std::vector< unsigned int > lexicographic_to_hierarchic
std::vector< unsigned int > hierarchic_to_lexicographic
virtual std::unique_ptr< TensorPolynomialsBase< dim > > clone() const override
const AnisotropicPolynomials< dim > polynomial_space
std::vector< Point< dim > > get_polynomial_support_points() const
const Point< dim > & point(const unsigned int i) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 > > &points)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1672
STL namespace.