25 #ifndef OPM_POLYNOMIAL_UTILS_HH
26 #define OPM_POLYNOMIAL_UTILS_HH
48 template <
class Scalar,
class SolContainer>
54 if (
std::abs(Toolbox::value(a)) < 1e-30)
77 template <
class Scalar,
class SolContainer>
86 if (
std::abs(Toolbox::value(a)) < 1e-30)
90 Scalar Delta = b*b - 4*a*c;
95 sol[0] = (- b + Delta)/(2*a);
96 sol[1] = (- b - Delta)/(2*a);
100 std::swap(sol[0], sol[1]);
105 template <
class Scalar,
class SolContainer>
106 void invertCubicPolynomialPostProcess_(SolContainer &sol,
113 typedef MathToolbox<Scalar> Toolbox;
117 for (
int i = 0; i < numSol; ++i) {
119 Scalar fOld = d + x*(c + x*(b + x*a));
121 Scalar fPrime = c + x*(2*b + x*3*a);
122 if (
std::abs(Toolbox::value(fPrime)) < 1e-30)
126 Scalar fNew = d + x*(c + x*(b + x*a));
150 template <
class Scalar,
class SolContainer>
160 if (
std::abs(Toolbox::value(a)) < 1e-30)
170 Scalar p = c - b*b/3;
171 Scalar q = d + (2*b*b*b - 9*b*c)/27;
173 if (
std::abs(Toolbox::value(p)) > 1e-30 &&
std::abs(Toolbox::value(q)) > 1e-30) {
204 Scalar wDisc = q*q/4 + p*p*p/27;
213 sol[0] = u - p/(3*u) - b/3;
217 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
221 Scalar uCubedRe = - q/2;
266 for (
int i = 0; i < 3; ++i) {
273 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
276 std::sort(sol, sol + 3);
283 else if (
std::abs(Toolbox::value(p)) < 1e-30 &&
std::abs(Toolbox::value(q)) < 1e-30) {
285 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
288 else if (
std::abs(Toolbox::value(p)) < 1e-30 &&
std::abs(Toolbox::value(q)) > 1e-30) {
300 assert(
std::abs(Toolbox::value(p)) > 1e-30 &&
std::abs(Toolbox::value(q)) > 1e-30);
Definition: Air_Mesitylene.hpp:31
Evaluation< Scalar, VarSetTag, numVars > sqrt(const Evaluation< Scalar, VarSetTag, numVars > &x)
Definition: Math.hpp:278
Evaluation< Scalar, VarSetTag, numVars > cos(const Evaluation< Scalar, VarSetTag, numVars > &x)
Definition: Math.hpp:248
Evaluation< Scalar, VarSetTag, numVars > atan2(const Evaluation< Scalar, VarSetTag, numVars > &x, const Evaluation< Scalar, VarSetTag, numVars > &y)
Definition: Math.hpp:198
int invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:151
int invertLinearPolynomial(SolContainer &sol, Scalar a, Scalar b)
Invert a linear polynomial analytically.
Definition: PolynomialUtils.hpp:49
Evaluation< Scalar, VarSetTag, numVars > abs(const Evaluation< Scalar, VarSetTag, numVars > &)
Definition: Math.hpp:41
Evaluation< Scalar, VarSetTag, numVars > pow(const Evaluation< Scalar, VarSetTag, numVars > &base, Scalar exp)
Definition: Math.hpp:312
int invertQuadraticPolynomial(SolContainer &sol, Scalar a, Scalar b, Scalar c)
Invert a quadratic polynomial analytically.
Definition: PolynomialUtils.hpp:78