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

Use Utilities::fixed_power() where possible. #16441

Merged
merged 2 commits into from
Jan 16, 2024
Merged

Use Utilities::fixed_power() where possible. #16441

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bfe7245
Use Utilities::fixed_power() where possible.
bangerth Jan 10, 2024
907d0f0
Use std::hypot() instead of std::sqrt(sum of squares).
bangerth Jan 10, 2024
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
12 changes: 7 additions & 5 deletions examples/step-25/step-25.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -336,11 +336,12 @@ namespace Step25
// First we assemble the Jacobian matrix $F'_h(U^{n,l})$, where $U^{n,l}$
// is stored in the vector <code>solution</code> for convenience.
system_matrix.copy_from(mass_matrix);
system_matrix.add(std::pow(time_step * theta, 2), laplace_matrix);
system_matrix.add(Utilities::fixed_power<2>(time_step * theta),
laplace_matrix);

SparseMatrix<double> tmp_matrix(sparsity_pattern);
compute_nl_matrix(old_solution, solution, tmp_matrix);
system_matrix.add(std::pow(time_step * theta, 2), tmp_matrix);
system_matrix.add(Utilities::fixed_power<2>(time_step * theta), tmp_matrix);

// Next we compute the right-hand side vector. This is just the
// combination of matrix-vector products implied by the description of
Expand All @@ -351,17 +352,18 @@ namespace Step25

mass_matrix.vmult(system_rhs, solution);
laplace_matrix.vmult(tmp_vector, solution);
system_rhs.add(std::pow(time_step * theta, 2), tmp_vector);
system_rhs.add(Utilities::fixed_power<2>(time_step * theta), tmp_vector);

mass_matrix.vmult(tmp_vector, old_solution);
system_rhs.add(-1.0, tmp_vector);
laplace_matrix.vmult(tmp_vector, old_solution);
system_rhs.add(std::pow(time_step, 2) * theta * (1 - theta), tmp_vector);
system_rhs.add(Utilities::fixed_power<2>(time_step) * theta * (1 - theta),
tmp_vector);

system_rhs.add(-time_step, M_x_velocity);

compute_nl_term(old_solution, solution, tmp_vector);
system_rhs.add(std::pow(time_step, 2) * theta, tmp_vector);
system_rhs.add(Utilities::fixed_power<2>(time_step) * theta, tmp_vector);

system_rhs *= -1.;
}
Expand Down
2 changes: 1 addition & 1 deletion examples/step-43/step-43.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -300,7 +300,7 @@ namespace Step43

const double numerator =
2.0 * S * temp - S * S * (2.0 * S - 2.0 * viscosity * (1 - S));
const double denominator = std::pow(temp, 2.0);
const double denominator = Utilities::fixed_power<2>(temp);

const double F_prime = numerator / denominator;

Expand Down
2 changes: 1 addition & 1 deletion examples/step-44/step-44.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1862,7 +1862,7 @@ namespace Step44
const double det_F_qp = lqph[q_point]->get_det_F();
const double J_tilde_qp = lqph[q_point]->get_J_tilde();
const double the_error_qp_squared =
std::pow((det_F_qp - J_tilde_qp), 2);
Utilities::fixed_power<2>((det_F_qp - J_tilde_qp));
const double JxW = fe_values.JxW(q_point);

dil_L2_error += the_error_qp_squared * JxW;
Expand Down
2 changes: 1 addition & 1 deletion examples/step-47/step-47.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ namespace Step47
const unsigned int /*component*/ = 0) const override

{
return 4 * std::pow(PI, 4.0) * std::sin(PI * p[0]) *
return 4 * Utilities::fixed_power<4>(PI) * std::sin(PI * p[0]) *
std::sin(PI * p[1]);
}
};
Expand Down
6 changes: 3 additions & 3 deletions examples/step-53/step-53.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,9 +263,9 @@ namespace Step53
const double th = std::atan2(R * x(2), b * p);
const double phi = std::atan2(x(1), x(0));
const double theta =
std::atan2(x(2) + ep * ep * b * std::pow(std::sin(th), 3),
(p -
(ellipticity * ellipticity * R * std::pow(std::cos(th), 3))));
std::atan2(x(2) + ep * ep * b * Utilities::fixed_power<3>(std::sin(th)),
(p - (ellipticity * ellipticity * R *
Utilities::fixed_power<3>(std::cos(th)))));
const double R_bar =
R / (std::sqrt(1 - ellipticity * ellipticity * std::sin(theta) *
std::sin(theta)));
Expand Down
6 changes: 4 additions & 2 deletions examples/step-55/step-55.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,13 +210,15 @@ namespace Step55
std::exp(R_x * (-2 * std::sqrt(25.0 + 4 * pi2) + 10.0)) -
0.4 * pi2 * std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
std::cos(2 * R_y * pi) +
0.1 * std::pow(-std::sqrt(25.0 + 4 * pi2) + 5.0, 2) *
0.1 *
Utilities::fixed_power<2>(-std::sqrt(25.0 + 4 * pi2) + 5.0) *
std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
std::cos(2 * R_y * pi);
values[1] = 0.2 * pi * (-std::sqrt(25.0 + 4 * pi2) + 5.0) *
std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
std::sin(2 * R_y * pi) -
0.05 * std::pow(-std::sqrt(25.0 + 4 * pi2) + 5.0, 3) *
0.05 *
Utilities::fixed_power<3>(-std::sqrt(25.0 + 4 * pi2) + 5.0) *
std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
std::sin(2 * R_y * pi) / pi;

