Use Utilities::fixed_power() in step-82.

diff --git a/examples/step-82/step-82.cc b/examples/step-82/step-82.cc
index f865d9dcdb43452e231e2aff0ed6a726a7b15346..e05666259a24fd9d6b6875c98b529a7e4470ec9a 100644 (file)
--- a/examples/step-82/step-82.cc
+++ b/examples/step-82/step-82.cc
@@ -169,26 +169,28 @@ namespace Step82
 
     if (dim == 2)
       {
-        return_value = 24.0 * std::pow(p(1) * (1.0 - p(1)), 2) +
-                       +24.0 * std::pow(p(0) * (1.0 - p(0)), 2) +
+        return_value = 24.0 * Utilities::fixed_power<2>(p(1) * (1.0 - p(1))) +
+                       +24.0 * Utilities::fixed_power<2>(p(0) * (1.0 - p(0))) +
                        2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
                          (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1));
       }
     else if (dim == 3)
       {
-        return_value =
-          24.0 * std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2) +
-          24.0 * std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2) +
-          24.0 * std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2) +
-          2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
-            (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
-            std::pow(p(2) * (1.0 - p(2)), 2) +
-          2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
-            (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
-            std::pow(p(1) * (1.0 - p(1)), 2) +
-          2.0 * (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
-            (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
-            std::pow(p(0) * (1.0 - p(0)), 2);
+        return_value = 24.0 * Utilities::fixed_power<2>(p(1) * (1.0 - p(1)) *
+                                                        p(2) * (1.0 - p(2))) +
+                       24.0 * Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) *
+                                                        p(2) * (1.0 - p(2))) +
+                       24.0 * Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) *
+                                                        p(1) * (1.0 - p(1))) +
+                       2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
+                         (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
+                         Utilities::fixed_power<2>(p(2) * (1.0 - p(2))) +
+                       2.0 * (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
+                         (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
+                         Utilities::fixed_power<2>(p(1) * (1.0 - p(1))) +
+                       2.0 * (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
+                         (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
+                         Utilities::fixed_power<2>(p(0) * (1.0 - p(0)));
       }
     else
       Assert(false, ExcNotImplemented());
@@ -230,13 +232,13 @@ namespace Step82
 
     if (dim == 2)
       {
-        return_value = std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2);
+        return_value =
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)));
       }
     else if (dim == 3)
       {
-        return_value = std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)) *
-                                  p(2) * (1.0 - p(2)),
-                                2);
+        return_value = Utilities::fixed_power<2>(
+          p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)));
       }
     else
       Assert(false, ExcNotImplemented());
@@ -256,23 +258,28 @@ namespace Step82
     if (dim == 2)
       {
         return_gradient[0] =
-          (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
-          std::pow(p(1) * (1.0 - p(1)), 2);
+          (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) +
+           4.0 * Utilities::fixed_power<3>(p(0))) *
+          Utilities::fixed_power<2>(p(1) * (1.0 - p(1)));
         return_gradient[1] =
-          (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
-          std::pow(p(0) * (1.0 - p(0)), 2);
+          (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) +
+           4.0 * Utilities::fixed_power<3>(p(1))) *
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)));
       }
     else if (dim == 3)
       {
         return_gradient[0] =
-          (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
-          std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2);
+          (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) +
+           4.0 * Utilities::fixed_power<3>(p(0))) *
+          Utilities::fixed_power<2>(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)));
         return_gradient[1] =
-          (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
-          std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2);
+          (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) +
+           4.0 * Utilities::fixed_power<3>(p(1))) *
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)));
         return_gradient[2] =
-          (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) *
-          std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2);
+          (2.0 * p(2) - 6.0 * Utilities::fixed_power<2>(p(2)) +
+           4.0 * Utilities::fixed_power<3>(p(2))) *
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)));
       }
     else
       Assert(false, ExcNotImplemented());
@@ -292,36 +299,44 @@ namespace Step82
     if (dim == 2)
       {
         return_hessian[0][0] = (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
-                               std::pow(p(1) * (1.0 - p(1)), 2);
+                               Utilities::fixed_power<2>(p(1) * (1.0 - p(1)));
         return_hessian[0][1] =
-          (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
-          (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3));
+          (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) +
+           4.0 * Utilities::fixed_power<3>(p(0))) *
+          (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) +
+           4.0 * Utilities::fixed_power<3>(p(1)));
         return_hessian[1][1] = (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
-                               std::pow(p(0) * (1.0 - p(0)), 2);
+                               Utilities::fixed_power<2>(p(0) * (1.0 - p(0)));
       }
     else if (dim == 3)
       {
         return_hessian[0][0] =
           (2.0 - 12.0 * p(0) + 12.0 * p(0) * p(0)) *
-          std::pow(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)), 2);
+          Utilities::fixed_power<2>(p(1) * (1.0 - p(1)) * p(2) * (1.0 - p(2)));
         return_hessian[0][1] =
