From d4cd6a535c99bdd9899ecec6f6f4c9d4ecd2a414 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 10 Jan 2024 09:29:31 -0700 Subject: [PATCH] Use Utilities::fixed_power() in step-82. --- examples/step-82/step-82.cc | 141 ++++++++++++++++++++---------------- 1 file changed, 79 insertions(+), 62 deletions(-) diff --git a/examples/step-82/step-82.cc b/examples/step-82/step-82.cc index f865d9dcdb..e05666259a 100644 --- 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 -- 2.39.5