Reference documentation for deal.II version Git 06d359936b 2020-09-19 18:16:22 -0400
\(\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 - 2019 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 
22 #include <boost/container/small_vector.hpp>
23 
24 #include <array>
25 #include <memory>
26 
28 
29 
30 
31 /* ------------------- TensorProductPolynomials -------------- */
32 
33 
34 namespace 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 
82 template <int dim, typename PolynomialType>
83 inline void
85  const unsigned int i,
86  std::array<unsigned int, dim> &indices) const
87 {
88  Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
90  internal::compute_tensor_index(index_map[i],
91  polynomials.size(),
92  polynomials.size(),
93  indices);
94 }
95 
96 
97 
98 template <>
99 inline void
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 
110 template <int dim, typename PolynomialType>
111 void
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 
128 template <>
129 void
131  std::ostream &) const
132 {
133  constexpr int dim = 0;
134  AssertThrow(dim > 0, ExcNotImplemented());
135 }
136 
137 
138 
139 template <int dim, typename PolynomialType>
140 void
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 
154 template <>
155 void
157  const std::vector<unsigned int> &)
158 {
159  constexpr int dim = 0;
160  AssertThrow(dim > 0, ExcNotImplemented());
161 }
162 
163 
164 
165 template <int dim, typename PolynomialType>
166 double
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 
185 template <>
186 double
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 
199 template <int dim, typename PolynomialType>
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  // co-ordinate direction
212  std::array<std::array<double, 2>, dim> 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 
236 template <>
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 
250 template <int dim, typename PolynomialType>
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  std::array<std::array<double, 3>, dim> 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 
295 template <>
298  const unsigned int,
299  const Point<0> &) const
300 {
301  return {};
302 }
303 
304 
305 
306 template <int dim, typename PolynomialType>
307 void
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;
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<std::array<std::array<double, 5>, dim>, 20>
355  values_1d(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  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 
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 
469 template <>
470 void
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 
485 template <int dim, typename PolynomialType>
486 std::unique_ptr<ScalarPolynomialsBase<dim>>
488 {
489  return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
490 }
491 
492 
493 
494 template <int dim, typename PolynomialType>
495 std::size_t
497 {
498  return (MemoryConsumption::memory_consumption(polynomials) +
500  MemoryConsumption::memory_consumption(index_map_inverse));
501 }
502 
503 
504 
505 /* ------------------- AnisotropicPolynomials -------------- */
506 
507 
508 template <int dim>
510  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
511  : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
512  , polynomials(pols)
513 {
514  Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
515  for (const auto &pols_d : pols)
516  {
517  (void)pols_d;
518  Assert(pols_d.size() > 0,
519  ExcMessage("The number of polynomials must be larger than zero "
520  "for all coordinate directions."));
521  }
522 }
523 
524 
525 
526 template <int dim>
527 void
529  const unsigned int i,
530  std::array<unsigned int, dim> &indices) const
531 {
532 #ifdef DEBUG
533  unsigned int n_poly = 1;
534  for (unsigned int d = 0; d < dim; ++d)
535  n_poly *= polynomials[d].size();
536  Assert(i < n_poly, ExcInternalError());
537 #endif
538 
539  if (dim == 0)
540  {
541  }
542  else if (dim == 1)
543  internal::compute_tensor_index(i,
544  polynomials[0].size(),
545  0 /*not used*/,
546  indices);
547  else
548  internal::compute_tensor_index(i,
549  polynomials[0].size(),
550  polynomials[1].size(),
551  indices);
552 }
553 
554 
555 
556 template <>
557 void
559  std::array<unsigned int, 0> &) const
560 {
561  constexpr int dim = 0;
562  AssertThrow(dim > 0, ExcNotImplemented());
563 }
564 
565 
566 
567 template <int dim>
568 double
570  const Point<dim> & p) const
571 {
572  std::array<unsigned int, dim> indices;
573  compute_index(i, indices);
574 
575  double value = 1.;
576  for (unsigned int d = 0; d < dim; ++d)
577  value *= polynomials[d][indices[d]].value(p(d));
578 
579  return value;
580 }
581 
582 
583 
584 template <>
585 double
587  const Point<0> &) const
588 {
589  constexpr int dim = 0;
590  AssertThrow(dim > 0, ExcNotImplemented());
591 
592  return {};
593 }
594 
595 
596 
597 template <int dim>
600  const Point<dim> & p) const
601 {
602  std::array<unsigned int, dim> indices;
603  compute_index(i, indices);
604 
605  // compute values and
606  // uni-directional derivatives at
607  // the given point in each
608  // co-ordinate direction
609  std::array<std::array<double, 2>, dim> v;
610  for (unsigned int d = 0; d < dim; ++d)
611  polynomials[d][indices[d]].value(p(d), 1, v[d].data());
612 
613  Tensor<1, dim> grad;
614  for (unsigned int d = 0; d < dim; ++d)
615  {
616  grad[d] = 1.;
617  for (unsigned int x = 0; x < dim; ++x)
618  grad[d] *= v[x][d == x];
619  }
620 
621  return grad;
622 }
623 
624 
625 
626 template <>
629  const Point<0> &) const
630 {
631  constexpr int dim = 0;
632  AssertThrow(dim > 0, ExcNotImplemented());
633 
634  return {};
635 }
636 
637 
638 
639 template <int dim>
642  const Point<dim> & p) const
643 {
644  std::array<unsigned int, dim> indices;
645  compute_index(i, indices);
646 
647  std::array<std::array<double, 3>, dim> v;
648  for (unsigned int d = 0; d < dim; ++d)
649  polynomials[d][indices[d]].value(p(d), 2, v[d].data());
650 
651  Tensor<2, dim> grad_grad;
652  for (unsigned int d1 = 0; d1 < dim; ++d1)
653  for (unsigned int d2 = 0; d2 < dim; ++d2)
654  {
655  grad_grad[d1][d2] = 1.;
656  for (unsigned int x = 0; x < dim; ++x)
657  {
658  unsigned int derivative = 0;
659  if (d1 == x || d2 == x)
660  {
661  if (d1 == d2)
662  derivative = 2;
663  else
664  derivative = 1;
665  }
666  grad_grad[d1][d2] *= v[x][derivative];
667  }
668  }
669 
670  return grad_grad;
671 }
672 
673 
674 
675 template <>
678  const Point<0> &) const
679 {
680  constexpr int dim = 0;
681  AssertThrow(dim > 0, ExcNotImplemented());
682 
683  return {};
684 }
685 
686 
687 
688 template <int dim>
689 void
691  const Point<dim> & p,
692  std::vector<double> & values,
693  std::vector<Tensor<1, dim>> &grads,
694  std::vector<Tensor<2, dim>> &grad_grads,
695  std::vector<Tensor<3, dim>> &third_derivatives,
696  std::vector<Tensor<4, dim>> &fourth_derivatives) const
697 {
698  Assert(values.size() == this->n() || values.size() == 0,
699  ExcDimensionMismatch2(values.size(), this->n(), 0));
700  Assert(grads.size() == this->n() || grads.size() == 0,
701  ExcDimensionMismatch2(grads.size(), this->n(), 0));
702  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
703  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
704  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
705  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
706  Assert(fourth_derivatives.size() == this->n() ||
707  fourth_derivatives.size() == 0,
708  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
709 
710  const bool update_values = (values.size() == this->n()),
711  update_grads = (grads.size() == this->n()),
712  update_grad_grads = (grad_grads.size() == this->n()),
713  update_3rd_derivatives = (third_derivatives.size() == this->n()),
714  update_4th_derivatives = (fourth_derivatives.size() == this->n());
715 
716  // check how many
717  // values/derivatives we have to
718  // compute
719  unsigned int n_values_and_derivatives = 0;
720  if (update_values)
721  n_values_and_derivatives = 1;
722  if (update_grads)
723  n_values_and_derivatives = 2;
724  if (update_grad_grads)
725  n_values_and_derivatives = 3;
727  n_values_and_derivatives = 4;
728  if (update_4th_derivatives)
729  n_values_and_derivatives = 5;
730 
731  // compute the values (and
732  // derivatives, if necessary) of
733  // all polynomials at this
734  // evaluation point
735  std::size_t max_n_polynomials = 0;
736  for (unsigned int d = 0; d < dim; ++d)
737  max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
738 
739  // 5 is enough to store values and derivatives in all supported cases
740  Table<2, std::array<double, 5>> v(dim, max_n_polynomials);
741  for (unsigned int d = 0; d < dim; ++d)
742  for (unsigned int i = 0; i < polynomials[d].size(); ++i)
743  polynomials[d][i].value(p(d),
744  n_values_and_derivatives - 1,
745  v(d, i).data());
746 
747  for (unsigned int i = 0; i < this->n(); ++i)
748  {
749  // first get the
750  // one-dimensional indices of
751  // this particular tensor
752  // product polynomial
753  std::array<unsigned int, dim> indices;
754  compute_index(i, indices);
755 
756  if (update_values)
757  {
758  values[i] = 1;
759  for (unsigned int x = 0; x < dim; ++x)
760  values[i] *= v(x, indices[x])[0];
761  }
762 
763  if (update_grads)
764  for (unsigned int d = 0; d < dim; ++d)
765  {
766  grads[i][d] = 1.;
767  for (unsigned int x = 0; x < dim; ++x)
768  grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0];
769  }
770 
771  if (update_grad_grads)
772  for (unsigned int d1 = 0; d1 < dim; ++d1)
773  for (unsigned int d2 = 0; d2 < dim; ++d2)
774  {
775  grad_grads[i][d1][d2] = 1.;
776  for (unsigned int x = 0; x < dim; ++x)
777  {
778  unsigned int derivative = 0;
779  if (d1 == x)
780  ++derivative;
781  if (d2 == x)
782  ++derivative;
783 
784  grad_grads[i][d1][d2] *= v(x, indices[x])[derivative];
785  }
786  }
787 
789  for (unsigned int d1 = 0; d1 < dim; ++d1)
790  for (unsigned int d2 = 0; d2 < dim; ++d2)
791  for (unsigned int d3 = 0; d3 < dim; ++d3)
792  {
793  third_derivatives[i][d1][d2][d3] = 1.;
794  for (unsigned int x = 0; x < dim; ++x)
795  {
796  unsigned int derivative = 0;
797  if (d1 == x)
798  ++derivative;
799  if (d2 == x)
800  ++derivative;
801  if (d3 == x)
802  ++derivative;
803 
804  third_derivatives[i][d1][d2][d3] *=
805  v(x, indices[x])[derivative];
806  }
807  }
808 
809  if (update_4th_derivatives)
810  for (unsigned int d1 = 0; d1 < dim; ++d1)
811  for (unsigned int d2 = 0; d2 < dim; ++d2)
812  for (unsigned int d3 = 0; d3 < dim; ++d3)
813  for (unsigned int d4 = 0; d4 < dim; ++d4)
814  {
815  fourth_derivatives[i][d1][d2][d3][d4] = 1.;
816  for (unsigned int x = 0; x < dim; ++x)
817  {
818  unsigned int derivative = 0;
819  if (d1 == x)
820  ++derivative;
821  if (d2 == x)
822  ++derivative;
823  if (d3 == x)
824  ++derivative;
825  if (d4 == x)
826  ++derivative;
827 
828  fourth_derivatives[i][d1][d2][d3][d4] *=
829  v(x, indices[x])[derivative];
830  }
831  }
832  }
833 }
834 
835 
836 
837 template <>
838 void
840  std::vector<double> &,
841  std::vector<Tensor<1, 0>> &,
842  std::vector<Tensor<2, 0>> &,
843  std::vector<Tensor<3, 0>> &,
844  std::vector<Tensor<4, 0>> &) const
845 {
846  constexpr int dim = 0;
847  AssertThrow(dim > 0, ExcNotImplemented());
848 }
849 
850 
851 
852 template <int dim>
853 unsigned int
855  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
856 {
857  unsigned int y = 1;
858  for (unsigned int d = 0; d < dim; ++d)
859  y *= pols[d].size();
860  return y;
861 }
862 
863 
864 
865 template <>
866 unsigned int
868  const std::vector<std::vector<Polynomials::Polynomial<double>>> &)
869 {
870  constexpr int dim = 0;
871  AssertThrow(dim > 0, ExcNotImplemented());
872 
873  return {};
874 }
875 
876 
877 
878 template <int dim>
879 std::unique_ptr<ScalarPolynomialsBase<dim>>
881 {
882  return std::make_unique<AnisotropicPolynomials<dim>>(*this);
883 }
884 
885 
886 
887 /* ------------------- explicit instantiations -------------- */
892 
893 template class TensorProductPolynomials<
894  1,
896 template class TensorProductPolynomials<
897  2,
898  Polynomials::PiecewisePolynomial<double>>;
899 template class TensorProductPolynomials<
900  3,
901  Polynomials::PiecewisePolynomial<double>>;
902 
903 template class AnisotropicPolynomials<0>;
904 template class AnisotropicPolynomials<1>;
905 template class AnisotropicPolynomials<2>;
906 template class AnisotropicPolynomials<3>;
907 
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)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
static ::ExceptionBase & ExcMessage(std::string arg1)
Third derivatives of shape functions.
#define Assert(cond, exc)
Definition: exceptions.h:1411
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:363
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_NAMESPACE_OPEN
Definition: config.h:362
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:687
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
T max(const T &t, const MPI_Comm &mpi_communicator)
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