Reference documentation for deal.II version GIT 77cf9aa859 2023-06-10 04:00:02+00:00
\(\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 
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;
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  // coordinate direction
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 
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<ndarray<double, dim, 5>, 20> values_1d(
355  n_polynomials);
356  if (n_values_and_derivatives == 1)
357  for (unsigned int i = 0; i < n_polynomials; ++i)
358  for (unsigned int d = 0; d < dim; ++d)
359  values_1d[i][d][0] = polynomials[i].value(p(d));
360  else
361  for (unsigned int i = 0; i < n_polynomials; ++i)
362  for (unsigned d = 0; d < dim; ++d)
363  polynomials[i].value(p(d),
364  n_values_and_derivatives,
365  values_1d[i][d].data());
366 
367  unsigned int indices[3];
368  unsigned int ind = 0;
369  for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2])
370  for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1);
371  ++indices[1])
372  if (n_values_and_derivatives == 1)
373  for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
374  {
375  double value = values_1d[indices[0]][0][0];
376  for (unsigned int d = 1; d < dim; ++d)
377  value *= values_1d[indices[d]][d][0];
378  values[index_map_inverse[ind]] = value;
379  }
380  else
381  for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
382  {
383  const unsigned int i = index_map_inverse[ind];
384 
385  if (update_values)
386  {
387  double value = values_1d[indices[0]][0][0];
388  for (unsigned int x = 1; x < dim; ++x)
389  value *= values_1d[indices[x]][x][0];
390  values[i] = value;
391  }
392 
393  if (update_grads)
394  for (unsigned int d = 0; d < dim; ++d)
395  {
396  double grad = 1.;
397  for (unsigned int x = 0; x < dim; ++x)
398  grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
399  grads[i][d] = grad;
400  }
401 
402  if (update_grad_grads)
403  for (unsigned int d1 = 0; d1 < dim; ++d1)
404  for (unsigned int d2 = 0; d2 < dim; ++d2)
405  {
406  double der2 = 1.;
407  for (unsigned int x = 0; x < dim; ++x)
408  {
409  unsigned int derivative = 0;
410  if (d1 == x)
411  ++derivative;
412  if (d2 == x)
413  ++derivative;
414 
415  der2 *= values_1d[indices[x]][x][derivative];
416  }
417  grad_grads[i][d1][d2] = der2;
418  }
419 
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 template <int dim, typename PolynomialType>
506 std::vector<PolynomialType>
508  const
509 {
510  return polynomials;
511 }
512 
513 
514 
515 /* ------------------- AnisotropicPolynomials -------------- */
516 
517 
518 template <int dim>
520  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
521  : ScalarPolynomialsBase<dim>(1, get_n_tensor_pols(pols))
522  , polynomials(pols)
523 {
524  Assert(pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
525  for (const auto &pols_d : pols)
526  {
527  (void)pols_d;
528  Assert(pols_d.size() > 0,
529  ExcMessage("The number of polynomials must be larger than zero "
530  "for all coordinate directions."));
531  }
532 }
533 
534 
535 
536 template <int dim>
537 void
539  const unsigned int i,
540  std::array<unsigned int, dim> &indices) const
541 {
542 #ifdef DEBUG
543  unsigned int n_poly = 1;
544  for (unsigned int d = 0; d < dim; ++d)
545  n_poly *= polynomials[d].size();
546  Assert(i < n_poly, ExcInternalError());
547 #endif
548 
549  if (dim == 0)
550  {}
551  else if (dim == 1)
552  internal::compute_tensor_index(i,
553  polynomials[0].size(),
554  0 /*not used*/,
555  indices);
556  else
557  internal::compute_tensor_index(i,
558  polynomials[0].size(),
559  polynomials[1].size(),
560  indices);
561 }
562 
563 
564 
565 template <>
566 void
568  std::array<unsigned int, 0> &) const
569 {
570  constexpr int dim = 0;
571  AssertThrow(dim > 0, ExcNotImplemented());
572 }
573 
574 
575 
576 template <int dim>
577 double
579  const Point<dim> & p) const
580 {
581  std::array<unsigned int, dim> indices;
582  compute_index(i, indices);
583 
584  double value = 1.;
585  for (unsigned int d = 0; d < dim; ++d)
586  value *= polynomials[d][indices[d]].value(p(d));
587 
588  return value;
589 }
590 
591 
592 
593 template <>
594 double
596  const Point<0> &) const
597 {
598  constexpr int dim = 0;
599  AssertThrow(dim > 0, ExcNotImplemented());
600 
601  return {};
602 }
603 
604 
605 
606 template <int dim>
609  const Point<dim> & p) const
610 {
611  std::array<unsigned int, dim> indices;
612  compute_index(i, indices);
613 
614  // compute values and
615  // uni-directional derivatives at
616  // the given point in each
617  // coordinate direction
619  for (unsigned int d = 0; d < dim; ++d)
620  polynomials[d][indices[d]].value(p(d), 1, v[d].data());
621 
622  Tensor<1, dim> grad;
623  for (unsigned int d = 0; d < dim; ++d)
624  {
625  grad[d] = 1.;
626  for (unsigned int x = 0; x < dim; ++x)
627  grad[d] *= v[x][d == x];
628  }
629 
630  return grad;
631 }
632 
633 
634 
635 template <>
638  const Point<0> &) const
639 {
640  constexpr int dim = 0;
641  AssertThrow(dim > 0, ExcNotImplemented());
642 
643  return {};
644 }
645 
646 
647 
648 template <int dim>
651  const Point<dim> & p) const
652 {
653  std::array<unsigned int, dim> indices;
654  compute_index(i, indices);
655 
657  for (unsigned int d = 0; d < dim; ++d)
658  polynomials[d][indices[d]].value(p(d), 2, v[d].data());
659 
660  Tensor<2, dim> grad_grad;
661  for (unsigned int d1 = 0; d1 < dim; ++d1)
662  for (unsigned int d2 = 0; d2 < dim; ++d2)
663  {
664  grad_grad[d1][d2] = 1.;
665  for (unsigned int x = 0; x < dim; ++x)
666  {
667  unsigned int derivative = 0;
668  if (d1 == x || d2 == x)
669  {
670  if (d1 == d2)
671  derivative = 2;
672  else
673  derivative = 1;
674  }
675  grad_grad[d1][d2] *= v[x][derivative];
676  }
677  }
678 
679  return grad_grad;
680 }
681 
682 
683 
684 template <>
687  const Point<0> &) const
688 {
689  constexpr int dim = 0;
690  AssertThrow(dim > 0, ExcNotImplemented());
691 
692  return {};
693 }
694 
695 
696 
697 template <int dim>
698 void
700  const Point<dim> & p,
701  std::vector<double> & values,
702  std::vector<Tensor<1, dim>> &grads,
703  std::vector<Tensor<2, dim>> &grad_grads,
704  std::vector<Tensor<3, dim>> &third_derivatives,
705  std::vector<Tensor<4, dim>> &fourth_derivatives) const
706 {
707  Assert(values.size() == this->n() || values.size() == 0,
708  ExcDimensionMismatch2(values.size(), this->n(), 0));
709  Assert(grads.size() == this->n() || grads.size() == 0,
710  ExcDimensionMismatch2(grads.size(), this->n(), 0));
711  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
712  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
713  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
714  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
715  Assert(fourth_derivatives.size() == this->n() ||
716  fourth_derivatives.size() == 0,
717  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
718 
719  const bool update_values = (values.size() == this->n()),
720  update_grads = (grads.size() == this->n()),
721  update_grad_grads = (grad_grads.size() == this->n()),
722  update_3rd_derivatives = (third_derivatives.size() == this->n()),
723  update_4th_derivatives = (fourth_derivatives.size() == this->n());
724 
725  // check how many
726  // values/derivatives we have to
727  // compute
728  unsigned int n_values_and_derivatives = 0;
729  if (update_values)
730  n_values_and_derivatives = 1;
731  if (update_grads)
732  n_values_and_derivatives = 2;
733  if (update_grad_grads)
734  n_values_and_derivatives = 3;
736  n_values_and_derivatives = 4;
737  if (update_4th_derivatives)
738  n_values_and_derivatives = 5;
739 
740  // compute the values (and
741  // derivatives, if necessary) of
742  // all polynomials at this
743  // evaluation point
744  std::size_t max_n_polynomials = 0;
745  for (unsigned int d = 0; d < dim; ++d)
746  max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size());
747 
748  // 5 is enough to store values and derivatives in all supported cases
749  Table<2, std::array<double, 5>> v(dim, max_n_polynomials);
750  for (unsigned int d = 0; d < dim; ++d)
751  for (unsigned int i = 0; i < polynomials[d].size(); ++i)
752  polynomials[d][i].value(p(d),
753  n_values_and_derivatives - 1,
754  v(d, i).data());
755 
756  for (unsigned int i = 0; i < this->n(); ++i)
757  {
758  // first get the
759  // one-dimensional indices of
760  // this particular tensor
761  // product polynomial
762  std::array<unsigned int, dim> indices;
763  compute_index(i, indices);
764 
765  if (update_values)
766  {
767  values[i] = 1;
768  for (unsigned int x = 0; x < dim; ++x)
769  values[i] *= v(x, indices[x])[0];
770  }
771 
772  if (update_grads)
773  for (unsigned int d = 0; d < dim; ++d)
774  {
775  grads[i][d] = 1.;
776  for (unsigned int x = 0; x < dim; ++x)
777  grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0];
778  }
779 
780  if (update_grad_grads)
781  for (unsigned int d1 = 0; d1 < dim; ++d1)
782  for (unsigned int d2 = 0; d2 < dim; ++d2)
783  {
784  grad_grads[i][d1][d2] = 1.;
785  for (unsigned int x = 0; x < dim; ++x)
786  {
787  unsigned int derivative = 0;
788  if (d1 == x)
789  ++derivative;
790  if (d2 == x)
791  ++derivative;
792 
793  grad_grads[i][d1][d2] *= v(x, indices[x])[derivative];
794  }
795  }
796 
798  for (unsigned int d1 = 0; d1 < dim; ++d1)
799  for (unsigned int d2 = 0; d2 < dim; ++d2)
800  for (unsigned int d3 = 0; d3 < dim; ++d3)
801  {
802  third_derivatives[i][d1][d2][d3] = 1.;
803  for (unsigned int x = 0; x < dim; ++x)
804  {
805  unsigned int derivative = 0;
806  if (d1 == x)
807  ++derivative;
808  if (d2 == x)
809  ++derivative;
810  if (d3 == x)
811  ++derivative;
812 
813  third_derivatives[i][d1][d2][d3] *=
814  v(x, indices[x])[derivative];
815  }
816  }
817 
818  if (update_4th_derivatives)
819  for (unsigned int d1 = 0; d1 < dim; ++d1)
820  for (unsigned int d2 = 0; d2 < dim; ++d2)
821  for (unsigned int d3 = 0; d3 < dim; ++d3)
822  for (unsigned int d4 = 0; d4 < dim; ++d4)
823  {
824  fourth_derivatives[i][d1][d2][d3][d4] = 1.;
825  for (unsigned int x = 0; x < dim; ++x)
826  {
827  unsigned int derivative = 0;
828  if (d1 == x)
829  ++derivative;
830  if (d2 == x)
831  ++derivative;
832  if (d3 == x)
833  ++derivative;
834  if (d4 == x)
835  ++derivative;
836 
837  fourth_derivatives[i][d1][d2][d3][d4] *=
838  v(x, indices[x])[derivative];
839  }
840  }
841  }
842 }
843 
844 
845 
846 template <>
847 void
849  std::vector<double> &,
850  std::vector<Tensor<1, 0>> &,
851  std::vector<Tensor<2, 0>> &,
852  std::vector<Tensor<3, 0>> &,
853  std::vector<Tensor<4, 0>> &) const
854 {
855  constexpr int dim = 0;
856  AssertThrow(dim > 0, ExcNotImplemented());
857 }
858 
859 
860 
861 template <int dim>
862 unsigned int
864  const std::vector<std::vector<Polynomials::Polynomial<double>>> &pols)
865 {
866  unsigned int y = 1;
867  for (unsigned int d = 0; d < dim; ++d)
868  y *= pols[d].size();
869  return y;
870 }
871 
872 
873 
874 template <>
875 unsigned int
877  const std::vector<std::vector<Polynomials::Polynomial<double>>> &)
878 {
879  constexpr int dim = 0;
880  AssertThrow(dim > 0, ExcNotImplemented());
881 
882  return {};
883 }
884 
885 
886 
887 template <int dim>
888 std::unique_ptr<ScalarPolynomialsBase<dim>>
890 {
891  return std::make_unique<AnisotropicPolynomials<dim>>(*this);
892 }
893 
894 
895 
896 /* ------------------- explicit instantiations -------------- */
901 
902 template class TensorProductPolynomials<
903  1,
905 template class TensorProductPolynomials<
906  2,
908 template class TensorProductPolynomials<
909  3,
911 
912 template class AnisotropicPolynomials<0>;
913 template class AnisotropicPolynomials<1>;
914 template class AnisotropicPolynomials<2>;
915 template class AnisotropicPolynomials<3>;
916 
double compute_value(const unsigned int i, const Point< dim > &p) const override
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() 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
Definition: point.h:112
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
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
double compute_value(const unsigned int i, const Point< dim > &p) const override
virtual std::size_t memory_consumption() const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
std::vector< PolynomialType > get_underlying_polynomials() const
void set_numbering(const std::vector< unsigned int > &renumber)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:475
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:476
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108