Reference documentation for deal.II version GIT 8ad845fa53 2023-12-10 20:50:02+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\}}\)
tensor_product_polynomials.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 
19 #include <deal.II/base/table.h>
21 
22 #include <boost/container/small_vector.hpp>
23 
24 #include <array>
25 #include <memory>
26 
28 
29 
30 
31 /* ------------------- TensorProductPolynomials -------------- */
32 
33 
34 namespace internal
35 {
36  namespace
37  {
38  template <std::size_t dim>
39  inline void
40  compute_tensor_index(const unsigned int,
41  const unsigned int,
42  const unsigned int,
43  std::array<unsigned int, dim> &)
44  {
45  Assert(false, ExcNotImplemented());
46  }
47 
48  inline void
49  compute_tensor_index(const unsigned int n,
50  const unsigned int,
51  const unsigned int,
52  std::array<unsigned int, 1> &indices)
53  {
54  indices[0] = n;
55  }
56 
57  inline void
58  compute_tensor_index(const unsigned int n,
59  const unsigned int n_pols_0,
60  const unsigned int,
61  std::array<unsigned int, 2> &indices)
62  {
63  indices[0] = n % n_pols_0;
64  indices[1] = n / n_pols_0;
65  }
66 
67  inline void
68  compute_tensor_index(const unsigned int n,
69  const unsigned int n_pols_0,
70  const unsigned int n_pols_1,
71  std::array<unsigned int, 3> &indices)
72  {
73  indices[0] = n % n_pols_0;
74  indices[1] = (n / n_pols_0) % n_pols_1;
75  indices[2] = n / (n_pols_0 * n_pols_1);
76  }
77  } // namespace
78 } // namespace internal
79 
80 
81 
82 template <int dim, typename PolynomialType>
83 inline void
85  const unsigned int i,
86  std::array<unsigned int, dim> &indices) const
87 {
88  Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
90  internal::compute_tensor_index(index_map[i],
91  polynomials.size(),
92  polynomials.size(),
93  indices);
94 }
95 
96 
97 
98 template <>
99 inline void
101  const unsigned int,
102  std::array<unsigned int, 0> &) const
103 {
104  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
105 }
106 
107 
108 
109 template <int dim, typename PolynomialType>
110 void
112  std::ostream &out) const
113 {
114  std::array<unsigned int, dim> ix;
115  for (unsigned int i = 0; i < this->n(); ++i)
116  {
117  compute_index(i, ix);
118  out << i << "\t";
119  for (unsigned int d = 0; d < dim; ++d)
120  out << ix[d] << " ";
121  out << std::endl;
122  }
123 }
124 
125 
126 
127 template <>
128 void
130  std::ostream &) const
131 {
132  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
133 }
134 
135 
136 
137 template <int dim, typename PolynomialType>
138 void
140  const std::vector<unsigned int> &renumber)
141 {
142  Assert(renumber.size() == index_map.size(),
143  ExcDimensionMismatch(renumber.size(), index_map.size()));
144 
145  index_map = renumber;
146  for (unsigned int i = 0; i < index_map.size(); ++i)
147  index_map_inverse[index_map[i]] = i;
148 }
149 
150 
151 
152 template <>
153 void
155  const std::vector<unsigned int> &)
156 {
157  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
158 }
159 
160 
161 
162 template <int dim, typename PolynomialType>
163 double
165  const unsigned int i,
166  const Point<dim> &p) const
167 {
168  Assert(dim > 0, ExcNotImplemented());
169 
170  std::array<unsigned int, dim> indices;
171  compute_index(i, indices);
172 
173  double value = 1.;
174  for (unsigned int d = 0; d < dim; ++d)
175  value *= polynomials[indices[d]].value(p(d));
176 
177  return value;
178 }
179 
180 
181 
182 template <>
183 double
185  const unsigned int,
186  const Point<0> &) const
187 {
188  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
189 
190  return {};
191 }
192 
193 
194 
195 template <int dim, typename PolynomialType>
198  const unsigned int i,
199  const Point<dim> &p) const
200 {
201  std::array<unsigned int, dim> indices;
202  compute_index(i, indices);
203 
204  // compute values and
205  // uni-directional derivatives at
206  // the given point in each
207  // coordinate direction
209  {
210  std::vector<double> tmp(2);
211  for (unsigned int d = 0; d < dim; ++d)
212  {
213  polynomials[indices[d]].value(p(d), tmp);
214  v[d][0] = tmp[0];
215  v[d][1] = tmp[1];
216  }
217  }
218 
219  Tensor<1, dim> grad;
220  for (unsigned int d = 0; d < dim; ++d)
221  {
222  grad[d] = 1.;
223  for (unsigned int x = 0; x < dim; ++x)
224  grad[d] *= v[x][d == x];
225  }
226 
227  return grad;
228 }
229 
230 
231 
232 template <>
235  const unsigned int,
236  const Point<0> &) const
237 {
238  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
239 
240  return {};
241 }
242 
243 
244 
245 template <int dim, typename PolynomialType>
248  const unsigned int i,
249  const Point<dim> &p) const
250 {
251  std::array<unsigned int, dim> indices;
252  compute_index(i, indices);
253 
255  {
256  std::vector<double> tmp(3);
257  for (unsigned int d = 0; d < dim; ++d)
258  {
259  polynomials[indices[d]].value(p(d), tmp);
260  v[d][0] = tmp[0];
261  v[d][1] = tmp[1];
262  v[d][2] = tmp[2];
263  }
264  }
265 
266  Tensor<2, dim> grad_grad;
267  for (unsigned int d1 = 0; d1 < dim; ++d1)
268  for (unsigned int d2 = 0; d2 < dim; ++d2)
269  {
270  grad_grad[d1][d2] = 1.;
271  for (unsigned int x = 0; x < dim; ++x)
272  {
273  unsigned int derivative = 0;
274  if (d1 == x || d2 == x)
275  {
276  if (d1 == d2)
277  derivative = 2;
278  else
279  derivative = 1;
280  }
281  grad_grad[d1][d2] *= v[x][derivative];
282  }
283  }
284 
285  return grad_grad;
286 }
287 
288 
289 
290 template <>
293  const unsigned int,
294  const Point<0> &) const
295 {
296  return {};
297 }
298 
299 
300 
301 namespace internal
302 {
304  {
305  // This function computes the tensor product of some tabulated
306  // one-dimensional polynomials (also the anisotropic case is supported)
307  // with tensor product indices of all dimensions except the first one
308  // tabulated in the 'indices' array; the first dimension is manually
309  // iterated through because these are possibly performance-critical loops,
310  // so we want to avoid indirect addressing.
311  template <int dim, std::size_t dim1>
312  void
314  const unsigned int n_derivatives,
315  const boost::container::small_vector<::ndarray<double, 5, dim>, 10>
316  &values_1d,
317  const unsigned int size_x,
318  const boost::container::small_vector<std::array<unsigned int, dim1>, 64>
319  &indices,
320  const std::vector<unsigned int> &index_map,
321  std::vector<double> &values,
322  std::vector<Tensor<1, dim>> &grads,
323  std::vector<Tensor<2, dim>> &grad_grads,
324  std::vector<Tensor<3, dim>> &third_derivatives,
325  std::vector<Tensor<4, dim>> &fourth_derivatives)
326  {
327  const bool update_values = (values.size() == indices.size() * size_x);
328  const bool update_grads = (grads.size() == indices.size() * size_x);
329  const bool update_grad_grads =
330  (grad_grads.size() == indices.size() * size_x);
331  const bool update_3rd_derivatives =
332  (third_derivatives.size() == indices.size() * size_x);
333  const bool update_4th_derivatives =
334  (fourth_derivatives.size() == indices.size() * size_x);
335 
336  // For values, 1st and 2nd derivatives use a more lengthy code that
337  // minimizes the number of arithmetic operations and memory accesses
338  if (n_derivatives == 0)
339  for (unsigned int i = 0, i1 = 0; i1 < indices.size(); ++i1)
340  {
341  double value_outer = 1.;
342  if constexpr (dim > 1)
343  for (unsigned int d = 1; d < dim; ++d)
344  value_outer *= values_1d[indices[i1][d - 1]][0][d];
345  if (index_map.empty())
346  for (unsigned int ix = 0; ix < size_x; ++ix, ++i)
347  values[i] = value_outer * values_1d[ix][0][0];
348  else
349  for (unsigned int ix = 0; ix < size_x; ++ix, ++i)
350  values[index_map[i]] = value_outer * values_1d[ix][0][0];
351  }
352  else
353  for (unsigned int iy = 0, i1 = 0; i1 < indices.size(); ++i1)
354  {
355  // prepare parts of products in y (and z) directions
356  std::array<double, dim + (dim * (dim - 1)) / 2> value_outer;
357  value_outer[0] = 1.;
358  if constexpr (dim > 1)
359  {
360  for (unsigned int x = 1; x < dim; ++x)
361  value_outer[0] *= values_1d[indices[i1][x - 1]][0][x];
362  for (unsigned int d = 1; d < dim; ++d)
363  {
364  value_outer[d] = values_1d[indices[i1][d - 1]][1][d];
365  for (unsigned int x = 1; x < dim; ++x)
366  if (x != d)
367  value_outer[d] *= values_1d[indices[i1][x - 1]][0][x];
368  }
369  for (unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
370  for (unsigned int d2 = d1; d2 < dim; ++d2, ++count)
371  {
372  value_outer[count] = 1.;
373  for (unsigned int x = 1; x < dim; ++x)
374  {
375  unsigned int derivative = 0;
376  if (d1 == x)
377  ++derivative;
378  if (d2 == x)
379  ++derivative;
380 
381  value_outer[count] *=
382  values_1d[indices[i1][x - 1]][derivative][x];
383  }
384  }
385  }
386 
387  // now run the loop over x and multiply by the values/derivatives
388  // in x direction
389  for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
390  {
391  std::array<double, 3> val_x{{values_1d[ix][0][0],
392  values_1d[ix][1][0],
393  values_1d[ix][2][0]}};
394  const unsigned int index =
395  (index_map.empty() ? i : index_map[i]);
396 
397  if (update_values)
398  values[index] = value_outer[0] * val_x[0];
399 
400  if (update_grads)
401  {
402  grads[index][0] = value_outer[0] * val_x[1];
403  if constexpr (dim > 1)
404  for (unsigned int d = 1; d < dim; ++d)
405  grads[index][d] = value_outer[d] * val_x[0];
406  }
407 
408  if (update_grad_grads)
409  {
410  grad_grads[index][0][0] = value_outer[0] * val_x[2];
411  if constexpr (dim > 1)
412  {
413  for (unsigned int d = 1; d < dim; ++d)
414  grad_grads[index][0][d] = grad_grads[index][d][0] =
415  value_outer[d] * val_x[1];
416  for (unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
417  for (unsigned int d2 = d1; d2 < dim; ++d2, ++count)
418  grad_grads[index][d1][d2] =
419  grad_grads[index][d2][d1] =
420  value_outer[count] * val_x[0];
421  }
422  }
423  }
424 
425  // Use slower code for 3rd and 4th derivatives
427  for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
428  {
429  const unsigned int index =
430  (index_map.empty() ? i : index_map[i]);
431  std::array<unsigned int, dim> my_indices;
432  my_indices[0] = ix;
433  if constexpr (dim > 1)
434  for (unsigned int d = 1; d < dim; ++d)
435  my_indices[d] = indices[i1][d - 1];
436  for (unsigned int d1 = 0; d1 < dim; ++d1)
437  for (unsigned int d2 = 0; d2 < dim; ++d2)
438  for (unsigned int d3 = 0; d3 < dim; ++d3)
439  {
440  double der3 = 1.;
441  for (unsigned int x = 0; x < dim; ++x)
442  {
443  unsigned int derivative = 0;
444  if (d1 == x)
445  ++derivative;
446  if (d2 == x)
447  ++derivative;
448  if (d3 == x)
449  ++derivative;
450 
451  der3 *= values_1d[my_indices[x]][derivative][x];
452  }
453  third_derivatives[index][d1][d2][d3] = der3;
454  }
455  }
456 
457  if (update_4th_derivatives)
458  for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
459  {
460  const unsigned int index =
461  (index_map.empty() ? i : index_map[i]);
462  std::array<unsigned int, dim> my_indices;
463  my_indices[0] = ix;
464  if constexpr (dim > 1)
465  for (unsigned int d = 1; d < dim; ++d)
466  my_indices[d] = indices[i1][d - 1];
467  for (unsigned int d1 = 0; d1 < dim; ++d1)
468  for (unsigned int d2 = 0; d2 < dim; ++d2)
469  for (unsigned int d3 = 0; d3 < dim; ++d3)
470  for (unsigned int d4 = 0; d4 < dim; ++d4)
471  {
472  double der4 = 1.;
473  for (unsigned int x = 0; x < dim; ++x)
474  {
475  unsigned int derivative = 0;
476  if (d1 == x)
477  ++derivative;
478  if (d2 == x)
479  ++derivative;
480  if (d3 == x)
481  ++derivative;
482  if (d4 == x)
483  ++derivative;
484 
485  der4 *= values_1d[my_indices[x]][derivative][x];
486  }
487  fourth_derivatives[index][d1][d2][d3][d4] = der4;
488  }
489  }
490 
491  iy += size_x;
492  }
493  }
494  } // namespace TensorProductPolynomials
495 } // namespace internal
496 
497 
498 
499 template <int dim, typename PolynomialType>
500 void
502  const Point<dim> &p,
503  std::vector<double> &values,
504  std::vector<Tensor<1, dim>> &grads,
505  std::vector<Tensor<2, dim>> &grad_grads,
506  std::vector<Tensor<3, dim>> &third_derivatives,
507  std::vector<Tensor<4, dim>> &fourth_derivatives) const
508 {
509  Assert(dim <= 3, ExcNotImplemented());
510  Assert(values.size() == this->n() || values.empty(),
511  ExcDimensionMismatch2(values.size(), this->n(), 0));
512  Assert(grads.size() == this->n() || grads.empty(),
513  ExcDimensionMismatch2(grads.size(), this->n(), 0));
514  Assert(grad_grads.size() == this->n() || grad_grads.empty(),
515  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
516  Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
517  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
518  Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
519  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
520 
521  // check how many values/derivatives we have to compute
522  unsigned int n_derivatives = 0;
523  if (values.size() == this->n())
524  n_derivatives = 0;
525  if (grads.size() == this->n())
526  n_derivatives = 1;
527  if (grad_grads.size() == this->n())
528  n_derivatives = 2;
529  if (third_derivatives.size() == this->n())
530  n_derivatives = 3;
531  if (fourth_derivatives.size() == this->n())
532  n_derivatives = 4;
533 
534  // Compute the values (and derivatives, if necessary) of all 1d
535  // polynomials at this evaluation point. We can use the more optimized
536  // values_of_array function to compute 'dim' polynomials at once
537  const unsigned int n_polynomials = polynomials.size();
538  boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
539  n_polynomials);
540  if constexpr (std::is_same<PolynomialType,
542  {
543  std::array<double, dim> point_array;
544  for (unsigned int d = 0; d < dim; ++d)
545  point_array[d] = p[d];
546  for (unsigned int i = 0; i < n_polynomials; ++i)
547  polynomials[i].values_of_array(point_array,
548  n_derivatives,
549  values_1d[i].data());
550  }
551  else
552  for (unsigned int i = 0; i < n_polynomials; ++i)
553  for (unsigned int d = 0; d < dim; ++d)
554  {
555  std::array<double, 5> derivatives;
556  polynomials[i].value(p[d], n_derivatives, derivatives.data());
557  for (unsigned int j = 0; j <= n_derivatives; ++j)
558  values_1d[i][j][d] = derivatives[j];
559  }
560 
561  // Unroll the tensor product indices of all but the first dimension in
562  // arbitrary dimension
563  constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
564  boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
565  if constexpr (dim > 1)
566  for (unsigned int d = 1; d < dim; ++d)
567  {
568  const unsigned int size = indices.size();
569  for (unsigned int i = 1; i < n_polynomials; ++i)
570  for (unsigned int j = 0; j < size; ++j)
571  {
572  std::array<unsigned int, dim1> next_index = indices[j];
573  next_index[d - 1] = i;
574  indices.push_back(next_index);
575  }
576  }
577  AssertDimension(indices.size(), Utilities::pow(n_polynomials, dim - 1));
578 
579  internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
580  n_derivatives,
581  values_1d,
582  n_polynomials,
583  indices,
584  index_map_inverse,
585  values,
586  grads,
587  grad_grads,
588  third_derivatives,
589  fourth_derivatives);
590 }
591 
592 
593 
594 template <>
595 void
597  const Point<0> &,
598  std::vector<double> &,
599  std::vector<Tensor<1, 0>> &,
600  std::vector<Tensor<2, 0>> &,
601  std::vector<Tensor<3, 0>> &,
602  std::vector<Tensor<4, 0>> &) const
603 {
604  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
605 }
606 
607 
608 
609 template <int dim, typename PolynomialType>
610 std::unique_ptr<ScalarPolynomialsBase<dim>>
612 {
613  return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
614 }
615 
616 
617 
618 template <int dim, typename PolynomialType>
619 std::size_t
621 {
622  return (MemoryConsumption::memory_consumption(polynomials) +
624  MemoryConsumption::memory_consumption(index_map_inverse));
625 }
626 
627 
628 
629 template <int dim, typename PolynomialType>
630 std::vector<PolynomialType>
632  const
633 {
634  return polynomials;
635 }
636 
637 
638 
639 /* ------------------- AnisotropicPolynomials -------------- */
640 
641 
642 template <int dim>
644  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
645  : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
646  , polynomials(pols)
647 {
648  Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
649  for (const auto &pols_d : pols)
650  {
651  (void)pols_d;
652  Assert(pols_d.size() > 0,
653  ExcMessage("The number of polynomials must be larger than zero "
654  "for all coordinate directions."));
655  }
656 }
657 
658 
659 
660 template <int dim>
661 void
663  const unsigned int i,
664  std::array<unsigned int, dim> &indices) const
665 {
666 #ifdef DEBUG
667  unsigned int n_poly = 1;
668  for (unsigned int d = 0; d < dim; ++d)
669  n_poly *= polynomials[d].size();
670  Assert(i < n_poly, ExcInternalError());
671 #endif
672 
673  if (dim == 0)
674  {
675  }
676  else if (dim == 1)
677  internal::compute_tensor_index(i,
678  polynomials[0].size(),
679  0 /*not used*/,
680  indices);
681  else
682  internal::compute_tensor_index(i,
683  polynomials[0].size(),
684  polynomials[1].size(),
685  indices);
686 }
687 
688 
689 
690 template <>
691 void
693  std::array<unsigned int, 0> &) const
694 {
695  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
696 }
697 
698 
699 
700 template <int dim>
701 double
703  const Point<dim> &p) const
704 {
705  std::array<unsigned int, dim> indices;
706  compute_index(i, indices);
707 
708  double value = 1.;
709  for (unsigned int d = 0; d < dim; ++d)
710  value *= polynomials[d][indices[d]].value(p(d));
711 
712  return value;
713 }
714 
715 
716 
717 template <>
718 double
720  const Point<0> &) const
721 {
722  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
723 
724  return {};
725 }
726 
727 
728 
729 template <int dim>
732  const Point<dim> &p) const
733 {
734  std::array<unsigned int, dim> indices;
735  compute_index(i, indices);
736 
737  // compute values and
738  // uni-directional derivatives at
739  // the given point in each
740  // coordinate direction
742  for (unsigned int d = 0; d < dim; ++d)
743  polynomials[d][indices[d]].value(p(d), 1, v[d].data());
744 
745  Tensor<1, dim> grad;
746  for (unsigned int d = 0; d < dim; ++d)
747  {
748  grad[d] = 1.;
749  for (unsigned int x = 0; x < dim; ++x)
750  grad[d] *= v[x][d == x];
751  }
752 
753  return grad;
754 }
755 
756 
757 
758 template <>
761  const Point<0> &) const
762 {
763  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
764 
765  return {};
766 }
767 
768 
769 
770 template <int dim>
773  const Point<dim> &p) const
774 {
775  std::array<unsigned int, dim> indices;
776  compute_index(i, indices);
777 
779  for (unsigned int d = 0; d < dim; ++d)
780  polynomials[d][indices[d]].value(p(d), 2, v[d].data());
781 
782  Tensor<2, dim> grad_grad;
783  for (unsigned int d1 = 0; d1 < dim; ++d1)
784  for (unsigned int d2 = 0; d2 < dim; ++d2)
785  {
786  grad_grad[d1][d2] = 1.;
787  for (unsigned int x = 0; x < dim; ++x)
788  {
789  unsigned int derivative = 0;
790  if (d1 == x || d2 == x)
791  {
792  if (d1 == d2)
793  derivative = 2;
794  else
795  derivative = 1;
796  }
797  grad_grad[d1][d2] *= v[x][derivative];
798  }
799  }
800 
801  return grad_grad;
802 }
803 
804 
805 
806 template <>
809  const Point<0> &) const
810 {
811  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
812 
813  return {};
814 }
815 
816 
817 
818 template <int dim>
819 void
821  const Point<dim> &p,
822  std::vector<double> &values,
823  std::vector<Tensor<1, dim>> &grads,
824  std::vector<Tensor<2, dim>> &grad_grads,
825  std::vector<Tensor<3, dim>> &third_derivatives,
826  std::vector<Tensor<4, dim>> &fourth_derivatives) const
827 {
828  Assert(values.size() == this->n() || values.empty(),
829  ExcDimensionMismatch2(values.size(), this->n(), 0));
830  Assert(grads.size() == this->n() || grads.empty(),
831  ExcDimensionMismatch2(grads.size(), this->n(), 0));
832  Assert(grad_grads.size() == this->n() || grad_grads.empty(),
833  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
834  Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
835  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
836  Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
837  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
838 
839  // check how many values/derivatives we have to compute
840  unsigned int n_derivatives = 0;
841  if (values.size() == this->n())
842  n_derivatives = 0;
843  if (grads.size() == this->n())
844  n_derivatives = 1;
845  if (grad_grads.size() == this->n())
846  n_derivatives = 2;
847  if (third_derivatives.size() == this->n())
848  n_derivatives = 3;
849  if (fourth_derivatives.size() == this->n())
850  n_derivatives = 4;
851 
852  // compute the values (and derivatives, if necessary) of all polynomials at
853  // this evaluation point
854  std::size_t max_n_polynomials = 0;
855  for (unsigned int d = 0; d < dim; ++d)
856  max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
857 
858  // 5 is enough to store values and derivatives in all supported cases
859  boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
860  max_n_polynomials);
861  if (n_derivatives == 0)
862  for (unsigned int d = 0; d < dim; ++d)
863  for (unsigned int i = 0; i < polynomials[d].size(); ++i)
864  values_1d[i][0][d] = polynomials[d][i].value(p[d]);
865  else
866  for (unsigned int d = 0; d < dim; ++d)
867  for (unsigned int i = 0; i < polynomials[d].size(); ++i)
868  {
869  // The isotropic tensor product function wants us to use a different
870  // innermost index, so we cannot pass the values_1d array into the
871  // function directly
872  std::array<double, 5> derivatives;
873  polynomials[d][i].value(p[d], n_derivatives, derivatives.data());
874  for (unsigned int j = 0; j <= n_derivatives; ++j)
875  values_1d[i][j][d] = derivatives[j];
876  }
877 
878  // Unroll the tensor product indices in arbitrary dimension
879  constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
880  boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
881  for (unsigned int d = 1; d < dim; ++d)
882  {
883  const unsigned int size = indices.size();
884  for (unsigned int i = 1; i < polynomials[d].size(); ++i)
885  for (unsigned int j = 0; j < size; ++j)
886  {
887  std::array<unsigned int, dim1> next_index = indices[j];
888  next_index[d - 1] = i;
889  indices.push_back(next_index);
890  }
891  }
892 
893  internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
894  n_derivatives,
895  values_1d,
896  polynomials[0].size(),
897  indices,
898  {},
899  values,
900  grads,
901  grad_grads,
902  third_derivatives,
903  fourth_derivatives);
904 }
905 
906 
907 
908 template <>
909 void
911  std::vector<double> &,
912  std::vector<Tensor<1, 0>> &,
913  std::vector<Tensor<2, 0>> &,
914  std::vector<Tensor<3, 0>> &,
915  std::vector<Tensor<4, 0>> &) const
916 {
917  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
918 }
919 
920 
921 
922 template <int dim>
923 unsigned int
925  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
926 {
927  unsigned int y = 1;
928  for (unsigned int d = 0; d < dim; ++d)
929  y *= pols[d].size();
930  return y;
931 }
932 
933 
934 
935 template <>
936 unsigned int
938  const std::vector<std::vector<Polynomials::Polynomial<double>>> &)
939 {
940  AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
941 
942  return {};
943 }
944 
945 
946 
947 template <int dim>
948 std::unique_ptr<ScalarPolynomialsBase<dim>>
950 {
951  return std::make_unique<AnisotropicPolynomials<dim>>(*this);
952 }
953 
954 
955 
956 /* ------------------- explicit instantiations -------------- */
961 
962 template class TensorProductPolynomials<
963  1,
965 template class TensorProductPolynomials<
966  2,
968 template class TensorProductPolynomials<
969  3,
971 
972 template class AnisotropicPolynomials<0>;
973 template class AnisotropicPolynomials<1>;
974 template class AnisotropicPolynomials<2>;
975 template class AnisotropicPolynomials<3>;
976 
double compute_value(const unsigned int i, const Point< dim > &p) const override
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
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
Definition: point.h:112
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
void output_indices(std::ostream &out) const
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
double compute_value(const unsigned int i, const Point< dim > &p) const override
virtual std::size_t memory_consumption() const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
std::vector< PolynomialType > get_underlying_polynomials() const
void set_numbering(const std::vector< unsigned int > &renumber)
Tensor< 2, dim > compute_grad_grad(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
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
void evaluate_tensor_product(const unsigned int n_derivatives, const boost::container::small_vector<::ndarray< double, 5, dim >, 10 > &values_1d, const unsigned int size_x, const boost::container::small_vector< std::array< unsigned int, dim1 >, 64 > &indices, const std::vector< unsigned int > &index_map, 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)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108