From: Martin Kronbichler Date: Tue, 22 Aug 2023 16:23:56 +0000 (+0200) Subject: Use common code for tensor product polynomial evaluation X-Git-Tag: relicensing~547^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f39d9dd455d606e632a0788a94ab4228e7610257;p=dealii.git Use common code for tensor product polynomial evaluation --- diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc index b6b913f2f2..1f50ad2506 100644 --- a/source/base/tensor_product_polynomials.cc +++ b/source/base/tensor_product_polynomials.cc @@ -303,104 +303,72 @@ TensorProductPolynomials<0, Polynomials::Polynomial>::compute_grad_grad( -template -void -TensorProductPolynomials::evaluate( - const Point &p, - std::vector &values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const +namespace internal { - Assert(dim <= 3, ExcNotImplemented()); - Assert(values.size() == this->n() || values.empty(), - ExcDimensionMismatch2(values.size(), this->n(), 0)); - Assert(grads.size() == this->n() || grads.empty(), - ExcDimensionMismatch2(grads.size(), this->n(), 0)); - Assert(grad_grads.size() == this->n() || grad_grads.empty(), - ExcDimensionMismatch2(grad_grads.size(), this->n(), 0)); - Assert(third_derivatives.size() == this->n() || third_derivatives.empty(), - ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0)); - Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(), - ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0)); - - const bool update_values = (values.size() == this->n()), - update_grads = (grads.size() == this->n()), - update_grad_grads = (grad_grads.size() == this->n()), - update_3rd_derivatives = (third_derivatives.size() == this->n()), - update_4th_derivatives = (fourth_derivatives.size() == this->n()); - - // check how many values/derivatives we have to compute - unsigned int n_values_and_derivatives = 0; - if (update_values) - n_values_and_derivatives = 1; - if (update_grads) - n_values_and_derivatives = 2; - if (update_grad_grads) - n_values_and_derivatives = 3; - if (update_3rd_derivatives) - n_values_and_derivatives = 4; - if (update_4th_derivatives) - n_values_and_derivatives = 5; + namespace TensorProductPolynomials + { + // This function computes the tensor product of some tabulated + // one-dimensional polynomials (also the anisotropic case is supported) + // with tensor product indices as tabulated in the 'indices' array + template + void + evaluate_tensor_product( + const unsigned int n_values_and_derivatives, + const boost::container::small_vector, 10> + &values_1d, + const boost::container::small_vector, 125> + & indices, + const std::vector &index_map, + std::vector & values, + std::vector> & grads, + std::vector> & grad_grads, + std::vector> & third_derivatives, + std::vector> & fourth_derivatives) + { + const bool update_values = (values.size() == indices.size()), + update_grads = (grads.size() == indices.size()), + update_grad_grads = (grad_grads.size() == indices.size()), + update_3rd_derivatives = + (third_derivatives.size() == indices.size()), + update_4th_derivatives = + (fourth_derivatives.size() == indices.size()); - // Compute the values (and derivatives, if necessary) of all 1d polynomials - // at this evaluation point. We need to compute dim*n_polynomials - // evaluations, involving an evaluation of each polynomial for each - // coordinate direction. Once we have those values, we perform the - // multiplications for the tensor product in the arbitrary dimension. - const unsigned int n_polynomials = polynomials.size(); - boost::container::small_vector, 20> values_1d( - n_polynomials); - if (n_values_and_derivatives == 1) - for (unsigned int i = 0; i < n_polynomials; ++i) - for (unsigned int d = 0; d < dim; ++d) - values_1d[i][d][0] = polynomials[i].value(p(d)); - else - for (unsigned int i = 0; i < n_polynomials; ++i) - for (unsigned d = 0; d < dim; ++d) - polynomials[i].value(p(d), - n_values_and_derivatives, - values_1d[i][d].data()); - - unsigned int indices[3]; - unsigned int ind = 0; - for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2]) - for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1); - ++indices[1]) if (n_values_and_derivatives == 1) - for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind) + for (unsigned int i = 0; i < indices.size(); ++i) { - double value = values_1d[indices[0]][0][0]; + const std::array my_indices = + indices[index_map.empty() ? i : index_map[i]]; + double value = values_1d[my_indices[0]][0][0]; for (unsigned int d = 1; d < dim; ++d) - value *= values_1d[indices[d]][d][0]; - values[index_map_inverse[ind]] = value; + value *= values_1d[my_indices[d]][0][d]; + values[i] = value; } else - for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind) + for (unsigned int i = 0; i < indices.size(); ++i) { - const unsigned int i = index_map_inverse[ind]; - + const std::array my_indices = + indices[index_map.empty() ? i : index_map[i]]; if (update_values) { - double value = values_1d[indices[0]][0][0]; + double value = values_1d[my_indices[0]][0][0]; for (unsigned int x = 1; x < dim; ++x) - value *= values_1d[indices[x]][x][0]; + value *= values_1d[my_indices[x]][0][x]; values[i] = value; } if (update_grads) for (unsigned int d = 0; d < dim; ++d) { - double grad = 1.; + double grad = values_1d[my_indices[d]][1][d]; for (unsigned int x = 0; x < dim; ++x) - grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0]; + if (x != d) + grad *= values_1d[my_indices[x]][0][x]; grads[i][d] = grad; } if (update_grad_grads) for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) + for (unsigned int d2 = d1; d2 < dim; ++d2) { double der2 = 1.; for (unsigned int x = 0; x < dim; ++x) @@ -411,9 +379,10 @@ TensorProductPolynomials::evaluate( if (d2 == x) ++derivative; - der2 *= values_1d[indices[x]][x][derivative]; + der2 *= values_1d[my_indices[x]][derivative][x]; } grad_grads[i][d1][d2] = der2; + grad_grads[i][d2][d1] = der2; } if (update_3rd_derivatives) @@ -432,7 +401,7 @@ TensorProductPolynomials::evaluate( if (d3 == x) ++derivative; - der3 *= values_1d[indices[x]][x][derivative]; + der3 *= values_1d[my_indices[x]][derivative][x]; } third_derivatives[i][d1][d2][d3] = der3; } @@ -456,11 +425,106 @@ TensorProductPolynomials::evaluate( if (d4 == x) ++derivative; - der4 *= values_1d[indices[x]][x][derivative]; + der4 *= values_1d[my_indices[x]][derivative][x]; } fourth_derivatives[i][d1][d2][d3][d4] = der4; } } + } + } // namespace TensorProductPolynomials +} // namespace internal + + +template +void +TensorProductPolynomials::evaluate( + const Point & p, + std::vector & values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const +{ + Assert(dim <= 3, ExcNotImplemented()); + Assert(values.size() == this->n() || values.empty(), + ExcDimensionMismatch2(values.size(), this->n(), 0)); + Assert(grads.size() == this->n() || grads.empty(), + ExcDimensionMismatch2(grads.size(), this->n(), 0)); + Assert(grad_grads.size() == this->n() || grad_grads.empty(), + ExcDimensionMismatch2(grad_grads.size(), this->n(), 0)); + Assert(third_derivatives.size() == this->n() || third_derivatives.empty(), + ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0)); + Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(), + ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0)); + + // check how many values/derivatives we have to compute + unsigned int n_values_and_derivatives = 0; + if (values.size() == this->n()) + n_values_and_derivatives = 1; + if (grads.size() == this->n()) + n_values_and_derivatives = 2; + if (grad_grads.size() == this->n()) + n_values_and_derivatives = 3; + if (third_derivatives.size() == this->n()) + n_values_and_derivatives = 4; + if (fourth_derivatives.size() == this->n()) + n_values_and_derivatives = 5; + + // Compute the values (and derivatives, if necessary) of all 1d polynomials + // at this evaluation point. We can use the more optimized values_of_array + // function to compute 'dim' polynomials at once + const unsigned int n_polynomials = polynomials.size(); + boost::container::small_vector, 10> values_1d( + n_polynomials); + if constexpr (std::is_same>::value) + { + std::array point_array; + for (unsigned int d = 0; d < dim; ++d) + point_array[d] = p[d]; + for (unsigned int i = 0; i < n_polynomials; ++i) + polynomials[i].values_of_array(point_array, + n_values_and_derivatives, + values_1d[i].data()); + } + else + for (unsigned int i = 0; i < n_polynomials; ++i) + for (unsigned int d = 0; d < dim; ++d) + { + std::array derivatives; + polynomials[i].value(p[d], + n_values_and_derivatives, + derivatives.data()); + for (unsigned int j = 0; j < n_values_and_derivatives; ++j) + values_1d[i][j][d] = derivatives[j]; + } + + // Unroll the tensor product indices in arbitrary dimension + boost::container::small_vector, 125> indices(1); + indices.reserve(Utilities::pow(n_polynomials, dim)); + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int size = indices.size(); + for (unsigned int i = 1; i < n_polynomials; ++i) + for (unsigned int j = 0; j < size; ++j) + { + std::array next_index = indices[j]; + next_index[d] = i; + indices.push_back(next_index); + } + } + AssertDimension(indices.size(), Utilities::pow(n_polynomials, dim)); + + internal::TensorProductPolynomials::evaluate_tensor_product( + n_values_and_derivatives, + values_1d, + indices, + index_map, + values, + grads, + grad_grads, + third_derivatives, + fourth_derivatives); } @@ -715,129 +779,73 @@ AnisotropicPolynomials::evaluate( Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(), ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0)); - const bool update_values = (values.size() == this->n()), - update_grads = (grads.size() == this->n()), - update_grad_grads = (grad_grads.size() == this->n()), - update_3rd_derivatives = (third_derivatives.size() == this->n()), - update_4th_derivatives = (fourth_derivatives.size() == this->n()); - - // check how many - // values/derivatives we have to - // compute + // check how many values/derivatives we have to compute unsigned int n_values_and_derivatives = 0; - if (update_values) + if (values.size() == this->n()) n_values_and_derivatives = 1; - if (update_grads) + if (grads.size() == this->n()) n_values_and_derivatives = 2; - if (update_grad_grads) + if (grad_grads.size() == this->n()) n_values_and_derivatives = 3; - if (update_3rd_derivatives) + if (third_derivatives.size() == this->n()) n_values_and_derivatives = 4; - if (update_4th_derivatives) + if (fourth_derivatives.size() == this->n()) n_values_and_derivatives = 5; - // compute the values (and - // derivatives, if necessary) of - // all polynomials at this - // evaluation point + // compute the values (and derivatives, if necessary) of all polynomials at + // this evaluation point std::size_t max_n_polynomials = 0; for (unsigned int d = 0; d < dim; ++d) max_n_polynomials = std::max(max_n_polynomials, polynomials[d].size()); // 5 is enough to store values and derivatives in all supported cases - Table<2, std::array> v(dim, max_n_polynomials); - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int i = 0; i < polynomials[d].size(); ++i) - polynomials[d][i].value(p(d), - n_values_and_derivatives - 1, - v(d, i).data()); - - for (unsigned int i = 0; i < this->n(); ++i) - { - // first get the - // one-dimensional indices of - // this particular tensor - // product polynomial - std::array indices; - compute_index(i, indices); - - if (update_values) + boost::container::small_vector, 10> values_1d( + max_n_polynomials); + if (n_values_and_derivatives == 1) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int i = 0; i < polynomials[d].size(); ++i) + values_1d[i][0][d] = polynomials[d][i].value(p[d]); + else + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int i = 0; i < polynomials[d].size(); ++i) { - values[i] = 1; - for (unsigned int x = 0; x < dim; ++x) - values[i] *= v(x, indices[x])[0]; + // The isotropic tensor product function wants us to use a different + // innermost index, so we cannot pass the values_1d array into the + // function directly + std::array derivatives; + polynomials[d][i].value(p[d], + n_values_and_derivatives, + derivatives.data()); + for (unsigned int j = 0; j < n_values_and_derivatives; ++j) + values_1d[i][j][d] = derivatives[j]; } - if (update_grads) - for (unsigned int d = 0; d < dim; ++d) + // Unroll the tensor product indices in arbitrary dimension + boost::container::small_vector, 125> indices(1); + indices.reserve(this->n()); + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int size = indices.size(); + for (unsigned int i = 1; i < polynomials[d].size(); ++i) + for (unsigned int j = 0; j < size; ++j) { - grads[i][d] = 1.; - for (unsigned int x = 0; x < dim; ++x) - grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0]; + std::array next_index = indices[j]; + next_index[d] = i; + indices.push_back(next_index); } - - if (update_grad_grads) - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - { - grad_grads[i][d1][d2] = 1.; - for (unsigned int x = 0; x < dim; ++x) - { - unsigned int derivative = 0; - if (d1 == x) - ++derivative; - if (d2 == x) - ++derivative; - - grad_grads[i][d1][d2] *= v(x, indices[x])[derivative]; - } - } - - if (update_3rd_derivatives) - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - for (unsigned int d3 = 0; d3 < dim; ++d3) - { - third_derivatives[i][d1][d2][d3] = 1.; - for (unsigned int x = 0; x < dim; ++x) - { - unsigned int derivative = 0; - if (d1 == x) - ++derivative; - if (d2 == x) - ++derivative; - if (d3 == x) - ++derivative; - - third_derivatives[i][d1][d2][d3] *= - v(x, indices[x])[derivative]; - } - } - - if (update_4th_derivatives) - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - for (unsigned int d3 = 0; d3 < dim; ++d3) - for (unsigned int d4 = 0; d4 < dim; ++d4) - { - fourth_derivatives[i][d1][d2][d3][d4] = 1.; - for (unsigned int x = 0; x < dim; ++x) - { - unsigned int derivative = 0; - if (d1 == x) - ++derivative; - if (d2 == x) - ++derivative; - if (d3 == x) - ++derivative; - if (d4 == x) - ++derivative; - - fourth_derivatives[i][d1][d2][d3][d4] *= - v(x, indices[x])[derivative]; - } - } } + AssertDimension(indices.size(), this->n()); + + internal::TensorProductPolynomials::evaluate_tensor_product( + n_values_and_derivatives, + values_1d, + indices, + {}, + values, + grads, + grad_grads, + third_derivatives, + fourth_derivatives); }