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