Reference documentation for deal.II version 9.2.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\}}\)
polynomial_space.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2020 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/point.h>
27 #include <deal.II/base/tensor.h>
28 
29 #include <vector>
30 
32 
98 template <int dim>
100 {
101 public:
106  static const unsigned int dimension = dim;
107 
115  template <class Pol>
116  PolynomialSpace(const std::vector<Pol> &pols);
117 
121  template <class StreamType>
122  void
123  output_indices(StreamType &out) const;
124 
129  void
130  set_numbering(const std::vector<unsigned int> &renumber);
131 
145  void
146  evaluate(const Point<dim> & unit_point,
147  std::vector<double> & values,
148  std::vector<Tensor<1, dim>> &grads,
149  std::vector<Tensor<2, dim>> &grad_grads,
150  std::vector<Tensor<3, dim>> &third_derivatives,
151  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
152 
159  double
160  compute_value(const unsigned int i, const Point<dim> &p) const;
161 
170  template <int order>
172  compute_derivative(const unsigned int i, const Point<dim> &p) const;
173 
181  compute_grad(const unsigned int i, const Point<dim> &p) const;
182 
190  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
191 
198  static unsigned int
199  n_polynomials(const unsigned int n);
200 
204  std::string
205  name() const override;
206 
210  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
211  clone() const override;
212 
213 protected:
222  std::array<unsigned int, dim>
223  compute_index(const unsigned int n) const;
224 
225 private:
229  const std::vector<Polynomials::Polynomial<double>> polynomials;
230 
234  std::vector<unsigned int> index_map;
235 
239  std::vector<unsigned int> index_map_inverse;
240 };
241 
242 
243 /* -------------- declaration of explicit specializations --- */
244 
245 template <>
246 std::array<unsigned int, 1>
247 PolynomialSpace<1>::compute_index(const unsigned int n) const;
248 template <>
249 std::array<unsigned int, 2>
250 PolynomialSpace<2>::compute_index(const unsigned int n) const;
251 template <>
252 std::array<unsigned int, 3>
253 PolynomialSpace<3>::compute_index(const unsigned int n) const;
254 
255 
256 
257 /* -------------- inline and template functions ------------- */
258 
259 template <int dim>
260 template <class Pol>
261 PolynomialSpace<dim>::PolynomialSpace(const std::vector<Pol> &pols)
262  : ScalarPolynomialsBase<dim>(pols.size(), n_polynomials(pols.size()))
263  , polynomials(pols.begin(), pols.end())
264  , index_map(n_polynomials(pols.size()))
265  , index_map_inverse(n_polynomials(pols.size()))
266 {
267  // per default set this index map
268  // to identity. This map can be
269  // changed by the user through the
270  // set_numbering function
271  for (unsigned int i = 0; i < this->n(); ++i)
272  {
273  index_map[i] = i;
274  index_map_inverse[i] = i;
275  }
276 }
277 
278 
279 
280 template <int dim>
281 inline std::string
283 {
284  return "PolynomialSpace";
285 }
286 
287 
288 template <int dim>
289 template <class StreamType>
290 void
292 {
293  for (unsigned int i = 0; i < this->n(); ++i)
294  {
295  const std::array<unsigned int, dim> ix = compute_index(i);
296  out << i << "\t";
297  for (unsigned int d = 0; d < dim; ++d)
298  out << ix[d] << " ";
299  out << std::endl;
300  }
301 }
302 
303 template <int dim>
304 template <int order>
307  const Point<dim> & p) const
308 {
309  const std::array<unsigned int, dim> indices = compute_index(i);
310 
311  double v[dim][order + 1];
312  {
313  std::vector<double> tmp(order + 1);
314  for (unsigned int d = 0; d < dim; ++d)
315  {
316  polynomials[indices[d]].value(p(d), tmp);
317  for (unsigned int j = 0; j < order + 1; ++j)
318  v[d][j] = tmp[j];
319  }
320  }
321 
322  Tensor<order, dim> derivative;
323  switch (order)
324  {
325  case 1:
326  {
327  Tensor<1, dim> &derivative_1 =
328  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
329  for (unsigned int d = 0; d < dim; ++d)
330  {
331  derivative_1[d] = 1.;
332  for (unsigned int x = 0; x < dim; ++x)
333  {
334  unsigned int x_order = 0;
335  if (d == x)
336  ++x_order;
337 
338  derivative_1[d] *= v[x][x_order];
339  }
340  }
341 
342  return derivative;
343  }
344  case 2:
345  {
346  Tensor<2, dim> &derivative_2 =
347  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
348  for (unsigned int d1 = 0; d1 < dim; ++d1)
349  for (unsigned int d2 = 0; d2 < dim; ++d2)
350  {
351  derivative_2[d1][d2] = 1.;
352  for (unsigned int x = 0; x < dim; ++x)
353  {
354  unsigned int x_order = 0;
355  if (d1 == x)
356  ++x_order;
357  if (d2 == x)
358  ++x_order;
359 
360  derivative_2[d1][d2] *= v[x][x_order];
361  }
362  }
363 
364  return derivative;
365  }
366  case 3:
367  {
368  Tensor<3, dim> &derivative_3 =
369  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
370  for (unsigned int d1 = 0; d1 < dim; ++d1)
371  for (unsigned int d2 = 0; d2 < dim; ++d2)
372  for (unsigned int d3 = 0; d3 < dim; ++d3)
373  {
374  derivative_3[d1][d2][d3] = 1.;
375  for (unsigned int x = 0; x < dim; ++x)
376  {
377  unsigned int x_order = 0;
378  if (d1 == x)
379  ++x_order;
380  if (d2 == x)
381  ++x_order;
382  if (d3 == x)
383  ++x_order;
384 
385  derivative_3[d1][d2][d3] *= v[x][x_order];
386  }
387  }
388 
389  return derivative;
390  }
391  case 4:
392  {
393  Tensor<4, dim> &derivative_4 =
394  *reinterpret_cast<Tensor<4, dim> *>(&derivative);
395  for (unsigned int d1 = 0; d1 < dim; ++d1)
396  for (unsigned int d2 = 0; d2 < dim; ++d2)
397  for (unsigned int d3 = 0; d3 < dim; ++d3)
398  for (unsigned int d4 = 0; d4 < dim; ++d4)
399  {
400  derivative_4[d1][d2][d3][d4] = 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  if (d4 == x)
411  ++x_order;
412 
413  derivative_4[d1][d2][d3][d4] *= v[x][x_order];
414  }
415  }
416 
417  return derivative;
418  }
419  default:
420  {
421  Assert(false, ExcNotImplemented());
422  return derivative;
423  }
424  }
425 }
426 
427 
429 
430 #endif
void set_numbering(const std::vector< unsigned int > &renumber)
static const unsigned int dimension
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
std::vector< unsigned int > index_map_inverse
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
static unsigned int n_polynomials(const unsigned int n)
std::string name() const override
std::array< unsigned int, dim > compute_index(const unsigned int n) const
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
void output_indices(StreamType &out) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
PolynomialSpace(const std::vector< Pol > &pols)
static ::ExceptionBase & ExcNotImplemented()
std::vector< unsigned int > index_map
const std::vector< Polynomials::Polynomial< double > > polynomials