22 #include <boost/container/small_vector.hpp>
38 template <std::
size_t dim>
40 compute_tensor_index(
const unsigned int,
43 std::array<unsigned int, dim> &)
49 compute_tensor_index(
const unsigned int n,
52 std::array<unsigned int, 1> &indices)
58 compute_tensor_index(
const unsigned int n,
59 const unsigned int n_pols_0,
61 std::array<unsigned int, 2> &indices)
63 indices[0] = n % n_pols_0;
64 indices[1] = n / n_pols_0;
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)
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);
82 template <
int dim,
typename PolynomialType>
86 std::array<unsigned int, dim> &indices)
const
88 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
90 internal::compute_tensor_index(index_map[i],
102 std::array<unsigned int, 0> &)
const
109 template <
int dim,
typename PolynomialType>
112 std::ostream &out)
const
114 std::array<unsigned int, dim> ix;
115 for (
unsigned int i = 0; i < this->n(); ++i)
117 compute_index(i, ix);
119 for (
unsigned int d = 0;
d < dim; ++
d)
130 std::ostream &)
const
137 template <
int dim,
typename PolynomialType>
140 const std::vector<unsigned int> &renumber)
142 Assert(renumber.size() == index_map.size(),
145 index_map = renumber;
146 for (
unsigned int i = 0; i < index_map.size(); ++i)
147 index_map_inverse[index_map[i]] = i;
155 const std::vector<unsigned int> &)
162 template <
int dim,
typename PolynomialType>
165 const unsigned int i,
170 std::array<unsigned int, dim> indices;
171 compute_index(i, indices);
174 for (
unsigned int d = 0;
d < dim; ++
d)
175 value *= polynomials[indices[
d]].value(p(
d));
195 template <
int dim,
typename PolynomialType>
198 const unsigned int i,
201 std::array<unsigned int, dim> indices;
202 compute_index(i, indices);
210 std::vector<double> tmp(2);
211 for (
unsigned int d = 0;
d < dim; ++
d)
213 polynomials[indices[
d]].value(p(
d), tmp);
220 for (
unsigned int d = 0;
d < dim; ++
d)
223 for (
unsigned int x = 0; x < dim; ++x)
224 grad[
d] *= v[x][
d == x];
245 template <
int dim,
typename PolynomialType>
248 const unsigned int i,
251 std::array<unsigned int, dim> indices;
252 compute_index(i, indices);
256 std::vector<double> tmp(3);
257 for (
unsigned int d = 0;
d < dim; ++
d)
259 polynomials[indices[
d]].value(p(
d), tmp);
267 for (
unsigned int d1 = 0; d1 < dim; ++d1)
268 for (
unsigned int d2 = 0; d2 < dim; ++d2)
270 grad_grad[d1][d2] = 1.;
271 for (
unsigned int x = 0; x < dim; ++x)
273 unsigned int derivative = 0;
274 if (d1 == x || d2 == x)
281 grad_grad[d1][d2] *= v[x][derivative];
311 template <
int dim, std::
size_t dim1>
314 const unsigned int n_derivatives,
317 const unsigned int size_x,
318 const boost::container::small_vector<std::array<unsigned int, dim1>, 64>
320 const std::vector<unsigned int> &index_map,
321 std::vector<double> &
values,
328 const bool update_grads = (grads.size() == indices.size() * size_x);
329 const bool update_grad_grads =
330 (grad_grads.size() == indices.size() * size_x);
332 (third_derivatives.size() == indices.size() * size_x);
333 const bool update_4th_derivatives =
334 (fourth_derivatives.size() == indices.size() * size_x);
338 if (n_derivatives == 0)
339 for (
unsigned int i = 0, i1 = 0; i1 < indices.size(); ++i1)
341 double value_outer = 1.;
342 if constexpr (dim > 1)
343 for (
unsigned int d = 1;
d < dim; ++
d)
344 value_outer *= values_1d[indices[i1][
d - 1]][0][
d];
345 if (index_map.empty())
346 for (
unsigned int ix = 0; ix < size_x; ++ix, ++i)
347 values[i] = value_outer * values_1d[ix][0][0];
349 for (
unsigned int ix = 0; ix < size_x; ++ix, ++i)
350 values[index_map[i]] = value_outer * values_1d[ix][0][0];
353 for (
unsigned int iy = 0, i1 = 0; i1 < indices.size(); ++i1)
356 std::array<
double, dim + (dim * (dim - 1)) / 2> value_outer;
358 if constexpr (dim > 1)
360 for (
unsigned int x = 1; x < dim; ++x)
361 value_outer[0] *= values_1d[indices[i1][x - 1]][0][x];
362 for (
unsigned int d = 1;
d < dim; ++
d)
364 value_outer[
d] = values_1d[indices[i1][
d - 1]][1][
d];
365 for (
unsigned int x = 1; x < dim; ++x)
367 value_outer[
d] *= values_1d[indices[i1][x - 1]][0][x];
369 for (
unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
370 for (
unsigned int d2 = d1; d2 < dim; ++d2, ++count)
372 value_outer[count] = 1.;
373 for (
unsigned int x = 1; x < dim; ++x)
375 unsigned int derivative = 0;
381 value_outer[count] *=
382 values_1d[indices[i1][x - 1]][derivative][x];
389 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
391 std::array<double, 3> val_x{{values_1d[ix][0][0],
393 values_1d[ix][2][0]}};
394 const unsigned int index =
395 (index_map.empty() ? i : index_map[i]);
402 grads[
index][0] = value_outer[0] * val_x[1];
403 if constexpr (dim > 1)
404 for (
unsigned int d = 1;
d < dim; ++
d)
405 grads[
index][
d] = value_outer[
d] * val_x[0];
408 if (update_grad_grads)
410 grad_grads[
index][0][0] = value_outer[0] * val_x[2];
411 if constexpr (dim > 1)
413 for (
unsigned int d = 1;
d < dim; ++
d)
415 value_outer[
d] * val_x[1];
416 for (
unsigned int d1 = 1, count = dim; d1 < dim; ++d1)
417 for (
unsigned int d2 = d1; d2 < dim; ++d2, ++count)
418 grad_grads[
index][d1][d2] =
419 grad_grads[
index][d2][d1] =
420 value_outer[count] * val_x[0];
427 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
429 const unsigned int index =
430 (index_map.empty() ? i : index_map[i]);
431 std::array<unsigned int, dim> my_indices;
433 if constexpr (dim > 1)
434 for (
unsigned int d = 1;
d < dim; ++
d)
435 my_indices[
d] = indices[i1][
d - 1];
436 for (
unsigned int d1 = 0; d1 < dim; ++d1)
437 for (
unsigned int d2 = 0; d2 < dim; ++d2)
438 for (
unsigned int d3 = 0; d3 < dim; ++d3)
441 for (
unsigned int x = 0; x < dim; ++x)
443 unsigned int derivative = 0;
451 der3 *= values_1d[my_indices[x]][derivative][x];
453 third_derivatives[
index][d1][d2][d3] = der3;
457 if (update_4th_derivatives)
458 for (
unsigned int ix = 0, i = iy; ix < size_x; ++ix, ++i)
460 const unsigned int index =
461 (index_map.empty() ? i : index_map[i]);
462 std::array<unsigned int, dim> my_indices;
464 if constexpr (dim > 1)
465 for (
unsigned int d = 1;
d < dim; ++
d)
466 my_indices[
d] = indices[i1][
d - 1];
467 for (
unsigned int d1 = 0; d1 < dim; ++d1)
468 for (
unsigned int d2 = 0; d2 < dim; ++d2)
469 for (
unsigned int d3 = 0; d3 < dim; ++d3)
470 for (
unsigned int d4 = 0; d4 < dim; ++d4)
473 for (
unsigned int x = 0; x < dim; ++x)
475 unsigned int derivative = 0;
485 der4 *= values_1d[my_indices[x]][derivative][x];
487 fourth_derivatives[
index][d1][d2][d3][d4] = der4;
499 template <
int dim,
typename PolynomialType>
503 std::vector<double> &
values,
512 Assert(grads.size() == this->n() || grads.empty(),
514 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
516 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
518 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
522 unsigned int n_derivatives = 0;
523 if (
values.size() == this->n())
525 if (grads.size() == this->n())
527 if (grad_grads.size() == this->n())
529 if (third_derivatives.size() == this->n())
531 if (fourth_derivatives.size() == this->n())
537 const unsigned int n_polynomials = polynomials.size();
538 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
540 if constexpr (std::is_same<PolynomialType,
543 std::array<double, dim> point_array;
544 for (
unsigned int d = 0;
d < dim; ++
d)
545 point_array[
d] = p[
d];
546 for (
unsigned int i = 0; i < n_polynomials; ++i)
547 polynomials[i].values_of_array(point_array,
549 values_1d[i].data());
552 for (
unsigned int i = 0; i < n_polynomials; ++i)
553 for (
unsigned int d = 0;
d < dim; ++
d)
555 std::array<double, 5> derivatives;
556 polynomials[i].value(p[
d], n_derivatives, derivatives.data());
557 for (
unsigned int j = 0; j <= n_derivatives; ++j)
558 values_1d[i][j][
d] = derivatives[j];
563 constexpr
unsigned int dim1 = dim > 1 ? dim - 1 : 1;
564 boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
565 if constexpr (dim > 1)
566 for (
unsigned int d = 1;
d < dim; ++
d)
568 const unsigned int size = indices.size();
569 for (
unsigned int i = 1; i < n_polynomials; ++i)
570 for (
unsigned int j = 0; j < size; ++j)
572 std::array<unsigned int, dim1> next_index = indices[j];
573 next_index[
d - 1] = i;
574 indices.push_back(next_index);
579 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
598 std::vector<double> &,
609 template <
int dim,
typename PolynomialType>
610 std::unique_ptr<ScalarPolynomialsBase<dim>>
613 return std::make_unique<TensorProductPolynomials<dim, PolynomialType>>(*this);
618 template <
int dim,
typename PolynomialType>
629 template <
int dim,
typename PolynomialType>
630 std::vector<PolynomialType>
649 for (
const auto &pols_d : pols)
653 ExcMessage(
"The number of polynomials must be larger than zero "
654 "for all coordinate directions."));
663 const unsigned int i,
664 std::array<unsigned int, dim> &indices)
const
667 unsigned int n_poly = 1;
668 for (
unsigned int d = 0;
d < dim; ++
d)
669 n_poly *= polynomials[
d].size();
677 internal::compute_tensor_index(i,
678 polynomials[0].size(),
682 internal::compute_tensor_index(i,
683 polynomials[0].size(),
684 polynomials[1].size(),
693 std::array<unsigned int, 0> &)
const
705 std::array<unsigned int, dim> indices;
706 compute_index(i, indices);
709 for (
unsigned int d = 0;
d < dim; ++
d)
710 value *= polynomials[
d][indices[
d]].value(p(
d));
734 std::array<unsigned int, dim> indices;
735 compute_index(i, indices);
742 for (
unsigned int d = 0;
d < dim; ++
d)
743 polynomials[
d][indices[
d]].value(p(
d), 1, v[
d].data());
746 for (
unsigned int d = 0;
d < dim; ++
d)
749 for (
unsigned int x = 0; x < dim; ++x)
750 grad[
d] *= v[x][
d == x];
775 std::array<unsigned int, dim> indices;
776 compute_index(i, indices);
779 for (
unsigned int d = 0;
d < dim; ++
d)
780 polynomials[
d][indices[
d]].value(p(
d), 2, v[
d].data());
783 for (
unsigned int d1 = 0; d1 < dim; ++d1)
784 for (
unsigned int d2 = 0; d2 < dim; ++d2)
786 grad_grad[d1][d2] = 1.;
787 for (
unsigned int x = 0; x < dim; ++x)
789 unsigned int derivative = 0;
790 if (d1 == x || d2 == x)
797 grad_grad[d1][d2] *= v[x][derivative];
822 std::vector<double> &
values,
830 Assert(grads.size() == this->n() || grads.empty(),
832 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
834 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
836 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
840 unsigned int n_derivatives = 0;
841 if (
values.size() == this->n())
843 if (grads.size() == this->n())
845 if (grad_grads.size() == this->n())
847 if (third_derivatives.size() == this->n())
849 if (fourth_derivatives.size() == this->n())
854 std::size_t max_n_polynomials = 0;
855 for (
unsigned int d = 0;
d < dim; ++
d)
856 max_n_polynomials =
std::max(max_n_polynomials, polynomials[
d].size());
859 boost::container::small_vector<ndarray<double, 5, dim>, 10> values_1d(
861 if (n_derivatives == 0)
862 for (
unsigned int d = 0;
d < dim; ++
d)
863 for (
unsigned int i = 0; i < polynomials[
d].size(); ++i)
864 values_1d[i][0][
d] = polynomials[
d][i].value(p[
d]);
866 for (
unsigned int d = 0;
d < dim; ++
d)
867 for (
unsigned int i = 0; i < polynomials[
d].size(); ++i)
872 std::array<double, 5> derivatives;
873 polynomials[
d][i].value(p[
d], n_derivatives, derivatives.data());
874 for (
unsigned int j = 0; j <= n_derivatives; ++j)
875 values_1d[i][j][
d] = derivatives[j];
879 constexpr
unsigned int dim1 = dim > 1 ? dim - 1 : 1;
880 boost::container::small_vector<std::array<unsigned int, dim1>, 64> indices(1);
881 for (
unsigned int d = 1;
d < dim; ++
d)
883 const unsigned int size = indices.size();
884 for (
unsigned int i = 1; i < polynomials[
d].size(); ++i)
885 for (
unsigned int j = 0; j < size; ++j)
887 std::array<unsigned int, dim1> next_index = indices[j];
888 next_index[
d - 1] = i;
889 indices.push_back(next_index);
893 internal::TensorProductPolynomials::evaluate_tensor_product<dim>(
896 polynomials[0].size(),
911 std::vector<double> &,
928 for (
unsigned int d = 0;
d < dim; ++
d)
948 std::unique_ptr<ScalarPolynomialsBase<dim>>
951 return std::make_unique<AnisotropicPolynomials<dim>>(*this);
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
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 & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
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)
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
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)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray