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 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
89 internal::compute_tensor_index(index_map[i],
101 std::array<unsigned int, 0> &)
const
108template <
int dim,
typename PolynomialType>
111 std::ostream &out)
const
113 std::array<unsigned int, dim> ix;
114 for (
unsigned int i = 0; i < this->n(); ++i)
116 compute_index(i, ix);
118 for (
unsigned int d = 0; d < dim; ++d)
129 std::ostream &)
const
136template <
int dim,
typename PolynomialType>
139 const std::vector<unsigned int> &renumber)
141 Assert(renumber.size() == index_map.size(),
144 index_map = renumber;
145 for (
unsigned int i = 0; i < index_map.size(); ++i)
146 index_map_inverse[index_map[i]] = i;
154 const std::vector<unsigned int> &)
161template <
int dim,
typename PolynomialType>
164 const unsigned int i,
169 std::array<unsigned int, dim> indices;
170 compute_index(i, indices);
173 for (
unsigned int d = 0; d < dim; ++d)
174 value *= polynomials[indices[d]].value(p[d]);
194template <
int dim,
typename PolynomialType>
197 const unsigned int i,
200 std::array<unsigned int, dim> indices;
201 compute_index(i, indices);
209 std::vector<double> tmp(2);
210 for (
unsigned int d = 0; d < dim; ++d)
212 polynomials[indices[d]].value(p[d], tmp);
219 for (
unsigned int d = 0;
d < dim; ++
d)
222 for (
unsigned int x = 0; x < dim; ++x)
223 grad[
d] *= v[x][
d == x];
244template <
int dim,
typename PolynomialType>
247 const unsigned int i,
250 std::array<unsigned int, dim> indices;
251 compute_index(i, indices);
255 std::vector<double> tmp(3);
256 for (
unsigned int d = 0; d < dim; ++d)
258 polynomials[indices[d]].value(p[d], tmp);
266 for (
unsigned int d1 = 0; d1 < dim; ++d1)
267 for (
unsigned int d2 = 0; d2 < dim; ++d2)
269 grad_grad[d1][d2] = 1.;
270 for (
unsigned int x = 0; x < dim; ++x)
272 unsigned int derivative = 0;
273 if (d1 == x || d2 == x)
280 grad_grad[d1][d2] *= v[x][derivative];
310 template <
int dim, std::
size_t dim1>
313 const unsigned int n_derivatives,
316 const unsigned int size_x,
317 const boost::container::small_vector<std::array<unsigned int, dim1>, 64>
319 const std::vector<unsigned int> &index_map,
320 std::vector<double> &values,
326 const bool update_values = (values.size() == indices.size() * size_x);
327 const bool update_grads = (grads.size() == indices.size() * size_x);
328 const bool update_grad_grads =
329 (grad_grads.size() == indices.size() * size_x);
331 (third_derivatives.size() == indices.size() * size_x);
332 const bool update_4th_derivatives =
333 (fourth_derivatives.size() == indices.size() * size_x);
337 if (n_derivatives == 0)
338 for (
unsigned int i = 0, i1 = 0; i1 < indices.size(); ++i1)
340 double value_outer = 1.;
341 if constexpr (dim > 1)
342 for (
unsigned int d = 1; d < dim; ++d)
343 value_outer *= values_1d[indices[i1][d - 1]][0][d];
344 if (index_map.empty())
345 for (
unsigned int ix = 0; ix < size_x; ++ix, ++i)
346 values[i] = value_outer * values_1d[ix][0][0];
348 for (
unsigned int ix = 0; ix < size_x; ++ix, ++i)
349 values[index_map[i]] = value_outer * values_1d[ix][0][0];
352 for (
unsigned int iy = 0, i1 = 0; i1 < indices.size(); ++i1)
355 std::array<double, dim + (dim * (dim - 1)) / 2> value_outer;
357 if constexpr (dim > 1)
359 for (
unsigned int x = 1; x < dim; ++x)
360 value_outer[0] *= values_1d[indices[i1][x - 1]][0][x];
361 for (
unsigned int d = 1; d < dim; ++d)
363 value_outer[d] = values_1d[indices[i1][d - 1]][1][d];
364 for (
unsigned int x = 1; x < dim; ++x)
366 value_outer[d] *= values_1d[indices[i1][x - 1]][0][x];
368 for (
unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
369 for (
unsigned int d2 = d1; d2 < dim; ++d2, ++count)
371 value_outer[count] = 1.;
372 for (
unsigned int x = 1; x < dim; ++x)
374 unsigned int derivative = 0;
380 value_outer[count] *=
381 values_1d[indices[i1][x - 1]][derivative][x];
388 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
390 std::array<double, 3> val_x{{values_1d[ix][0][0],
392 values_1d[ix][2][0]}};
393 const unsigned int index =
394 (index_map.empty() ? i : index_map[i]);
397 values[
index] = value_outer[0] * val_x[0];
401 grads[
index][0] = value_outer[0] * val_x[1];
402 if constexpr (dim > 1)
403 for (
unsigned int d = 1; d < dim; ++d)
404 grads[
index][d] = value_outer[d] * val_x[0];
407 if (update_grad_grads)
409 grad_grads[
index][0][0] = value_outer[0] * val_x[2];
410 if constexpr (dim > 1)
412 for (
unsigned int d = 1; d < dim; ++d)
413 grad_grads[
index][0][d] = grad_grads[
index][d][0] =
414 value_outer[d] * val_x[1];
415 for (
unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
416 for (
unsigned int d2 = d1; d2 < dim; ++d2, ++count)
417 grad_grads[
index][d1][d2] =
418 grad_grads[
index][d2][d1] =
419 value_outer[count] * val_x[0];
426 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
428 const unsigned int index =
429 (index_map.empty() ? i : index_map[i]);
430 std::array<unsigned int, dim> my_indices;
432 if constexpr (dim > 1)
433 for (
unsigned int d = 1; d < dim; ++d)
434 my_indices[d] = indices[i1][d - 1];
435 for (
unsigned int d1 = 0; d1 < dim; ++d1)
436 for (
unsigned int d2 = 0; d2 < dim; ++d2)
437 for (
unsigned int d3 = 0; d3 < dim; ++d3)
440 for (
unsigned int x = 0; x < dim; ++x)
442 unsigned int derivative = 0;
450 der3 *= values_1d[my_indices[x]][derivative][x];
452 third_derivatives[
index][d1][d2][d3] = der3;
456 if (update_4th_derivatives)
457 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
459 const unsigned int index =
460 (index_map.empty() ? i : index_map[i]);
461 std::array<unsigned int, dim> my_indices;
463 if constexpr (dim > 1)
464 for (
unsigned int d = 1; d < dim; ++d)
465 my_indices[d] = indices[i1][d - 1];
466 for (
unsigned int d1 = 0; d1 < dim; ++d1)
467 for (
unsigned int d2 = 0; d2 < dim; ++d2)
468 for (
unsigned int d3 = 0; d3 < dim; ++d3)
469 for (
unsigned int d4 = 0; d4 < dim; ++d4)
472 for (
unsigned int x = 0; x < dim; ++x)
474 unsigned int derivative = 0;
484 der4 *= values_1d[my_indices[x]][derivative][x];
486 fourth_derivatives[
index][d1][d2][d3][d4] = der4;
498template <
int dim,
typename PolynomialType>
502 std::vector<double> &values,
509 Assert(values.size() == this->n() || values.empty(),
511 Assert(grads.size() == this->n() || grads.empty(),
513 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
515 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
517 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
521 unsigned int n_derivatives = 0;
522 if (values.size() == this->n())
524 if (grads.size() == this->n())
526 if (grad_grads.size() == this->n())
528 if (third_derivatives.size() == this->n())
530 if (fourth_derivatives.size() == this->n())
536 const unsigned int n_polynomials = polynomials.size();
537 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
539 if constexpr (std::is_same<PolynomialType,
542 std::array<double, dim> point_array;
543 for (
unsigned int d = 0; d < dim; ++d)
544 point_array[d] = p[d];
545 for (
unsigned int i = 0; i < n_polynomials; ++i)
546 polynomials[i].values_of_array(point_array,
548 values_1d[i].data());
551 for (
unsigned int i = 0; i < n_polynomials; ++i)
552 for (
unsigned int d = 0; d < dim; ++d)
554 std::array<double, 5> derivatives;
555 polynomials[i].value(p[d], n_derivatives, derivatives.data());
556 for (
unsigned int j = 0; j <= n_derivatives; ++j)
557 values_1d[i][j][d] = derivatives[j];
562 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
563 boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
564 if constexpr (dim > 1)
565 for (
unsigned int d = 1; d < dim; ++d)
567 const unsigned int size = indices.size();
568 for (
unsigned int i = 1; i < n_polynomials; ++i)
569 for (
unsigned int j = 0; j < size; ++j)
571 std::array<unsigned int, dim1> next_index = indices[j];
572 next_index[d - 1] = i;
573 indices.push_back(next_index);
578 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
597 std::vector<double> &,
608template <
int dim,
typename PolynomialType>
609std::unique_ptr<ScalarPolynomialsBase<dim>>
612 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
617template <
int dim,
typename PolynomialType>
628template <
int dim,
typename PolynomialType>
629std::vector<PolynomialType>
648 for (
const auto &pols_d : pols)
652 ExcMessage(
"The number of polynomials must be larger than zero "
653 "for all coordinate directions."));
662 const unsigned int i,
663 std::array<unsigned int, dim> &indices)
const
666 unsigned int n_poly = 1;
667 for (
unsigned int d = 0; d < dim; ++d)
668 n_poly *= polynomials[d].size();
676 internal::compute_tensor_index(i,
677 polynomials[0].size(),
681 internal::compute_tensor_index(i,
682 polynomials[0].size(),
683 polynomials[1].size(),
692 std::array<unsigned int, 0> &)
const
704 std::array<unsigned int, dim> indices;
705 compute_index(i, indices);
708 for (
unsigned int d = 0; d < dim; ++d)
709 value *= polynomials[d][indices[d]].value(p[d]);
733 std::array<unsigned int, dim> indices;
734 compute_index(i, indices);
741 for (
unsigned int d = 0; d < dim; ++d)
742 polynomials[d][indices[d]].value(p[d], 1, v[d].data());
745 for (
unsigned int d = 0; d < dim; ++d)
748 for (
unsigned int x = 0; x < dim; ++x)
749 grad[d] *= v[x][d == x];
774 std::array<unsigned int, dim> indices;
775 compute_index(i, indices);
778 for (
unsigned int d = 0; d < dim; ++d)
779 polynomials[d][indices[d]].value(p[d], 2, v[d].data());
782 for (
unsigned int d1 = 0; d1 < dim; ++d1)
783 for (
unsigned int d2 = 0; d2 < dim; ++d2)
785 grad_grad[d1][d2] = 1.;
786 for (
unsigned int x = 0; x < dim; ++x)
788 unsigned int derivative = 0;
789 if (d1 == x || d2 == x)
796 grad_grad[d1][d2] *= v[x][derivative];
821 std::vector<double> &values,
827 Assert(values.size() == this->n() || values.empty(),
829 Assert(grads.size() == this->n() || grads.empty(),
831 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
833 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
835 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
839 unsigned int n_derivatives = 0;
840 if (values.size() == this->n())
842 if (grads.size() == this->n())
844 if (grad_grads.size() == this->n())
846 if (third_derivatives.size() == this->n())
848 if (fourth_derivatives.size() == this->n())
853 std::size_t max_n_polynomials = 0;
854 for (
unsigned int d = 0; d < dim; ++d)
855 max_n_polynomials =
std::max(max_n_polynomials, polynomials[d].size());
858 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
860 if (n_derivatives == 0)
861 for (
unsigned int d = 0; d < dim; ++d)
862 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
863 values_1d[i][0][d] = polynomials[d][i].value(p[d]);
865 for (
unsigned int d = 0; d < dim; ++d)
866 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
871 std::array<double, 5> derivatives;
872 polynomials[d][i].value(p[d], n_derivatives, derivatives.data());
873 for (
unsigned int j = 0; j <= n_derivatives; ++j)
874 values_1d[i][j][d] = derivatives[j];
878 constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1;
879 boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
880 for (
unsigned int d = 1; d < dim; ++d)
882 const unsigned int size = indices.size();
883 for (
unsigned int i = 1; i < polynomials[d].size(); ++i)
884 for (
unsigned int j = 0; j < size; ++j)
886 std::array<unsigned int, dim1> next_index = indices[j];
887 next_index[d - 1] = i;
888 indices.push_back(next_index);
892 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
895 polynomials[0].size(),
910 std::vector<double> &,
927 for (
unsigned int d = 0; d < dim; ++d)
947std::unique_ptr<ScalarPolynomialsBase<dim>>
950 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_NAMESPACE_CLOSE
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.
#define DEAL_II_NOT_IMPLEMENTED()
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