27 #ifndef OPM_POLYNOMIAL_UTILS_HH
28 #define OPM_POLYNOMIAL_UTILS_HH
50 template <
class Scalar,
class SolContainer>
51 unsigned invertLinearPolynomial(SolContainer& sol,
55 if (std::abs(Opm::scalarValue(a)) < 1e-30)
78 template <
class Scalar,
class SolContainer>
79 unsigned invertQuadraticPolynomial(SolContainer& sol,
85 if (std::abs(Opm::scalarValue(a)) < 1e-30)
86 return invertLinearPolynomial(sol, b, c);
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>
148 unsigned invertCubicPolynomial(SolContainer* sol,
155 if (std::abs(Opm::scalarValue(a)) < 1e-30)
156 return invertQuadraticPolynomial(sol, b, c, d);
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;