00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00028 #ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
00029 #define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
00030
00031 #include <opm/common/Valgrind.hpp>
00032 #include <opm/common/Exceptions.hpp>
00033 #include <opm/common/ErrorMacros.hpp>
00034 #include <opm/common/Unused.hpp>
00035 #include <opm/material/common/MathToolbox.hpp>
00036
00037 #include <iostream>
00038 #include <vector>
00039 #include <limits>
00040 #include <tuple>
00041
00042 #include <assert.h>
00043
00044 namespace Opm {
00045
00054 template <class Scalar>
00055 class UniformXTabulated2DFunction
00056 {
00057 typedef std::tuple<Scalar, Scalar, Scalar> SamplePoint;
00058
00059 public:
00060 UniformXTabulated2DFunction()
00061 { }
00062
00066 Scalar xMin() const
00067 { return xPos_.front(); }
00068
00072 Scalar xMax() const
00073 { return xPos_.back(); }
00074
00078 Scalar xAt(size_t i) const
00079 { return xPos_[i]; }
00080
00084 Scalar yAt(size_t i, size_t j) const
00085 { return std::get<1>(samples_[i][j]); }
00086
00090 Scalar valueAt(size_t i, size_t j) const
00091 { return std::get<2>(samples_[i][j]); }
00092
00096 size_t numX() const
00097 { return xPos_.size(); }
00098
00102 Scalar yMin(size_t i) const
00103 { return std::get<1>(samples_.at(i).front()); }
00104
00108 Scalar yMax(size_t i) const
00109 { return std::get<1>(samples_.at(i).back()); }
00110
00114 size_t numY(size_t i) const
00115 { return samples_.at(i).size(); }
00116
00120 Scalar iToX(size_t i) const
00121 {
00122 assert(0 <= i && i < numX());
00123
00124 return xPos_.at(i);
00125 }
00126
00130 Scalar jToY(size_t i, size_t j) const
00131 {
00132 assert(0 <= i && i < numX());
00133 assert(0 <= j && size_t(j) < samples_[i].size());
00134
00135 return std::get<1>(samples_.at(i).at(j));
00136 }
00137
00146 template <class Evaluation>
00147 Evaluation xToI(const Evaluation& x, bool extrapolate OPM_OPTIM_UNUSED = false) const
00148 {
00149 assert(extrapolate || (xMin() <= x && x <= xMax()));
00150
00151
00152 assert(xPos_.size() >= 2);
00153
00154 size_t segmentIdx;
00155 if (x <= xPos_[1])
00156 segmentIdx = 0;
00157 else if (x >= xPos_[xPos_.size() - 2])
00158 segmentIdx = xPos_.size() - 2;
00159 else {
00160
00161 segmentIdx = 1;
00162 size_t upperIdx = xPos_.size() - 2;
00163 while (segmentIdx + 1 < upperIdx) {
00164 size_t pivotIdx = (segmentIdx + upperIdx) / 2;
00165 if (x < xPos_[pivotIdx])
00166 upperIdx = pivotIdx;
00167 else
00168 segmentIdx = pivotIdx;
00169 }
00170
00171 assert(xPos_[segmentIdx] <= x);
00172 assert(x <= xPos_[segmentIdx + 1]);
00173 }
00174
00175 Scalar x1 = xPos_[segmentIdx];
00176 Scalar x2 = xPos_[segmentIdx + 1];
00177 return Scalar(segmentIdx) + (x - x1)/(x2 - x1);
00178 }
00179
00188 template <class Evaluation>
00189 Evaluation yToJ(size_t i, const Evaluation& y, bool extrapolate OPM_OPTIM_UNUSED = false) const
00190 {
00191 assert(0 <= i && i < numX());
00192 const auto& colSamplePoints = samples_.at(i);
00193
00194 assert(extrapolate || (yMin(i) <= y && y <= yMax(i)));
00195
00196 Scalar y1;
00197 Scalar y2;
00198
00199
00200 size_t lowerIdx = 0;
00201 size_t upperIdx = colSamplePoints.size() - 1;
00202 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
00203 while (lowerIdx + 1 < upperIdx) {
00204 if (y < std::get<1>(colSamplePoints[pivotIdx]))
00205 upperIdx = pivotIdx;
00206 else
00207 lowerIdx = pivotIdx;
00208 pivotIdx = (lowerIdx + upperIdx) / 2;
00209 }
00210
00211 y1 = std::get<1>(colSamplePoints[lowerIdx]);
00212 y2 = std::get<1>(colSamplePoints[lowerIdx + 1]);
00213
00214 assert(y1 <= y || (extrapolate && lowerIdx == 0));
00215 assert(y <= y2 || (extrapolate && lowerIdx == colSamplePoints.size() - 2));
00216
00217 return Scalar(lowerIdx) + (y - y1)/(y2 - y1);
00218 }
00219
00223 bool applies(Scalar x, Scalar y) const
00224 {
00225 if (x < xMin() || xMax() < x)
00226 return false;
00227
00228 Scalar i = xToI(x, false);
00229 const auto& col1SamplePoints = samples_.at(unsigned(i));
00230 const auto& col2SamplePoints = samples_.at(unsigned(i));
00231 Scalar alpha = i - static_cast<int>(i);
00232
00233 Scalar minY =
00234 alpha*std::get<1>(col1SamplePoints.front()) +
00235 (1 - alpha)*std::get<1>(col2SamplePoints.front());
00236
00237 Scalar maxY =
00238 alpha*std::get<1>(col1SamplePoints.back()) +
00239 (1 - alpha)*std::get<1>(col2SamplePoints.back());
00240
00241 return minY <= y && y <= maxY;
00242 }
00249 template <class Evaluation>
00250 Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
00251 {
00252 #ifndef NDEBUG
00253 if (!extrapolate && !applies(Opm::scalarValue(x), Opm::scalarValue(y))) {
00254 OPM_THROW(NumericalProblem,
00255 "Attempt to get undefined table value (" << x << ", " << y << ")");
00256 };
00257 #endif
00258
00259
00260
00261 Evaluation alpha = xToI(x, extrapolate);
00262 size_t i =
00263 static_cast<size_t>(std::max(0, std::min(static_cast<int>(numX()) - 2,
00264 static_cast<int>(Opm::scalarValue(alpha)))));
00265 alpha -= i;
00266
00267 Evaluation beta1;
00268 Evaluation beta2;
00269
00270 beta1 = yToJ(i, y, extrapolate);
00271 beta2 = yToJ(i + 1, y, extrapolate);
00272
00273 size_t j1 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i)) - 2,
00274 static_cast<int>(Opm::scalarValue(beta1)))));
00275 size_t j2 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i + 1)) - 2,
00276 static_cast<int>(Opm::scalarValue(beta2)))));
00277
00278 beta1 -= j1;
00279 beta2 -= j2;
00280
00281
00282 Evaluation s1, s2;
00283 s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
00284 s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
00285
00286 Valgrind::CheckDefined(s1);
00287 Valgrind::CheckDefined(s2);
00288
00289
00290 Evaluation result;
00291 result = s1*(1.0 - alpha) + s2*alpha;
00292 Valgrind::CheckDefined(result);
00293
00294 return result;
00295 }
00296
00302 size_t appendXPos(Scalar nextX)
00303 {
00304 if (xPos_.empty() || xPos_.back() < nextX) {
00305 xPos_.push_back(nextX);
00306 samples_.resize(xPos_.size());
00307 return xPos_.size() - 1;
00308 }
00309 else if (xPos_.front() > nextX) {
00310
00311 xPos_.insert(xPos_.begin(), nextX);
00312 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
00313 return 0;
00314 }
00315 OPM_THROW(std::invalid_argument,
00316 "Sampling points should be specified either monotonically "
00317 "ascending or descending.");
00318 }
00319
00325 size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
00326 {
00327 assert(0 <= i && i < numX());
00328
00329 Scalar x = iToX(i);
00330 if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
00331 samples_[i].push_back(SamplePoint(x, y, value));
00332 return samples_[i].size() - 1;
00333 }
00334 else if (std::get<1>(samples_[i].front()) > y) {
00335
00336 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
00337 return 0;
00338 }
00339
00340 OPM_THROW(std::invalid_argument,
00341 "Sampling points must be specified in either monotonically "
00342 "ascending or descending order.");
00343 }
00344
00351 void print(std::ostream& os = std::cout) const
00352 {
00353 Scalar x0 = xMin();
00354 Scalar x1 = xMax();
00355 int m = numX();
00356
00357 Scalar y0 = 1e30;
00358 Scalar y1 = -1e30;
00359 int n = 0;
00360 for (int i = 0; i < m; ++ i) {
00361 y0 = std::min(y0, yMin(i));
00362 y1 = std::max(y1, yMax(i));
00363 n = std::max(n, numY(i));
00364 }
00365
00366 m *= 3;
00367 n *= 3;
00368 for (int i = 0; i <= m; ++i) {
00369 Scalar x = x0 + (x1 - x0)*i/m;
00370 for (int j = 0; j <= n; ++j) {
00371 Scalar y = y0 + (y1 - y0)*j/n;
00372 os << x << " " << y << " " << eval(x, y) << "\n";
00373 }
00374 os << "\n";
00375 }
00376 }
00377
00378 private:
00379
00380
00381
00382 std::vector<std::vector<SamplePoint> > samples_;
00383
00384
00385 std::vector<Scalar> xPos_;
00386 };
00387 }
00388
00389 #endif