deal.II version GIT relicensing-2901-g19332422bd 2025-03-23 19:50:00+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 if constexpr (dim == 0)
88 {
89 (void)i;
90 (void)indices;
92 }
93 else
94 {
95 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
97 internal::compute_tensor_index(index_map[i],
98 polynomials.size(),
99 polynomials.size(),
100 indices);
101 }
102}
103
105
106template <int dim, typename PolynomialType>
107void
109 std::ostream &out) const
110{
111 if constexpr (dim == 0)
112 {
113 (void)out;
115 }
116 else
117 {
118 std::array<unsigned int, dim> ix;
119 for (unsigned int i = 0; i < this->n(); ++i)
120 {
121 compute_index(i, ix);
122 out << i << "\t";
123 for (unsigned int d = 0; d < dim; ++d)
124 out << ix[d] << " ";
125 out << std::endl;
126 }
127 }
128}
129
130
132template <int dim>
133inline const std::vector<unsigned int> &
135{
136 return index_map;
137}
138
139
140
141template <int dim>
142inline const std::vector<unsigned int> &
144{
145 return index_map_inverse;
146}
147
148
149
150template <int dim>
151void
153 const std::vector<unsigned int> &renumber)
154{
155 Assert(renumber.size() == index_map.size(),
156 ExcDimensionMismatch(renumber.size(), index_map.size()));
157
158 index_map = renumber;
159 for (unsigned int i = 0; i < index_map.size(); ++i)
160 index_map_inverse[index_map[i]] = i;
161}
162
163
164
165template <>
166void
167AnisotropicPolynomials<0>::set_numbering(const std::vector<unsigned int> &)
168{
169 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
170}
171
172
173
174template <int dim, typename PolynomialType>
175void
177 const std::vector<unsigned int> &renumber)
178{
179 Assert(renumber.size() == index_map.size(),
180 ExcDimensionMismatch(renumber.size(), index_map.size()));
181
182 index_map = renumber;
183 for (unsigned int i = 0; i < index_map.size(); ++i)
184 index_map_inverse[index_map[i]] = i;
185}
186
187
188
189template <>
190void
192 const std::vector<unsigned int> &)
193{
194 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
195}
196
197
198
199template <int dim, typename PolynomialType>
200double
202 const unsigned int i,
203 const Point<dim> &p) const
204{
205 if constexpr (dim == 0)
206 {
207 (void)i;
208 (void)p;
210 return 0;
211 }
212 else
213 {
214 std::array<unsigned int, dim> indices;
215 compute_index(i, indices);
216
217 double value = 1.;
218 for (unsigned int d = 0; d < dim; ++d)
219 value *= polynomials[indices[d]].value(p[d]);
220
221 return value;
222 }
223}
224
225
226
227template <int dim, typename PolynomialType>
230 const unsigned int i,
231 const Point<dim> &p) const
232{
233 if constexpr (dim == 0)
234 {
235 (void)i;
236 (void)p;
238 return {};
240 else
241 {
242 std::array<unsigned int, dim> indices;
243 compute_index(i, indices);
244
245 // compute values and
246 // uni-directional derivatives at
247 // the given point in each
248 // coordinate direction
250 {
251 std::vector<double> tmp(2);
252 for (unsigned int d = 0; d < dim; ++d)
253 {
254 polynomials[indices[d]].value(p[d], tmp);
255 v[d][0] = tmp[0];
256 v[d][1] = tmp[1];
257 }
258 }
259
260 Tensor<1, dim> grad;
261 for (unsigned int d = 0; d < dim; ++d)
262 {
263 grad[d] = 1.;
264 for (unsigned int x = 0; x < dim; ++x)
265 grad[d] *= v[x][d == x];
266 }
267
268 return grad;
269 }
270}
271
272
273
274template <int dim, typename PolynomialType>
277 const unsigned int i,
278 const Point<dim> &p) const
279{
280 if constexpr (dim == 0)
281 {
282 (void)i;
283 (void)p;
285 return {};
286 }
287 else
288 {
289 std::array<unsigned int, dim> indices;
290 compute_index(i, indices);
291
293 {
294 std::vector<double> tmp(3);
295 for (unsigned int d = 0; d < dim; ++d)
296 {
297 polynomials[indices[d]].value(p[d], tmp);
298 v[d][0] = tmp[0];
299 v[d][1] = tmp[1];
300 v[d][2] = tmp[2];
301 }
302 }
303
304 Tensor<2, dim> grad_grad;
305 for (unsigned int d1 = 0; d1 < dim; ++d1)
306 for (unsigned int d2 = 0; d2 < dim; ++d2)
307 {
308 grad_grad[d1][d2] = 1.;
309 for (unsigned int x = 0; x < dim; ++x)
310 {
311 unsigned int derivative = 0;
312 if (d1 == x || d2 == x)
313 {
314 if (d1 == d2)
315 derivative = 2;
316 else
317 derivative = 1;
318 }
319 grad_grad[d1][d2] *= v[x][derivative];
320 }
321 }
322
323 return grad_grad;
324 }
325}
326
327
328
329namespace internal
330{
332 {
333 // This function computes the tensor product of some tabulated
334 // one-dimensional polynomials (also the anisotropic case is supported)
335 // with tensor product indices of all dimensions except the first one
336 // tabulated in the 'indices' array; the first dimension is manually
337 // iterated through because these are possibly performance-critical loops,
338 // so we want to avoid indirect addressing.
339 template <int dim, std::size_t dim1>
340 void
342 const unsigned int n_derivatives,
343 const boost::container::small_vector<::ndarray<double, 5, dim>, 10>
344 &values_1d,
345 const unsigned int size_x,
346 const boost::container::small_vector<std::array<unsigned int, dim1>, 64>
347 &indices,
348 const std::vector<unsigned int> &index_map,
349 std::vector<double> &values,
350 std::vector<Tensor<1, dim>> &grads,
351 std::vector<Tensor<2, dim>> &grad_grads,
352 std::vector<Tensor<3, dim>> &third_derivatives,
353 std::vector<Tensor<4, dim>> &fourth_derivatives)
354 {
355 const bool update_values = (values.size() == indices.size() * size_x);
356 const bool update_grads = (grads.size() == indices.size() * size_x);
357 const bool update_grad_grads =
358 (grad_grads.size() == indices.size() * size_x);
359 const bool update_3rd_derivatives =
360 (third_derivatives.size() == indices.size() * size_x);
361 const bool update_4th_derivatives =
362 (fourth_derivatives.size() == indices.size() * size_x);
363
364 // For values, 1st and 2nd derivatives use a more lengthy code that
365 // minimizes the number of arithmetic operations and memory accesses
366 if (n_derivatives == 0)
367 for (unsigned int i = 0, i1 = 0; i1 < indices.size(); ++i1)
368 {
369 double value_outer = 1.;
370 if constexpr (dim > 1)
371 for (unsigned int d = 1; d < dim; ++d)
372 value_outer *= values_1d[indices[i1][d - 1]][0][d];
373 if (index_map.empty())
374 for (unsigned int ix = 0; ix < size_x; ++ix, ++i)
375 values[i] = value_outer * values_1d[ix][0][0];
376 else
377 for (unsigned int ix = 0; ix < size_x; ++ix, ++i)
378 values[index_map[i]] = value_outer * values_1d[ix][0][0];
379 }
380 else
381 for (unsigned int iy = 0, i1 = 0; i1 < indices.size(); ++i1)
382 {
383 // prepare parts of products in y (and z) directions
384 std::array<double, dim + (dim * (dim - 1)) / 2> value_outer;
385 value_outer[0] = 1.;
386 if constexpr (dim > 1)
387 {
388 for (unsigned int x = 1; x < dim; ++x)
389 value_outer[0] *= values_1d[indices[i1][x - 1]][0][x];
390 for (unsigned int d = 1; d < dim; ++d)
391 {
392 value_outer[d] = values_1d[indices[i1][d - 1]][1][d];
393 for (unsigned int x = 1; x < dim; ++x)
394 if (x != d)
395 value_outer[d] *= values_1d[indices[i1][x - 1]][0][x];
396 }
397 for (unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
398 for (unsigned int d2 = d1; d2 < dim; ++d2, ++count)
399 {
400 value_outer[count] = 1.;
401 for (unsigned int x = 1; x < dim; ++x)
402 {
403 unsigned int derivative = 0;
404 if (d1 == x)
405 ++derivative;
406 if (d2 == x)
407 ++derivative;
408
409 value_outer[count] *=
410 values_1d[indices[i1][x - 1]][derivative][x];
411 }
412 }
413 }
414
415 // now run the loop over x and multiply by the values/derivatives
416 // in x direction
417 for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
418 {
419 std::array<double, 3> val_x{{values_1d[ix][0][0],
420 values_1d[ix][1][0],
421 values_1d[ix][2][0]}};
422 const unsigned int index =
423 (index_map.empty() ? i : index_map[i]);
424
425 if (update_values)
426 values[index] = value_outer[0] * val_x[0];
427
428 if (update_grads)
429 {
430 grads[index][0] = value_outer[0] * val_x[1];
431 if constexpr (dim > 1)
432 for (unsigned int d = 1; d < dim; ++d)
433 grads[index][d] = value_outer[d] * val_x[0];
434 }
435
436 if (update_grad_grads)
437 {
438 grad_grads[index][0][0] = value_outer[0] * val_x[2];
439 if constexpr (dim > 1)
440 {
441 for (unsigned int d = 1; d < dim; ++d)
442 grad_grads[index][0][d] = grad_grads[index][d][0] =
443 value_outer[d] * val_x[1];
444 for (unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
445 for (unsigned int d2 = d1; d2 < dim; ++d2, ++count)
446 grad_grads[index][d1][d2] =
447 grad_grads[index][d2][d1] =
448 value_outer[count] * val_x[0];
449 }
450 }
451 }
452
453 // Use slower code for 3rd and 4th derivatives
455 for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
456 {
457 const unsigned int index =
458 (index_map.empty() ? i : index_map[i]);
459 std::array<unsigned int, dim> my_indices;
460 my_indices[0] = ix;
461 if constexpr (dim > 1)
462 for (unsigned int d = 1; d < dim; ++d)
463 my_indices[d] = indices[i1][d - 1];
464 for (unsigned int d1 = 0; d1 < dim; ++d1)
465 for (unsigned int d2 = 0; d2 < dim; ++d2)
466 for (unsigned int d3 = 0; d3 < dim; ++d3)
467 {
468 double der3 = 1.;
469 for (unsigned int x = 0; x < dim; ++x)
470 {
471 unsigned int derivative = 0;
472 if (d1 == x)
473 ++derivative;
474 if (d2 == x)
475 ++derivative;
476 if (d3 == x)
477 ++derivative;
478
479 der3 *= values_1d[my_indices[x]][derivative][x];
480 }
481 third_derivatives[index][d1][d2][d3] = der3;
482 }
483 }
484
485 if (update_4th_derivatives)
486 for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
487 {
488 const unsigned int index =
489 (index_map.empty() ? i : index_map[i]);
490 std::array<unsigned int, dim> my_indices;
491 my_indices[0] = ix;
492 if constexpr (dim > 1)
493 for (unsigned int d = 1; d < dim; ++d)
494 my_indices[d] = indices[i1][d - 1];
495 for (unsigned int d1 = 0; d1 < dim; ++d1)
496 for (unsigned int d2 = 0; d2 < dim; ++d2)
497 for (unsigned int d3 = 0; d3 < dim; ++d3)
498 for (unsigned int d4 = 0; d4 < dim; ++d4)
499 {
500 double der4 = 1.;
501 for (unsigned int x = 0; x < dim; ++x)
502 {
503 unsigned int derivative = 0;
504 if (d1 == x)
505 ++derivative;
506 if (d2 == x)
507 ++derivative;
508 if (d3 == x)
509 ++derivative;
510 if (d4 == x)
511 ++derivative;
512
513 der4 *= values_1d[my_indices[x]][derivative][x];
514 }
515 fourth_derivatives[index][d1][d2][d3][d4] = der4;
516 }
517 }
518
519 iy += size_x;
520 }
521 }
522 } // namespace TensorProductPolynomials
523} // namespace internal
524
525
526
527template <int dim, typename PolynomialType>
528void
530 const Point<dim> &p,
531 std::vector<double> &values,
532 std::vector<Tensor<1, dim>> &grads,
533 std::vector<Tensor<2, dim>> &grad_grads,
534 std::vector<Tensor<3, dim>> &third_derivatives,
535 std::vector<Tensor<4, dim>> &fourth_derivatives) const
536{
537 if constexpr (dim == 0)
538 {
539 (void)p;
540 (void)values;
541 (void)grads;
542 (void)grad_grads;
543 (void)third_derivatives;
544 (void)fourth_derivatives;
546 }
547 else
548 {
549 Assert(dim <= 3, ExcNotImplemented());
550 Assert(values.size() == this->n() || values.empty(),
551 ExcDimensionMismatch2(values.size(), this->n(), 0));
552 Assert(grads.size() == this->n() || grads.empty(),
553 ExcDimensionMismatch2(grads.size(), this->n(), 0));
554 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
555 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
556 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
557 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
558 Assert(fourth_derivatives.size() == this->n() ||
559 fourth_derivatives.empty(),
560 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
561
562 // check how many values/derivatives we have to compute
563 unsigned int n_derivatives = 0;
564 if (values.size() == this->n())
565 n_derivatives = 0;
566 if (grads.size() == this->n())
567 n_derivatives = 1;
568 if (grad_grads.size() == this->n())
569 n_derivatives = 2;
570 if (third_derivatives.size() == this->n())
571 n_derivatives = 3;
572 if (fourth_derivatives.size() == this->n())
573 n_derivatives = 4;
574
575 // Compute the values (and derivatives, if necessary) of all 1d
576 // polynomials at this evaluation point. We can use the more optimized
577 // values_of_array function to compute 'dim' polynomials at once
578 const unsigned int n_polynomials = polynomials.size();
579 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
580 n_polynomials);
581 if constexpr (std::is_same_v<PolynomialType,
583 {
584 std::array<double, dim> point_array;
585 for (unsigned int d = 0; d < dim; ++d)
586 point_array[d] = p[d];
587 for (unsigned int i = 0; i < n_polynomials; ++i)
588 polynomials[i].values_of_array(point_array,
589 n_derivatives,
590 values_1d[i].data());
591 }
592 else
593 for (unsigned int i = 0; i < n_polynomials; ++i)
594 for (unsigned int d = 0; d < dim; ++d)
595 {
596 std::array<double, 5> derivatives;
597 polynomials[i].value(p[d], n_derivatives, derivatives.data());
598 for (unsigned int j = 0; j <= n_derivatives; ++j)
599 values_1d[i][j][d] = derivatives[j];
600 }
601
602 // Unroll the tensor product indices of all but the first dimension in
603 // arbitrary dimension
604 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
605 boost::container::small_vector<std::array<unsigned int, dim1>, 64>
606 indices(1);
607 if constexpr (dim > 1)
608 for (unsigned int d = 1; d < dim; ++d)
609 {
610 const unsigned int size = indices.size();
611 for (unsigned int i = 1; i < n_polynomials; ++i)
612 for (unsigned int j = 0; j < size; ++j)
613 {
614 std::array<unsigned int, dim1> next_index = indices[j];
615 next_index[d - 1] = i;
616 indices.push_back(next_index);
617 }
618 }
619 AssertDimension(indices.size(), Utilities::pow(n_polynomials, dim - 1));
620
621 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
622 n_derivatives,
623 values_1d,
624 n_polynomials,
625 indices,
626 index_map_inverse,
627 values,
628 grads,
629 grad_grads,
630 third_derivatives,
631 fourth_derivatives);
632 }
633}
634
635
636
637template <int dim, typename PolynomialType>
638std::unique_ptr<ScalarPolynomialsBase<dim>>
640{
641 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
642}
643
644
645
646template <int dim, typename PolynomialType>
647std::size_t
654
655
656
657template <int dim, typename PolynomialType>
658std::vector<PolynomialType>
664
665
666
667/* ------------------- AnisotropicPolynomials -------------- */
668
669
670template <int dim>
672 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
673 : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
674 , polynomials(pols)
675 , index_map(this->n())
676 , index_map_inverse(this->n())
677{
678 Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
679 for (const auto &pols_d : pols)
680 {
681 (void)pols_d;
682 Assert(pols_d.size() > 0,
683 ExcMessage("The number of polynomials must be larger than zero "
684 "for all coordinate directions."));
685 }
686
687 // per default set this index map to identity. This map can be changed by
688 // the user through the set_numbering() function
689 for (unsigned int i = 0; i < this->n(); ++i)
690 {
691 index_map[i] = i;
692 index_map_inverse[i] = i;
693 }
694}
695
696
697
698template <int dim>
699void
701 const unsigned int i,
702 std::array<unsigned int, dim> &indices) const
703{
704 if constexpr (dim == 0)
705 {
706 (void)i;
707 (void)indices;
709 }
710 else
711 {
712 if constexpr (running_in_debug_mode())
713 {
714 unsigned int n_poly = 1;
715 for (unsigned int d = 0; d < dim; ++d)
716 n_poly *= polynomials[d].size();
717 Assert(i < n_poly, ExcInternalError());
718 }
719
720 if (dim == 0)
721 {
722 }
723 else if (dim == 1)
724 internal::compute_tensor_index(index_map[i],
725 polynomials[0].size(),
726 0 /*not used*/,
727 indices);
728 else
729 internal::compute_tensor_index(index_map[i],
730 polynomials[0].size(),
731 polynomials[1].size(),
732 indices);
733 }
734}
735
736
737
738template <int dim>
739double
741 const Point<dim> &p) const
742{
743 if constexpr (dim == 0)
744 {
745 (void)i;
746 (void)p;
748 return {};
749 }
750 else
751 {
752 std::array<unsigned int, dim> indices;
753 compute_index(i, indices);
754
755 double value = 1.;
756 for (unsigned int d = 0; d < dim; ++d)
757 value *= polynomials[d][indices[d]].value(p[d]);
758
759 return value;
760 }
761}
762
763
764
765template <int dim>
768 const Point<dim> &p) const
769{
770 if constexpr (dim == 0)
771 {
772 (void)i;
773 (void)p;
775 return {};
776 }
777 else
778 {
779 std::array<unsigned int, dim> indices;
780 compute_index(i, indices);
781
782 // compute values and
783 // uni-directional derivatives at
784 // the given point in each
785 // coordinate direction
787 for (unsigned int d = 0; d < dim; ++d)
788 polynomials[d][indices[d]].value(p[d], 1, v[d].data());
789
790 Tensor<1, dim> grad;
791 for (unsigned int d = 0; d < dim; ++d)
792 {
793 grad[d] = 1.;
794 for (unsigned int x = 0; x < dim; ++x)
795 grad[d] *= v[x][d == x];
796 }
797
798 return grad;
799 }
800}
801
802
803
804template <int dim>
807 const Point<dim> &p) const
808{
809 if constexpr (dim == 0)
810 {
811 (void)i;
812 (void)p;
814 return {};
815 }
816 else
817 {
818 std::array<unsigned int, dim> indices;
819 compute_index(i, indices);
820
822 for (unsigned int d = 0; d < dim; ++d)
823 polynomials[d][indices[d]].value(p[d], 2, v[d].data());
824
825 Tensor<2, dim> grad_grad;
826 for (unsigned int d1 = 0; d1 < dim; ++d1)
827 for (unsigned int d2 = 0; d2 < dim; ++d2)
828 {
829 grad_grad[d1][d2] = 1.;
830 for (unsigned int x = 0; x < dim; ++x)
831 {
832 unsigned int derivative = 0;
833 if (d1 == x || d2 == x)
834 {
835 if (d1 == d2)
836 derivative = 2;
837 else
838 derivative = 1;
839 }
840 grad_grad[d1][d2] *= v[x][derivative];
841 }
842 }
843
844 return grad_grad;
845 }
846}
847
848
849
850template <int dim>
851void
853 const Point<dim> &p,
854 std::vector<double> &values,
855 std::vector<Tensor<1, dim>> &grads,
856 std::vector<Tensor<2, dim>> &grad_grads,
857 std::vector<Tensor<3, dim>> &third_derivatives,
858 std::vector<Tensor<4, dim>> &fourth_derivatives) const
859{
860 if constexpr (dim == 0)
861 {
862 (void)p;
863 (void)values;
864 (void)grads;
865 (void)grad_grads;
866 (void)third_derivatives;
867 (void)fourth_derivatives;
869 }
870 else
871 {
872 Assert(values.size() == this->n() || values.empty(),
873 ExcDimensionMismatch2(values.size(), this->n(), 0));
874 Assert(grads.size() == this->n() || grads.empty(),
875 ExcDimensionMismatch2(grads.size(), this->n(), 0));
876 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
877 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
878 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
879 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
880 Assert(fourth_derivatives.size() == this->n() ||
881 fourth_derivatives.empty(),
882 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
883
884 // check how many values/derivatives we have to compute
885 unsigned int n_derivatives = 0;
886 if (values.size() == this->n())
887 n_derivatives = 0;
888 if (grads.size() == this->n())
889 n_derivatives = 1;
890 if (grad_grads.size() == this->n())
891 n_derivatives = 2;
892 if (third_derivatives.size() == this->n())
893 n_derivatives = 3;
894 if (fourth_derivatives.size() == this->n())
895 n_derivatives = 4;
896
897 // compute the values (and derivatives, if necessary) of all polynomials
898 // at this evaluation point
899 std::size_t max_n_polynomials = 0;
900 for (unsigned int d = 0; d < dim; ++d)
901 max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
902
903 // 5 is enough to store values and derivatives in all supported cases
904 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
905 max_n_polynomials);
906 if (n_derivatives == 0)
907 for (unsigned int d = 0; d < dim; ++d)
908 for (unsigned int i = 0; i < polynomials[d].size(); ++i)
909 values_1d[i][0][d] = polynomials[d][i].value(p[d]);
910 else
911 for (unsigned int d = 0; d < dim; ++d)
912 for (unsigned int i = 0; i < polynomials[d].size(); ++i)
913 {
914 // The isotropic tensor product function wants us to use a
915 // different innermost index, so we cannot pass the values_1d
916 // array into the function directly
917 std::array<double, 5> derivatives;
918 polynomials[d][i].value(p[d], n_derivatives, derivatives.data());
919 for (unsigned int j = 0; j <= n_derivatives; ++j)
920 values_1d[i][j][d] = derivatives[j];
921 }
922
923 // Unroll the tensor product indices in arbitrary dimension
924 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
925 boost::container::small_vector<std::array<unsigned int, dim1>, 64>
926 indices(1);
927 for (unsigned int d = 1; d < dim; ++d)
928 {
929 const unsigned int size = indices.size();
930 for (unsigned int i = 1; i < polynomials[d].size(); ++i)
931 for (unsigned int j = 0; j < size; ++j)
932 {
933 std::array<unsigned int, dim1> next_index = indices[j];
934 next_index[d - 1] = i;
935 indices.push_back(next_index);
936 }
937 }
938
939 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
940 n_derivatives,
941 values_1d,
942 polynomials[0].size(),
943 indices,
944 index_map_inverse,
945 values,
946 grads,
947 grad_grads,
948 third_derivatives,
949 fourth_derivatives);
950 }
951}
952
953
954
955template <int dim>
956unsigned int
958 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
959{
960 if constexpr (dim == 0)
961 {
962 (void)pols;
964 return {};
965 }
966 else
967 {
968 unsigned int y = 1;
969 for (unsigned int d = 0; d < dim; ++d)
970 y *= pols[d].size();
971 return y;
972 }
973}
974
975
976
977template <int dim>
978std::unique_ptr<ScalarPolynomialsBase<dim>>
980{
981 return std::make_unique<AnisotropicPolynomials<dim>>(*this);
982}
983
984
985
986/* ------------------- explicit instantiations -------------- */
991
992template class TensorProductPolynomials<
993 1,
995template class TensorProductPolynomials<
996 2,
998template class TensorProductPolynomials<
999 3,
1001
1002template class AnisotropicPolynomials<0>;
1003template class AnisotropicPolynomials<1>;
1004template class AnisotropicPolynomials<2>;
1005template class AnisotropicPolynomials<3>;
1006
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
void set_numbering(const std::vector< unsigned int > &renumber)
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)
std::vector< unsigned int > index_map
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
const std::vector< unsigned int > & get_numbering() const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map_inverse
const std::vector< unsigned int > & get_numbering_inverse() const
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Definition point.h:113
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:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
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.
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
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:967
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