From: Wolfgang Bangerth Date: Wed, 10 Jan 2024 19:48:37 +0000 (-0700) Subject: Address more places where we call std::pow. X-Git-Tag: relicensing~178^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16452%2Fhead;p=dealii.git Address more places where we call std::pow. --- diff --git a/source/base/function_lib.cc b/source/base/function_lib.cc index 221ceaaa35..7e11ff326a 100644 --- a/source/base/function_lib.cc +++ b/source/base/function_lib.cc @@ -1181,7 +1181,7 @@ namespace Functions const double phi = std::atan2(y, -x) + numbers::PI; const double r_squared = x * x + y * y; - return std::pow(r_squared, 1. / 3.) * std::sin(2. / 3. * phi); + return std::cbrt(r_squared) * std::sin(2. / 3. * phi); } @@ -1206,7 +1206,7 @@ namespace Functions const double phi = std::atan2(y, -x) + numbers::PI; const double r_squared = x * x + y * y; - values[i] = std::pow(r_squared, 1. / 3.) * std::sin(2. / 3. * phi); + values[i] = std::cbrt(r_squared) * std::sin(2. / 3. * phi); } } } @@ -1235,8 +1235,7 @@ namespace Functions const double phi = std::atan2(y, -x) + numbers::PI; const double r_squared = x * x + y * y; - values[i](0) = - std::pow(r_squared, 1. / 3.) * std::sin(2. / 3. * phi); + values[i](0) = std::cbrt(r_squared) * std::sin(2. / 3. * phi); } } } @@ -2965,10 +2964,12 @@ namespace Functions const double pi_y = numbers::PI * point(1); const double pi_t = numbers::PI / T * this->get_time(); - values[0] = -2 * std::cos(pi_t) * std::pow(std::sin(pi_x), 2) * - std::sin(pi_y) * std::cos(pi_y); - values[1] = +2 * std::cos(pi_t) * std::pow(std::sin(pi_y), 2) * - std::sin(pi_x) * std::cos(pi_x); + values[0] = -2 * std::cos(pi_t) * + Utilities::fixed_power<2>(std::sin(pi_x)) * std::sin(pi_y) * + std::cos(pi_y); + values[1] = +2 * std::cos(pi_t) * + Utilities::fixed_power<2>(std::sin(pi_y)) * std::sin(pi_x) * + std::cos(pi_x); if (dim == 3) values[2] = 0; diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index fd1207cff0..8e640c3463 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -127,8 +127,8 @@ namespace internal std::vector w(q); const long double factor = - std::pow(2., alpha + beta + 1) * gamma(alpha + q) * gamma(beta + q) / - ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1)); + Utilities::pow(2, alpha + beta + 1) * gamma(alpha + q) * + gamma(beta + q) / ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1)); for (unsigned int i = 0; i < q; ++i) { const long double s = @@ -1168,22 +1168,19 @@ QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity) // We need to check if the singularity is at the boundary of the interval. if (std::abs(eta_star) <= tol) { - gamma_bar = - std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) + - std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) + - eta_bar; + gamma_bar = std::cbrt(eta_bar * eta_star + std::abs(eta_star)) + + std::cbrt(eta_bar * eta_star - std::abs(eta_star)) + eta_bar; } else { - gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) / - std::abs(eta_bar * eta_star + std::abs(eta_star)) * - std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)), - 1.0 / 3.0) + - (eta_bar * eta_star - std::abs(eta_star)) / - std::abs(eta_bar * eta_star - std::abs(eta_star)) * - std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)), - 1.0 / 3.0) + - eta_bar; + gamma_bar = + (eta_bar * eta_star + std::abs(eta_star)) / + std::abs(eta_bar * eta_star + std::abs(eta_star)) * + std::cbrt(std::abs(eta_bar * eta_star + std::abs(eta_star))) + + (eta_bar * eta_star - std::abs(eta_star)) / + std::abs(eta_bar * eta_star - std::abs(eta_star)) * + std::cbrt(std::abs(eta_bar * eta_star - std::abs(eta_star))) + + eta_bar; } for (unsigned int q = 0; q < quadrature_points.size(); ++q) { diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index 794f0d0e1f..4c5e51ada7 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -484,21 +484,17 @@ FE_Q_Hierarchical::build_dofs_cell( // factor == k * (k-1) * ... * (k-j+1) / j! = k! / (k-j)! / j! if (c == 0) { - dofs_subcell[c](j, k) = - ((k + j) % 2 == 0) ? - std::pow(.5, static_cast(k)) * factor : - -std::pow(.5, static_cast(k)) * factor; - dofs_cell[c](j, k) = - std::pow(2., static_cast(j)) * factor; + dofs_subcell[c](j, k) = ((k + j) % 2 == 0) ? + Utilities::pow(.5, k) * factor : + -Utilities::pow(.5, k) * factor; + dofs_cell[c](j, k) = Utilities::pow(2, j) * factor; } else { - dofs_subcell[c](j, k) = - std::pow(.5, static_cast(k)) * factor; - dofs_cell[c](j, k) = - ((k + j) % 2 == 0) ? - std::pow(2., static_cast(j)) * factor : - -std::pow(2., static_cast(j)) * factor; + dofs_subcell[c](j, k) = Utilities::pow(.5, k) * factor; + dofs_cell[c](j, k) = ((k + j) % 2 == 0) ? + Utilities::pow(2, j) * factor : + -Utilities::pow(2, j) * factor; } } } @@ -1082,7 +1078,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( for (unsigned int i = 2; i < this->n_dofs_per_face(face_no); ++i) { - interpolation_matrix(i, i) = std::pow(0.5, i); + interpolation_matrix(i, i) = Utilities::pow(0.5, i); factorial_i *= i; int factorial_j = factorial_i; int factorial_ij = 1; @@ -1096,12 +1092,12 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( if (((i + j) & 1) != 0u) interpolation_matrix(i, j) = - -1.0 * std::pow(0.5, j) * factorial_j / + -1.0 * Utilities::pow(0.5, j) * factorial_j / (factorial_i * factorial_ij); else interpolation_matrix(i, j) = - std::pow(0.5, j) * factorial_j / + Utilities::pow(0.5, j) * factorial_j / (factorial_i * factorial_ij); } } @@ -1128,7 +1124,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( for (unsigned int i = 2; i < this->n_dofs_per_face(face_no); ++i) { - interpolation_matrix(i, i) = std::pow(0.5, i); + interpolation_matrix(i, i) = Utilities::pow(0.5, i); factorial_i *= i; int factorial_j = factorial_i; int factorial_ij = 1; @@ -1140,7 +1136,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( factorial_ij *= j - i; factorial_j *= j; interpolation_matrix(i, j) = - std::pow(0.5, j) * factorial_j / + Utilities::pow(0.5, j) * factorial_j / (factorial_i * factorial_ij); } } @@ -1194,8 +1190,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( for (unsigned int i = 2; i <= this->degree; ++i) { - double tmp = std::pow(0.5, i); - interpolation_matrix(i + 2, i + 2) = tmp; + double tmp = Utilities::pow(0.5, i); + interpolation_matrix(i + 2, i + 2) = tmp; interpolation_matrix(i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; tmp *= 0.5; @@ -1227,7 +1223,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix(i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = - std::pow(0.5, i + j); + Utilities::pow(0.5, i + j); factorial_k *= j; int factorial_kl = 1; int factorial_l = factorial_k; @@ -1241,14 +1237,14 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = - -1.0 * std::pow(0.5, i + k) * factorial_l / - (factorial_k * factorial_kl); + -1.0 * Utilities::pow(0.5, i + k) * + factorial_l / (factorial_k * factorial_kl); else interpolation_matrix( i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = - std::pow(0.5, i + k) * factorial_l / + Utilities::pow(0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } @@ -1264,8 +1260,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( if (((i + j) & 1) != 0u) { - tmp = -1.0 * std::pow(0.5, j) * factorial_j / - (factorial_i * factorial_ij); + tmp = -1.0 * Utilities::pow(0.5, j) * + factorial_j / (factorial_i * factorial_ij); interpolation_matrix(i + 2, j + 2) = tmp; interpolation_matrix(i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; @@ -1276,7 +1272,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = - tmp * std::pow(0.5, k); + tmp * Utilities::pow(0.5, k); factorial_k *= k; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1292,7 +1288,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - -1.0 * tmp * std::pow(0.5, l) * + -1.0 * tmp * Utilities::pow(0.5, l) * factorial_l / (factorial_k * factorial_kl); @@ -1300,7 +1296,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - tmp * std::pow(0.5, l) * factorial_l / + tmp * Utilities::pow(0.5, l) * + factorial_l / (factorial_k * factorial_kl); } } @@ -1330,7 +1327,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( } else { - tmp = std::pow(0.5, j) * factorial_j / + tmp = Utilities::pow(0.5, j) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix(i + 2, j + 2) = tmp; interpolation_matrix(i + 2 * source_fe.degree, @@ -1342,7 +1339,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = - tmp * std::pow(0.5, k); + tmp * Utilities::pow(0.5, k); factorial_k *= k; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1358,7 +1355,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - -1.0 * tmp * std::pow(0.5, l) * + -1.0 * tmp * Utilities::pow(0.5, l) * factorial_l / (factorial_k * factorial_kl); @@ -1366,7 +1363,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - tmp * std::pow(0.5, l) * factorial_l / + tmp * Utilities::pow(0.5, l) * + factorial_l / (factorial_k * factorial_kl); } } @@ -1440,8 +1438,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( for (unsigned int i = 2; i <= this->degree; ++i) { - double tmp = std::pow(0.5, i + 1); - interpolation_matrix(i + 2, i + 2) = tmp; + double tmp = Utilities::pow(0.5, i + 1); + interpolation_matrix(i + 2, i + 2) = tmp; interpolation_matrix(i + 2, i + this->degree + 1) = tmp; interpolation_matrix(i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp; @@ -1473,7 +1471,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( { factorial_ij *= j - i; factorial_j *= j; - tmp = std::pow(0.5, j) * factorial_j / + tmp = Utilities::pow(0.5, j) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix(i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; @@ -1484,7 +1482,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = - tmp * std::pow(0.5, k); + tmp * Utilities::pow(0.5, k); factorial_k *= k; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1499,7 +1497,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - -1.0 * tmp * std::pow(0.5, l) * + -1.0 * tmp * Utilities::pow(0.5, l) * factorial_l / (factorial_k * factorial_kl); @@ -1507,7 +1505,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - tmp * std::pow(0.5, l) * factorial_l / + tmp * Utilities::pow(0.5, l) * + factorial_l / (factorial_k * factorial_kl); } } @@ -1555,7 +1554,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix(i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = - std::pow(0.5, i + j); + Utilities::pow(0.5, i + j); factorial_k *= j; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1569,14 +1568,14 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = - -1.0 * std::pow(0.5, i + k) * factorial_l / - (factorial_k * factorial_kl); + -1.0 * Utilities::pow(0.5, i + k) * + factorial_l / (factorial_k * factorial_kl); else interpolation_matrix( i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = - std::pow(0.5, i + k) * factorial_l / + Utilities::pow(0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } @@ -1624,8 +1623,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( for (unsigned int i = 2; i <= this->degree; ++i) { - double tmp = std::pow(0.5, i); - interpolation_matrix(i + 2, i + 2) = tmp; + double tmp = Utilities::pow(0.5, i); + interpolation_matrix(i + 2, i + 2) = tmp; interpolation_matrix(i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; tmp *= 0.5; @@ -1657,7 +1656,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix(i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = - std::pow(0.5, i + j); + Utilities::pow(0.5, i + j); factorial_k *= j; int factorial_kl = 1; int factorial_l = factorial_k; @@ -1669,7 +1668,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = - std::pow(0.5, i + k) * factorial_l / + Utilities::pow(0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } @@ -1682,7 +1681,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( { factorial_ij *= j - i; factorial_j *= j; - tmp = std::pow(0.5, j) * factorial_j / + tmp = Utilities::pow(0.5, j) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix(i + 2, j + 2) = tmp; tmp *= -1.0; @@ -1718,7 +1717,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = - tmp * std::pow(0.5, k); + tmp * Utilities::pow(0.5, k); factorial_k *= k; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1731,7 +1730,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - tmp * std::pow(0.5, l) * factorial_l / + tmp * Utilities::pow(0.5, l) * factorial_l / (factorial_k * factorial_kl); } } @@ -1790,8 +1789,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( for (unsigned int i = 2; i <= this->degree; ++i) { - double tmp = std::pow(0.5, i + 1); - interpolation_matrix(i + 2, i + 2) = tmp; + double tmp = Utilities::pow(0.5, i + 1); + interpolation_matrix(i + 2, i + 2) = tmp; interpolation_matrix(i + 2, i + this->degree + 1) = tmp; interpolation_matrix(i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; @@ -1822,7 +1821,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix(i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = - std::pow(0.5, i + j); + Utilities::pow(0.5, i + j); factorial_k *= j; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1834,7 +1833,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = - std::pow(0.5, i + k) * factorial_l / + Utilities::pow(0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } @@ -1847,7 +1846,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( { factorial_ij *= j - i; factorial_j *= j; - tmp = std::pow(0.5, j + 1) * factorial_j / + tmp = Utilities::pow(0.5, j + 1) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix(i + 2, j + 2) = tmp; interpolation_matrix(i + 2, j + this->degree + 1) = @@ -1868,7 +1867,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = - tmp * std::pow(0.5, k); + tmp * Utilities::pow(0.5, k); factorial_k *= k; int factorial_l = factorial_k; int factorial_kl = 1; @@ -1881,7 +1880,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix( i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = - tmp * std::pow(0.5, l) * factorial_l / + tmp * Utilities::pow(0.5, l) * factorial_l / (factorial_k * factorial_kl); } } diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index f598404ba2..03a7078410 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -665,7 +665,7 @@ namespace GridGenerator const double x = i * 1 / (1.0 * number_points - 1); const double y_t = 5 * t * - (0.2969 * std::pow(x, 0.5) - 0.126 * x - + (0.2969 * std::sqrt(x) - 0.126 * x - 0.3516 * Utilities::fixed_power<2>(x) + 0.2843 * Utilities::fixed_power<3>(x) - 0.1036 * Utilities::fixed_power<4>( @@ -696,7 +696,7 @@ namespace GridGenerator const double y_t = 5 * t * - (0.2969 * std::pow(x, 0.5) - 0.126 * x - + (0.2969 * std::sqrt(x) - 0.126 * x - 0.3516 * Utilities::fixed_power<2>(x) + 0.2843 * Utilities::fixed_power<3>(x) - 0.1036 * Utilities::fixed_power<4>( diff --git a/source/grid/grid_refinement.cc b/source/grid/grid_refinement.cc index 9d5fd9d402..eae7ca90e5 100644 --- a/source/grid/grid_refinement.cc +++ b/source/grid/grid_refinement.cc @@ -472,14 +472,15 @@ GridRefinement::refine_and_coarsen_optimize(Triangulation &tria, double min_cost = std::numeric_limits::max(); std::size_t min_arg = 0; + const double reduction_factor = (1. - std::pow(2., -1. * order)); for (std::size_t M = 0; M < criteria.size(); ++M) { - expected_error_reduction += - (1 - std::pow(2., -1. * order)) * criteria(cell_indices[M]); + expected_error_reduction += reduction_factor * criteria(cell_indices[M]); - const double cost = std::pow(((std::pow(2., dim) - 1) * (1 + M) + N), - static_cast(order) / dim) * - (original_error - expected_error_reduction); + const double cost = + std::pow(((Utilities::fixed_power(2) - 1) * (1 + M) + N), + static_cast(order) / dim) * + (original_error - expected_error_reduction); if (cost <= min_cost) { min_cost = cost; diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index c1d33574e7..dad68daa4e 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2427,10 +2427,10 @@ namespace GridTools double objective = 0; for (unsigned int c = 0; c < object->n_children(); ++c) for (const unsigned int i : object->child(c)->vertex_indices()) - objective += - (child_alternating_forms[c][i] - - average_parent_alternating_form / std::pow(2., 1. * structdim)) - .norm_square(); + objective += (child_alternating_forms[c][i] - + average_parent_alternating_form / + Utilities::fixed_power(2)) + .norm_square(); return objective; } diff --git a/source/grid/tria.cc b/source/grid/tria.cc index f1e48f669c..752e9b4ddd 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -1577,7 +1577,8 @@ namespace GeometryInfo::alternating_form_at_vertices(vertices, determinants); for (const unsigned int i : GeometryInfo::vertex_indices()) - if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim)) + if (determinants[i] <= + 1e-9 * Utilities::fixed_power(cell->diameter())) { distorted_cells.distorted_cells.push_back(cell); break; @@ -1612,7 +1613,7 @@ namespace for (const unsigned int i : GeometryInfo::vertex_indices()) if (determinants[i] <= - 1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim)) + 1e-9 * Utilities::fixed_power(cell->child(c)->diameter())) return true; } diff --git a/source/hp/refinement.cc b/source/hp/refinement.cc index 36f6851f34..fa3cec1e8d 100644 --- a/source/hp/refinement.cc +++ b/source/hp/refinement.cc @@ -633,16 +633,26 @@ namespace hp // step 1: exponential decay with p-adaptation if (cell->future_fe_index_set()) { - predicted_errors[cell->active_cell_index()] *= - std::pow(gamma_p, - int(future_fe_degree) - int(cell->get_fe().degree)); + if (future_fe_degree > cell->get_fe().degree) + predicted_errors[cell->active_cell_index()] *= + Utilities::pow(gamma_p, + future_fe_degree - cell->get_fe().degree); + else if (future_fe_degree < cell->get_fe().degree) + predicted_errors[cell->active_cell_index()] /= + Utilities::pow(gamma_p, + cell->get_fe().degree - future_fe_degree); + else + { + // The two degrees are the same; we do not need to + // adapt the predicted error + } } // step 2: algebraic decay with h-adaptation if (cell->refine_flag_set()) { predicted_errors[cell->active_cell_index()] *= - (gamma_h * std::pow(.5, future_fe_degree)); + (gamma_h * Utilities::pow(.5, future_fe_degree)); // predicted error will be split on children cells // after adaptation via CellDataTransfer @@ -650,7 +660,7 @@ namespace hp else if (cell->coarsen_flag_set()) { predicted_errors[cell->active_cell_index()] /= - (gamma_h * std::pow(.5, future_fe_degree)); + (gamma_h * Utilities::pow(.5, future_fe_degree)); // predicted error will be summed up on parent cell // after adaptation via CellDataTransfer