Expand Down
11 changes: 6 additions & 5 deletions examples/step-62/step-62.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -402,10 +402,11 @@ namespace step62
std::abs(p[1] - force_center[1]) < max_force_width_y / 2)
{
return max_force_amplitude *
std::exp(-(std::pow(p[0] - force_center[0], 2) /
(2 * std::pow(force_sigma_x, 2)) +
std::pow(p[1] - force_center[1], 2) /
(2 * std::pow(force_sigma_y, 2))));
std::exp(
-(Utilities::fixed_power<2>(p[0] - force_center[0]) /
(2 * Utilities::fixed_power<2>(force_sigma_x)) +
Utilities::fixed_power<2>(p[1] - force_center[1]) /
(2 * Utilities::fixed_power<2>(force_sigma_y))));
}
else
{
Expand Down Expand Up @@ -967,7 +968,7 @@ namespace step62
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
std::complex<double> matrix_sum = 0;
matrix_sum += -std::pow(omega, 2) *
matrix_sum += -Utilities::fixed_power<2>(omega) *
quadrature_data.mass_coefficient[i][j];
matrix_sum += quadrature_data.stiffness_coefficient[i][j];
cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;
Expand Down
2 changes: 1 addition & 1 deletion examples/step-71/step-71.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2319,7 +2319,7 @@ namespace Step71
// The first derivative of the saturation function, noting that
// $\frac{d \tanh(x)}{dx} = \text{sech}^{2}(x)$.
const double dtanh_two_h_dot_h_div_h_sat_squ =
std::pow(1.0 / std::cosh(two_h_dot_h_div_h_sat_squ), 2.0);
Utilities::fixed_power<2>(1.0 / std::cosh(two_h_dot_h_div_h_sat_squ));
const Tensor<1, dim> dtwo_h_dot_h_div_h_sat_squ_dH =
2.0 * 2.0 / (this->get_mu_e_h_sat() * this->get_mu_e_h_sat()) * H;

Expand Down
2 changes: 1 addition & 1 deletion examples/step-85/step-85.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -650,7 +650,7 @@ namespace Step85
const double error_at_point =
solution_values.at(q) - analytical_solution.value(point);
error_L2_squared +=
std::pow(error_at_point, 2) * fe_values->JxW(q);
Utilities::fixed_power<2>(error_at_point) * fe_values->JxW(q);
}
}
}
Expand Down
7 changes: 3 additions & 4 deletions source/base/data_out_base.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4166,10 +4166,9 @@ namespace DataOutBase
h1(0) * h2(1) - h1(1) * h2(0);

// normalize Vector
double norm =
std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
std::pow(nrml[i * d1 + j * d2](1), 2.) +
std::pow(nrml[i * d1 + j * d2](2), 2.));
Comment on lines -4170 to -4172
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@drwells I think this is an example where std::pow and Utilities::fixed_power will be different. std::pow(x,2.) does not expand to x*x as far as I remember trying. You need to have std::pow(x,2) with an integer. I think this PR makes this a lot clearer. I prefer expliciting the fixed_power.
<

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, you could alternatively also use std::hypot (works in 3D since C++17).

double norm = std::hypot(nrml[i * d1 + j * d2](0),
nrml[i * d1 + j * d2](1),
nrml[i * d1 + j * d2](2));

if (nrml[i * d1 + j * d2](1) < 0)
norm *= -1.;
Expand Down
10 changes: 6 additions & 4 deletions source/base/function_lib.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2965,10 +2965,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
14 changes: 7 additions & 7 deletions source/grid/grid_out.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2455,9 +2455,9 @@ GridOut::write_svg(const Triangulation<2, 2> &tria, std::ostream &out) const
}

const double distance_to_camera =
std::sqrt(std::pow(point[0] - camera_position[0], 2.) +
std::pow(point[1] - camera_position[1], 2.) +
std::pow(point[2] - camera_position[2], 2.));
std::hypot(point[0] - camera_position[0],
point[1] - camera_position[1],
point[2] - camera_position[2]);
const double distance_factor =
distance_to_camera / (2. * std::max(x_dimension, y_dimension));

Expand Down Expand Up @@ -2625,10 +2625,10 @@ GridOut::write_svg(const Triangulation<2, 2> &tria, std::ostream &out) const

if (svg_flags.label_boundary_id)
{
const double distance_to_camera = std::sqrt(
std::pow(point[0] - camera_position[0], 2.) +
std::pow(point[1] - camera_position[1], 2.) +
std::pow(point[2] - camera_position[2], 2.));
Comment on lines -2628 to -2631
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here

const double distance_to_camera =
std::hypot(point[0] - camera_position[0],
point[1] - camera_position[1],
point[2] - camera_position[2]);
const double distance_factor =
distance_to_camera /
(2. * std::max(x_dimension, y_dimension));
Expand Down