-          (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
-          (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
-          std::pow(p(2) * (1.0 - p(2)), 2);
+          (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) +
+           4.0 * Utilities::fixed_power<3>(p(0))) *
+          (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) +
+           4.0 * Utilities::fixed_power<3>(p(1))) *
+          Utilities::fixed_power<2>(p(2) * (1.0 - p(2)));
         return_hessian[0][2] =
-          (2.0 * p(0) - 6.0 * std::pow(p(0), 2) + 4.0 * std::pow(p(0), 3)) *
-          (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) *
-          std::pow(p(1) * (1.0 - p(1)), 2);
+          (2.0 * p(0) - 6.0 * Utilities::fixed_power<2>(p(0)) +
+           4.0 * Utilities::fixed_power<3>(p(0))) *
+          (2.0 * p(2) - 6.0 * Utilities::fixed_power<2>(p(2)) +
+           4.0 * Utilities::fixed_power<3>(p(2))) *
+          Utilities::fixed_power<2>(p(1) * (1.0 - p(1)));
         return_hessian[1][1] =
           (2.0 - 12.0 * p(1) + 12.0 * p(1) * p(1)) *
-          std::pow(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)), 2);
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(2) * (1.0 - p(2)));
         return_hessian[1][2] =
-          (2.0 * p(1) - 6.0 * std::pow(p(1), 2) + 4.0 * std::pow(p(1), 3)) *
-          (2.0 * p(2) - 6.0 * std::pow(p(2), 2) + 4.0 * std::pow(p(2), 3)) *
-          std::pow(p(0) * (1.0 - p(0)), 2);
+          (2.0 * p(1) - 6.0 * Utilities::fixed_power<2>(p(1)) +
+           4.0 * Utilities::fixed_power<3>(p(1))) *
+          (2.0 * p(2) - 6.0 * Utilities::fixed_power<2>(p(2)) +
+           4.0 * Utilities::fixed_power<3>(p(2))) *
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)));
         return_hessian[2][2] =
           (2.0 - 12.0 * p(2) + 12.0 * p(2) * p(2)) *
-          std::pow(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)), 2);
+          Utilities::fixed_power<2>(p(0) * (1.0 - p(0)) * p(1) * (1.0 - p(1)));
       }
     else
       Assert(false, ExcNotImplemented());
@@ -711,7 +726,7 @@ namespace Step82
 
             const double mesh_inv = 1.0 / face->diameter(); // h_e^{-1}
             const double mesh3_inv =
-              1.0 / std::pow(face->diameter(), 3); // h_e^{-3}
+              1.0 / Utilities::fixed_power<3>(face->diameter()); // h_e^{-3}
 
             fe_face.reinit(cell, face_no);
 
@@ -971,9 +986,9 @@ namespace Step82
                          solution_gradients_cell[q])
                           .norm_square() *
                         dx;
-            error_L2 += std::pow(u_exact.value(fe_values.quadrature_point(q)) -
-                                   solution_values_cell[q],
-                                 2) *
+            error_L2 += Utilities::fixed_power<2>(
+                          u_exact.value(fe_values.quadrature_point(q)) -
+                          solution_values_cell[q]) *
                         dx;
           } // for quadrature points
 
@@ -985,7 +1000,7 @@ namespace Step82
 
             const double mesh_inv = 1.0 / face->diameter(); // h^{-1}
             const double mesh3_inv =
-              1.0 / std::pow(face->diameter(), 3); // h^{-3}
+              1.0 / Utilities::fixed_power<3>(face->diameter()); // h^{-3}
 
             fe_face.reinit(cell, face_no);
 
@@ -1008,10 +1023,12 @@ namespace Step82
                       (u_exact_grad_q - solution_gradients[q]).norm_square() *
                       dx;
                     error_H2 += mesh3_inv *
-                                std::pow(u_exact_q - solution_values[q], 2) *
+                                Utilities::fixed_power<2>(u_exact_q -
+                                                          solution_values[q]) *
                                 dx;
                     error_H1 += mesh_inv *
-                                std::pow(u_exact_q - solution_values[q], 2) *
+                                Utilities::fixed_power<2>(u_exact_q -
+                                                          solution_values[q]) *
                                 dx;
                   }
               }
@@ -1052,16 +1069,16 @@ namespace Step82
                           (solution_gradients_neigh[q] - solution_gradients[q])
                             .norm_square() *
                           dx;
-                        error_H2 += mesh3_inv *
-                                    std::pow(solution_values_neigh[q] -
-                                               solution_values[q],
-                                             2) *
-                                    dx;
-                        error_H1 += mesh_inv *
-                                    std::pow(solution_values_neigh[q] -
-                                               solution_values[q],
-                                             2) *
-                                    dx;
+                        error_H2 +=
+                          mesh3_inv *
+                          Utilities::fixed_power<2>(solution_values_neigh[q] -
+                                                    solution_values[q]) *
+                          dx;
+                        error_H1 +=
+                          mesh_inv *
+                          Utilities::fixed_power<2>(solution_values_neigh[q] -
+                                                    solution_values[q]) *
+                          dx;
                       }
                   } // face not visited yet