00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_POLYNOMIAL_UTILS_HH
00028 #define OPM_POLYNOMIAL_UTILS_HH
00029
00030 #include <cmath>
00031 #include <algorithm>
00032
00033 #include <opm/material/common/MathToolbox.hpp>
00034
00035 namespace Opm {
00050 template <class Scalar, class SolContainer>
00051 unsigned invertLinearPolynomial(SolContainer& sol,
00052 Scalar a,
00053 Scalar b)
00054 {
00055 if (std::abs(Opm::scalarValue(a)) < 1e-30)
00056 return 0;
00057
00058 sol[0] = -b/a;
00059 return 1;
00060 }
00061
00078 template <class Scalar, class SolContainer>
00079 unsigned invertQuadraticPolynomial(SolContainer& sol,
00080 Scalar a,
00081 Scalar b,
00082 Scalar c)
00083 {
00084
00085 if (std::abs(Opm::scalarValue(a)) < 1e-30)
00086 return invertLinearPolynomial(sol, b, c);
00087
00088
00089 Scalar Delta = b*b - 4*a*c;
00090 if (Delta < 0)
00091 return 0;
00092
00093 Delta = Opm::sqrt(Delta);
00094 sol[0] = (- b + Delta)/(2*a);
00095 sol[1] = (- b - Delta)/(2*a);
00096
00097
00098 if (sol[0] > sol[1])
00099 std::swap(sol[0], sol[1]);
00100 return 2;
00101 }
00102
00104 template <class Scalar, class SolContainer>
00105 void invertCubicPolynomialPostProcess_(SolContainer& sol,
00106 int numSol,
00107 Scalar a,
00108 Scalar b,
00109 Scalar c,
00110 Scalar d)
00111 {
00112
00113
00114 for (int i = 0; i < numSol; ++i) {
00115 Scalar x = sol[i];
00116 Scalar fOld = d + x*(c + x*(b + x*a));
00117
00118 Scalar fPrime = c + x*(2*b + x*3*a);
00119 if (std::abs(Opm::scalarValue(fPrime)) < 1e-30)
00120 continue;
00121 x -= fOld/fPrime;
00122
00123 Scalar fNew = d + x*(c + x*(b + x*a));
00124 if (std::abs(Opm::scalarValue(fNew)) < std::abs(Opm::scalarValue(fOld)))
00125 sol[i] = x;
00126 }
00127 }
00129
00147 template <class Scalar, class SolContainer>
00148 unsigned invertCubicPolynomial(SolContainer* sol,
00149 Scalar a,
00150 Scalar b,
00151 Scalar c,
00152 Scalar d)
00153 {
00154
00155 if (std::abs(Opm::scalarValue(a)) < 1e-30)
00156 return invertQuadraticPolynomial(sol, b, c, d);
00157
00158
00159 b /= a;
00160 c /= a;
00161 d /= a;
00162 a = 1;
00163
00164
00165 Scalar p = c - b*b/3;
00166 Scalar q = d + (2*b*b*b - 9*b*c)/27;
00167
00168 if (std::abs(Opm::scalarValue(p)) > 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30) {
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 Scalar wDisc = q*q/4 + p*p*p/27;
00200 if (wDisc >= 0) {
00201
00202 Scalar u = - q/2 + Opm::sqrt(wDisc);
00203 if (u < 0) u = - Opm::pow(-u, 1.0/3);
00204 else u = Opm::pow(u, 1.0/3);
00205
00206
00207
00208 sol[0] = u - p/(3*u) - b/3;
00209
00210
00211
00212 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
00213 return 1;
00214 }
00215 else {
00216 Scalar uCubedRe = - q/2;
00217 Scalar uCubedIm = Opm::sqrt(-wDisc);
00218
00219 Scalar uAbs = Opm::pow(Opm::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
00220 Scalar phi = Opm::atan2(uCubedIm, uCubedRe)/3;
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 for (int i = 0; i < 3; ++i) {
00262 sol[i] = Opm::cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
00263 phi += 2*M_PI/3;
00264 }
00265
00266
00267
00268 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
00269
00270
00271 std::sort(sol, sol + 3);
00272
00273 return 3;
00274 }
00275 }
00276
00277
00278 else if (std::abs(Opm::scalarValue(p)) < 1e-30 && std::abs(Opm::scalarValue(q)) < 1e-30) {
00279
00280 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
00281 return 3;
00282 }
00283 else if (std::abs(Opm::scalarValue(p)) < 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30) {
00284
00285
00286
00287 Scalar t;
00288 if (-q > 0) t = Opm::pow(-q, 1./3);
00289 else t = - Opm::pow(q, 1./3);
00290 sol[0] = t - b/3;
00291
00292 return 1;
00293 }
00294
00295 assert(std::abs(Opm::scalarValue(p)) > 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30);
00296
00297
00298
00299
00300 if (p > 0) {
00301 sol[0] = 0.0 - b/3;
00302 return 1;
00303 }
00304
00305
00306 sol[0] = -Opm::sqrt(-p) - b/3;;
00307 sol[1] = 0.0 - b/3;
00308 sol[2] = Opm::sqrt(-p) - b/3;
00309
00310 return 3;
00311 }
00312 }
00313
00314 #endif