From d6bcb6880d149af2c89d72dbe4e90f9fefe392ee Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 13 Aug 2024 16:04:07 -0600 Subject: [PATCH] TensorProductPolynomialsBubbles: remove some reinterpret_casts. --- .../base/tensor_product_polynomials_bubbles.h | 211 +++++++++--------- 1 file changed, 102 insertions(+), 109 deletions(-) diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h index f2a8396978..96083ac2a1 100644 --- a/include/deal.II/base/tensor_product_polynomials_bubbles.h +++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h @@ -326,136 +326,129 @@ TensorProductPolynomialsBubbles::compute_derivative( if (i < max_q_indices) return tensor_polys.template compute_derivative(i, p); - const unsigned int comp = i - tensor_polys.n(); + [[maybe_unused]] const unsigned int comp = i - tensor_polys.n(); - Tensor derivative; - switch (order) + if constexpr (order == 1) { - case 1: + Tensor<1, dim> derivative; + for (unsigned int d = 0; d < dim; ++d) { - Tensor<1, dim> &derivative_1 = - *reinterpret_cast *>(&derivative); - - for (unsigned int d = 0; d < dim; ++d) - { - derivative_1[d] = 1.; - // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p) - for (unsigned j = 0; j < dim; ++j) - derivative_1[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) - derivative_1[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; - derivative_1[comp] += tmp; - } - - return derivative; + derivative[d] = 1.; + // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p) + for (unsigned j = 0; j < dim; ++j) + derivative[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) + derivative[d] *= 2 * p[comp] - 1; } - case 2: + + if (q_degree >= 2) { - Tensor<2, dim> &derivative_2 = - *reinterpret_cast *>(&derivative); + // 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; + derivative[comp] += tmp; + } - double v[dim + 1][3]; - { - for (unsigned int c = 0; c < dim; ++c) - { - v[c][0] = 4 * p[c] * (1 - p[c]); - v[c][1] = 4 * (1 - 2 * p[c]); - v[c][2] = -8; - } + return derivative; + } + else if constexpr (order == 2) + { + Tensor<2, dim> derivative; - double tmp = 1.; - for (unsigned int i = 0; i < q_degree - 1; ++i) - tmp *= 2 * p[comp] - 1; - v[dim][0] = tmp; + double v[dim + 1][3]; + { + for (unsigned int c = 0; c < dim; ++c) + { + v[c][0] = 4 * p[c] * (1 - p[c]); + v[c][1] = 4 * (1 - 2 * p[c]); + v[c][2] = -8; + } - if (q_degree >= 2) - { - double tmp = 2 * (q_degree - 1); - for (unsigned int i = 0; i < q_degree - 2; ++i) - tmp *= 2 * p[comp] - 1; - v[dim][1] = tmp; - } - else - v[dim][1] = 0.; + double tmp = 1.; + for (unsigned int i = 0; i < q_degree - 1; ++i) + tmp *= 2 * p[comp] - 1; + v[dim][0] = tmp; - if (q_degree >= 3) - { - double tmp = 4 * (q_degree - 2) * (q_degree - 1); - for (unsigned int i = 0; i < q_degree - 3; ++i) - tmp *= 2 * p[comp] - 1; - v[dim][2] = tmp; - } - else - v[dim][2] = 0.; + if (q_degree >= 2) + { + double tmp = 2 * (q_degree - 1); + for (unsigned int i = 0; i < q_degree - 2; ++i) + tmp *= 2 * p[comp] - 1; + v[dim][1] = tmp; } + else + v[dim][1] = 0.; - // calculate (\partial_j \partial_k \psi) * monomial - Tensor<2, dim> grad_grad_1; - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) + if (q_degree >= 3) + { + double tmp = 4 * (q_degree - 2) * (q_degree - 1); + for (unsigned int i = 0; i < q_degree - 3; ++i) + tmp *= 2 * p[comp] - 1; + v[dim][2] = tmp; + } + else + v[dim][2] = 0.; + } + + // calculate (\partial_j \partial_k \psi) * monomial + Tensor<2, dim> grad_grad_1; + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + { + grad_grad_1[d1][d2] = v[dim][0]; + for (unsigned int x = 0; x < dim; ++x) { - grad_grad_1[d1][d2] = v[dim][0]; - for (unsigned int x = 0; x < dim; ++x) + unsigned int derivative = 0; + if (d1 == x || d2 == x) { - unsigned int derivative = 0; - if (d1 == x || d2 == x) - { - if (d1 == d2) - derivative = 2; - else - derivative = 1; - } - grad_grad_1[d1][d2] *= v[x][derivative]; + if (d1 == d2) + derivative = 2; + else + derivative = 1; } + grad_grad_1[d1][d2] *= v[x][derivative]; } + } - // calculate (\partial_j \psi) *(\partial_k monomial) - // and (\partial_k \psi) *(\partial_j monomial) - Tensor<2, dim> grad_grad_2; - Tensor<2, dim> grad_grad_3; - for (unsigned int d = 0; d < dim; ++d) + // calculate (\partial_j \psi) *(\partial_k monomial) + // and (\partial_k \psi) *(\partial_j monomial) + Tensor<2, dim> grad_grad_2; + Tensor<2, dim> grad_grad_3; + for (unsigned int d = 0; d < dim; ++d) + { + grad_grad_2[d][comp] = v[dim][1]; + grad_grad_3[comp][d] = v[dim][1]; + for (unsigned int x = 0; x < dim; ++x) { - grad_grad_2[d][comp] = v[dim][1]; - grad_grad_3[comp][d] = v[dim][1]; - for (unsigned int x = 0; x < dim; ++x) - { - grad_grad_2[d][comp] *= v[x][d == x]; - grad_grad_3[comp][d] *= v[x][d == x]; - } + grad_grad_2[d][comp] *= v[x][d == x]; + grad_grad_3[comp][d] *= v[x][d == x]; } + } - // calculate \psi *(\partial j \partial_k monomial) and sum - double psi_value = 1.; - for (unsigned int x = 0; x < dim; ++x) - psi_value *= v[x][0]; + // calculate \psi *(\partial j \partial_k monomial) and sum + double psi_value = 1.; + for (unsigned int x = 0; x < dim; ++x) + psi_value *= v[x][0]; - for (unsigned int d1 = 0; d1 < dim; ++d1) - for (unsigned int d2 = 0; d2 < dim; ++d2) - derivative_2[d1][d2] = - grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2]; - derivative_2[comp][comp] += psi_value * v[dim][2]; + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + derivative[d1][d2] = + grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2]; + derivative[comp][comp] += psi_value * v[dim][2]; - return derivative; - } - default: - { - DEAL_II_NOT_IMPLEMENTED(); - return derivative; - } + return derivative; + } + else + { + DEAL_II_NOT_IMPLEMENTED(); + return {}; } } -- 2.39.5