Reference documentation for deal.II version Git 73c87d96ef 2021-11-30 22:54:44 +0100
\(\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 
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 
36 namespace 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 
84 template <int dim, typename PolynomialType>
85 inline 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 
98 
99 
100 template <>
101 inline void
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 
112 template <int dim, typename PolynomialType>
113 void
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 
130 template <>
131 void
133  std::ostream &) const
134 {
135  constexpr int dim = 0;
136  AssertThrow(dim > 0, ExcNotImplemented());
137 }
138 
139 
140 
141 template <int dim, typename PolynomialType>
142 void
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 
156 template <>
157 void
159  const std::vector<unsigned int> &)
160 {
161  constexpr int dim = 0;
162  AssertThrow(dim > 0, ExcNotImplemented());
163 }
164 
165 
166 
167 template <int dim, typename PolynomialType>
168 double
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 
187 template <>
188 double
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 
201 template <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 
238 template <>
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 
252 template <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 
297 template <>
300  const unsigned int,
301  const Point<0> &) const
302 {
303  return {};
304 }
305 
306 
307 
308 template <int dim, typename PolynomialType>
309 void
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 
471 template <>
472 void
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 
487 template <int dim, typename PolynomialType>
488 std::unique_ptr<ScalarPolynomialsBase<dim>>
490 {
491  return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
492 }
493 
494 
495 
496 template <int dim, typename PolynomialType>
497 std::size_t
499 {
500  return (MemoryConsumption::memory_consumption(polynomials) +
502  MemoryConsumption::memory_consumption(index_map_inverse));
503 }
504 
505 
506 
507 template <int dim, typename PolynomialType>
508 std::vector<PolynomialType>
510  const
511 {
512  return polynomials;
513 }
514 
515 
516 
517 /* ------------------- AnisotropicPolynomials -------------- */
518 
519 
520 template <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 
538 template <int dim>
539 void
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 
567 template <>
568 void
570  std::array<unsigned int, 0> &) const
571 {
572  constexpr int dim = 0;
573  AssertThrow(dim > 0, ExcNotImplemented());
574 }
575 
576 
577 
578 template <int dim>
579 double
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 
595 template <>
596 double
598  const Point<0> &) const
599 {
600  constexpr int dim = 0;
601  AssertThrow(dim > 0, ExcNotImplemented());
602 
603  return {};
604 }
605 
606 
607 
608 template <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 
637 template <>
640  const Point<0> &) const
641 {
642  constexpr int dim = 0;
643  AssertThrow(dim > 0, ExcNotImplemented());
644 
645  return {};
646 }
647 
648 
649 
650 template <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 
686 template <>
689  const Point<0> &) const
690 {
691  constexpr int dim = 0;
692  AssertThrow(dim > 0, ExcNotImplemented());
693 
694  return {};
695 }
696 
697 
698 
699 template <int dim>
700 void
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 
848 template <>
849 void
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 
863 template <int dim>
864 unsigned 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 
876 template <>
877 unsigned 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 
889 template <int dim>
890 std::unique_ptr<ScalarPolynomialsBase<dim>>
892 {
893  return std::make_unique<AnisotropicPolynomials<dim>>(*this);
894 }
895 
896 
897 
898 /* ------------------- explicit instantiations -------------- */
903 
904 template class TensorProductPolynomials<
905  1,
907 template class TensorProductPolynomials<
908  2,
909  Polynomials::PiecewisePolynomial<double>>;
910 template class TensorProductPolynomials<
911  3,
912  Polynomials::PiecewisePolynomial<double>>;
913 
914 template class AnisotropicPolynomials<0>;
915 template class AnisotropicPolynomials<1>;
916 template class AnisotropicPolynomials<2>;
917 template class AnisotropicPolynomials<3>;
918 
Shape function values.
double compute_value(const unsigned int i, const Point< dim > &p) const override
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
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
std::vector< PolynomialType > get_underlying_polynomials() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:414
static ::ExceptionBase & ExcMessage(std::string arg1)
Third derivatives of shape functions.
#define Assert(cond, exc)
Definition: exceptions.h:1461
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
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
virtual std::size_t memory_consumption() const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
const std::vector< std::vector< Polynomials::Polynomial< double > > > polynomials
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
double compute_value(const unsigned int i, const Point< dim > &p) const override
unsigned int compute_index()
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:452
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
static ::ExceptionBase & ExcInternalError()
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override