deal.II version GIT relicensing-1991-gec5bc16adb 2024-10-14 08:00: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
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
24#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 if constexpr (order == 1)
349 {
350 Tensor<1, dim> derivative;
351 for (unsigned int d = 0; d < dim; ++d)
352 {
353 derivative[d] = 1.;
354 for (unsigned int x = 0; x < dim; ++x)
355 {
356 unsigned int x_order = 0;
357 if (d == x)
358 ++x_order;
359
360 derivative[d] *= v[x][x_order];
361 }
362 }
363
364 return derivative;
365 }
366 else if constexpr (order == 2)
367 {
368 Tensor<2, dim> derivative;
369 for (unsigned int d1 = 0; d1 < dim; ++d1)
370 for (unsigned int d2 = 0; d2 < dim; ++d2)
371 {
372 derivative[d1][d2] = 1.;
373 for (unsigned int x = 0; x < dim; ++x)
374 {
375 unsigned int x_order = 0;
376 if (d1 == x)
377 ++x_order;
378 if (d2 == x)
379 ++x_order;
380
381 derivative[d1][d2] *= v[x][x_order];
382 }
383 }
384
385 return derivative;
386 }
387 else if constexpr (order == 3)
388 {
389 Tensor<3, dim> derivative;
390 for (unsigned int d1 = 0; d1 < dim; ++d1)
391 for (unsigned int d2 = 0; d2 < dim; ++d2)
392 for (unsigned int d3 = 0; d3 < dim; ++d3)
393 {
394 derivative[d1][d2][d3] = 1.;
395 for (unsigned int x = 0; x < dim; ++x)
396 {
397 unsigned int x_order = 0;
398 if (d1 == x)
399 ++x_order;
400 if (d2 == x)
401 ++x_order;
402 if (d3 == x)
403 ++x_order;
404
405 derivative[d1][d2][d3] *= v[x][x_order];
406 }
407 }
408
409 return derivative;
410 }
411 else if constexpr (order == 4)
412 {
413 Tensor<4, dim> derivative;
414 for (unsigned int d1 = 0; d1 < dim; ++d1)
415 for (unsigned int d2 = 0; d2 < dim; ++d2)
416 for (unsigned int d3 = 0; d3 < dim; ++d3)
417 for (unsigned int d4 = 0; d4 < dim; ++d4)
418 {
419 derivative[d1][d2][d3][d4] = 1.;
420 for (unsigned int x = 0; x < dim; ++x)
421 {
422 unsigned int x_order = 0;
423 if (d1 == x)
424 ++x_order;
425 if (d2 == x)
426 ++x_order;
427 if (d3 == x)
428 ++x_order;
429 if (d4 == x)
430 ++x_order;
431
432 derivative[d1][d2][d3][d4] *= v[x][x_order];
433 }
434 }
435
436 return derivative;
437 }
438 else
439 {
441 return {};
442 }
443}
444
445
446
447template <int dim>
448inline Tensor<1, dim>
450 const Point<dim> &p) const
451{
452 return compute_derivative<1>(i, p);
453}
454
455
456
457template <int dim>
458inline Tensor<2, dim>
460 const Point<dim> &p) const
461{
462 return compute_derivative<2>(i, p);
463}
464
465
466
467template <int dim>
468inline Tensor<3, dim>
470 const Point<dim> &p) const
471{
472 return compute_derivative<3>(i, p);
473}
474
475
476
477template <int dim>
478inline Tensor<4, dim>
480 const Point<dim> &p) const
481{
482 return compute_derivative<4>(i, p);
483}
484
486
487#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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#define DEAL_II_NOT_IMPLEMENTED()
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107