23#include <boost/container/small_vector.hpp>
40 template <std::
size_t dim>
42 compute_tensor_index(
const unsigned int,
45 std::array<unsigned int, dim> &)
51 compute_tensor_index(
const unsigned int n,
54 std::array<unsigned int, 1> &indices)
60 compute_tensor_index(
const unsigned int n,
61 const unsigned int n_pols_0,
63 std::array<unsigned int, 2> &indices)
65 indices[0] = n % n_pols_0;
66 indices[1] = n / n_pols_0;
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)
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);
84template <
int dim,
typename PolynomialType>
88 std::array<unsigned int, dim> &indices)
const
90 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
92 internal::compute_tensor_index(index_map[i],
104 std::array<unsigned int, 0> &)
const
106 constexpr int dim = 0;
112template <
int dim,
typename PolynomialType>
115 std::ostream &out)
const
117 std::array<unsigned int, dim> ix;
118 for (
unsigned int i = 0; i < this->n(); ++i)
120 compute_index(i, ix);
122 for (
unsigned int d = 0; d < dim; ++d)
133 std::ostream &)
const
135 constexpr int dim = 0;
141template <
int dim,
typename PolynomialType>
144 const std::vector<unsigned int> &renumber)
146 Assert(renumber.size() == index_map.size(),
149 index_map = renumber;
150 for (
unsigned int i = 0; i < index_map.size(); ++i)
151 index_map_inverse[index_map[i]] = i;
159 const std::vector<unsigned int> &)
161 constexpr int dim = 0;
167template <
int dim,
typename PolynomialType>
170 const unsigned int i,
175 std::array<unsigned int, dim> indices;
176 compute_index(i, indices);
179 for (
unsigned int d = 0; d < dim; ++d)
180 value *= polynomials[indices[d]].value(p(d));
193 constexpr int dim = 0;
201template <
int dim,
typename PolynomialType>
204 const unsigned int i,
207 std::array<unsigned int, dim> indices;
208 compute_index(i, indices);
216 std::vector<double> tmp(2);
217 for (
unsigned int d = 0; d < dim; ++d)
219 polynomials[indices[d]].value(p(d), tmp);
226 for (
unsigned int d = 0; d < dim; ++d)
229 for (
unsigned int x = 0; x < dim; ++x)
230 grad[d] *= v[x][d == x];
244 constexpr int dim = 0;
252template <
int dim,
typename PolynomialType>
255 const unsigned int i,
258 std::array<unsigned int, dim> indices;
259 compute_index(i, indices);
263 std::vector<double> tmp(3);
264 for (
unsigned int d = 0; d < dim; ++d)
266 polynomials[indices[d]].value(p(d), tmp);
274 for (
unsigned int d1 = 0; d1 < dim; ++d1)
275 for (
unsigned int d2 = 0; d2 < dim; ++d2)
277 grad_grad[d1][d2] = 1.;
278 for (
unsigned int x = 0; x < dim; ++x)
280 unsigned int derivative = 0;
281 if (d1 == x || d2 == x)
288 grad_grad[d1][d2] *= v[x][derivative];
308template <
int dim,
typename PolynomialType>
312 std::vector<double> & values,
319 Assert(values.size() == this->n() || values.size() == 0,
321 Assert(grads.size() == this->n() || grads.size() == 0,
323 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
325 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
327 Assert(fourth_derivatives.size() == this->n() ||
328 fourth_derivatives.size() == 0,
332 update_grads = (grads.size() == this->n()),
333 update_grad_grads = (grad_grads.size() == this->n()),
335 update_4th_derivatives = (fourth_derivatives.size() == this->n());
338 unsigned int n_values_and_derivatives = 0;
340 n_values_and_derivatives = 1;
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;
355 const unsigned int n_polynomials = polynomials.size();
356 boost::container::small_vector<ndarray<double, dim, 5>, 20> values_1d(
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));
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());
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);
374 if (n_values_and_derivatives == 1)
375 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
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;
383 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
385 const unsigned int i = index_map_inverse[ind];
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];
396 for (
unsigned int d = 0; d < dim; ++d)
399 for (
unsigned int x = 0; x < dim; ++x)
400 grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
404 if (update_grad_grads)
405 for (
unsigned int d1 = 0; d1 < dim; ++d1)
406 for (
unsigned int d2 = 0; d2 < dim; ++d2)
409 for (
unsigned int x = 0; x < dim; ++x)
411 unsigned int derivative = 0;
417 der2 *= values_1d[indices[x]][x][derivative];
419 grad_grads[i][d1][d2] = der2;
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)
428 for (
unsigned int x = 0; x < dim; ++x)
430 unsigned int derivative = 0;
438 der3 *= values_1d[indices[x]][x][derivative];
440 third_derivatives[i][d1][d2][d3] = der3;
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)
450 for (
unsigned int x = 0; x < dim; ++x)
452 unsigned int derivative = 0;
462 der4 *= values_1d[indices[x]][x][derivative];
464 fourth_derivatives[i][d1][d2][d3][d4] = der4;
475 std::vector<double> &,
481 constexpr int dim = 0;
487template <
int dim,
typename PolynomialType>
488std::unique_ptr<ScalarPolynomialsBase<dim>>
491 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
496template <
int dim,
typename PolynomialType>
507template <
int dim,
typename PolynomialType>
508std::vector<PolynomialType>
527 for (
const auto &pols_d : pols)
531 ExcMessage(
"The number of polynomials must be larger than zero "
532 "for all coordinate directions."));
541 const unsigned int i,
542 std::array<unsigned int, dim> &indices)
const
545 unsigned int n_poly = 1;
546 for (
unsigned int d = 0; d < dim; ++d)
547 n_poly *= polynomials[d].size();
554 internal::compute_tensor_index(i,
555 polynomials[0].size(),
559 internal::compute_tensor_index(i,
560 polynomials[0].size(),
561 polynomials[1].size(),
570 std::array<unsigned int, 0> &)
const
572 constexpr int dim = 0;
583 std::array<unsigned int, dim> indices;
584 compute_index(i, indices);
587 for (
unsigned int d = 0; d < dim; ++d)
588 value *= polynomials[d][indices[d]].value(p(d));
600 constexpr int dim = 0;
613 std::array<unsigned int, dim> indices;
614 compute_index(i, indices);
621 for (
unsigned int d = 0; d < dim; ++d)
622 polynomials[d][indices[d]].value(p(d), 1, v[d].data());
625 for (
unsigned int d = 0; d < dim; ++d)
628 for (
unsigned int x = 0; x < dim; ++x)
629 grad[d] *= v[x][d == x];
642 constexpr int dim = 0;
655 std::array<unsigned int, dim> indices;
656 compute_index(i, indices);
659 for (
unsigned int d = 0; d < dim; ++d)
660 polynomials[d][indices[d]].value(p(d), 2, v[d].data());
663 for (
unsigned int d1 = 0; d1 < dim; ++d1)
664 for (
unsigned int d2 = 0; d2 < dim; ++d2)
666 grad_grad[d1][d2] = 1.;
667 for (
unsigned int x = 0; x < dim; ++x)
669 unsigned int derivative = 0;
670 if (d1 == x || d2 == x)
677 grad_grad[d1][d2] *= v[x][derivative];
691 constexpr int dim = 0;
703 std::vector<double> & values,
709 Assert(values.size() == this->n() || values.size() == 0,
711 Assert(grads.size() == this->n() || grads.size() == 0,
713 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
715 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
717 Assert(fourth_derivatives.size() == this->n() ||
718 fourth_derivatives.size() == 0,
722 update_grads = (grads.size() == this->n()),
723 update_grad_grads = (grad_grads.size() == this->n()),
725 update_4th_derivatives = (fourth_derivatives.size() == this->n());
730 unsigned int n_values_and_derivatives = 0;
732 n_values_and_derivatives = 1;
734 n_values_and_derivatives = 2;
735 if (update_grad_grads)
736 n_values_and_derivatives = 3;
738 n_values_and_derivatives = 4;
739 if (update_4th_derivatives)
740 n_values_and_derivatives = 5;
746 std::size_t max_n_polynomials = 0;
747 for (
unsigned int d = 0; d < dim; ++d)
748 max_n_polynomials =
std::max(max_n_polynomials, polynomials[d].size());
752 for (
unsigned int d = 0; d < dim; ++d)
753 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
754 polynomials[d][i].value(p(d),
755 n_values_and_derivatives - 1,
758 for (
unsigned int i = 0; i < this->n(); ++i)
764 std::array<unsigned int, dim> indices;
765 compute_index(i, indices);
770 for (
unsigned int x = 0; x < dim; ++x)
771 values[i] *= v(x, indices[x])[0];
775 for (
unsigned int d = 0; d < dim; ++d)
778 for (
unsigned int x = 0; x < dim; ++x)
779 grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0];
782 if (update_grad_grads)
783 for (
unsigned int d1 = 0; d1 < dim; ++d1)
784 for (
unsigned int d2 = 0; d2 < dim; ++d2)
786 grad_grads[i][d1][d2] = 1.;
787 for (
unsigned int x = 0; x < dim; ++x)
789 unsigned int derivative = 0;
795 grad_grads[i][d1][d2] *= v(x, indices[x])[derivative];
800 for (
unsigned int d1 = 0; d1 < dim; ++d1)
801 for (
unsigned int d2 = 0; d2 < dim; ++d2)
802 for (
unsigned int d3 = 0; d3 < dim; ++d3)
804 third_derivatives[i][d1][d2][d3] = 1.;
805 for (
unsigned int x = 0; x < dim; ++x)
807 unsigned int derivative = 0;
815 third_derivatives[i][d1][d2][d3] *=
816 v(x, indices[x])[derivative];
820 if (update_4th_derivatives)
821 for (
unsigned int d1 = 0; d1 < dim; ++d1)
822 for (
unsigned int d2 = 0; d2 < dim; ++d2)
823 for (
unsigned int d3 = 0; d3 < dim; ++d3)
824 for (
unsigned int d4 = 0; d4 < dim; ++d4)
826 fourth_derivatives[i][d1][d2][d3][d4] = 1.;
827 for (
unsigned int x = 0; x < dim; ++x)
829 unsigned int derivative = 0;
839 fourth_derivatives[i][d1][d2][d3][d4] *=
840 v(x, indices[x])[derivative];
851 std::vector<double> &,
857 constexpr int dim = 0;
869 for (
unsigned int d = 0; d < dim; ++d)
881 constexpr int dim = 0;
890std::unique_ptr<ScalarPolynomialsBase<dim>>
893 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
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)
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
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
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
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray