Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tensor_product_polynomials.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 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 #include <deal.II/base/exceptions.h>
17 #include <deal.II/base/polynomials_piecewise.h>
18 #include <deal.II/base/tensor_product_polynomials.h>
19 
20 #include <boost/container/small_vector.hpp>
21 
22 #include <array>
23 
24 DEAL_II_NAMESPACE_OPEN
25 
26 
27 
28 /* ------------------- TensorProductPolynomials -------------- */
29 
30 
31 namespace internal
32 {
33  namespace
34  {
35  template <int dim>
36  inline void
37  compute_tensor_index(const unsigned int,
38  const unsigned int,
39  const unsigned int,
40  unsigned int (&)[dim])
41  {
42  Assert(false, ExcNotImplemented());
43  }
44 
45  inline void
46  compute_tensor_index(const unsigned int n,
47  const unsigned int,
48  const unsigned int,
49  unsigned int (&indices)[1])
50  {
51  indices[0] = n;
52  }
53 
54  inline void
55  compute_tensor_index(const unsigned int n,
56  const unsigned int n_pols_0,
57  const unsigned int,
58  unsigned int (&indices)[2])
59  {
60  indices[0] = n % n_pols_0;
61  indices[1] = n / n_pols_0;
62  }
63 
64  inline void
65  compute_tensor_index(const unsigned int n,
66  const unsigned int n_pols_0,
67  const unsigned int n_pols_1,
68  unsigned int (&indices)[3])
69  {
70  indices[0] = n % n_pols_0;
71  indices[1] = (n / n_pols_0) % n_pols_1;
72  indices[2] = n / (n_pols_0 * n_pols_1);
73  }
74  } // namespace
75 } // namespace internal
76 
77 
78 
79 template <int dim, typename PolynomialType>
80 inline void
82  const unsigned int i,
83  unsigned int (&indices)[(dim > 0 ? dim : 1)]) const
84 {
85  Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
87  internal::compute_tensor_index(index_map[i],
88  polynomials.size(),
89  polynomials.size(),
90  indices);
91 }
92 
93 
94 
95 template <int dim, typename PolynomialType>
96 void
98  std::ostream &out) const
99 {
100  unsigned int ix[dim];
101  for (unsigned int i = 0; i < n_tensor_pols; ++i)
102  {
103  compute_index(i, ix);
104  out << i << "\t";
105  for (unsigned int d = 0; d < dim; ++d)
106  out << ix[d] << " ";
107  out << std::endl;
108  }
109 }
110 
111 
112 
113 template <int dim, typename PolynomialType>
114 void
116  const std::vector<unsigned int> &renumber)
117 {
118  Assert(renumber.size() == index_map.size(),
119  ExcDimensionMismatch(renumber.size(), index_map.size()));
120 
121  index_map = renumber;
122  for (unsigned int i = 0; i < index_map.size(); ++i)
123  index_map_inverse[index_map[i]] = i;
124 }
125 
126 
127 
128 template <>
129 double
131  const unsigned int,
132  const Point<0> &) const
133 {
134  Assert(false, ExcNotImplemented());
135  return 0;
136 }
137 
138 
139 
140 template <int dim, typename PolynomialType>
141 double
143  const unsigned int i,
144  const Point<dim> & p) const
145 {
146  Assert(dim > 0, ExcNotImplemented());
147 
148  unsigned int indices[dim];
149  compute_index(i, indices);
150 
151  double value = 1.;
152  for (unsigned int d = 0; d < dim; ++d)
153  value *= polynomials[indices[d]].value(p(d));
154 
155  return value;
156 }
157 
158 
159 
160 template <int dim, typename PolynomialType>
163  const unsigned int i,
164  const Point<dim> & p) const
165 {
166  unsigned int indices[dim];
167  compute_index(i, indices);
168 
169  // compute values and
170  // uni-directional derivatives at
171  // the given point in each
172  // co-ordinate direction
173  double v[dim][2];
174  {
175  std::vector<double> tmp(2);
176  for (unsigned int d = 0; d < dim; ++d)
177  {
178  polynomials[indices[d]].value(p(d), tmp);
179  v[d][0] = tmp[0];
180  v[d][1] = tmp[1];
181  }
182  }
183 
184  Tensor<1, dim> grad;
185  for (unsigned int d = 0; d < dim; ++d)
186  {
187  grad[d] = 1.;
188  for (unsigned int x = 0; x < dim; ++x)
189  grad[d] *= v[x][d == x];
190  }
191 
192  return grad;
193 }
194 
195 
196 
197 template <int dim, typename PolynomialType>
200  const unsigned int i,
201  const Point<dim> & p) const
202 {
203  unsigned int indices[dim];
204  compute_index(i, indices);
205 
206  double v[dim][3];
207  {
208  std::vector<double> tmp(3);
209  for (unsigned int d = 0; d < dim; ++d)
210  {
211  polynomials[indices[d]].value(p(d), tmp);
212  v[d][0] = tmp[0];
213  v[d][1] = tmp[1];
214  v[d][2] = tmp[2];
215  }
216  }
217 
218  Tensor<2, dim> grad_grad;
219  for (unsigned int d1 = 0; d1 < dim; ++d1)
220  for (unsigned int d2 = 0; d2 < dim; ++d2)
221  {
222  grad_grad[d1][d2] = 1.;
223  for (unsigned int x = 0; x < dim; ++x)
224  {
225  unsigned int derivative = 0;
226  if (d1 == x || d2 == x)
227  {
228  if (d1 == d2)
229  derivative = 2;
230  else
231  derivative = 1;
232  }
233  grad_grad[d1][d2] *= v[x][derivative];
234  }
235  }
236 
237  return grad_grad;
238 }
239 
240 
241 
242 template <int dim, typename PolynomialType>
243 void
245  const Point<dim> & p,
246  std::vector<double> & values,
247  std::vector<Tensor<1, dim>> &grads,
248  std::vector<Tensor<2, dim>> &grad_grads,
249  std::vector<Tensor<3, dim>> &third_derivatives,
250  std::vector<Tensor<4, dim>> &fourth_derivatives) const
251 {
252  Assert(dim <= 3, ExcNotImplemented());
253  Assert(values.size() == n_tensor_pols || values.size() == 0,
254  ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
255  Assert(grads.size() == n_tensor_pols || grads.size() == 0,
256  ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
257  Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
258  ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
259  Assert(third_derivatives.size() == n_tensor_pols ||
260  third_derivatives.size() == 0,
261  ExcDimensionMismatch2(third_derivatives.size(), n_tensor_pols, 0));
262  Assert(fourth_derivatives.size() == n_tensor_pols ||
263  fourth_derivatives.size() == 0,
264  ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));
265 
266  const bool update_values = (values.size() == n_tensor_pols),
267  update_grads = (grads.size() == n_tensor_pols),
268  update_grad_grads = (grad_grads.size() == n_tensor_pols),
270  (third_derivatives.size() == n_tensor_pols),
271  update_4th_derivatives =
272  (fourth_derivatives.size() == n_tensor_pols);
273 
274  // check how many values/derivatives we have to compute
275  unsigned int n_values_and_derivatives = 0;
276  if (update_values)
277  n_values_and_derivatives = 1;
278  if (update_grads)
279  n_values_and_derivatives = 2;
280  if (update_grad_grads)
281  n_values_and_derivatives = 3;
283  n_values_and_derivatives = 4;
284  if (update_4th_derivatives)
285  n_values_and_derivatives = 5;
286 
287  // Compute the values (and derivatives, if necessary) of all 1D polynomials
288  // at this evaluation point. We need to compute dim*n_polynomials
289  // evaluations, involving an evaluation of each polynomial for each
290  // coordinate direction. Once we have those values, we perform the
291  // multiplications for the tensor product in the arbitrary dimension.
292  const unsigned int n_polynomials = polynomials.size();
293  boost::container::small_vector<std::array<std::array<double, 5>, dim>, 20>
294  values_1d(n_polynomials);
295  if (n_values_and_derivatives == 1)
296  for (unsigned int i = 0; i < n_polynomials; ++i)
297  for (unsigned int d = 0; d < dim; ++d)
298  values_1d[i][d][0] = polynomials[i].value(p(d));
299  else
300  for (unsigned int i = 0; i < n_polynomials; ++i)
301  for (unsigned d = 0; d < dim; ++d)
302  polynomials[i].value(p(d),
303  n_values_and_derivatives,
304  values_1d[i][d].data());
305 
306  unsigned int indices[3];
307  unsigned int ind = 0;
308  for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2])
309  for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1);
310  ++indices[1])
311  if (n_values_and_derivatives == 1)
312  for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
313  {
314  double value = values_1d[indices[0]][0][0];
315  for (unsigned int d = 1; d < dim; ++d)
316  value *= values_1d[indices[d]][d][0];
317  values[index_map_inverse[ind]] = value;
318  }
319  else
320  for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
321  {
322  unsigned int i = index_map_inverse[ind];
323 
324  if (update_values)
325  {
326  double value = values_1d[indices[0]][0][0];
327  for (unsigned int x = 1; x < dim; ++x)
328  value *= values_1d[indices[x]][x][0];
329  values[i] = value;
330  }
331 
332  if (update_grads)
333  for (unsigned int d = 0; d < dim; ++d)
334  {
335  double grad = 1.;
336  for (unsigned int x = 0; x < dim; ++x)
337  grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
338  grads[i][d] = grad;
339  }
340 
341  if (update_grad_grads)
342  for (unsigned int d1 = 0; d1 < dim; ++d1)
343  for (unsigned int d2 = 0; d2 < dim; ++d2)
344  {
345  double der2 = 1.;
346  for (unsigned int x = 0; x < dim; ++x)
347  {
348  unsigned int derivative = 0;
349  if (d1 == x)
350  ++derivative;
351  if (d2 == x)
352  ++derivative;
353 
354  der2 *= values_1d[indices[x]][x][derivative];
355  }
356  grad_grads[i][d1][d2] = der2;
357  }
358 
360  for (unsigned int d1 = 0; d1 < dim; ++d1)
361  for (unsigned int d2 = 0; d2 < dim; ++d2)
362  for (unsigned int d3 = 0; d3 < dim; ++d3)
363  {
364  double der3 = 1.;
365  for (unsigned int x = 0; x < dim; ++x)
366  {
367  unsigned int derivative = 0;
368  if (d1 == x)
369  ++derivative;
370  if (d2 == x)
371  ++derivative;
372  if (d3 == x)
373  ++derivative;
374 
375  der3 *= values_1d[indices[x]][x][derivative];
376  }
377  third_derivatives[i][d1][d2][d3] = der3;
378  }
379 
380  if (update_4th_derivatives)
381  for (unsigned int d1 = 0; d1 < dim; ++d1)
382  for (unsigned int d2 = 0; d2 < dim; ++d2)
383  for (unsigned int d3 = 0; d3 < dim; ++d3)
384  for (unsigned int d4 = 0; d4 < dim; ++d4)
385  {
386  double der4 = 1.;
387  for (unsigned int x = 0; x < dim; ++x)
388  {
389  unsigned int derivative = 0;
390  if (d1 == x)
391  ++derivative;
392  if (d2 == x)
393  ++derivative;
394  if (d3 == x)
395  ++derivative;
396  if (d4 == x)
397  ++derivative;
398 
399  der4 *= values_1d[indices[x]][x][derivative];
400  }
401  fourth_derivatives[i][d1][d2][d3][d4] = der4;
402  }
403  }
404 }
405 
406 
407 
408 /* ------------------- AnisotropicPolynomials -------------- */
409 
410 
411 template <int dim>
413  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
414  : polynomials(pols)
415  , n_tensor_pols(get_n_tensor_pols(pols))
416 {
417  Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
418  for (unsigned int d = 0; d < dim; ++d)
419  Assert(pols[d].size() > 0,
420  ExcMessage("The number of polynomials must be larger than zero "
421  "for all coordinate directions."));
422 }
423 
424 
425 
426 template <int dim>
427 void
429  unsigned int (&indices)[dim]) const
430 {
431 #ifdef DEBUG
432  unsigned int n_poly = 1;
433  for (unsigned int d = 0; d < dim; ++d)
434  n_poly *= polynomials[d].size();
435  Assert(i < n_poly, ExcInternalError());
436 #endif
437 
438  if (dim == 1)
439  internal::compute_tensor_index(i,
440  polynomials[0].size(),
441  0 /*not used*/,
442  indices);
443  else
444  internal::compute_tensor_index(i,
445  polynomials[0].size(),
446  polynomials[1].size(),
447  indices);
448 }
449 
450 
451 
452 template <int dim>
453 double
455  const Point<dim> & p) const
456 {
457  unsigned int indices[dim];
458  compute_index(i, indices);
459 
460  double value = 1.;
461  for (unsigned int d = 0; d < dim; ++d)
462  value *= polynomials[d][indices[d]].value(p(d));
463 
464  return value;
465 }
466 
467 
468 template <int dim>
471  const Point<dim> & p) const
472 {
473  unsigned int indices[dim];
474  compute_index(i, indices);
475 
476  // compute values and
477  // uni-directional derivatives at
478  // the given point in each
479  // co-ordinate direction
480  std::vector<std::vector<double>> v(dim, std::vector<double>(2));
481  for (unsigned int d = 0; d < dim; ++d)
482  polynomials[d][indices[d]].value(p(d), v[d]);
483 
484  Tensor<1, dim> grad;
485  for (unsigned int d = 0; d < dim; ++d)
486  {
487  grad[d] = 1.;
488  for (unsigned int x = 0; x < dim; ++x)
489  grad[d] *= v[x][d == x];
490  }
491 
492  return grad;
493 }
494 
495 
496 template <int dim>
499  const Point<dim> & p) const
500 {
501  unsigned int indices[dim];
502  compute_index(i, indices);
503 
504  std::vector<std::vector<double>> v(dim, std::vector<double>(3));
505  for (unsigned int d = 0; d < dim; ++d)
506  polynomials[d][indices[d]].value(p(d), v[d]);
507 
508  Tensor<2, dim> grad_grad;
509  for (unsigned int d1 = 0; d1 < dim; ++d1)
510  for (unsigned int d2 = 0; d2 < dim; ++d2)
511  {
512  grad_grad[d1][d2] = 1.;
513  for (unsigned int x = 0; x < dim; ++x)
514  {
515  unsigned int derivative = 0;
516  if (d1 == x || d2 == x)
517  {
518  if (d1 == d2)
519  derivative = 2;
520  else
521  derivative = 1;
522  }
523  grad_grad[d1][d2] *= v[x][derivative];
524  }
525  }
526 
527  return grad_grad;
528 }
529 
530 
531 
532 template <int dim>
533 void
535  const Point<dim> & p,
536  std::vector<double> & values,
537  std::vector<Tensor<1, dim>> &grads,
538  std::vector<Tensor<2, dim>> &grad_grads,
539  std::vector<Tensor<3, dim>> &third_derivatives,
540  std::vector<Tensor<4, dim>> &fourth_derivatives) const
541 {
542  Assert(values.size() == n_tensor_pols || values.size() == 0,
543  ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
544  Assert(grads.size() == n_tensor_pols || grads.size() == 0,
545  ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
546  Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
547  ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
548  Assert(third_derivatives.size() == n_tensor_pols ||
549  third_derivatives.size() == 0,
550  ExcDimensionMismatch2(third_derivatives.size(), n_tensor_pols, 0));
551  Assert(fourth_derivatives.size() == n_tensor_pols ||
552  fourth_derivatives.size() == 0,
553  ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));
554 
555  const bool update_values = (values.size() == n_tensor_pols),
556  update_grads = (grads.size() == n_tensor_pols),
557  update_grad_grads = (grad_grads.size() == n_tensor_pols),
559  (third_derivatives.size() == n_tensor_pols),
560  update_4th_derivatives =
561  (fourth_derivatives.size() == n_tensor_pols);
562 
563  // check how many
564  // values/derivatives we have to
565  // compute
566  unsigned int n_values_and_derivatives = 0;
567  if (update_values)
568  n_values_and_derivatives = 1;
569  if (update_grads)
570  n_values_and_derivatives = 2;
571  if (update_grad_grads)
572  n_values_and_derivatives = 3;
574  n_values_and_derivatives = 4;
575  if (update_4th_derivatives)
576  n_values_and_derivatives = 5;
577 
578  // compute the values (and
579  // derivatives, if necessary) of
580  // all polynomials at this
581  // evaluation point
582  std::vector<std::vector<std::vector<double>>> v(dim);
583  for (unsigned int d = 0; d < dim; ++d)
584  {
585  v[d].resize(polynomials[d].size());
586  for (unsigned int i = 0; i < polynomials[d].size(); ++i)
587  {
588  v[d][i].resize(n_values_and_derivatives, 0.);
589  polynomials[d][i].value(p(d), v[d][i]);
590  }
591  }
592 
593  for (unsigned int i = 0; i < n_tensor_pols; ++i)
594  {
595  // first get the
596  // one-dimensional indices of
597  // this particular tensor
598  // product polynomial
599  unsigned int indices[dim];
600  compute_index(i, indices);
601 
602  if (update_values)
603  {
604  values[i] = 1;
605  for (unsigned int x = 0; x < dim; ++x)
606  values[i] *= v[x][indices[x]][0];
607  }
608 
609  if (update_grads)
610  for (unsigned int d = 0; d < dim; ++d)
611  {
612  grads[i][d] = 1.;
613  for (unsigned int x = 0; x < dim; ++x)
614  grads[i][d] *= v[x][indices[x]][d == x ? 1 : 0];
615  }
616 
617  if (update_grad_grads)
618  for (unsigned int d1 = 0; d1 < dim; ++d1)
619  for (unsigned int d2 = 0; d2 < dim; ++d2)
620  {
621  grad_grads[i][d1][d2] = 1.;
622  for (unsigned int x = 0; x < dim; ++x)
623  {
624  unsigned int derivative = 0;
625  if (d1 == x)
626  ++derivative;
627  if (d2 == x)
628  ++derivative;
629 
630  grad_grads[i][d1][d2] *= v[x][indices[x]][derivative];
631  }
632  }
633 
635  for (unsigned int d1 = 0; d1 < dim; ++d1)
636  for (unsigned int d2 = 0; d2 < dim; ++d2)
637  for (unsigned int d3 = 0; d3 < dim; ++d3)
638  {
639  third_derivatives[i][d1][d2][d3] = 1.;
640  for (unsigned int x = 0; x < dim; ++x)
641  {
642  unsigned int derivative = 0;
643  if (d1 == x)
644  ++derivative;
645  if (d2 == x)
646  ++derivative;
647  if (d3 == x)
648  ++derivative;
649 
650  third_derivatives[i][d1][d2][d3] *=
651  v[x][indices[x]][derivative];
652  }
653  }
654 
655  if (update_4th_derivatives)
656  for (unsigned int d1 = 0; d1 < dim; ++d1)
657  for (unsigned int d2 = 0; d2 < dim; ++d2)
658  for (unsigned int d3 = 0; d3 < dim; ++d3)
659  for (unsigned int d4 = 0; d4 < dim; ++d4)
660  {
661  fourth_derivatives[i][d1][d2][d3][d4] = 1.;
662  for (unsigned int x = 0; x < dim; ++x)
663  {
664  unsigned int derivative = 0;
665  if (d1 == x)
666  ++derivative;
667  if (d2 == x)
668  ++derivative;
669  if (d3 == x)
670  ++derivative;
671  if (d4 == x)
672  ++derivative;
673 
674  fourth_derivatives[i][d1][d2][d3][d4] *=
675  v[x][indices[x]][derivative];
676  }
677  }
678  }
679 }
680 
681 
682 
683 template <int dim>
684 unsigned int
686 {
687  return n_tensor_pols;
688 }
689 
690 
691 template <int dim>
692 unsigned int
694  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
695 {
696  unsigned int y = 1;
697  for (unsigned int d = 0; d < dim; ++d)
698  y *= pols[d].size();
699  return y;
700 }
701 
702 
703 
704 /* ------------------- explicit instantiations -------------- */
708 
709 template class TensorProductPolynomials<
710  1,
712 template class TensorProductPolynomials<
713  2,
715 template class TensorProductPolynomials<
716  3,
718 
719 template class AnisotropicPolynomials<1>;
720 template class AnisotropicPolynomials<2>;
721 template class AnisotropicPolynomials<3>;
722 
723 DEAL_II_NAMESPACE_CLOSE
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
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
Shape function values.
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
void output_indices(std::ostream &out) const
void set_numbering(const std::vector< unsigned int > &renumber)
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
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
Definition: point.h:110
void compute_index(const unsigned int i, unsigned int(&indices)[dim]) const
static ::ExceptionBase & ExcMessage(std::string arg1)
Third derivatives of shape functions.
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
double compute_value(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcNotImplemented()
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
void compute_index(const unsigned int i, unsigned int(&indices)[(dim > 0 ? dim :1)]) const
double compute_value(const unsigned int i, const Point< dim > &p) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcInternalError()