Reference documentation for deal.II version GIT 4efb66ecd0 2023-02-07 13:45: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\}}\)
polynomials_raviart_thomas.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 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 
20 
21 #include <iomanip>
22 #include <iostream>
23 #include <memory>
24 
25 
27 
28 
29 namespace
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 
63 template <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 
125 template <int dim>
127  : PolynomialsRaviartThomas(k + 1, k)
128 {}
129 
130 
131 
132 template <int dim>
133 void
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 
226 template <int dim>
227 std::string
229 {
230  return "RaviartThomas";
231 }
232 
233 
234 
235 template <int dim>
236 unsigned int
238 {
239  return n_polynomials(degree + 1, degree);
240 }
241 
242 
243 
244 template <int dim>
245 unsigned 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 
256 template <int dim>
257 std::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 
321 template <int dim>
322 std::unique_ptr<TensorPolynomialsBase<dim>>
324 {
325  return std::make_unique<PolynomialsRaviartThomas<dim>>(*this);
326 }
327 
328 
329 
330 template <int dim>
331 std::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 
361 template class PolynomialsRaviartThomas<1>;
362 template class PolynomialsRaviartThomas<2>;
363 template class PolynomialsRaviartThomas<3>;
364 
365 
Definition: point.h:111
static std::vector< unsigned int > get_lexicographic_numbering(const unsigned int normal_degree, const unsigned int tangential_degree)
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::array< std::vector< unsigned int >, dim > renumber_aniso
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
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1583
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
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)
Definition: polynomial.cc:702
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:450
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
Definition: utilities.h:1655