Reference documentation for deal.II version 9.6.0
\(\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>
137inline const std::vector<unsigned int> &
139{
140 return index_map;
141}
142
143
144
145template <int dim>
146inline const std::vector<unsigned int> &
148{
149 return index_map_inverse;
150}
152
153
154template <int dim>
155void
157 const std::vector<unsigned int> &renumber)
158{
159 Assert(renumber.size() == index_map.size(),
160 ExcDimensionMismatch(renumber.size(), index_map.size()));
161
162 index_map = renumber;
163 for (unsigned int i = 0; i < index_map.size(); ++i)
164 index_map_inverse[index_map[i]] = i;
165}
166
167
168
169template <>
170void
171AnisotropicPolynomials<0>::set_numbering(const std::vector<unsigned int> &)
172{
173 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
174}
175
176
177
178template <int dim, typename PolynomialType>
179void
181 const std::vector<unsigned int> &renumber)
182{
183 Assert(renumber.size() == index_map.size(),
184 ExcDimensionMismatch(renumber.size(), index_map.size()));
185
186 index_map = renumber;
187 for (unsigned int i = 0; i < index_map.size(); ++i)
188 index_map_inverse[index_map[i]] = i;
189}
190
191
192
193template <>
194void
196 const std::vector<unsigned int> &)
197{
198 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
199}
200
201
202
203template <int dim, typename PolynomialType>
204double
206 const unsigned int i,
207 const Point<dim> &p) const
208{
209 Assert(dim > 0, ExcNotImplemented());
210
211 std::array<unsigned int, dim> indices;
212 compute_index(i, indices);
213
214 double value = 1.;
215 for (unsigned int d = 0; d < dim; ++d)
216 value *= polynomials[indices[d]].value(p[d]);
217
218 return value;
219}
220
221
222
223template <>
224double
226 const unsigned int,
227 const Point<0> &) const
228{
229 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
230
231 return {};
232}
233
234
235
236template <int dim, typename PolynomialType>
239 const unsigned int i,
240 const Point<dim> &p) const
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
273template <>
276 const unsigned int,
277 const Point<0> &) const
278{
279 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
280
281 return {};
282}
283
284
285
286template <int dim, typename PolynomialType>
289 const unsigned int i,
290 const Point<dim> &p) const
291{
292 std::array<unsigned int, dim> indices;
293 compute_index(i, indices);
294
296 {
297 std::vector<double> tmp(3);
298 for (unsigned int d = 0; d < dim; ++d)
299 {
300 polynomials[indices[d]].value(p[d], tmp);
301 v[d][0] = tmp[0];
302 v[d][1] = tmp[1];
303 v[d][2] = tmp[2];
304 }
305 }
306
307 Tensor<2, dim> grad_grad;
308 for (unsigned int d1 = 0; d1 < dim; ++d1)
309 for (unsigned int d2 = 0; d2 < dim; ++d2)
310 {
311 grad_grad[d1][d2] = 1.;
312 for (unsigned int x = 0; x < dim; ++x)
313 {
314 unsigned int derivative = 0;
315 if (d1 == x || d2 == x)
316 {
317 if (d1 == d2)
318 derivative = 2;
319 else
320 derivative = 1;
321 }
322 grad_grad[d1][d2] *= v[x][derivative];
323 }
324 }
325
326 return grad_grad;
327}
328
329
330
331template <>
334 const unsigned int,
335 const Point<0> &) const
336{
337 return {};
338}
339
340
341
342namespace internal
343{
345 {
346 // This function computes the tensor product of some tabulated
347 // one-dimensional polynomials (also the anisotropic case is supported)
348 // with tensor product indices of all dimensions except the first one
349 // tabulated in the 'indices' array; the first dimension is manually
350 // iterated through because these are possibly performance-critical loops,
351 // so we want to avoid indirect addressing.
352 template <int dim, std::size_t dim1>
353 void
355 const unsigned int n_derivatives,
356 const boost::container::small_vector<::ndarray<double, 5, dim>, 10>
357 &values_1d,
358 const unsigned int size_x,
359 const boost::container::small_vector<std::array<unsigned int, dim1>, 64>
360 &indices,
361 const std::vector<unsigned int> &index_map,
362 std::vector<double> &values,
363 std::vector<Tensor<1, dim>> &grads,
364 std::vector<Tensor<2, dim>> &grad_grads,
365 std::vector<Tensor<3, dim>> &third_derivatives,
366 std::vector<Tensor<4, dim>> &fourth_derivatives)
367 {
368 const bool update_values = (values.size() == indices.size() * size_x);
369 const bool update_grads = (grads.size() == indices.size() * size_x);
370 const bool update_grad_grads =
371 (grad_grads.size() == indices.size() * size_x);
372 const bool update_3rd_derivatives =
373 (third_derivatives.size() == indices.size() * size_x);
374 const bool update_4th_derivatives =
375 (fourth_derivatives.size() == indices.size() * size_x);
376
377 // For values, 1st and 2nd derivatives use a more lengthy code that
378 // minimizes the number of arithmetic operations and memory accesses
379 if (n_derivatives == 0)
380 for (unsigned int i = 0, i1 = 0; i1 < indices.size(); ++i1)
381 {
382 double value_outer = 1.;
383 if constexpr (dim > 1)
384 for (unsigned int d = 1; d < dim; ++d)
385 value_outer *= values_1d[indices[i1][d - 1]][0][d];
386 if (index_map.empty())
387 for (unsigned int ix = 0; ix < size_x; ++ix, ++i)
388 values[i] = value_outer * values_1d[ix][0][0];
389 else
390 for (unsigned int ix = 0; ix < size_x; ++ix, ++i)
391 values[index_map[i]] = value_outer * values_1d[ix][0][0];
392 }
393 else
394 for (unsigned int iy = 0, i1 = 0; i1 < indices.size(); ++i1)
395 {
396 // prepare parts of products in y (and z) directions
397 std::array<double, dim + (dim * (dim - 1)) / 2> value_outer;
398 value_outer[0] = 1.;
399 if constexpr (dim > 1)
400 {
401 for (unsigned int x = 1; x < dim; ++x)
402 value_outer[0] *= values_1d[indices[i1][x - 1]][0][x];
403 for (unsigned int d = 1; d < dim; ++d)
404 {
405 value_outer[d] = values_1d[indices[i1][d - 1]][1][d];
406 for (unsigned int x = 1; x < dim; ++x)
407 if (x != d)
408 value_outer[d] *= values_1d[indices[i1][x - 1]][0][x];
409 }
410 for (unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
411 for (unsigned int d2 = d1; d2 < dim; ++d2, ++count)
412 {
413 value_outer[count] = 1.;
414 for (unsigned int x = 1; x < dim; ++x)
415 {
416 unsigned int derivative = 0;
417 if (d1 == x)
418 ++derivative;
419 if (d2 == x)
420 ++derivative;
421
422 value_outer[count] *=
423 values_1d[indices[i1][x - 1]][derivative][x];
424 }
425 }
426 }
427
428 // now run the loop over x and multiply by the values/derivatives
429 // in x direction
430 for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
431 {
432 std::array<double, 3> val_x{{values_1d[ix][0][0],
433 values_1d[ix][1][0],
434 values_1d[ix][2][0]}};
435 const unsigned int index =
436 (index_map.empty() ? i : index_map[i]);
437
438 if (update_values)
439 values[index] = value_outer[0] * val_x[0];
440
441 if (update_grads)
442 {
443 grads[index][0] = value_outer[0] * val_x[1];
444 if constexpr (dim > 1)
445 for (unsigned int d = 1; d < dim; ++d)
446 grads[index][d] = value_outer[d] * val_x[0];
447 }
448
449 if (update_grad_grads)
450 {
451 grad_grads[index][0][0] = value_outer[0] * val_x[2];
452 if constexpr (dim > 1)
453 {
454 for (unsigned int d = 1; d < dim; ++d)
455 grad_grads[index][0][d] = grad_grads[index][d][0] =
456 value_outer[d] * val_x[1];
457 for (unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
458 for (unsigned int d2 = d1; d2 < dim; ++d2, ++count)
459 grad_grads[index][d1][d2] =
460 grad_grads[index][d2][d1] =
461 value_outer[count] * val_x[0];
462 }
463 }
464 }
465
466 // Use slower code for 3rd and 4th derivatives
468 for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
469 {
470 const unsigned int index =
471 (index_map.empty() ? i : index_map[i]);
472 std::array<unsigned int, dim> my_indices;
473 my_indices[0] = ix;
474 if constexpr (dim > 1)
475 for (unsigned int d = 1; d < dim; ++d)
476 my_indices[d] = indices[i1][d - 1];
477 for (unsigned int d1 = 0; d1 < dim; ++d1)
478 for (unsigned int d2 = 0; d2 < dim; ++d2)
479 for (unsigned int d3 = 0; d3 < dim; ++d3)
480 {
481 double der3 = 1.;
482 for (unsigned int x = 0; x < dim; ++x)
483 {
484 unsigned int derivative = 0;
485 if (d1 == x)
486 ++derivative;
487 if (d2 == x)
488 ++derivative;
489 if (d3 == x)
490 ++derivative;
491
492 der3 *= values_1d[my_indices[x]][derivative][x];
493 }
494 third_derivatives[index][d1][d2][d3] = der3;
495 }
496 }
497
498 if (update_4th_derivatives)
499 for (unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
500 {
501 const unsigned int index =
502 (index_map.empty() ? i : index_map[i]);
503 std::array<unsigned int, dim> my_indices;
504 my_indices[0] = ix;
505 if constexpr (dim > 1)
506 for (unsigned int d = 1; d < dim; ++d)
507 my_indices[d] = indices[i1][d - 1];
508 for (unsigned int d1 = 0; d1 < dim; ++d1)
509 for (unsigned int d2 = 0; d2 < dim; ++d2)
510 for (unsigned int d3 = 0; d3 < dim; ++d3)
511 for (unsigned int d4 = 0; d4 < dim; ++d4)
512 {
513 double der4 = 1.;
514 for (unsigned int x = 0; x < dim; ++x)
515 {
516 unsigned int derivative = 0;
517 if (d1 == x)
518 ++derivative;
519 if (d2 == x)
520 ++derivative;
521 if (d3 == x)
522 ++derivative;
523 if (d4 == x)
524 ++derivative;
525
526 der4 *= values_1d[my_indices[x]][derivative][x];
527 }
528 fourth_derivatives[index][d1][d2][d3][d4] = der4;
529 }
530 }
531
532 iy += size_x;
533 }
534 }
535 } // namespace TensorProductPolynomials
536} // namespace internal
537
538
539
540template <int dim, typename PolynomialType>
541void
543 const Point<dim> &p,
544 std::vector<double> &values,
545 std::vector<Tensor<1, dim>> &grads,
546 std::vector<Tensor<2, dim>> &grad_grads,
547 std::vector<Tensor<3, dim>> &third_derivatives,
548 std::vector<Tensor<4, dim>> &fourth_derivatives) const
549{
550 Assert(dim <= 3, ExcNotImplemented());
551 Assert(values.size() == this->n() || values.empty(),
552 ExcDimensionMismatch2(values.size(), this->n(), 0));
553 Assert(grads.size() == this->n() || grads.empty(),
554 ExcDimensionMismatch2(grads.size(), this->n(), 0));
555 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
556 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
557 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
558 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
559 Assert(fourth_derivatives.size() == this->n() || 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<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> indices(1);
606 if constexpr (dim > 1)
607 for (unsigned int d = 1; d < dim; ++d)
608 {
609 const unsigned int size = indices.size();
610 for (unsigned int i = 1; i < n_polynomials; ++i)
611 for (unsigned int j = 0; j < size; ++j)
612 {
613 std::array<unsigned int, dim1> next_index = indices[j];
614 next_index[d - 1] = i;
615 indices.push_back(next_index);
616 }
617 }
618 AssertDimension(indices.size(), Utilities::pow(n_polynomials, dim - 1));
619
621 n_derivatives,
622 values_1d,
623 n_polynomials,
624 indices,
625 index_map_inverse,
626 values,
627 grads,
628 grad_grads,
629 third_derivatives,
630 fourth_derivatives);
631}
632
633
634
635template <>
636void
638 const Point<0> &,
639 std::vector<double> &,
640 std::vector<Tensor<1, 0>> &,
641 std::vector<Tensor<2, 0>> &,
642 std::vector<Tensor<3, 0>> &,
643 std::vector<Tensor<4, 0>> &) const
644{
645 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
646}
647
648
649
650template <int dim, typename PolynomialType>
651std::unique_ptr<ScalarPolynomialsBase<dim>>
653{
654 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
655}
656
657
658
659template <int dim, typename PolynomialType>
660std::size_t
667
668
669
670template <int dim, typename PolynomialType>
671std::vector<PolynomialType>
677
678
679
680/* ------------------- AnisotropicPolynomials -------------- */
681
682
683template <int dim>
685 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
686 : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
687 , polynomials(pols)
688 , index_map(this->n())
689 , index_map_inverse(this->n())
690{
691 Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
692 for (const auto &pols_d : pols)
693 {
694 (void)pols_d;
695 Assert(pols_d.size() > 0,
696 ExcMessage("The number of polynomials must be larger than zero "
697 "for all coordinate directions."));
698 }
699
700 // per default set this index map to identity. This map can be changed by
701 // the user through the set_numbering() function
702 for (unsigned int i = 0; i < this->n(); ++i)
703 {
704 index_map[i] = i;
705 index_map_inverse[i] = i;
706 }
707}
708
709
710
711template <int dim>
712void
714 const unsigned int i,
715 std::array<unsigned int, dim> &indices) const
716{
717#ifdef DEBUG
718 unsigned int n_poly = 1;
719 for (unsigned int d = 0; d < dim; ++d)
720 n_poly *= polynomials[d].size();
721 Assert(i < n_poly, ExcInternalError());
722#endif
723
724 if (dim == 0)
725 {
726 }
727 else if (dim == 1)
728 internal::compute_tensor_index(index_map[i],
729 polynomials[0].size(),
730 0 /*not used*/,
731 indices);
732 else
733 internal::compute_tensor_index(index_map[i],
734 polynomials[0].size(),
735 polynomials[1].size(),
736 indices);
737}
738
739
740
741template <>
742void
744 std::array<unsigned int, 0> &) const
745{
746 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
747}
748
749
750
751template <int dim>
752double
754 const Point<dim> &p) const
755{
756 std::array<unsigned int, dim> indices;
757 compute_index(i, indices);
758
759 double value = 1.;
760 for (unsigned int d = 0; d < dim; ++d)
761 value *= polynomials[d][indices[d]].value(p[d]);
762
763 return value;
764}
765
766
767
768template <>
769double
771 const Point<0> &) const
772{
773 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
774
775 return {};
776}
777
778
779
780template <int dim>
783 const Point<dim> &p) const
784{
785 std::array<unsigned int, dim> indices;
786 compute_index(i, indices);
787
788 // compute values and
789 // uni-directional derivatives at
790 // the given point in each
791 // coordinate direction
793 for (unsigned int d = 0; d < dim; ++d)
794 polynomials[d][indices[d]].value(p[d], 1, v[d].data());
795
796 Tensor<1, dim> grad;
797 for (unsigned int d = 0; d < dim; ++d)
798 {
799 grad[d] = 1.;
800 for (unsigned int x = 0; x < dim; ++x)
801 grad[d] *= v[x][d == x];
802 }
803
804 return grad;
805}
806
807
808
809template <>
812 const Point<0> &) const
813{
814 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
815
816 return {};
817}
818
819
820
821template <int dim>
824 const Point<dim> &p) const
825{
826 std::array<unsigned int, dim> indices;
827 compute_index(i, indices);
828
830 for (unsigned int d = 0; d < dim; ++d)
831 polynomials[d][indices[d]].value(p[d], 2, v[d].data());
832
833 Tensor<2, dim> grad_grad;
834 for (unsigned int d1 = 0; d1 < dim; ++d1)
835 for (unsigned int d2 = 0; d2 < dim; ++d2)
836 {
837 grad_grad[d1][d2] = 1.;
838 for (unsigned int x = 0; x < dim; ++x)
839 {
840 unsigned int derivative = 0;
841 if (d1 == x || d2 == x)
842 {
843 if (d1 == d2)
844 derivative = 2;
845 else
846 derivative = 1;
847 }
848 grad_grad[d1][d2] *= v[x][derivative];
849 }
850 }
851
852 return grad_grad;
853}
854
855
856
857template <>
860 const Point<0> &) const
861{
862 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
863
864 return {};
865}
866
867
868
869template <int dim>
870void
872 const Point<dim> &p,
873 std::vector<double> &values,
874 std::vector<Tensor<1, dim>> &grads,
875 std::vector<Tensor<2, dim>> &grad_grads,
876 std::vector<Tensor<3, dim>> &third_derivatives,
877 std::vector<Tensor<4, dim>> &fourth_derivatives) const
878{
879 Assert(values.size() == this->n() || values.empty(),
880 ExcDimensionMismatch2(values.size(), this->n(), 0));
881 Assert(grads.size() == this->n() || grads.empty(),
882 ExcDimensionMismatch2(grads.size(), this->n(), 0));
883 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
884 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
885 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
886 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
887 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
888 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
889
890 // check how many values/derivatives we have to compute
891 unsigned int n_derivatives = 0;
892 if (values.size() == this->n())
893 n_derivatives = 0;
894 if (grads.size() == this->n())
895 n_derivatives = 1;
896 if (grad_grads.size() == this->n())
897 n_derivatives = 2;
898 if (third_derivatives.size() == this->n())
899 n_derivatives = 3;
900 if (fourth_derivatives.size() == this->n())
901 n_derivatives = 4;
902
903 // compute the values (and derivatives, if necessary) of all polynomials at
904 // this evaluation point
905 std::size_t max_n_polynomials = 0;
906 for (unsigned int d = 0; d < dim; ++d)
907 max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
908
909 // 5 is enough to store values and derivatives in all supported cases
910 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
911 max_n_polynomials);
912 if (n_derivatives == 0)
913 for (unsigned int d = 0; d < dim; ++d)
914 for (unsigned int i = 0; i < polynomials[d].size(); ++i)
915 values_1d[i][0][d] = polynomials[d][i].value(p[d]);
916 else
917 for (unsigned int d = 0; d < dim; ++d)
918 for (unsigned int i = 0; i < polynomials[d].size(); ++i)
919 {
920 // The isotropic tensor product function wants us to use a different
921 // innermost index, so we cannot pass the values_1d array into the
922 // function directly
923 std::array<double, 5> derivatives;
924 polynomials[d][i].value(p[d], n_derivatives, derivatives.data());
925 for (unsigned int j = 0; j <= n_derivatives; ++j)
926 values_1d[i][j][d] = derivatives[j];
927 }
928
929 // Unroll the tensor product indices in arbitrary dimension
930 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
931 boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
932 for (unsigned int d = 1; d < dim; ++d)
933 {
934 const unsigned int size = indices.size();
935 for (unsigned int i = 1; i < polynomials[d].size(); ++i)
936 for (unsigned int j = 0; j < size; ++j)
937 {
938 std::array<unsigned int, dim1> next_index = indices[j];
939 next_index[d - 1] = i;
940 indices.push_back(next_index);
941 }
942 }
943
945 n_derivatives,
946 values_1d,
947 polynomials[0].size(),
948 indices,
949 index_map_inverse,
950 values,
951 grads,
952 grad_grads,
953 third_derivatives,
954 fourth_derivatives);
955}
956
957
958
959template <>
960void
962 std::vector<double> &,
963 std::vector<Tensor<1, 0>> &,
964 std::vector<Tensor<2, 0>> &,
965 std::vector<Tensor<3, 0>> &,
966 std::vector<Tensor<4, 0>> &) const
967{
968 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
969}
970
971
972
973template <int dim>
974unsigned int
976 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
977{
978 unsigned int y = 1;
979 for (unsigned int d = 0; d < dim; ++d)
980 y *= pols[d].size();
981 return y;
982}
983
984
985
986template <>
987unsigned int
989 const std::vector<std::vector<Polynomials::Polynomial<double>>> &)
990{
991 AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!"));
992
993 return {};
994}
995
996
997
998template <int dim>
999std::unique_ptr<ScalarPolynomialsBase<dim>>
1001{
1002 return std::make_unique<AnisotropicPolynomials<dim>>(*this);
1003}
1004
1005
1006
1007/* ------------------- explicit instantiations -------------- */
1012
1013template class TensorProductPolynomials<
1014 1,
1016template class TensorProductPolynomials<
1017 2,
1019template class TensorProductPolynomials<
1020 3,
1022
1023template class AnisotropicPolynomials<0>;
1024template class AnisotropicPolynomials<1>;
1025template class AnisotropicPolynomials<2>;
1026template class AnisotropicPolynomials<3>;
1027
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: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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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 fixed_power(const T t)
Definition utilities.h:942
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