Skip to content

Commit

Permalink
Use Utilities::fixed_power() in step-82.
Browse files Browse the repository at this point in the history
  • Loading branch information
@bangerth
bangerth committed Jan 10, 2024
1 parent df76733 commit b528e0d
Showing 1 changed file with 64 additions and 46 deletions.
110 changes: 64 additions & 46 deletions examples/step-82/step-82.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down Expand Up @@ -230,7 +232,8 @@ 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)
{
Expand All @@ -256,23 +259,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());
Expand All @@ -292,36 +300,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());
Expand Down Expand Up @@ -711,7 +727,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);

Expand Down Expand Up @@ -985,7 +1001,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);

Expand All @@ -1008,10 +1024,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;
}
}
Expand Down

0 comments on commit b528e0d

Please sign in to comment.