Reference documentation for deal.II version 9.3.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\}}\)
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  }
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 
568 template <>
569 void
571  std::array<unsigned int, 0> &) const
572 {
573  constexpr int dim = 0;
574  AssertThrow(dim > 0, ExcNotImplemented());
575 }
576 
577 
578 
579 template <int dim>
580 double
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 
596 template <>
597 double
599  const Point<0> &) const
600 {
601  constexpr int dim = 0;
602  AssertThrow(dim > 0, ExcNotImplemented());
603 
604  return {};
605 }
606 
607 
608 
609 template <int dim>
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
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 
638 template <>
641  const Point<0> &) const
642 {
643  constexpr int dim = 0;
644  AssertThrow(dim > 0, ExcNotImplemented());
645 
646  return {};
647 }
648 
649 
650 
651 template <int dim>
654  const Point<dim> & p) const
655 {
656  std::array<unsigned int, dim> indices;
657  compute_index(i, indices);
658 
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 
687 template <>
690  const Point<0> &) const
691 {
692  constexpr int dim = 0;
693  AssertThrow(dim > 0, ExcNotImplemented());
694 
695  return {};
696 }
697 
698 
699 
700 template <int dim>
701 void
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;
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 
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 
849 template <>
850 void
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 
864 template <int dim>
865 unsigned int
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 
877 template <>
878 unsigned int
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 
890 template <int dim>
891 std::unique_ptr<ScalarPolynomialsBase<dim>>
893 {
894  return std::make_unique<AnisotropicPolynomials<dim>>(*this);
895 }
896 
897 
898 
899 /* ------------------- explicit instantiations -------------- */
904 
905 template class TensorProductPolynomials<
906  1,
908 template class TensorProductPolynomials<
909  2,
910  Polynomials::PiecewisePolynomial<double>>;
911 template class TensorProductPolynomials<
912  3,
913  Polynomials::PiecewisePolynomial<double>>;
914 
915 template class AnisotropicPolynomials<0>;
916 template class AnisotropicPolynomials<1>;
917 template class AnisotropicPolynomials<2>;
918 template class AnisotropicPolynomials<3>;
919 
Shape function values.
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
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:1575
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:408
static ::ExceptionBase & ExcMessage(std::string arg1)
Third derivatives of shape functions.
#define Assert(cond, exc)
Definition: exceptions.h:1465
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:395
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:445
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
static ::ExceptionBase & ExcNotImplemented()
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