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