Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
polynomial_space.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2018 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 #include <deal.II/base/polynomial.h>
25 #include <deal.II/base/smartpointer.h>
26 #include <deal.II/base/tensor.h>
27 
28 #include <vector>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
97 template <int dim>
99 {
100 public:
105  static const 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  compute(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;
151 
158  double
159  compute_value(const unsigned int i, const Point<dim> &p) const;
160 
169  template <int order>
171  compute_derivative(const unsigned int i, const Point<dim> &p) const;
172 
180  compute_grad(const unsigned int i, const Point<dim> &p) const;
181 
189  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
190 
197  unsigned int
198  n() const;
199 
206  unsigned int
207  degree() const;
208 
216  static unsigned int
217  compute_n_pols(const unsigned int n);
218 
219 protected:
228  std::array<unsigned int, dim>
229  compute_index(const unsigned int n) const;
230 
231 private:
235  const std::vector<Polynomials::Polynomial<double>> polynomials;
236 
240  const unsigned int n_pols;
241 
245  std::vector<unsigned int> index_map;
246 
250  std::vector<unsigned int> index_map_inverse;
251 };
252 
253 
254 /* -------------- declaration of explicit specializations --- */
255 
256 template <>
257 std::array<unsigned int, 1>
258 PolynomialSpace<1>::compute_index(const unsigned int n) const;
259 template <>
260 std::array<unsigned int, 2>
261 PolynomialSpace<2>::compute_index(const unsigned int n) const;
262 template <>
263 std::array<unsigned int, 3>
264 PolynomialSpace<3>::compute_index(const unsigned int n) const;
265 
266 
267 
268 /* -------------- inline and template functions ------------- */
269 
270 template <int dim>
271 template <class Pol>
272 PolynomialSpace<dim>::PolynomialSpace(const std::vector<Pol> &pols)
273  : polynomials(pols.begin(), pols.end())
274  , n_pols(compute_n_pols(polynomials.size()))
275  , index_map(n_pols)
276  , index_map_inverse(n_pols)
277 {
278  // per default set this index map
279  // to identity. This map can be
280  // changed by the user through the
281  // set_numbering function
282  for (unsigned int i = 0; i < n_pols; ++i)
283  {
284  index_map[i] = i;
285  index_map_inverse[i] = i;
286  }
287 }
288 
289 
290 template <int dim>
291 inline unsigned int
293 {
294  return n_pols;
295 }
296 
297 
298 
299 template <int dim>
300 inline unsigned int
302 {
303  return polynomials.size();
304 }
305 
306 
307 template <int dim>
308 template <class StreamType>
309 void
311 {
312  for (unsigned int i = 0; i < n_pols; ++i)
313  {
314  const std::array<unsigned int, dim> ix = compute_index(i);
315  out << i << "\t";
316  for (unsigned int d = 0; d < dim; ++d)
317  out << ix[d] << " ";
318  out << std::endl;
319  }
320 }
321 
322 template <int dim>
323 template <int order>
326  const Point<dim> & p) const
327 {
328  const std::array<unsigned int, dim> indices = compute_index(i);
329 
330  double v[dim][order + 1];
331  {
332  std::vector<double> tmp(order + 1);
333  for (unsigned int d = 0; d < dim; ++d)
334  {
335  polynomials[indices[d]].value(p(d), tmp);
336  for (unsigned int j = 0; j < order + 1; ++j)
337  v[d][j] = tmp[j];
338  }
339  }
340 
341  Tensor<order, dim> derivative;
342  switch (order)
343  {
344  case 1:
345  {
346  Tensor<1, dim> &derivative_1 =
347  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
348  for (unsigned int d = 0; d < dim; ++d)
349  {
350  derivative_1[d] = 1.;
351  for (unsigned int x = 0; x < dim; ++x)
352  {
353  unsigned int x_order = 0;
354  if (d == x)
355  ++x_order;
356 
357  derivative_1[d] *= v[x][x_order];
358  }
359  }
360 
361  return derivative;
362  }
363  case 2:
364  {
365  Tensor<2, dim> &derivative_2 =
366  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
367  for (unsigned int d1 = 0; d1 < dim; ++d1)
368  for (unsigned int d2 = 0; d2 < dim; ++d2)
369  {
370  derivative_2[d1][d2] = 1.;
371  for (unsigned int x = 0; x < dim; ++x)
372  {
373  unsigned int x_order = 0;
374  if (d1 == x)
375  ++x_order;
376  if (d2 == x)
377  ++x_order;
378 
379  derivative_2[d1][d2] *= v[x][x_order];
380  }
381  }
382 
383  return derivative;
384  }
385  case 3:
386  {
387  Tensor<3, dim> &derivative_3 =
388  *reinterpret_cast<Tensor<3, dim> *>(&derivative);
389  for (unsigned int d1 = 0; d1 < dim; ++d1)
390  for (unsigned int d2 = 0; d2 < dim; ++d2)
391  for (unsigned int d3 = 0; d3 < dim; ++d3)
392  {
393  derivative_3[d1][d2][d3] = 1.;
394  for (unsigned int x = 0; x < dim; ++x)
395  {
396  unsigned int x_order = 0;
397  if (d1 == x)
398  ++x_order;
399  if (d2 == x)
400  ++x_order;
401  if (d3 == x)
402  ++x_order;
403 
404  derivative_3[d1][d2][d3] *= v[x][x_order];
405  }
406  }
407 
408  return derivative;
409  }
410  case 4:
411  {
412  Tensor<4, dim> &derivative_4 =
413  *reinterpret_cast<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_4[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_4[d1][d2][d3][d4] *= v[x][x_order];
433  }
434  }
435 
436  return derivative;
437  }
438  default:
439  {
440  Assert(false, ExcNotImplemented());
441  return derivative;
442  }
443  }
444 }
445 
446 
447 DEAL_II_NAMESPACE_CLOSE
448 
449 #endif
void set_numbering(const std::vector< unsigned int > &renumber)
static const unsigned int dimension
std::vector< unsigned int > index_map_inverse
#define Assert(cond, exc)
Definition: exceptions.h:1407
unsigned int n() const
std::array< unsigned int, dim > compute_index(const unsigned int n) const
const unsigned int n_pols
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
void output_indices(StreamType &out) const
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
static unsigned int compute_n_pols(const unsigned int n)
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
unsigned int degree() const
PolynomialSpace(const std::vector< Pol > &pols)
static ::ExceptionBase & ExcNotImplemented()
std::vector< unsigned int > index_map
const std::vector< Polynomials::Polynomial< double > > polynomials
void compute(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