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