27 #ifndef OPM_POLYNOMIAL_UTILS_HH 28 #define OPM_POLYNOMIAL_UTILS_HH 50 template <
class Scalar,
class SolContainer>
55 if (std::abs(Opm::scalarValue(a)) < 1e-30)
78 template <
class Scalar,
class SolContainer>
85 if (std::abs(Opm::scalarValue(a)) < 1e-30)
89 Scalar Delta = b*b - 4*a*c;
93 Delta = Opm::sqrt(Delta);
94 sol[0] = (- b + Delta)/(2*a);
95 sol[1] = (- b - Delta)/(2*a);
99 std::swap(sol[0], sol[1]);
104 template <
class Scalar,
class SolContainer>
105 void invertCubicPolynomialPostProcess_(SolContainer& sol,
114 for (
int i = 0; i < numSol; ++i) {
116 Scalar fOld = d + x*(c + x*(b + x*a));
118 Scalar fPrime = c + x*(2*b + x*3*a);
119 if (std::abs(Opm::scalarValue(fPrime)) < 1e-30)
123 Scalar fNew = d + x*(c + x*(b + x*a));
124 if (std::abs(Opm::scalarValue(fNew)) < std::abs(Opm::scalarValue(fOld)))
147 template <
class Scalar,
class SolContainer>
155 if (std::abs(Opm::scalarValue(a)) < 1e-30)
165 Scalar p = c - b*b/3;
166 Scalar q = d + (2*b*b*b - 9*b*c)/27;
168 if (std::abs(Opm::scalarValue(p)) > 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30) {
199 Scalar wDisc = q*q/4 + p*p*p/27;
202 Scalar u = - q/2 + Opm::sqrt(wDisc);
203 if (u < 0) u = - Opm::pow(-u, 1.0/3);
204 else u = Opm::pow(u, 1.0/3);
208 sol[0] = u - p/(3*u) - b/3;
212 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
216 Scalar uCubedRe = - q/2;
217 Scalar uCubedIm = Opm::sqrt(-wDisc);
219 Scalar uAbs = Opm::pow(Opm::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
220 Scalar phi = Opm::atan2(uCubedIm, uCubedRe)/3;
261 for (
int i = 0; i < 3; ++i) {
262 sol[i] = Opm::cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
268 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
271 std::sort(sol, sol + 3);
278 else if (std::abs(Opm::scalarValue(p)) < 1e-30 && std::abs(Opm::scalarValue(q)) < 1e-30) {
280 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
283 else if (std::abs(Opm::scalarValue(p)) < 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30) {
288 if (-q > 0) t = Opm::pow(-q, 1./3);
289 else t = - Opm::pow(q, 1./3);
295 assert(std::abs(Opm::scalarValue(p)) > 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30);
306 sol[0] = -Opm::sqrt(-p) - b/3;;
308 sol[2] = Opm::sqrt(-p) - b/3;
Definition: Air_Mesitylene.hpp:33
unsigned invertLinearPolynomial(SolContainer &sol, Scalar a, Scalar b)
Invert a linear polynomial analytically.
Definition: PolynomialUtils.hpp:51
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
unsigned invertQuadraticPolynomial(SolContainer &sol, Scalar a, Scalar b, Scalar c)
Invert a quadratic polynomial analytically.
Definition: PolynomialUtils.hpp:79