From 23c29c83152b6c10c05fe9da2830974c406e8688 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 24 Feb 2025 14:00:11 -0700 Subject: [PATCH] Avoid explicit specializations. --- source/base/tensor_product_polynomials.cc | 791 ++++++++---------- .../tensor_product_polynomials_bubbles.cc | 83 +- .../base/tensor_product_polynomials_const.cc | 34 +- 3 files changed, 426 insertions(+), 482 deletions(-) diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc index cd809a5c6f..4116f11bc5 100644 --- a/source/base/tensor_product_polynomials.cc +++ b/source/base/tensor_product_polynomials.cc @@ -84,23 +84,19 @@ TensorProductPolynomials::compute_index( const unsigned int i, std::array &indices) const { - Assert(i < Utilities::fixed_power(polynomials.size()), - ExcInternalError()); - internal::compute_tensor_index(index_map[i], - polynomials.size(), - polynomials.size(), - indices); -} - - - -template <> -inline void -TensorProductPolynomials<0, Polynomials::Polynomial>::compute_index( - const unsigned int, - std::array &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + } + else + { + Assert(i < Utilities::fixed_power(polynomials.size()), + ExcInternalError()); + internal::compute_tensor_index(index_map[i], + polynomials.size(), + polynomials.size(), + indices); + } } @@ -110,25 +106,22 @@ void TensorProductPolynomials::output_indices( std::ostream &out) const { - std::array ix; - for (unsigned int i = 0; i < this->n(); ++i) + if constexpr (dim == 0) { - compute_index(i, ix); - out << i << "\t"; - for (unsigned int d = 0; d < dim; ++d) - out << ix[d] << " "; - out << std::endl; + DEAL_II_NOT_IMPLEMENTED(); + } + else + { + std::array ix; + for (unsigned int i = 0; i < this->n(); ++i) + { + compute_index(i, ix); + out << i << "\t"; + for (unsigned int d = 0; d < dim; ++d) + out << ix[d] << " "; + out << std::endl; + } } -} - - - -template <> -void -TensorProductPolynomials<0, Polynomials::Polynomial>::output_indices( - std::ostream &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); } @@ -206,29 +199,22 @@ TensorProductPolynomials::compute_value( const unsigned int i, const Point &p) const { - Assert(dim > 0, ExcNotImplemented()); - - std::array indices; - compute_index(i, indices); - - double value = 1.; - for (unsigned int d = 0; d < dim; ++d) - value *= polynomials[indices[d]].value(p[d]); - - return value; -} - - + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return 0; + } + else + { + std::array indices; + compute_index(i, indices); -template <> -double -TensorProductPolynomials<0, Polynomials::Polynomial>::compute_value( - const unsigned int, - const Point<0> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + double value = 1.; + for (unsigned int d = 0; d < dim; ++d) + value *= polynomials[indices[d]].value(p[d]); - return {}; + return value; + } } @@ -239,46 +225,41 @@ TensorProductPolynomials::compute_grad( const unsigned int i, const Point &p) const { - std::array indices; - compute_index(i, indices); - - // compute values and - // uni-directional derivatives at - // the given point in each - // coordinate direction - ndarray v; - { - std::vector tmp(2); - for (unsigned int d = 0; d < dim; ++d) - { - polynomials[indices[d]].value(p[d], tmp); - v[d][0] = tmp[0]; - v[d][1] = tmp[1]; - } - } - - Tensor<1, dim> grad; - for (unsigned int d = 0; d < dim; ++d) + if constexpr (dim == 0) { - grad[d] = 1.; - for (unsigned int x = 0; x < dim; ++x) - grad[d] *= v[x][d == x]; + DEAL_II_NOT_IMPLEMENTED(); + return {}; } + else + { + std::array indices; + compute_index(i, indices); + + // compute values and + // uni-directional derivatives at + // the given point in each + // coordinate direction + ndarray v; + { + std::vector tmp(2); + for (unsigned int d = 0; d < dim; ++d) + { + polynomials[indices[d]].value(p[d], tmp); + v[d][0] = tmp[0]; + v[d][1] = tmp[1]; + } + } - return grad; -} - - - -template <> -Tensor<1, 0> -TensorProductPolynomials<0, Polynomials::Polynomial>::compute_grad( - const unsigned int, - const Point<0> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + Tensor<1, dim> grad; + for (unsigned int d = 0; d < dim; ++d) + { + grad[d] = 1.; + for (unsigned int x = 0; x < dim; ++x) + grad[d] *= v[x][d == x]; + } - return {}; + return grad; + } } @@ -289,52 +270,49 @@ TensorProductPolynomials::compute_grad_grad( const unsigned int i, const Point &p) const { - std::array indices; - compute_index(i, indices); + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return {}; + } + else + { + std::array indices; + compute_index(i, indices); - ndarray v; - { - std::vector tmp(3); - for (unsigned int d = 0; d < dim; ++d) + ndarray v; { - polynomials[indices[d]].value(p[d], tmp); - v[d][0] = tmp[0]; - v[d][1] = tmp[1]; - v[d][2] = tmp[2]; + std::vector tmp(3); + for (unsigned int d = 0; d < dim; ++d) + { + polynomials[indices[d]].value(p[d], tmp); + v[d][0] = tmp[0]; + v[d][1] = tmp[1]; + v[d][2] = tmp[2]; + } } - } - Tensor<2, dim> grad_grad; - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - { - grad_grad[d1][d2] = 1.; - for (unsigned int x = 0; x < dim; ++x) + Tensor<2, dim> grad_grad; + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) { - unsigned int derivative = 0; - if (d1 == x || d2 == x) + grad_grad[d1][d2] = 1.; + for (unsigned int x = 0; x < dim; ++x) { - if (d1 == d2) - derivative = 2; - else - derivative = 1; + unsigned int derivative = 0; + if (d1 == x || d2 == x) + { + if (d1 == d2) + derivative = 2; + else + derivative = 1; + } + grad_grad[d1][d2] *= v[x][derivative]; } - grad_grad[d1][d2] *= v[x][derivative]; } - } - - return grad_grad; -} - - -template <> -Tensor<2, 0> -TensorProductPolynomials<0, Polynomials::Polynomial>::compute_grad_grad( - const unsigned int, - const Point<0> &) const -{ - return {}; + return grad_grad; + } } @@ -547,102 +525,97 @@ TensorProductPolynomials::evaluate( 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_derivatives = 0; - if (values.size() == this->n()) - n_derivatives = 0; - if (grads.size() == this->n()) - n_derivatives = 1; - if (grad_grads.size() == this->n()) - n_derivatives = 2; - if (third_derivatives.size() == this->n()) - n_derivatives = 3; - if (fourth_derivatives.size() == this->n()) - n_derivatives = 4; - - // 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) + if constexpr (dim == 0) { - 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_derivatives, - values_1d[i].data()); + DEAL_II_NOT_IMPLEMENTED(); } else - for (unsigned int i = 0; i < n_polynomials; ++i) - for (unsigned int d = 0; d < dim; ++d) + { + 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_derivatives = 0; + if (values.size() == this->n()) + n_derivatives = 0; + if (grads.size() == this->n()) + n_derivatives = 1; + if (grad_grads.size() == this->n()) + n_derivatives = 2; + if (third_derivatives.size() == this->n()) + n_derivatives = 3; + if (fourth_derivatives.size() == this->n()) + n_derivatives = 4; + + // 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< + PolynomialType, + dealii::Polynomials::Polynomial>::value) { - std::array derivatives; - polynomials[i].value(p[d], n_derivatives, derivatives.data()); - for (unsigned int j = 0; j <= n_derivatives; ++j) - values_1d[i][j][d] = derivatives[j]; + 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_derivatives, + values_1d[i].data()); } - - // Unroll the tensor product indices of all but the first dimension in - // arbitrary dimension - constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1; - boost::container::small_vector, 64> indices(1); - if constexpr (dim > 1) - for (unsigned int d = 1; 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) + else + for (unsigned int i = 0; i < n_polynomials; ++i) + for (unsigned int d = 0; d < dim; ++d) { - std::array next_index = indices[j]; - next_index[d - 1] = i; - indices.push_back(next_index); + std::array derivatives; + polynomials[i].value(p[d], n_derivatives, derivatives.data()); + for (unsigned int j = 0; j <= n_derivatives; ++j) + values_1d[i][j][d] = derivatives[j]; } - } - AssertDimension(indices.size(), Utilities::pow(n_polynomials, dim - 1)); - - internal::TensorProductPolynomials::evaluate_tensor_product( - n_derivatives, - values_1d, - n_polynomials, - indices, - index_map_inverse, - values, - grads, - grad_grads, - third_derivatives, - fourth_derivatives); -} - - -template <> -void -TensorProductPolynomials<0, Polynomials::Polynomial>::evaluate( - const Point<0> &, - std::vector &, - std::vector> &, - std::vector> &, - std::vector> &, - std::vector> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + // Unroll the tensor product indices of all but the first dimension in + // arbitrary dimension + constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1; + boost::container::small_vector, 64> + indices(1); + if constexpr (dim > 1) + for (unsigned int d = 1; 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 - 1] = i; + indices.push_back(next_index); + } + } + AssertDimension(indices.size(), Utilities::pow(n_polynomials, dim - 1)); + + internal::TensorProductPolynomials::evaluate_tensor_product( + n_derivatives, + values_1d, + n_polynomials, + indices, + index_map_inverse, + values, + grads, + grad_grads, + third_derivatives, + fourth_derivatives); + } } @@ -714,36 +687,33 @@ AnisotropicPolynomials::compute_index( const unsigned int i, std::array &indices) const { -#ifdef DEBUG - unsigned int n_poly = 1; - for (unsigned int d = 0; d < dim; ++d) - n_poly *= polynomials[d].size(); - Assert(i < n_poly, ExcInternalError()); -#endif - - if (dim == 0) + if constexpr (dim == 0) { + DEAL_II_NOT_IMPLEMENTED(); } - else if (dim == 1) - internal::compute_tensor_index(index_map[i], - polynomials[0].size(), - 0 /*not used*/, - indices); else - internal::compute_tensor_index(index_map[i], - polynomials[0].size(), - polynomials[1].size(), - indices); -} - - + { +#ifdef DEBUG + unsigned int n_poly = 1; + for (unsigned int d = 0; d < dim; ++d) + n_poly *= polynomials[d].size(); + Assert(i < n_poly, ExcInternalError()); +#endif -template <> -void -AnisotropicPolynomials<0>::compute_index(const unsigned int, - std::array &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + if (dim == 0) + { + } + else if (dim == 1) + internal::compute_tensor_index(index_map[i], + polynomials[0].size(), + 0 /*not used*/, + indices); + else + internal::compute_tensor_index(index_map[i], + polynomials[0].size(), + polynomials[1].size(), + indices); + } } @@ -753,26 +723,22 @@ double AnisotropicPolynomials::compute_value(const unsigned int i, const Point &p) const { - std::array indices; - compute_index(i, indices); - - double value = 1.; - for (unsigned int d = 0; d < dim; ++d) - value *= polynomials[d][indices[d]].value(p[d]); - - return value; -} - - + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return {}; + } + else + { + std::array indices; + compute_index(i, indices); -template <> -double -AnisotropicPolynomials<0>::compute_value(const unsigned int, - const Point<0> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + double value = 1.; + for (unsigned int d = 0; d < dim; ++d) + value *= polynomials[d][indices[d]].value(p[d]); - return {}; + return value; + } } @@ -782,38 +748,34 @@ Tensor<1, dim> AnisotropicPolynomials::compute_grad(const unsigned int i, const Point &p) const { - std::array indices; - compute_index(i, indices); - - // compute values and - // uni-directional derivatives at - // the given point in each - // coordinate direction - ndarray v; - for (unsigned int d = 0; d < dim; ++d) - polynomials[d][indices[d]].value(p[d], 1, v[d].data()); - - Tensor<1, dim> grad; - for (unsigned int d = 0; d < dim; ++d) + if constexpr (dim == 0) { - grad[d] = 1.; - for (unsigned int x = 0; x < dim; ++x) - grad[d] *= v[x][d == x]; + DEAL_II_NOT_IMPLEMENTED(); + return {}; } + else + { + std::array indices; + compute_index(i, indices); + + // compute values and + // uni-directional derivatives at + // the given point in each + // coordinate direction + ndarray v; + for (unsigned int d = 0; d < dim; ++d) + polynomials[d][indices[d]].value(p[d], 1, v[d].data()); - return grad; -} - - - -template <> -Tensor<1, 0> -AnisotropicPolynomials<0>::compute_grad(const unsigned int, - const Point<0> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + Tensor<1, dim> grad; + for (unsigned int d = 0; d < dim; ++d) + { + grad[d] = 1.; + for (unsigned int x = 0; x < dim; ++x) + grad[d] *= v[x][d == x]; + } - return {}; + return grad; + } } @@ -823,45 +785,41 @@ Tensor<2, dim> AnisotropicPolynomials::compute_grad_grad(const unsigned int i, const Point &p) const { - std::array indices; - compute_index(i, indices); + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return {}; + } + else + { + std::array indices; + compute_index(i, indices); - ndarray v; - for (unsigned int d = 0; d < dim; ++d) - polynomials[d][indices[d]].value(p[d], 2, v[d].data()); + ndarray v; + for (unsigned int d = 0; d < dim; ++d) + polynomials[d][indices[d]].value(p[d], 2, v[d].data()); - Tensor<2, dim> grad_grad; - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - { - grad_grad[d1][d2] = 1.; - for (unsigned int x = 0; x < dim; ++x) + Tensor<2, dim> grad_grad; + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) { - unsigned int derivative = 0; - if (d1 == x || d2 == x) + grad_grad[d1][d2] = 1.; + for (unsigned int x = 0; x < dim; ++x) { - if (d1 == d2) - derivative = 2; - else - derivative = 1; + unsigned int derivative = 0; + if (d1 == x || d2 == x) + { + if (d1 == d2) + derivative = 2; + else + derivative = 1; + } + grad_grad[d1][d2] *= v[x][derivative]; } - grad_grad[d1][d2] *= v[x][derivative]; } - } - - return grad_grad; -} - - -template <> -Tensor<2, 0> -AnisotropicPolynomials<0>::compute_grad_grad(const unsigned int, - const Point<0> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); - - return {}; + return grad_grad; + } } @@ -876,96 +834,91 @@ AnisotropicPolynomials::evaluate( std::vector> &third_derivatives, std::vector> &fourth_derivatives) const { - 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_derivatives = 0; - if (values.size() == this->n()) - n_derivatives = 0; - if (grads.size() == this->n()) - n_derivatives = 1; - if (grad_grads.size() == this->n()) - n_derivatives = 2; - if (third_derivatives.size() == this->n()) - n_derivatives = 3; - if (fourth_derivatives.size() == this->n()) - n_derivatives = 4; - - // 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 - boost::container::small_vector, 10> values_1d( - max_n_polynomials); - if (n_derivatives == 0) - 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) - { - // 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_derivatives, derivatives.data()); - for (unsigned int j = 0; j <= n_derivatives; ++j) - values_1d[i][j][d] = derivatives[j]; - } - - // Unroll the tensor product indices in arbitrary dimension - constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1; - boost::container::small_vector, 64> indices(1); - for (unsigned int d = 1; d < dim; ++d) + if constexpr (dim == 0) { - const unsigned int size = indices.size(); - for (unsigned int i = 1; i < polynomials[d].size(); ++i) - for (unsigned int j = 0; j < size; ++j) - { - std::array next_index = indices[j]; - next_index[d - 1] = i; - indices.push_back(next_index); - } + DEAL_II_NOT_IMPLEMENTED(); } + else + { + 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_derivatives = 0; + if (values.size() == this->n()) + n_derivatives = 0; + if (grads.size() == this->n()) + n_derivatives = 1; + if (grad_grads.size() == this->n()) + n_derivatives = 2; + if (third_derivatives.size() == this->n()) + n_derivatives = 3; + if (fourth_derivatives.size() == this->n()) + n_derivatives = 4; + + // 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()); - internal::TensorProductPolynomials::evaluate_tensor_product( - n_derivatives, - values_1d, - polynomials[0].size(), - indices, - index_map_inverse, - values, - grads, - grad_grads, - third_derivatives, - fourth_derivatives); -} - + // 5 is enough to store values and derivatives in all supported cases + boost::container::small_vector, 10> values_1d( + max_n_polynomials); + if (n_derivatives == 0) + 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) + { + // 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_derivatives, derivatives.data()); + for (unsigned int j = 0; j <= n_derivatives; ++j) + values_1d[i][j][d] = derivatives[j]; + } + // Unroll the tensor product indices in arbitrary dimension + constexpr unsigned int dim1 = dim > 1 ? dim - 1 : 1; + boost::container::small_vector, 64> + indices(1); + for (unsigned int d = 1; 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) + { + std::array next_index = indices[j]; + next_index[d - 1] = i; + indices.push_back(next_index); + } + } -template <> -void -AnisotropicPolynomials<0>::evaluate(const Point<0> &, - std::vector &, - std::vector> &, - std::vector> &, - std::vector> &, - std::vector> &) const -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); + internal::TensorProductPolynomials::evaluate_tensor_product( + n_derivatives, + values_1d, + polynomials[0].size(), + indices, + index_map_inverse, + values, + grads, + grad_grads, + third_derivatives, + fourth_derivatives); + } } @@ -975,22 +928,18 @@ unsigned int AnisotropicPolynomials::get_n_tensor_pols( const std::vector>> &pols) { - unsigned int y = 1; - for (unsigned int d = 0; d < dim; ++d) - y *= pols[d].size(); - return y; -} - - - -template <> -unsigned int -AnisotropicPolynomials<0>::get_n_tensor_pols( - const std::vector>> &) -{ - AssertThrow(false, ExcNotImplemented("This function does not work in 0-d!")); - - return {}; + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return {}; + } + else + { + unsigned int y = 1; + for (unsigned int d = 0; d < dim; ++d) + y *= pols[d].size(); + return y; + } } diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc index 0f4503345e..fab074977e 100644 --- a/source/base/tensor_product_polynomials_bubbles.cc +++ b/source/base/tensor_product_polynomials_bubbles.cc @@ -93,59 +93,56 @@ TensorProductPolynomialsBubbles::compute_value(const unsigned int i, -template <> -double -TensorProductPolynomialsBubbles<0>::compute_value(const unsigned int, - const Point<0> &) const -{ - DEAL_II_NOT_IMPLEMENTED(); - return 0.; -} - - - template Tensor<1, dim> TensorProductPolynomialsBubbles::compute_grad(const unsigned int i, const Point &p) const { - const unsigned int q_degree = tensor_polys.polynomials.size() - 1; - const unsigned int max_q_indices = tensor_polys.n(); - Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim), - ExcInternalError()); + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return 0.; + } + else + { + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; + const unsigned int max_q_indices = tensor_polys.n(); + Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim), + ExcInternalError()); - // treat the regular basis functions - if (i < max_q_indices) - return tensor_polys.compute_grad(i, p); + // treat the regular basis functions + if (i < max_q_indices) + return tensor_polys.compute_grad(i, p); - const unsigned int comp = i - tensor_polys.n(); - Tensor<1, dim> grad; + const unsigned int comp = i - tensor_polys.n(); + Tensor<1, dim> grad; - for (unsigned int d = 0; d < dim; ++d) - { - grad[d] = 1.; - // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p) - for (unsigned j = 0; j < dim; ++j) - grad[d] *= (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j])); - // and multiply with (2*x_i-1)^{r-1} - for (unsigned int i = 0; i < q_degree - 1; ++i) - grad[d] *= 2 * p[comp] - 1; - } + for (unsigned int d = 0; d < dim; ++d) + { + grad[d] = 1.; + // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p) + for (unsigned j = 0; j < dim; ++j) + grad[d] *= (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j])); + // and multiply with (2*x_i-1)^{r-1} + for (unsigned int i = 0; i < q_degree - 1; ++i) + grad[d] *= 2 * p[comp] - 1; + } - if (q_degree >= 2) - { - // add \prod_{i=1}^d 4*(x_i(1-x_i))(p) - double value = 1.; - for (unsigned int j = 0; j < dim; ++j) - value *= 4 * p[j] * (1 - p[j]); - // and multiply with grad(2*x_i-1)^{r-1} - double tmp = value * 2 * (q_degree - 1); - for (unsigned int i = 0; i < q_degree - 2; ++i) - tmp *= 2 * p[comp] - 1; - grad[comp] += tmp; - } + if (q_degree >= 2) + { + // add \prod_{i=1}^d 4*(x_i(1-x_i))(p) + double value = 1.; + for (unsigned int j = 0; j < dim; ++j) + value *= 4 * p[j] * (1 - p[j]); + // and multiply with grad(2*x_i-1)^{r-1} + double tmp = value * 2 * (q_degree - 1); + for (unsigned int i = 0; i < q_degree - 2; ++i) + tmp *= 2 * p[comp] - 1; + grad[comp] += tmp; + } - return grad; + return grad; + } } diff --git a/source/base/tensor_product_polynomials_const.cc b/source/base/tensor_product_polynomials_const.cc index 667072c4be..0e5d32ad55 100644 --- a/source/base/tensor_product_polynomials_const.cc +++ b/source/base/tensor_product_polynomials_const.cc @@ -81,30 +81,28 @@ TensorProductPolynomialsConst::compute_value(const unsigned int i, -template <> -double -TensorProductPolynomialsConst<0>::compute_value(const unsigned int, - const Point<0> &) const -{ - DEAL_II_NOT_IMPLEMENTED(); - return 0.; -} - - template Tensor<1, dim> TensorProductPolynomialsConst::compute_grad(const unsigned int i, const Point &p) const { - const unsigned int max_indices = tensor_polys.n(); - Assert(i <= max_indices, ExcInternalError()); - - // treat the regular basis functions - if (i < max_indices) - return tensor_polys.compute_grad(i, p); + if constexpr (dim == 0) + { + DEAL_II_NOT_IMPLEMENTED(); + return 0.; + } else - // this is for the constant function - return Tensor<1, dim>(); + { + const unsigned int max_indices = tensor_polys.n(); + Assert(i <= max_indices, ExcInternalError()); + + // treat the regular basis functions + if (i < max_indices) + return tensor_polys.compute_grad(i, p); + else + // this is for the constant function + return Tensor<1, dim>(); + } } template -- 2.39.5