Merge pull request #16441 from bangerth/pow-2

Use Utilities::fixed_power() where possible.

examples/step-25/step-25.cc

@@ -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
@@ -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.;
   }

examples/step-43/step-43.cc

@@ -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;
 

examples/step-44/step-44.cc

@@ -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;

examples/step-47/step-47.cc

@@ -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]);
       }
     };

examples/step-53/step-53.cc

@@ -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)));

examples/step-55/step-55.cc

@@ -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;
 

examples/step-62/step-62.cc

@@ -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
           {
@@ -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;

examples/step-71/step-71.cc

@@ -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;
 

examples/step-85/step-85.cc

@@ -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);
               }
           }
       }

source/base/data_out_base.cc

@@ -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.));
+                        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.;

source/grid/grid_out.cc

@@ -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));
 
@@ -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.));
+                          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));