Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+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\}}\)
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 
23 #include <deal.II/base/ndarray.h>
24 #include <deal.II/base/point.h>
28 #include <deal.II/base/tensor.h>
29 
30 #include <vector>
31 
33 
97 template <int dim>
99 {
100 public:
105  static constexpr unsigned int dimension = dim;
106 
114  template <class Pol>
115  PolynomialSpace(const std::vector<Pol> &pols);
116 
120  template <typename 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 
240 protected:
249  std::array<unsigned int, dim>
250  compute_index(const unsigned int n) const;
251 
252 private:
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 
272 template <>
273 std::array<unsigned int, 1>
274 PolynomialSpace<1>::compute_index(const unsigned int n) const;
275 template <>
276 std::array<unsigned int, 2>
277 PolynomialSpace<2>::compute_index(const unsigned int n) const;
278 template <>
279 std::array<unsigned int, 3>
280 PolynomialSpace<3>::compute_index(const unsigned int n) const;
281 
282 
283 
284 /* -------------- inline and template functions ------------- */
285 
286 template <int dim>
287 template <class Pol>
288 PolynomialSpace<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 
307 template <int dim>
308 inline std::string
310 {
311  return "PolynomialSpace";
312 }
313 
314 
315 template <int dim>
316 template <typename StreamType>
317 void
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 
330 template <int dim>
331 template <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 
456 template <int dim>
457 inline Tensor<1, dim>
459  const Point<dim> &p) const
460 {
461  return compute_derivative<1>(i, p);
462 }
463 
464 
465 
466 template <int dim>
467 inline Tensor<2, dim>
469  const Point<dim> &p) const
470 {
471  return compute_derivative<2>(i, p);
472 }
473 
474 
475 
476 template <int dim>
477 inline Tensor<3, dim>
479  const Point<dim> &p) const
480 {
481  return compute_derivative<3>(i, p);
482 }
483 
484 
485 
486 template <int dim>
487 inline 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)
std::array< unsigned int, dim > compute_index(const unsigned int n) const
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::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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108