#include <deal.II/base/logstream.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/std_cxx17/optional.h>
+#include <deal.II/base/utilities.h>
#include <deal.II/numerics/history.h>
const NumberType r1 = f2 - f1 - g1 * x2_shift;
const NumberType r2 = f3 - f1 - g1 * x3_shift;
const NumberType denom =
- std::pow(x2_shift * x3_shift, 2) * (x2_shift - x3_shift);
+ Utilities::fixed_power<2>(x2_shift * x3_shift) * (x2_shift - x3_shift);
if (denom == 0.)
return {};
- const NumberType A =
- (r1 * std::pow(x3_shift, 2) - r2 * std::pow(x2_shift, 2)) / denom;
- const NumberType B =
- (r2 * std::pow(x2_shift, 3) - r1 * std::pow(x3_shift, 3)) / denom;
+ const NumberType A = (r1 * Utilities::fixed_power<2>(x3_shift) -
+ r2 * Utilities::fixed_power<2>(x2_shift)) /
+ denom;
+ const NumberType B = (r2 * Utilities::fixed_power<3>(x2_shift) -
+ r1 * Utilities::fixed_power<3>(x3_shift)) /
+ denom;
const NumberType &C = g1;
// now get the minimizer:
if (d >= r)
return Tensor<1, dim>();
const double e = -d * d / (r - d) / (r + d);
- return ((e < -50) ? Point<dim>() :
- (p - this->center) / d *
- (-2.0 * r * r / std::pow(-r * r + d * d, 2.0) * d *
- std::exp(e)) *
- this->rescaling);
+ return ((e < -50) ?
+ Point<dim>() :
+ (p - this->center) / d *
+ (-2.0 * r * r / Utilities::fixed_power<2>(-r * r + d * d) * d *
+ std::exp(e)) *
+ this->rescaling);
}
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function_signed_distance.h>
+#include <deal.II/base/utilities.h>
#include <algorithm>
const SymmetricTensor<2, dim> hess =
unit_symmetric_tensor<dim>() / distance -
symmetrize(outer_product(center_to_point, center_to_point)) /
- std::pow(distance, 3);
+ Utilities::fixed_power<3>(distance);
return hess;
}
{
double val = 0.0;
for (unsigned int d = 0; d < dim; ++d)
- val += std::pow((point[d] - center[d]) / radii[d], 2);
+ val += Utilities::fixed_power<2>((point[d] - center[d]) / radii[d]);
return val - 1.0;
}
do
{
// compute the ellipse evolute (center of curvature) for the current t
- const double ex = (a * a - b * b) * std::pow(std::cos(t), 3) / a;
- const double ey = (b * b - a * a) * std::pow(std::sin(t), 3) / b;
+ const double ex =
+ (a * a - b * b) * Utilities::fixed_power<3>(std::cos(t)) / a;
+ const double ey =
+ (b * b - a * a) * Utilities::fixed_power<3>(std::sin(t)) / b;
// compute distances from current point on ellipse to its evolute
const double rx = x - ex;
const double ry = y - ey;
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/polynomial.h>
#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_values.h>
for (unsigned int q = 0; q < quadrature_points.size(); ++q)
{
double gamma = quadrature_points[q][0] * 2 - 1;
- double eta = (std::pow(gamma - gamma_bar, 3.0) +
+ double eta = (Utilities::fixed_power<3>(gamma - gamma_bar) +
gamma_bar * (gamma_bar * gamma_bar + 3)) /
(1 + 3 * gamma_bar * gamma_bar);
const double y_t =
5 * t *
(0.2969 * std::pow(x, 0.5) - 0.126 * x -
- 0.3516 * std::pow(x, 2) + 0.2843 * std::pow(x, 3) -
- 0.1036 * std::pow(x, 4)); // half thickness at a position x
+ 0.3516 * Utilities::fixed_power<2>(x) +
+ 0.2843 * Utilities::fixed_power<3>(x) -
+ 0.1036 * Utilities::fixed_power<4>(
+ x)); // half thickness at a position x
if (is_upper)
naca_points.emplace_back(x, +y_t);
const double x = i * 1 / (1.0 * number_points - 1);
const double y_c =
- (x <= p) ? m / std::pow(p, 2) * (2 * p * x - std::pow(x, 2)) :
- m / std::pow(1 - p, 2) *
- ((1 - 2 * p) + 2 * p * x - std::pow(x, 2));
+ (x <= p) ?
+ m / Utilities::fixed_power<2>(p) *
+ (2 * p * x - Utilities::fixed_power<2>(x)) :
+ m / Utilities::fixed_power<2>(1 - p) *
+ ((1 - 2 * p) + 2 * p * x - Utilities::fixed_power<2>(x));
- const double dy_c = (x <= p) ?
- 2 * m / std::pow(p, 2) * (p - x) :
- 2 * m / std::pow(1 - p, 2) * (p - x);
+ const double dy_c =
+ (x <= p) ? 2 * m / Utilities::fixed_power<2>(p) * (p - x) :
+ 2 * m / Utilities::fixed_power<2>(1 - p) * (p - x);
const double y_t =
5 * t *
(0.2969 * std::pow(x, 0.5) - 0.126 * x -
- 0.3516 * std::pow(x, 2) + 0.2843 * std::pow(x, 3) -
- 0.1036 * std::pow(x, 4)); // half thickness at a position x
+ 0.3516 * Utilities::fixed_power<2>(x) +
+ 0.2843 * Utilities::fixed_power<3>(x) -
+ 0.1036 * Utilities::fixed_power<4>(
+ x)); // half thickness at a position x
const double theta = std::atan(dy_c);
double y = p(2);
double phi = std::atan2(y, x);
double theta = std::atan2(z, std::sqrt(x * x + y * y) - R);
- double w = std::sqrt(std::pow(y - std::sin(phi) * R, 2.0) +
- std::pow(x - std::cos(phi) * R, 2.0) + z * z) /
- r;
+ double w =
+ std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * R) +
+ Utilities::fixed_power<2>(x - std::cos(phi) * R) + z * z) /
+ r;
return {phi, theta, w};
}
// ---------------------------------------------------------------------
#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
#include <deal.II/base/work_stream.h>
#include <deal.II/dofs/dof_accessor.h>
s[2][2] * s[2][2] + 2 * (ss01 + ss02 + ss12)) /
2.;
const double J3 =
- (std::pow(s[0][0], 3) + std::pow(s[1][1], 3) + std::pow(s[2][2], 3) +
- 3. * s[0][0] * (ss01 + ss02) + 3. * s[1][1] * (ss01 + ss12) +
- 3. * s[2][2] * (ss02 + ss12) + 6. * s[0][1] * s[0][2] * s[1][2]) /
+ (Utilities::fixed_power<3>(s[0][0]) +
+ Utilities::fixed_power<3>(s[1][1]) +
+ Utilities::fixed_power<3>(s[2][2]) + 3. * s[0][0] * (ss01 + ss02) +
+ 3. * s[1][1] * (ss01 + ss12) + 3. * s[2][2] * (ss02 + ss12) +
+ 6. * s[0][1] * s[0][2] * s[1][2]) /
3.;
const double R = std::sqrt(4. * J2 / 3.);