Skip to content

Commit

Permalink
Merge pull request #16452 from bangerth/pow-5
Browse files Browse the repository at this point in the history 
Address more places where we call std::pow.
  • Loading branch information
@blaisb
blaisb committed Jan 11, 2024
2 parents bfc8f4b + c739370 commit 253e516
Show file tree
Hide file tree
Showing 8 changed files with 106 additions and 97 deletions.
17 changes: 9 additions & 8 deletions source/base/function_lib.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}


Expand All @@ -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);
}
}
}
Expand Down Expand Up @@ -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);
}
}
}
Expand Down Expand Up @@ -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;
Expand Down
27 changes: 12 additions & 15 deletions source/base/quadrature_lib.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,8 @@ namespace internal
std::vector<long double> 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 =
Expand Down Expand Up @@ -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)
{
Expand Down
111 changes: 55 additions & 56 deletions source/fe/fe_q_hierarchical.cc
View file

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions source/grid/grid_generator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>(
Expand Down Expand Up @@ -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>(
Expand Down
11 changes: 6 additions & 5 deletions source/grid/grid_refinement.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -472,14 +472,15 @@ GridRefinement::refine_and_coarsen_optimize(Triangulation<dim, spacedim> &tria,
double min_cost = std::numeric_limits<double>::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<double>(order) / dim) *
(original_error - expected_error_reduction);
const double cost =
std::pow(((Utilities::fixed_power<dim>(2) - 1) * (1 + M) + N),
static_cast<double>(order) / dim) *
(original_error - expected_error_reduction);
if (cost <= min_cost)
{
min_cost = cost;
Expand Down
8 changes: 4 additions & 4 deletions source/grid/grid_tools.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<structdim>(2))
.norm_square();

return objective;
}
Expand Down
5 changes: 3 additions & 2 deletions source/grid/tria.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1577,7 +1577,8 @@ namespace
GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);

for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
if (determinants[i] <= 1e-9 * std::pow(cell->diameter(), 1. * dim))
if (determinants[i] <=
1e-9 * Utilities::fixed_power<dim>(cell->diameter()))
{
distorted_cells.distorted_cells.push_back(cell);
break;
Expand Down Expand Up @@ -1612,7 +1613,7 @@ namespace

for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
if (determinants[i] <=
1e-9 * std::pow(cell->child(c)->diameter(), 1. * dim))
1e-9 * Utilities::fixed_power<dim>(cell->child(c)->diameter()))
return true;
}

Expand Down
20 changes: 15 additions & 5 deletions source/hp/refinement.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -633,24 +633,34 @@ 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
}
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
Expand Down

0 comments on commit 253e516

Please sign in to comment.