Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Address more places where we call std::pow. #16452

Merged
merged 1 commit into from
Jan 11, 2024
+106 −97
Merged

Address more places where we call std::pow. #16452

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
17 changes: 9 additions & 8 deletions source/base/function_lib.cc
View file
Open in desktop
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
Open in desktop
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
Open in desktop

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions source/grid/grid_generator.cc
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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