Reference documentation for deal.II version 9.3.3
\(\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\}}\)
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
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
22DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
23#include <boost/container/small_vector.hpp>
24DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
25
26#include <array>
27#include <memory>
28
29DEAL_II_NAMESPACE_OPEN
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
86TensorProductPolynomials<dim, PolynomialType>::compute_index(
87 const unsigned int i,
88 std::array<unsigned int, dim> &indices) const
89{
90 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
91 ExcInternalError());
92 internal::compute_tensor_index(index_map[i],
93 polynomials.size(),
94 polynomials.size(),
95 indices);
96}
97
98
99
100template <>
101inline void
102TensorProductPolynomials<0, Polynomials::Polynomial<double>>::compute_index(
103 const unsigned int,
104 std::array<unsigned int, 0> &) const
105{
106 constexpr int dim = 0;
107 AssertThrow(dim > 0, ExcNotImplemented());
108}
109
110
111
112template <int dim, typename PolynomialType>
113void
114TensorProductPolynomials<dim, PolynomialType>::output_indices(
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
132TensorProductPolynomials<0, Polynomials::Polynomial<double>>::output_indices(
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
143TensorProductPolynomials<dim, PolynomialType>::set_numbering(
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;
152}
153
154
155
156template <>
157void
158TensorProductPolynomials<0, Polynomials::Polynomial<double>>::set_numbering(
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
169TensorProductPolynomials<dim, PolynomialType>::compute_value(
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
189TensorProductPolynomials<0, Polynomials::Polynomial<double>>::compute_value(
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>
202Tensor<1, dim>
203TensorProductPolynomials<dim, PolynomialType>::compute_grad(
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
214 ndarray<double, dim, 2> v;
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 <>
239Tensor<1, 0>
240TensorProductPolynomials<0, Polynomials::Polynomial<double>>::compute_grad(
241 const unsigned int,
242 const Point<0> &) const
243{
244 constexpr int dim = 0;
245 AssertThrow(dim > 0, ExcNotImplemented());
246
247 return {};
248}
249
250
251
252template <int dim, typename PolynomialType>
253Tensor<2, dim>
254TensorProductPolynomials<dim, PolynomialType>::compute_grad_grad(
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
261 ndarray<double, dim, 3> v;
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 <>
298Tensor<2, 0>
299TensorProductPolynomials<0, Polynomials::Polynomial<double>>::compute_grad_grad(
300 const unsigned int,
301 const Point<0> &) const
302{
303 return {};
304}
305
306
307
308template <int dim, typename PolynomialType>
309void
310TensorProductPolynomials<dim, PolynomialType>::evaluate(
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;
345 if (update_3rd_derivatives)
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
422 if (update_3rd_derivatives)
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
473TensorProductPolynomials<0, Polynomials::Polynomial<double>>::evaluate(
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>>
489TensorProductPolynomials<dim, PolynomialType>::clone() const
490{
491 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
492}
493
494
495
496template <int dim, typename PolynomialType>
497std::size_t
498TensorProductPolynomials<dim, PolynomialType>::memory_consumption() const
499{
500 return (MemoryConsumption::memory_consumption(polynomials) +
501 MemoryConsumption::memory_consumption(index_map) +
502 MemoryConsumption::memory_consumption(index_map_inverse));
503}
504
505
506
507template <int dim, typename PolynomialType>
508std::vector<PolynomialType>
509TensorProductPolynomials<dim, PolynomialType>::get_underlying_polynomials()
510 const
511{
512 return polynomials;
513}
514
515
516
517/* ------------------- AnisotropicPolynomials -------------- */
518
519
520template <int dim>
521AnisotropicPolynomials<dim>::AnisotropicPolynomials(
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
540AnisotropicPolynomials<dim>::compute_index(
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 }
554 else if (dim == 1)
555 internal::compute_tensor_index(i,
556 polynomials[0].size(),
557 0 /*not used*/,
558 indices);
559 else
560 internal::compute_tensor_index(i,
561 polynomials[0].size(),
562 polynomials[1].size(),
563 indices);
564}
565
566
567
568template <>
569void
570AnisotropicPolynomials<0>::compute_index(const unsigned int,
571 std::array<unsigned int, 0> &) const
572{
573 constexpr int dim = 0;
574 AssertThrow(dim > 0, ExcNotImplemented());
575}
576
577
578
579template <int dim>
580double
581AnisotropicPolynomials<dim>::compute_value(const unsigned int i,
582 const Point<dim> & p) const
583{
584 std::array<unsigned int, dim> indices;
585 compute_index(i, indices);
586
587 double value = 1.;
588 for (unsigned int d = 0; d < dim; ++d)
589 value *= polynomials[d][indices[d]].value(p(d));
590
591 return value;
592}
593
594
595
596template <>
597double
598AnisotropicPolynomials<0>::compute_value(const unsigned int,
599 const Point<0> &) const
600{
601 constexpr int dim = 0;
602 AssertThrow(dim > 0, ExcNotImplemented());
603
604 return {};
605}
606
607
608
609template <int dim>
610Tensor<1, dim>
611AnisotropicPolynomials<dim>::compute_grad(const unsigned int i,
612 const Point<dim> & p) const
613{
614 std::array<unsigned int, dim> indices;
615 compute_index(i, indices);
616
617 // compute values and
618 // uni-directional derivatives at
619 // the given point in each
620 // coordinate direction
621 ndarray<double, dim, 2> v;
622 for (unsigned int d = 0; d < dim; ++d)
623 polynomials[d][indices[d]].value(p(d), 1, v[d].data());
624
625 Tensor<1, dim> grad;
626 for (unsigned int d = 0; d < dim; ++d)
627 {
628 grad[d] = 1.;
629 for (unsigned int x = 0; x < dim; ++x)
630 grad[d] *= v[x][d == x];
631 }
632
633 return grad;
634}
635
636
637
638template <>
639Tensor<1, 0>
640AnisotropicPolynomials<0>::compute_grad(const unsigned int,
641 const Point<0> &) const
642{
643 constexpr int dim = 0;
644 AssertThrow(dim > 0, ExcNotImplemented());
645
646 return {};
647}
648
649
650
651template <int dim>
652Tensor<2, dim>
653AnisotropicPolynomials<dim>::compute_grad_grad(const unsigned int i,
654 const Point<dim> & p) const
655{
656 std::array<unsigned int, dim> indices;
657 compute_index(i, indices);
658
659 ndarray<double, dim, 3> v;
660 for (unsigned int d = 0; d < dim; ++d)
661 polynomials[d][indices[d]].value(p(d), 2, v[d].data());
662
663 Tensor<2, dim> grad_grad;
664 for (unsigned int d1 = 0; d1 < dim; ++d1)
665 for (unsigned int d2 = 0; d2 < dim; ++d2)
666 {
667 grad_grad[d1][d2] = 1.;
668 for (unsigned int x = 0; x < dim; ++x)
669 {
670 unsigned int derivative = 0;
671 if (d1 == x || d2 == x)
672 {
673 if (d1 == d2)
674 derivative = 2;
675 else
676 derivative = 1;
677 }
678 grad_grad[d1][d2] *= v[x][derivative];
679 }
680 }
681
682 return grad_grad;
683}
684
685
686
687template <>
688Tensor<2, 0>
689AnisotropicPolynomials<0>::compute_grad_grad(const unsigned int,
690 const Point<0> &) const
691{
692 constexpr int dim = 0;
693 AssertThrow(dim > 0, ExcNotImplemented());
694
695 return {};
696}
697
698
699
700template <int dim>
701void
702AnisotropicPolynomials<dim>::evaluate(
703 const Point<dim> & p,
704 std::vector<double> & values,
705 std::vector<Tensor<1, dim>> &grads,
706 std::vector<Tensor<2, dim>> &grad_grads,
707 std::vector<Tensor<3, dim>> &third_derivatives,
708 std::vector<Tensor<4, dim>> &fourth_derivatives) const
709{
710 Assert(values.size() == this->n() || values.size() == 0,
711 ExcDimensionMismatch2(values.size(), this->n(), 0));
712 Assert(grads.size() == this->n() || grads.size() == 0,
713 ExcDimensionMismatch2(grads.size(), this->n(), 0));
714 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
715 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
716 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
717 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
718 Assert(fourth_derivatives.size() == this->n() ||
719 fourth_derivatives.size() == 0,
720 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
721
722 const bool update_values = (values.size() == this->n()),
723 update_grads = (grads.size() == this->n()),
724 update_grad_grads = (grad_grads.size() == this->n()),
725 update_3rd_derivatives = (third_derivatives.size() == this->n()),
726 update_4th_derivatives = (fourth_derivatives.size() == this->n());
727
728 // check how many
729 // values/derivatives we have to
730 // compute
731 unsigned int n_values_and_derivatives = 0;
732 if (update_values)
733 n_values_and_derivatives = 1;
734 if (update_grads)
735 n_values_and_derivatives = 2;
736 if (update_grad_grads)
737 n_values_and_derivatives = 3;
738 if (update_3rd_derivatives)
739 n_values_and_derivatives = 4;
740 if (update_4th_derivatives)
741 n_values_and_derivatives = 5;
742
743 // compute the values (and
744 // derivatives, if necessary) of
745 // all polynomials at this
746 // evaluation point
747 std::size_t max_n_polynomials = 0;
748 for (unsigned int d = 0; d < dim; ++d)
749 max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
750
751 // 5 is enough to store values and derivatives in all supported cases
752 Table<2, std::array<double, 5>> v(dim, max_n_polynomials);
753 for (unsigned int d = 0; d < dim; ++d)
754 for (unsigned int i = 0; i < polynomials[d].size(); ++i)
755 polynomials[d][i].value(p(d),
756 n_values_and_derivatives - 1,
757 v(d, i).data());
758
759 for (unsigned int i = 0; i < this->n(); ++i)
760 {
761 // first get the
762 // one-dimensional indices of
763 // this particular tensor
764 // product polynomial
765 std::array<unsigned int, dim> indices;
766 compute_index(i, indices);
767
768 if (update_values)
769 {
770 values[i] = 1;
771 for (unsigned int x = 0; x < dim; ++x)
772 values[i] *= v(x, indices[x])[0];
773 }
774
775 if (update_grads)
776 for (unsigned int d = 0; d < dim; ++d)
777 {
778 grads[i][d] = 1.;
779 for (unsigned int x = 0; x < dim; ++x)
780 grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0];
781 }
782
783 if (update_grad_grads)
784 for (unsigned int d1 = 0; d1 < dim; ++d1)
785 for (unsigned int d2 = 0; d2 < dim; ++d2)
786 {
787 grad_grads[i][d1][d2] = 1.;
788 for (unsigned int x = 0; x < dim; ++x)
789 {
790 unsigned int derivative = 0;
791 if (d1 == x)
792 ++derivative;
793 if (d2 == x)
794 ++derivative;
795
796 grad_grads[i][d1][d2] *= v(x, indices[x])[derivative];
797 }
798 }
799
800 if (update_3rd_derivatives)
801 for (unsigned int d1 = 0; d1 < dim; ++d1)
802 for (unsigned int d2 = 0; d2 < dim; ++d2)
803 for (unsigned int d3 = 0; d3 < dim; ++d3)
804 {
805 third_derivatives[i][d1][d2][d3] = 1.;
806 for (unsigned int x = 0; x < dim; ++x)
807 {
808 unsigned int derivative = 0;
809 if (d1 == x)
810 ++derivative;
811 if (d2 == x)
812 ++derivative;
813 if (d3 == x)
814 ++derivative;
815
816 third_derivatives[i][d1][d2][d3] *=
817 v(x, indices[x])[derivative];
818 }
819 }
820
821 if (update_4th_derivatives)
822 for (unsigned int d1 = 0; d1 < dim; ++d1)
823 for (unsigned int d2 = 0; d2 < dim; ++d2)
824 for (unsigned int d3 = 0; d3 < dim; ++d3)
825 for (unsigned int d4 = 0; d4 < dim; ++d4)
826 {
827 fourth_derivatives[i][d1][d2][d3][d4] = 1.;
828 for (unsigned int x = 0; x < dim; ++x)
829 {
830 unsigned int derivative = 0;
831 if (d1 == x)
832 ++derivative;
833 if (d2 == x)
834 ++derivative;
835 if (d3 == x)
836 ++derivative;
837 if (d4 == x)
838 ++derivative;
839
840 fourth_derivatives[i][d1][d2][d3][d4] *=
841 v(x, indices[x])[derivative];
842 }
843 }
844 }
845}
846
847
848
849template <>
850void
851AnisotropicPolynomials<0>::evaluate(const Point<0> &,
852 std::vector<double> &,
853 std::vector<Tensor<1, 0>> &,
854 std::vector<Tensor<2, 0>> &,
855 std::vector<Tensor<3, 0>> &,
856 std::vector<Tensor<4, 0>> &) const
857{
858 constexpr int dim = 0;
859 AssertThrow(dim > 0, ExcNotImplemented());
860}
861
862
863
864template <int dim>
865unsigned int
866AnisotropicPolynomials<dim>::get_n_tensor_pols(
867 const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
868{
869 unsigned int y = 1;
870 for (unsigned int d = 0; d < dim; ++d)
871 y *= pols[d].size();
872 return y;
873}
874
875
876
877template <>
878unsigned int
879AnisotropicPolynomials<0>::get_n_tensor_pols(
880 const std::vector<std::vector<Polynomials::Polynomial<double>>> &)
881{
882 constexpr int dim = 0;
883 AssertThrow(dim > 0, ExcNotImplemented());
884
885 return {};
886}
887
888
889
890template <int dim>
891std::unique_ptr<ScalarPolynomialsBase<dim>>
892AnisotropicPolynomials<dim>::clone() const
893{
894 return std::make_unique<AnisotropicPolynomials<dim>>(*this);
895}
896
897
898
899/* ------------------- explicit instantiations -------------- */
900template class TensorProductPolynomials<0, Polynomials::Polynomial<double>>;
901template class TensorProductPolynomials<1, Polynomials::Polynomial<double>>;
902template class TensorProductPolynomials<2, Polynomials::Polynomial<double>>;
903template class TensorProductPolynomials<3, Polynomials::Polynomial<double>>;
904
905template class TensorProductPolynomials<
906 1,
907 Polynomials::PiecewisePolynomial<double>>;
908template class TensorProductPolynomials<
909 2,
910 Polynomials::PiecewisePolynomial<double>>;
911template class TensorProductPolynomials<
912 3,
913 Polynomials::PiecewisePolynomial<double>>;
914
915template class AnisotropicPolynomials<0>;
916template class AnisotropicPolynomials<1>;
917template class AnisotropicPolynomials<2>;
918template class AnisotropicPolynomials<3>;
919
920DEAL_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:702
AnisotropicPolynomials::compute_value
double compute_value(const unsigned int i, const Point< dim > &p) const override
Definition: tensor_product_polynomials.cc:581
AnisotropicPolynomials::AnisotropicPolynomials
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double > > > &base_polynomials)
Definition: tensor_product_polynomials.cc:521
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:866
AnisotropicPolynomials::compute_index
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Definition: tensor_product_polynomials.cc:540
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:653
AnisotropicPolynomials::compute_grad
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
Definition: tensor_product_polynomials.cc:611
AnisotropicPolynomials::clone
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Definition: tensor_product_polynomials.cc:892
Point
Definition: point.h:111
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:310
TensorProductPolynomials::output_indices
void output_indices(std::ostream &out) const
Definition: tensor_product_polynomials.cc:114
TensorProductPolynomials::compute_index
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Definition: tensor_product_polynomials.cc:86
TensorProductPolynomials::compute_value
double compute_value(const unsigned int i, const Point< dim > &p) const override
Definition: tensor_product_polynomials.cc:169
TensorProductPolynomials::memory_consumption
virtual std::size_t memory_consumption() const override
Definition: tensor_product_polynomials.cc:498
TensorProductPolynomials::compute_grad
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
Definition: tensor_product_polynomials.cc:203
TensorProductPolynomials::clone
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Definition: tensor_product_polynomials.cc:489
TensorProductPolynomials::get_underlying_polynomials
std::vector< PolynomialType > get_underlying_polynomials() const
Definition: tensor_product_polynomials.cc:509
TensorProductPolynomials::set_numbering
void set_numbering(const std::vector< unsigned int > &renumber)
Definition: tensor_product_polynomials.cc:143
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:254
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:416
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:454
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
update_3rd_derivatives
@ update_3rd_derivatives
Third derivatives of shape functions.
Definition: fe_update_flags.h:96
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
StandardExceptions::ExcDimensionMismatch2
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int 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:1575
memory_consumption.h
CUDAWrappers::internal::compute_index
unsigned int compute_index()
Definition: cuda_fe_evaluation.h:49
MemoryConsumption::memory_consumption
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Definition: memory_consumption.h:266
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5484
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.4