Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
tensor_product_polynomials.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2021 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
23#include <boost/container/small_vector.hpp>
25
26#include <array>
27#include <memory>
28
30
31
32
33/* ------------------- TensorProductPolynomials -------------- */
34
35
36namespace internal
37{
38 namespace
39 {
40 template <std::size_t dim>
41 inline void
42 compute_tensor_index(const unsigned int,
43 const unsigned int,
44 const unsigned int,
45 std::array<unsigned int, dim> &)
46 {
47 Assert(false, ExcNotImplemented());
48 }
49
50 inline void
51 compute_tensor_index(const unsigned int n,
52 const unsigned int,
53 const unsigned int,
54 std::array<unsigned int, 1> &indices)
55 {
56 indices[0] = n;
57 }
58
59 inline void
60 compute_tensor_index(const unsigned int n,
61 const unsigned int n_pols_0,
62 const unsigned int,
63 std::array<unsigned int, 2> &indices)
64 {
65 indices[0] = n % n_pols_0;
66 indices[1] = n / n_pols_0;
67 }
68
69 inline void
70 compute_tensor_index(const unsigned int n,
71 const unsigned int n_pols_0,
72 const unsigned int n_pols_1,
73 std::array<unsigned int, 3> &indices)
74 {
75 indices[0] = n % n_pols_0;
76 indices[1] = (n / n_pols_0) % n_pols_1;
77 indices[2] = n / (n_pols_0 * n_pols_1);
78 }
79 } // namespace
80} // namespace internal
81
82
83
84template <int dim, typename PolynomialType>
85inline void
87 const unsigned int i,
88 std::array<unsigned int, dim> &indices) const
89{
90 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
92 internal::compute_tensor_index(index_map[i],
93 polynomials.size(),
94 polynomials.size(),
95 indices);
96}
97
99
100template <>
101inline void
103 const unsigned int,
104 std::array<unsigned int, 0> &) const
106 constexpr int dim = 0;
107 AssertThrow(dim > 0, ExcNotImplemented());
108}
109
110
111
112template <int dim, typename PolynomialType>
113void
115 std::ostream &out) const
116{
117 std::array<unsigned int, dim> ix;
118 for (unsigned int i = 0; i < this->n(); ++i)
119 {
120 compute_index(i, ix);
121 out << i << "\t";
122 for (unsigned int d = 0; d < dim; ++d)
123 out << ix[d] << " ";
124 out << std::endl;
125 }
126}
127
128
129
130template <>
131void
133 std::ostream &) const
134{
135 constexpr int dim = 0;
136 AssertThrow(dim > 0, ExcNotImplemented());
137}
138
139
140
141template <int dim, typename PolynomialType>
142void
144 const std::vector<unsigned int> &renumber)
145{
146 Assert(renumber.size() == index_map.size(),
147 ExcDimensionMismatch(renumber.size(), index_map.size()));
148
149 index_map = renumber;
150 for (unsigned int i = 0; i < index_map.size(); ++i)
151 index_map_inverse[index_map[i]] = i;
153
154
155
156template <>
157void
159 const std::vector<unsigned int> &)
160{
161 constexpr int dim = 0;
162 AssertThrow(dim > 0, ExcNotImplemented());
163}
164
165
166
167template <int dim, typename PolynomialType>
168double
170 const unsigned int i,
171 const Point<dim> & p) const
172{
173 Assert(dim > 0, ExcNotImplemented());
174
175 std::array<unsigned int, dim> indices;
176 compute_index(i, indices);
177
178 double value = 1.;
179 for (unsigned int d = 0; d < dim; ++d)
180 value *= polynomials[indices[d]].value(p(d));
181
182 return value;
183}
184
185
186
187template <>
188double
190 const unsigned int,
191 const Point<0> &) const
192{
193 constexpr int dim = 0;
194 AssertThrow(dim > 0, ExcNotImplemented());
195
196 return {};
197}
198
199
200
201template <int dim, typename PolynomialType>
204 const unsigned int i,
205 const Point<dim> & p) const
206{
207 std::array<unsigned int, dim> indices;
208 compute_index(i, indices);
209
210 // compute values and
211 // uni-directional derivatives at
212 // the given point in each
213 // coordinate direction
215 {
216 std::vector<double> tmp(2);
217 for (unsigned int d = 0; d < dim; ++d)
218 {
219 polynomials[indices[d]].value(p(d), tmp);
220 v[d][0] = tmp[0];
221 v[d][1] = tmp[1];
222 }
223 }
224
225 Tensor<1, dim> grad;
226 for (unsigned int d = 0; d < dim; ++d)
227 {
228 grad[d] = 1.;
229 for (unsigned int x = 0; x < dim; ++x)
230 grad[d] *= v[x][d == x];
231 }
232
233 return grad;
234}
235
236
237
238template <>
241 const unsigned int,
242 const Point<0> &) const
243{
244 constexpr int dim = 0;
245 AssertThrow(dim > 0, ExcNotImplemented());
247 return {};
248}
249
250
251
252template <int dim, typename PolynomialType>
255 const unsigned int i,
256 const Point<dim> & p) const
257{
258 std::array<unsigned int, dim> indices;
259 compute_index(i, indices);
260
262 {
263 std::vector<double> tmp(3);
264 for (unsigned int d = 0; d < dim; ++d)
265 {
266 polynomials[indices[d]].value(p(d), tmp);
267 v[d][0] = tmp[0];
268 v[d][1] = tmp[1];
269 v[d][2] = tmp[2];
270 }
271 }
272
273 Tensor<2, dim> grad_grad;
274 for (unsigned int d1 = 0; d1 < dim; ++d1)
275 for (unsigned int d2 = 0; d2 < dim; ++d2)
276 {
277 grad_grad[d1][d2] = 1.;
278 for (unsigned int x = 0; x < dim; ++x)
279 {
280 unsigned int derivative = 0;
281 if (d1 == x || d2 == x)
282 {
283 if (d1 == d2)
284 derivative = 2;
285 else
286 derivative = 1;
287 }
288 grad_grad[d1][d2] *= v[x][derivative];
289 }
290 }
291
292 return grad_grad;
293}
294
295
296
297template <>
300 const unsigned int,
301 const Point<0> &) const
302{
303 return {};
304}
305
306
307
308template <int dim, typename PolynomialType>
309void
311 const Point<dim> & p,
312 std::vector<double> & values,
313 std::vector<Tensor<1, dim>> &grads,
314 std::vector<Tensor<2, dim>> &grad_grads,
315 std::vector<Tensor<3, dim>> &third_derivatives,
316 std::vector<Tensor<4, dim>> &fourth_derivatives) const
317{
318 Assert(dim <= 3, ExcNotImplemented());
319 Assert(values.size() == this->n() || values.size() == 0,
320 ExcDimensionMismatch2(values.size(), this->n(), 0));
321 Assert(grads.size() == this->n() || grads.size() == 0,
322 ExcDimensionMismatch2(grads.size(), this->n(), 0));
323 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
324 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
325 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
326 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
327 Assert(fourth_derivatives.size() == this->n() ||
328 fourth_derivatives.size() == 0,
329 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
330
331 const bool update_values = (values.size() == this->n()),
332 update_grads = (grads.size() == this->n()),
333 update_grad_grads = (grad_grads.size() == this->n()),
334 update_3rd_derivatives = (third_derivatives.size() == this->n()),
335 update_4th_derivatives = (fourth_derivatives.size() == this->n());
336
337 // check how many values/derivatives we have to compute
338 unsigned int n_values_and_derivatives = 0;
339 if (update_values)
340 n_values_and_derivatives = 1;
341 if (update_grads)
342 n_values_and_derivatives = 2;
343 if (update_grad_grads)
344 n_values_and_derivatives = 3;
346 n_values_and_derivatives = 4;
347 if (update_4th_derivatives)
348 n_values_and_derivatives = 5;
349
350 // Compute the values (and derivatives, if necessary) of all 1D polynomials
351 // at this evaluation point. We need to compute dim*n_polynomials
352 // evaluations, involving an evaluation of each polynomial for each
353 // coordinate direction. Once we have those values, we perform the
354 // multiplications for the tensor product in the arbitrary dimension.
355 const unsigned int n_polynomials = polynomials.size();
356 boost::container::small_vector<ndarray<double, dim, 5>, 20> values_1d(
357 n_polynomials);
358 if (n_values_and_derivatives == 1)
359 for (unsigned int i = 0; i < n_polynomials; ++i)
360 for (unsigned int d = 0; d < dim; ++d)
361 values_1d[i][d][0] = polynomials[i].value(p(d));
362 else
363 for (unsigned int i = 0; i < n_polynomials; ++i)
364 for (unsigned d = 0; d < dim; ++d)
365 polynomials[i].value(p(d),
366 n_values_and_derivatives,
367 values_1d[i][d].data());
368
369 unsigned int indices[3];
370 unsigned int ind = 0;
371 for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2])
372 for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1);
373 ++indices[1])
374 if (n_values_and_derivatives == 1)
375 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
376 {
377 double value = values_1d[indices[0]][0][0];
378 for (unsigned int d = 1; d < dim; ++d)
379 value *= values_1d[indices[d]][d][0];
380 values[index_map_inverse[ind]] = value;
381 }
382 else
383 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
384 {
385 const unsigned int i = index_map_inverse[ind];
386
387 if (update_values)
388 {
389 double value = values_1d[indices[0]][0][0];
390 for (unsigned int x = 1; x < dim; ++x)
391 value *= values_1d[indices[x]][x][0];
392 values[i] = value;
393 }
394
395 if (update_grads)
396 for (unsigned int d = 0; d < dim; ++d)
397 {
398 double grad = 1.;
399 for (unsigned int x = 0; x < dim; ++x)
400 grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
401 grads[i][d] = grad;
402 }
403
404 if (update_grad_grads)
405 for (unsigned int d1 = 0; d1 < dim; ++d1)
406 for (unsigned int d2 = 0; d2 < dim; ++d2)
407 {
408 double der2 = 1.;
409 for (unsigned int x = 0; x < dim; ++x)
410 {
411 unsigned int derivative = 0;
412 if (d1 == x)
413 ++derivative;
414 if (d2 == x)
415 ++derivative;
416
417 der2 *= values_1d[indices[x]][x][derivative];
418 }
419 grad_grads[i][d1][d2] = der2;
420 }
421
423 for (unsigned int d1 = 0; d1 < dim; ++d1)
424 for (unsigned int d2 = 0; d2 < dim; ++d2)
425 for (unsigned int d3 = 0; d3 < dim; ++d3)
426 {
427 double der3 = 1.;
428 for (unsigned int x = 0; x < dim; ++x)
429 {
430 unsigned int derivative = 0;
431 if (d1 == x)
432 ++derivative;
433 if (d2 == x)
434 ++derivative;
435 if (d3 == x)
436 ++derivative;
437
438 der3 *= values_1d[indices[x]][x][derivative];
439 }
440 third_derivatives[i][d1][d2][d3] = der3;
441 }
442
443 if (update_4th_derivatives)
444 for (unsigned int d1 = 0; d1 < dim; ++d1)
445 for (unsigned int d2 = 0; d2 < dim; ++d2)
446 for (unsigned int d3 = 0; d3 < dim; ++d3)
447 for (unsigned int d4 = 0; d4 < dim; ++d4)
448 {
449 double der4 = 1.;
450 for (unsigned int x = 0; x < dim; ++x)
451 {
452 unsigned int derivative = 0;
453 if (d1 == x)
454 ++derivative;
455 if (d2 == x)
456 ++derivative;
457 if (d3 == x)
458 ++derivative;
459 if (d4 == x)
460 ++derivative;
461
462 der4 *= values_1d[indices[x]][x][derivative];
463 }
464 fourth_derivatives[i][d1][d2][d3][d4] = der4;
465 }
466 }
467}
468
469
470
471template <>
472void
474 const Point<0> &,
475 std::vector<double> &,
476 std::vector<Tensor<1, 0>> &,
477 std::vector<Tensor<2, 0>> &,
478 std::vector<Tensor<3, 0>> &,
479 std::vector<Tensor<4, 0>> &) const
480{
481 constexpr int dim = 0;
482 AssertThrow(dim > 0, ExcNotImplemented());
483}
484
485
486
487template <int dim, typename PolynomialType>
488std::unique_ptr<ScalarPolynomialsBase<dim>>
490{
491 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
492}
493
494
495
496template <int dim, typename PolynomialType>
497std::size_t
499{
500 return (MemoryConsumption::memory_consumption(polynomials) +
502 MemoryConsumption::memory_consumption(index_map_inverse));
503}
504
505
506
507template <int dim, typename PolynomialType>
508std::vector<PolynomialType>
510 const
511{
512 return polynomials;
513}
514
515
516
517/* ------------------- AnisotropicPolynomials -------------- */
518
519
520template <int dim>
522 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
523 : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
524 , polynomials(pols)
525{
526 Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
527 for (const auto &pols_d : pols)
528 {
529 (void)pols_d;
530 Assert(pols_d.size() > 0,
531 ExcMessage("The number of polynomials must be larger than zero "
532 "for all coordinate directions."));
533 }
534}
535
536
537
538template <int dim>
539void
541 const unsigned int i,
542 std::array<unsigned int, dim> &indices) const
543{
544#ifdef DEBUG
545 unsigned int n_poly = 1;
546 for (unsigned int d = 0; d < dim; ++d)
547 n_poly *= polynomials[d].size();
548 Assert(i < n_poly, ExcInternalError());
549#endif
550
551 if (dim == 0)
552 {}
553 else if (dim == 1)
554 internal::compute_tensor_index(i,
555 polynomials[0].size(),
556 0 /*not used*/,
557 indices);
558 else
559 internal::compute_tensor_index(i,
560 polynomials[0].size(),
561 polynomials[1].size(),
562 indices);
563}
564
565
566
567template <>
568void
570 std::array<unsigned int, 0> &) const
571{
572 constexpr int dim = 0;
573 AssertThrow(dim > 0, ExcNotImplemented());
574}
575
576
577
578template <int dim>
579double
581 const Point<dim> & p) const
582{
583 std::array<unsigned int, dim> indices;
584 compute_index(i, indices);
585
586 double value = 1.;
587 for (unsigned int d = 0; d < dim; ++d)
588 value *= polynomials[d][indices[d]].value(p(d));
589
590 return value;
591}
592
593
594
595template <>
596double
598 const Point<0> &) const
599{
600 constexpr int dim = 0;
601 AssertThrow(dim > 0, ExcNotImplemented());
602
603 return {};
604}
605
606
607
608template <int dim>
611 const Point<dim> & p) const
612{
613 std::array<unsigned int, dim> indices;
614 compute_index(i, indices);
615
616 // compute values and
617 // uni-directional derivatives at
618 // the given point in each
619 // coordinate direction
621 for (unsigned int d = 0; d < dim; ++d)
622 polynomials[d][indices[d]].value(p(d), 1, v[d].data());
623
624 Tensor<1, dim> grad;
625 for (unsigned int d = 0; d < dim; ++d)
626 {
627 grad[d] = 1.;
628 for (unsigned int x = 0; x < dim; ++x)
629 grad[d] *= v[x][d == x];
630 }
631
632 return grad;
633}
634
635
636
637template <>
640 const Point<0> &) const
641{
642 constexpr int dim = 0;
643 AssertThrow(dim > 0, ExcNotImplemented());
644
645 return {};
646}
647
648
649
650template <int dim>
653 const Point<dim> & p) const
654{
655 std::array<unsigned int, dim> indices;
656 compute_index(i, indices);
657
659 for (unsigned int d = 0; d < dim; ++d)
660 polynomials[d][indices[d]].value(p(d), 2, v[d].data());
661
662 Tensor<2, dim> grad_grad;
663 for (unsigned int d1 = 0; d1 < dim; ++d1)
664 for (unsigned int d2 = 0; d2 < dim; ++d2)
665 {
666 grad_grad[d1][d2] = 1.;
667 for (unsigned int x = 0; x < dim; ++x)
668 {
669 unsigned int derivative = 0;
670 if (d1 == x || d2 == x)
671 {
672 if (d1 == d2)
673 derivative = 2;
674 else
675 derivative = 1;
676 }
677 grad_grad[d1][d2] *= v[x][derivative];
678 }
679 }
680
681 return grad_grad;
682}
683
684
685
686template <>
689 const Point<0> &) const
690{
691 constexpr int dim = 0;
692 AssertThrow(dim > 0, ExcNotImplemented());
693
694 return {};
695}
696
697
698
699template <int dim>
700void
702 const Point<dim> & p,
703 std::vector<double> & values,
704 std::vector<Tensor<1, dim>> &grads,
705 std::vector<Tensor<2, dim>> &grad_grads,
706 std::vector<Tensor<3, dim>> &third_derivatives,
707 std::vector<Tensor<4, dim>> &fourth_derivatives) const
708{
709 Assert(values.size() == this->n() || values.size() == 0,
710 ExcDimensionMismatch2(values.size(), this->n(), 0));
711 Assert(grads.size() == this->n() || grads.size() == 0,
712 ExcDimensionMismatch2(grads.size(), this->n(), 0));
713 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
714 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
715 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
716 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
717 Assert(fourth_derivatives.size() == this->n() ||
718 fourth_derivatives.size() == 0,
719 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
720
721 const bool update_values = (values.size() == this->n()),
722 update_grads = (grads.size() == this->n()),
723 update_grad_grads = (grad_grads.size() == this->n()),
724 update_3rd_derivatives = (third_derivatives.size() == this->n()),
725 update_4th_derivatives = (fourth_derivatives.size() == this->n());
726
727 // check how many
728 // values/derivatives we have to
729 // compute
730 unsigned int n_values_and_derivatives = 0;
731 if (update_values)
732 n_values_and_derivatives = 1;
733 if (update_grads)
734 n_values_and_derivatives = 2;
735 if (update_grad_grads)
736 n_values_and_derivatives = 3;
738 n_values_and_derivatives = 4;
739 if (update_4th_derivatives)
740 n_values_and_derivatives = 5;
741
742 // compute the values (and
743 // derivatives, if necessary) of
744 // all polynomials at this
745 // evaluation point
746 std::size_t max_n_polynomials = 0;
747 for (unsigned int d = 0; d < dim; ++d)
748 max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
749
750 // 5 is enough to store values and derivatives in all supported cases
751 Table<2, std::array<double, 5>> v(dim, max_n_polynomials);
752 for (unsigned int d = 0; d < dim; ++d)
753 for (unsigned int i = 0; i < polynomials[d].size(); ++i)
754 polynomials[d][i].value(p(d),
755 n_values_and_derivatives - 1,
756 v(d, i).data());
757
758 for (unsigned int i = 0; i < this->n(); ++i)
759 {
760 // first get the
761 // one-dimensional indices of
762 // this particular tensor
763 // product polynomial
764 std::array<unsigned int, dim> indices;
765 compute_index(i, indices);
766
767 if (update_values)
768 {
769 values[i] = 1;
770 for (unsigned int x = 0; x < dim; ++x)
771 values[i] *= v(x, indices[x])[0];
772 }
773
774 if (update_grads)
775 for (unsigned int d = 0; d < dim; ++d)
776 {
777 grads[i][d] = 1.;
778 for (unsigned int x = 0; x < dim; ++x)
779 grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0];
780 }
781
782 if (update_grad_grads)
783 for (unsigned int d1 = 0; d1 < dim; ++d1)
784 for (unsigned int d2 = 0; d2 < dim; ++d2)
785 {
786 grad_grads[i][d1][d2] = 1.;
787 for (unsigned int x = 0; x < dim; ++x)
788 {
789 unsigned int derivative = 0;
790 if (d1 == x)
791 ++derivative;
792 if (d2 == x)
793 ++derivative;
794
795 grad_grads[i][d1][d2] *= v(x, indices[x])[derivative];
796 }
797 }
798
800 for (unsigned int d1 = 0; d1 < dim; ++d1)
801 for (unsigned int d2 = 0; d2 < dim; ++d2)
802 for (unsigned int d3 = 0; d3 < dim; ++d3)
803 {
804 third_derivatives[i][d1][d2][d3] = 1.;
805 for (unsigned int x = 0; x < dim; ++x)
806 {
807 unsigned int derivative = 0;
808 if (d1 == x)
809 ++derivative;
810 if (d2 == x)
811 ++derivative;
812 if (d3 == x)
813 ++derivative;
814
815 third_derivatives[i][d1][d2][d3] *=
816 v(x, indices[x])[derivative];
817 }
818 }
819
820 if (update_4th_derivatives)
821 for (unsigned int d1 = 0; d1 < dim; ++d1)
822 for (unsigned int d2 = 0; d2 < dim; ++d2)
823 for (unsigned int d3 = 0; d3 < dim; ++d3)
824 for (unsigned int d4 = 0; d4 < dim; ++d4)
825 {
826 fourth_derivatives[i][d1][d2][d3][d4] = 1.;
827 for (unsigned int x = 0; x < dim; ++x)
828 {
829 unsigned int derivative = 0;
830 if (d1 == x)
831 ++derivative;
832 if (d2 == x)
833 ++derivative;
834 if (d3 == x)
835 ++derivative;
836 if (d4 == x)
837 ++derivative;
838
839 fourth_derivatives[i][d1][d2][d3][d4] *=
840 v(x, indices[x])[derivative];
841 }
842 }
843 }
844}
845
846
847
848template <>
849void
851 std::vector<double> &,
852 std::vector<Tensor<1, 0>> &,
853 std::vector<Tensor<2, 0>> &,
854 std::vector<Tensor<3, 0>> &,
855 std::vector<Tensor<4, 0>> &) const
856{
857 constexpr int dim = 0;
858 AssertThrow(dim > 0, ExcNotImplemented());
859}
860
861
862
863template <int dim>
864unsigned int
866 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
867{
868 unsigned int y = 1;
869 for (unsigned int d = 0; d < dim; ++d)
870 y *= pols[d].size();
871 return y;
872}
873
874
875
876template <>
877unsigned int
879 const std::vector<std::vector<Polynomials::Polynomial<double>>> &)
880{
881 constexpr int dim = 0;
882 AssertThrow(dim > 0, ExcNotImplemented());
883
884 return {};
885}
886
887
888
889template <int dim>
890std::unique_ptr<ScalarPolynomialsBase<dim>>
892{
893 return std::make_unique<AnisotropicPolynomials<dim>>(*this);
894}
895
896
897
898/* ------------------- explicit instantiations -------------- */
903
904template class TensorProductPolynomials<
905 1,
907template class TensorProductPolynomials<
908 2,
910template class TensorProductPolynomials<
911 3,
913
914template class AnisotropicPolynomials<0>;
915template class AnisotropicPolynomials<1>;
916template class AnisotropicPolynomials<2>;
917template class AnisotropicPolynomials<3>;
918
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
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108