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