Reference documentation for deal.II version GIT relicensing-437-g81ec864850 2024-04-19 07:30: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\}}\)
Loading...
Searching...
No Matches
polynomial_space.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2002 - 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#ifndef dealii_polynomial_space_h
16#define dealii_polynomial_space_h
17
18
19#include <deal.II/base/config.h>
20
23#include <deal.II/base/point.h>
27#include <deal.II/base/tensor.h>
28
29#include <vector>
30
32
96template <int dim>
98{
99public:
104 static constexpr unsigned int dimension = dim;
105
113 template <class Pol>
114 PolynomialSpace(const std::vector<Pol> &pols);
115
119 template <typename StreamType>
120 void
121 output_indices(StreamType &out) const;
122
127 void
128 set_numbering(const std::vector<unsigned int> &renumber);
129
143 void
144 evaluate(const Point<dim> &unit_point,
145 std::vector<double> &values,
146 std::vector<Tensor<1, dim>> &grads,
147 std::vector<Tensor<2, dim>> &grad_grads,
148 std::vector<Tensor<3, dim>> &third_derivatives,
149 std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
150
157 double
158 compute_value(const unsigned int i, const Point<dim> &p) const override;
159
168 template <int order>
170 compute_derivative(const unsigned int i, const Point<dim> &p) const;
171
175 virtual Tensor<1, dim>
176 compute_1st_derivative(const unsigned int i,
177 const Point<dim> &p) const override;
178
182 virtual Tensor<2, dim>
183 compute_2nd_derivative(const unsigned int i,
184 const Point<dim> &p) const override;
185
189 virtual Tensor<3, dim>
190 compute_3rd_derivative(const unsigned int i,
191 const Point<dim> &p) const override;
192
196 virtual Tensor<4, dim>
197 compute_4th_derivative(const unsigned int i,
198 const Point<dim> &p) const override;
199
207 compute_grad(const unsigned int i, const Point<dim> &p) const override;
208
216 compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
217
224 static unsigned int
225 n_polynomials(const unsigned int n);
226
230 std::string
231 name() const override;
232
236 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
237 clone() const override;
238
239protected:
248 std::array<unsigned int, dim>
249 compute_index(const unsigned int n) const;
250
251private:
255 const std::vector<Polynomials::Polynomial<double>> polynomials;
256
260 std::vector<unsigned int> index_map;
261
265 std::vector<unsigned int> index_map_inverse;
266};
267
268
269/* -------------- declaration of explicit specializations --- */
270
271template <>
272std::array<unsigned int, 1>
273PolynomialSpace<1>::compute_index(const unsigned int n) const;
274template <>
275std::array<unsigned int, 2>
276PolynomialSpace<2>::compute_index(const unsigned int n) const;
277template <>
278std::array<unsigned int, 3>
279PolynomialSpace<3>::compute_index(const unsigned int n) const;
280
281
282
283/* -------------- inline and template functions ------------- */
284
285template <int dim>
286template <class Pol>
287PolynomialSpace<dim>::PolynomialSpace(const std::vector<Pol> &pols)
288 : ScalarPolynomialsBase<dim>(pols.size(), n_polynomials(pols.size()))
289 , polynomials(pols.begin(), pols.end())
290 , index_map(n_polynomials(pols.size()))
291 , index_map_inverse(n_polynomials(pols.size()))
292{
293 // per default set this index map
294 // to identity. This map can be
295 // changed by the user through the
296 // set_numbering function
297 for (unsigned int i = 0; i < this->n(); ++i)
298 {
299 index_map[i] = i;
300 index_map_inverse[i] = i;
301 }
302}
303
304
305
306template <int dim>
307inline std::string
309{
310 return "PolynomialSpace";
311}
312
313
314template <int dim>
315template <typename StreamType>
316void
318{
319 for (unsigned int i = 0; i < this->n(); ++i)
320 {
321 const std::array<unsigned int, dim> ix = compute_index(i);
322 out << i << "\t";
323 for (unsigned int d = 0; d < dim; ++d)
324 out << ix[d] << ' ';
325 out << std::endl;
326 }
327}
328
329template <int dim>
330template <int order>
333 const Point<dim> &p) const
334{
335 const std::array<unsigned int, dim> indices = compute_index(i);
336
338 {
339 std::vector<double> tmp(order + 1);
340 for (unsigned int d = 0; d < dim; ++d)
341 {
342 polynomials[indices[d]].value(p[d], tmp);
343 for (unsigned int j = 0; j < order + 1; ++j)
344 v[d][j] = tmp[j];
345 }
346 }
347
348 Tensor<order, dim> derivative;
349 switch (order)
350 {
351 case 1:
352 {
353 Tensor<1, dim> &derivative_1 =
354 *reinterpret_cast<Tensor<1, dim> *>(&derivative);
355 for (unsigned int d = 0; d < dim; ++d)
356 {
357 derivative_1[d] = 1.;
358 for (unsigned int x = 0; x < dim; ++x)
359 {
360 unsigned int x_order = 0;
361 if (d == x)
362 ++x_order;
363
364 derivative_1[d] *= v[x][x_order];
365 }
366 }
367
368 return derivative;
369 }
370 case 2:
371 {
372 Tensor<2, dim> &derivative_2 =
373 *reinterpret_cast<Tensor<2, dim> *>(&derivative);
374 for (unsigned int d1 = 0; d1 < dim; ++d1)
375 for (unsigned int d2 = 0; d2 < dim; ++d2)
376 {
377 derivative_2[d1][d2] = 1.;
378 for (unsigned int x = 0; x < dim; ++x)
379 {
380 unsigned int x_order = 0;
381 if (d1 == x)
382 ++x_order;
383 if (d2 == x)
384 ++x_order;
385
386 derivative_2[d1][d2] *= v[x][x_order];
387 }
388 }
389
390 return derivative;
391 }
392 case 3:
393 {
394 Tensor<3, dim> &derivative_3 =
395 *reinterpret_cast<Tensor<3, dim> *>(&derivative);
396 for (unsigned int d1 = 0; d1 < dim; ++d1)
397 for (unsigned int d2 = 0; d2 < dim; ++d2)
398 for (unsigned int d3 = 0; d3 < dim; ++d3)
399 {
400 derivative_3[d1][d2][d3] = 1.;
401 for (unsigned int x = 0; x < dim; ++x)
402 {
403 unsigned int x_order = 0;
404 if (d1 == x)
405 ++x_order;
406 if (d2 == x)
407 ++x_order;
408 if (d3 == x)
409 ++x_order;
410
411 derivative_3[d1][d2][d3] *= v[x][x_order];
412 }
413 }
414
415 return derivative;
416 }
417 case 4:
418 {
419 Tensor<4, dim> &derivative_4 =
420 *reinterpret_cast<Tensor<4, dim> *>(&derivative);
421 for (unsigned int d1 = 0; d1 < dim; ++d1)
422 for (unsigned int d2 = 0; d2 < dim; ++d2)
423 for (unsigned int d3 = 0; d3 < dim; ++d3)
424 for (unsigned int d4 = 0; d4 < dim; ++d4)
425 {
426 derivative_4[d1][d2][d3][d4] = 1.;
427 for (unsigned int x = 0; x < dim; ++x)
428 {
429 unsigned int x_order = 0;
430 if (d1 == x)
431 ++x_order;
432 if (d2 == x)
433 ++x_order;
434 if (d3 == x)
435 ++x_order;
436 if (d4 == x)
437 ++x_order;
438
439 derivative_4[d1][d2][d3][d4] *= v[x][x_order];
440 }
441 }
442
443 return derivative;
444 }
445 default:
446 {
448 return derivative;
449 }
450 }
451}
452
453
454
455template <int dim>
456inline Tensor<1, dim>
458 const Point<dim> &p) const
459{
460 return compute_derivative<1>(i, p);
461}
462
463
464
465template <int dim>
466inline Tensor<2, dim>
468 const Point<dim> &p) const
469{
470 return compute_derivative<2>(i, p);
471}
472
473
474
475template <int dim>
476inline Tensor<3, dim>
478 const Point<dim> &p) const
479{
480 return compute_derivative<3>(i, p);
481}
482
483
484
485template <int dim>
486inline Tensor<4, dim>
488 const Point<dim> &p) const
489{
490 return compute_derivative<4>(i, p);
491}
492
494
495#endif
Definition point.h:111
const std::vector< Polynomials::Polynomial< double > > polynomials
void output_indices(StreamType &out) const
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
PolynomialSpace(const std::vector< Pol > &pols)
double compute_value(const unsigned int i, const Point< dim > &p) const override
static constexpr unsigned int dimension
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
static unsigned int n_polynomials(const unsigned int n)
std::array< unsigned int, dim > compute_index(const unsigned int n) const
std::string name() const override
virtual Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim > > &grads, std::vector< Tensor< 2, dim > > &grad_grads, std::vector< Tensor< 3, dim > > &third_derivatives, std::vector< Tensor< 4, dim > > &fourth_derivatives) const override
std::vector< unsigned int > index_map_inverse
virtual Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DEAL_II_NOT_IMPLEMENTED()
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107