21#include <boost/container/small_vector.hpp>
37 template <std::
size_t dim>
39 compute_tensor_index(
const unsigned int,
42 std::array<unsigned int, dim> &)
48 compute_tensor_index(
const unsigned int n,
51 std::array<unsigned int, 1> &indices)
57 compute_tensor_index(
const unsigned int n,
58 const unsigned int n_pols_0,
60 std::array<unsigned int, 2> &indices)
62 indices[0] = n % n_pols_0;
63 indices[1] = n / n_pols_0;
67 compute_tensor_index(
const unsigned int n,
68 const unsigned int n_pols_0,
69 const unsigned int n_pols_1,
70 std::array<unsigned int, 3> &indices)
72 indices[0] = n % n_pols_0;
73 indices[1] = (n / n_pols_0) % n_pols_1;
74 indices[2] = n / (n_pols_0 * n_pols_1);
81template <
int dim,
typename PolynomialType>
85 std::array<unsigned int, dim> &indices)
const
87 if constexpr (dim == 0)
95 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
97 internal::compute_tensor_index(index_map[i],
106template <
int dim,
typename PolynomialType>
109 std::ostream &out)
const
111 if constexpr (dim == 0)
118 std::array<unsigned int, dim> ix;
119 for (
unsigned int i = 0; i < this->n(); ++i)
121 compute_index(i, ix);
123 for (
unsigned int d = 0; d < dim; ++d)
133inline const std::vector<unsigned int> &
142inline const std::vector<unsigned int> &
145 return index_map_inverse;
153 const std::vector<unsigned int> &renumber)
155 Assert(renumber.size() == index_map.size(),
158 index_map = renumber;
159 for (
unsigned int i = 0; i < index_map.size(); ++i)
160 index_map_inverse[index_map[i]] = i;
174template <
int dim,
typename PolynomialType>
177 const std::vector<unsigned int> &renumber)
179 Assert(renumber.size() == index_map.size(),
182 index_map = renumber;
183 for (
unsigned int i = 0; i < index_map.size(); ++i)
184 index_map_inverse[index_map[i]] = i;
192 const std::vector<unsigned int> &)
199template <
int dim,
typename PolynomialType>
202 const unsigned int i,
205 if constexpr (dim == 0)
214 std::array<unsigned int, dim> indices;
215 compute_index(i, indices);
218 for (
unsigned int d = 0; d < dim; ++d)
219 value *= polynomials[indices[d]].value(p[d]);
227template <
int dim,
typename PolynomialType>
230 const unsigned int i,
233 if constexpr (dim == 0)
242 std::array<unsigned int, dim> indices;
243 compute_index(i, indices);
251 std::vector<double> tmp(2);
252 for (
unsigned int d = 0; d < dim; ++d)
254 polynomials[indices[d]].value(p[d], tmp);
261 for (
unsigned int d = 0;
d < dim; ++
d)
264 for (
unsigned int x = 0; x < dim; ++x)
265 grad[d] *= v[x][d == x];
274template <
int dim,
typename PolynomialType>
277 const unsigned int i,
280 if constexpr (dim == 0)
289 std::array<unsigned int, dim> indices;
290 compute_index(i, indices);
294 std::vector<double> tmp(3);
295 for (
unsigned int d = 0; d < dim; ++d)
297 polynomials[indices[d]].value(p[d], tmp);
305 for (
unsigned int d1 = 0; d1 < dim; ++d1)
306 for (
unsigned int d2 = 0; d2 < dim; ++d2)
308 grad_grad[d1][d2] = 1.;
309 for (
unsigned int x = 0; x < dim; ++x)
311 unsigned int derivative = 0;
312 if (d1 == x || d2 == x)
319 grad_grad[d1][d2] *= v[x][derivative];
339 template <
int dim, std::
size_t dim1>
342 const unsigned int n_derivatives,
345 const unsigned int size_x,
346 const boost::container::small_vector<std::array<unsigned int, dim1>, 64>
348 const std::vector<unsigned int> &index_map,
349 std::vector<double> &values,
355 const bool update_values = (values.size() == indices.size() * size_x);
356 const bool update_grads = (grads.size() == indices.size() * size_x);
357 const bool update_grad_grads =
358 (grad_grads.size() == indices.size() * size_x);
360 (third_derivatives.size() == indices.size() * size_x);
361 const bool update_4th_derivatives =
362 (fourth_derivatives.size() == indices.size() * size_x);
366 if (n_derivatives == 0)
367 for (
unsigned int i = 0, i1 = 0; i1 < indices.size(); ++i1)
369 double value_outer = 1.;
370 if constexpr (dim > 1)
371 for (
unsigned int d = 1; d < dim; ++d)
372 value_outer *= values_1d[indices[i1][d - 1]][0][d];
373 if (index_map.empty())
374 for (
unsigned int ix = 0; ix < size_x; ++ix, ++i)
375 values[i] = value_outer * values_1d[ix][0][0];
377 for (
unsigned int ix = 0; ix < size_x; ++ix, ++i)
378 values[index_map[i]] = value_outer * values_1d[ix][0][0];
381 for (
unsigned int iy = 0, i1 = 0; i1 < indices.size(); ++i1)
384 std::array<double, dim + (dim * (dim - 1)) / 2> value_outer;
386 if constexpr (dim > 1)
388 for (
unsigned int x = 1; x < dim; ++x)
389 value_outer[0] *= values_1d[indices[i1][x - 1]][0][x];
390 for (
unsigned int d = 1; d < dim; ++d)
392 value_outer[d] = values_1d[indices[i1][d - 1]][1][d];
393 for (
unsigned int x = 1; x < dim; ++x)
395 value_outer[d] *= values_1d[indices[i1][x - 1]][0][x];
397 for (
unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
398 for (
unsigned int d2 = d1; d2 < dim; ++d2, ++count)
400 value_outer[count] = 1.;
401 for (
unsigned int x = 1; x < dim; ++x)
403 unsigned int derivative = 0;
409 value_outer[count] *=
410 values_1d[indices[i1][x - 1]][derivative][x];
417 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
419 std::array<double, 3> val_x{{values_1d[ix][0][0],
421 values_1d[ix][2][0]}};
422 const unsigned int index =
423 (index_map.empty() ? i : index_map[i]);
426 values[
index] = value_outer[0] * val_x[0];
430 grads[
index][0] = value_outer[0] * val_x[1];
431 if constexpr (dim > 1)
432 for (
unsigned int d = 1; d < dim; ++d)
433 grads[
index][d] = value_outer[d] * val_x[0];
436 if (update_grad_grads)
438 grad_grads[
index][0][0] = value_outer[0] * val_x[2];
439 if constexpr (dim > 1)
441 for (
unsigned int d = 1; d < dim; ++d)
442 grad_grads[
index][0][d] = grad_grads[
index][d][0] =
443 value_outer[d] * val_x[1];
444 for (
unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
445 for (
unsigned int d2 = d1; d2 < dim; ++d2, ++count)
446 grad_grads[
index][d1][d2] =
447 grad_grads[
index][d2][d1] =
448 value_outer[count] * val_x[0];
455 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
457 const unsigned int index =
458 (index_map.empty() ? i : index_map[i]);
459 std::array<unsigned int, dim> my_indices;
461 if constexpr (dim > 1)
462 for (
unsigned int d = 1; d < dim; ++d)
463 my_indices[d] = indices[i1][d - 1];
464 for (
unsigned int d1 = 0; d1 < dim; ++d1)
465 for (
unsigned int d2 = 0; d2 < dim; ++d2)
466 for (
unsigned int d3 = 0; d3 < dim; ++d3)
469 for (
unsigned int x = 0; x < dim; ++x)
471 unsigned int derivative = 0;
479 der3 *= values_1d[my_indices[x]][derivative][x];
481 third_derivatives[
index][d1][d2][d3] = der3;
485 if (update_4th_derivatives)
486 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
488 const unsigned int index =
489 (index_map.empty() ? i : index_map[i]);
490 std::array<unsigned int, dim> my_indices;
492 if constexpr (dim > 1)
493 for (
unsigned int d = 1; d < dim; ++d)
494 my_indices[d] = indices[i1][d - 1];
495 for (
unsigned int d1 = 0; d1 < dim; ++d1)
496 for (
unsigned int d2 = 0; d2 < dim; ++d2)
497 for (
unsigned int d3 = 0; d3 < dim; ++d3)
498 for (
unsigned int d4 = 0; d4 < dim; ++d4)
501 for (
unsigned int x = 0; x < dim; ++x)
503 unsigned int derivative = 0;
513 der4 *= values_1d[my_indices[x]][derivative][x];
515 fourth_derivatives[
index][d1][d2][d3][d4] = der4;
527template <
int dim,
typename PolynomialType>
531 std::vector<double> &values,
537 if constexpr (dim == 0)
543 (void)third_derivatives;
544 (void)fourth_derivatives;
550 Assert(values.size() == this->n() || values.empty(),
552 Assert(grads.size() == this->n() || grads.empty(),
554 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
556 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
558 Assert(fourth_derivatives.size() == this->n() ||
559 fourth_derivatives.empty(),
563 unsigned int n_derivatives = 0;
564 if (values.size() == this->n())
566 if (grads.size() == this->n())
568 if (grad_grads.size() == this->n())
570 if (third_derivatives.size() == this->n())
572 if (fourth_derivatives.size() == this->n())
578 const unsigned int n_polynomials = polynomials.size();
579 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
581 if constexpr (std::is_same_v<PolynomialType,
584 std::array<double, dim> point_array;
585 for (
unsigned int d = 0; d < dim; ++d)
586 point_array[d] = p[d];
587 for (
unsigned int i = 0; i < n_polynomials; ++i)
588 polynomials[i].values_of_array(point_array,
590 values_1d[i].
data());
593 for (
unsigned int i = 0; i < n_polynomials; ++i)
594 for (
unsigned int d = 0; d < dim; ++d)
596 std::array<double, 5> derivatives;
597 polynomials[i].value(p[d], n_derivatives, derivatives.data());
598 for (
unsigned int j = 0; j <= n_derivatives; ++j)
599 values_1d[i][j][d] = derivatives[j];
604 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
605 boost::container::small_vector<std::array<unsigned int, dim1>, 64>
607 if constexpr (dim > 1)
608 for (
unsigned int d = 1; d < dim; ++d)
610 const unsigned int size = indices.size();
611 for (
unsigned int i = 1; i < n_polynomials; ++i)
612 for (
unsigned int j = 0; j <
size; ++j)
614 std::array<unsigned int, dim1> next_index = indices[j];
615 next_index[d - 1] = i;
616 indices.push_back(next_index);
621 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
637template <
int dim,
typename PolynomialType>
638std::unique_ptr<ScalarPolynomialsBase<dim>>
641 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
646template <
int dim,
typename PolynomialType>
657template <
int dim,
typename PolynomialType>
658std::vector<PolynomialType>
675 , index_map(this->n())
676 , index_map_inverse(this->n())
679 for (
const auto &pols_d : pols)
683 ExcMessage(
"The number of polynomials must be larger than zero "
684 "for all coordinate directions."));
689 for (
unsigned int i = 0; i < this->
n(); ++i)
701 const unsigned int i,
702 std::array<unsigned int, dim> &indices)
const
704 if constexpr (dim == 0)
714 unsigned int n_poly = 1;
715 for (
unsigned int d = 0; d < dim; ++d)
716 n_poly *= polynomials[d].
size();
724 internal::compute_tensor_index(index_map[i],
725 polynomials[0].
size(),
729 internal::compute_tensor_index(index_map[i],
730 polynomials[0].
size(),
731 polynomials[1].
size(),
743 if constexpr (dim == 0)
752 std::array<unsigned int, dim> indices;
753 compute_index(i, indices);
756 for (
unsigned int d = 0; d < dim; ++d)
757 value *= polynomials[d][indices[d]].value(p[d]);
770 if constexpr (dim == 0)
779 std::array<unsigned int, dim> indices;
780 compute_index(i, indices);
787 for (
unsigned int d = 0; d < dim; ++d)
788 polynomials[d][indices[d]].value(p[d], 1, v[d].
data());
791 for (
unsigned int d = 0; d < dim; ++d)
794 for (
unsigned int x = 0; x < dim; ++x)
795 grad[d] *= v[x][d == x];
809 if constexpr (dim == 0)
818 std::array<unsigned int, dim> indices;
819 compute_index(i, indices);
822 for (
unsigned int d = 0; d < dim; ++d)
823 polynomials[d][indices[d]].value(p[d], 2, v[d].
data());
826 for (
unsigned int d1 = 0; d1 < dim; ++d1)
827 for (
unsigned int d2 = 0; d2 < dim; ++d2)
829 grad_grad[d1][d2] = 1.;
830 for (
unsigned int x = 0; x < dim; ++x)
832 unsigned int derivative = 0;
833 if (d1 == x || d2 == x)
840 grad_grad[d1][d2] *= v[x][derivative];
854 std::vector<double> &values,
860 if constexpr (dim == 0)
866 (void)third_derivatives;
867 (void)fourth_derivatives;
872 Assert(values.size() == this->n() || values.empty(),
874 Assert(grads.size() == this->n() || grads.empty(),
876 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
878 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
880 Assert(fourth_derivatives.size() == this->n() ||
881 fourth_derivatives.empty(),
885 unsigned int n_derivatives = 0;
886 if (values.size() == this->n())
888 if (grads.size() == this->n())
890 if (grad_grads.size() == this->n())
892 if (third_derivatives.size() == this->n())
894 if (fourth_derivatives.size() == this->n())
899 std::size_t max_n_polynomials = 0;
900 for (
unsigned int d = 0; d < dim; ++d)
901 max_n_polynomials =
std::max(max_n_polynomials, polynomials[d].
size());
904 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
906 if (n_derivatives == 0)
907 for (
unsigned int d = 0; d < dim; ++d)
908 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
909 values_1d[i][0][d] = polynomials[d][i].value(p[d]);
911 for (
unsigned int d = 0; d < dim; ++d)
912 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
917 std::array<double, 5> derivatives;
918 polynomials[d][i].value(p[d], n_derivatives, derivatives.data());
919 for (
unsigned int j = 0; j <= n_derivatives; ++j)
920 values_1d[i][j][d] = derivatives[j];
924 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
925 boost::container::small_vector<std::array<unsigned int, dim1>, 64>
927 for (
unsigned int d = 1; d < dim; ++d)
929 const unsigned int size = indices.size();
930 for (
unsigned int i = 1; i < polynomials[d].size(); ++i)
931 for (
unsigned int j = 0; j <
size; ++j)
933 std::array<unsigned int, dim1> next_index = indices[j];
934 next_index[d - 1] = i;
935 indices.push_back(next_index);
939 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
942 polynomials[0].
size(),
960 if constexpr (dim == 0)
969 for (
unsigned int d = 0; d < dim; ++d)
978std::unique_ptr<ScalarPolynomialsBase<dim>>
981 return std::make_unique<AnisotropicPolynomials<dim>>(*this);
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
double compute_value(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double > > > &base_polynomials)
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double > > > &pols)
std::vector< unsigned int > index_map
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
const std::vector< unsigned int > & get_numbering() const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
std::vector< unsigned int > index_map_inverse
const std::vector< unsigned int > & get_numbering_inverse() const
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
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
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
std::vector< index_type > data
std::enable_if_t< std::is_fundamental_v< T >, 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)
constexpr T pow(const T base, const int iexp)
void evaluate_tensor_product(const unsigned int n_derivatives, const boost::container::small_vector<::ndarray< double, 5, dim >, 10 > &values_1d, const unsigned int size_x, const boost::container::small_vector< std::array< unsigned int, dim1 >, 64 > &indices, const std::vector< unsigned int > &index_map, 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)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray