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_TABULATED_1D_FUNCTION_HPP
00028 #define OPM_TABULATED_1D_FUNCTION_HPP
00029
00030 #include <opm/material/densead/Math.hpp>
00031 #include <opm/common/ErrorMacros.hpp>
00032 #include <opm/common/Exceptions.hpp>
00033 #include <opm/common/Unused.hpp>
00034
00035 #include <algorithm>
00036 #include <cassert>
00037 #include <iostream>
00038 #include <tuple>
00039 #include <vector>
00040
00041 namespace Opm {
00046 template <class Scalar>
00047 class Tabulated1DFunction
00048 {
00049 public:
00055 Tabulated1DFunction()
00056 {}
00057
00065 template <class ScalarArrayX, class ScalarArrayY>
00066 Tabulated1DFunction(size_t nSamples,
00067 const ScalarArrayX& x,
00068 const ScalarArrayY& y,
00069 bool sortInputs = true)
00070 { this->setXYArrays(nSamples, x, y, sortInputs); }
00071
00080 template <class ScalarContainer>
00081 Tabulated1DFunction(const ScalarContainer& x,
00082 const ScalarContainer& y,
00083 bool sortInputs = true)
00084 { this->setXYContainers(x, y, sortInputs); }
00085
00092 template <class PointContainer>
00093 Tabulated1DFunction(const PointContainer& points,
00094 bool sortInputs = true)
00095 { this->setContainerOfTuples(points, sortInputs); }
00096
00102 template <class ScalarArrayX, class ScalarArrayY>
00103 void setXYArrays(size_t nSamples,
00104 const ScalarArrayX& x,
00105 const ScalarArrayY& y,
00106 bool sortInputs = true)
00107 {
00108 assert(nSamples > 1);
00109
00110 resizeArrays_(nSamples);
00111 for (size_t i = 0; i < nSamples; ++i) {
00112 xValues_[i] = x[i];
00113 yValues_[i] = y[i];
00114 }
00115
00116 if (sortInputs)
00117 sortInput_();
00118 else if (xValues_[0] > xValues_[numSamples() - 1])
00119 reverseSamplingPoints_();
00120 }
00121
00127 template <class ScalarContainerX, class ScalarContainerY>
00128 void setXYContainers(const ScalarContainerX& x,
00129 const ScalarContainerY& y,
00130 bool sortInputs = true)
00131 {
00132 assert(x.size() == y.size());
00133 assert(x.size() > 1);
00134
00135 resizeArrays_(x.size());
00136 std::copy(x.begin(), x.end(), xValues_.begin());
00137 std::copy(y.begin(), y.end(), yValues_.begin());
00138
00139 if (sortInputs)
00140 sortInput_();
00141 else if (xValues_[0] > xValues_[numSamples() - 1])
00142 reverseSamplingPoints_();
00143 }
00144
00148 template <class PointArray>
00149 void setArrayOfPoints(size_t nSamples,
00150 const PointArray& points,
00151 bool sortInputs = true)
00152 {
00153
00154
00155 assert(nSamples > 1);
00156
00157 resizeArrays_(nSamples);
00158 for (size_t i = 0; i < nSamples; ++i) {
00159 xValues_[i] = points[i][0];
00160 yValues_[i] = points[i][1];
00161 }
00162
00163 if (sortInputs)
00164 sortInput_();
00165 else if (xValues_[0] > xValues_[numSamples() - 1])
00166 reverseSamplingPoints_();
00167 }
00168
00183 template <class XYContainer>
00184 void setContainerOfTuples(const XYContainer& points,
00185 bool sortInputs = true)
00186 {
00187
00188
00189 assert(points.size() > 1);
00190
00191 resizeArrays_(points.size());
00192 typename XYContainer::const_iterator it = points.begin();
00193 typename XYContainer::const_iterator endIt = points.end();
00194 for (unsigned i = 0; it != endIt; ++i, ++it) {
00195 xValues_[i] = std::get<0>(*it);
00196 yValues_[i] = std::get<1>(*it);
00197 }
00198
00199 if (sortInputs)
00200 sortInput_();
00201 else if (xValues_[0] > xValues_[numSamples() - 1])
00202 reverseSamplingPoints_();
00203 }
00204
00208 size_t numSamples() const
00209 { return xValues_.size(); }
00210
00214 Scalar xMin() const
00215 { return xValues_[0]; }
00216
00220 Scalar xMax() const
00221 { return xValues_[numSamples() - 1]; }
00222
00226 Scalar xAt(size_t i) const
00227 { return xValues_[i]; }
00228
00232 Scalar valueAt(size_t i) const
00233 { return yValues_[i]; }
00234
00238 template <class Evaluation>
00239 bool applies(const Evaluation& x) const
00240 { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
00241
00251 template <class Evaluation>
00252 Evaluation eval(const Evaluation& x, bool extrapolate = false) const
00253 {
00254 if (!extrapolate && !applies(x))
00255 OPM_THROW(Opm::NumericalProblem,
00256 "Tried to evaluate a tabulated function outside of its range");
00257
00258 size_t segIdx;
00259
00260 if (extrapolate && x < xValues_.front())
00261 segIdx = 0;
00262 else if (extrapolate && x > xValues_.back())
00263 segIdx = numSamples() - 2;
00264 else
00265 segIdx = findSegmentIndex_(Opm::scalarValue(x));
00266
00267 Scalar x0 = xValues_[segIdx];
00268 Scalar x1 = xValues_[segIdx + 1];
00269
00270 Scalar y0 = yValues_[segIdx];
00271 Scalar y1 = yValues_[segIdx + 1];
00272
00273 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
00274 }
00275
00287 template <class Evaluation>
00288 Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
00289 {
00290 if (!extrapolate && !applies(x)) {
00291 OPM_THROW(Opm::NumericalProblem,
00292 "Tried to evaluate a derivative of a tabulated"
00293 " function outside of its range");
00294 }
00295
00296 unsigned segIdx = findSegmentIndex_(Opm::scalarValue(x));
00297 return evalDerivative_(x, segIdx);
00298 }
00299
00314 template <class Evaluation>
00315 Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
00316 {
00317 if (!extrapolate && !applies(x)) {
00318 OPM_THROW(Opm::NumericalProblem,
00319 "Tried to evaluate a second derivative of a tabulated "
00320 " function outside of its range");
00321 }
00322
00323 return 0.0;
00324 }
00325
00340 template <class Evaluation>
00341 Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
00342 {
00343 if (!extrapolate && !applies(x)) {
00344 OPM_THROW(Opm::NumericalProblem,
00345 "Tried to evaluate a third derivative of a tabulated "
00346 " function outside of its range");
00347 }
00348
00349 return 0.0;
00350 }
00351
00360 int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED = false) const
00361 {
00362 assert(x0 != x1);
00363
00364
00365 if (x0 > x1)
00366 std::swap(x0, x1);
00367
00368 assert(x0 < x1);
00369
00370 int r = 3;
00371 if (x0 < xMin()) {
00372 assert(extrapolate);
00373
00374 x0 = xMin();
00375 };
00376
00377 size_t i = findSegmentIndex_(x0);
00378 if (xValues_[i + 1] >= x1) {
00379
00380
00381 updateMonotonicity_(i, r);
00382 return r;
00383 }
00384
00385
00386
00387 updateMonotonicity_(i, r);
00388 ++ i;
00389
00390
00391
00392 size_t iEnd = findSegmentIndex_(x1);
00393 for (; i < iEnd - 1; ++i) {
00394 updateMonotonicity_(i, r);
00395 if (!r)
00396 return 0;
00397 }
00398
00399
00400
00401
00402 if (x1 > xMax()) {
00403 assert(extrapolate);
00404
00405 Scalar m = evalDerivative_(xMax(), numSamples() - 2);
00406 if (m < 0)
00407 return (r < 0 || r==3)?-1:0;
00408 else if (m > 0)
00409 return (r > 0 || r==3)?1:0;
00410
00411 return r;
00412 }
00413
00414
00415 updateMonotonicity_(iEnd, r);
00416
00417 return r;
00418 }
00419
00424 int monotonic() const
00425 { return monotonic(xMin(), xMax()); }
00426
00443 void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os = std::cout) const
00444 {
00445 Scalar x0 = std::min(xi0, xi1);
00446 Scalar x1 = std::max(xi0, xi1);
00447 const int n = numSamples() - 1;
00448 for (int i = 0; i <= k; ++i) {
00449 double x = i*(x1 - x0)/k + x0;
00450 double y;
00451 double dy_dx;
00452 if (!applies(x)) {
00453 if (x < xValues_[0]) {
00454 dy_dx = evalDerivative(xValues_[0]);
00455 y = (x - xValues_[0])*dy_dx + yValues_[0];
00456 }
00457 else if (x > xValues_[n]) {
00458 dy_dx = evalDerivative(xValues_[n]);
00459 y = (x - xValues_[n])*dy_dx + yValues_[n];
00460 }
00461 else {
00462 OPM_THROW(std::runtime_error,
00463 "The sampling points given to a function must be sorted by their x value!");
00464 }
00465 }
00466 else {
00467 y = eval(x);
00468 dy_dx = evalDerivative(x);
00469 }
00470
00471 os << x << " " << y << " " << dy_dx << "\n";
00472 }
00473 }
00474
00475 private:
00476 size_t findSegmentIndex_(Scalar x) const
00477 {
00478
00479 assert(xValues_.size() >= 2);
00480
00481 if (x <= xValues_[1])
00482 return 0;
00483 else if (x >= xValues_[xValues_.size() - 2])
00484 return xValues_.size() - 2;
00485 else {
00486
00487 size_t segmentIdx = 1;
00488 size_t upperIdx = xValues_.size() - 2;
00489 while (segmentIdx + 1 < upperIdx) {
00490 size_t pivotIdx = (segmentIdx + upperIdx) / 2;
00491 if (x < xValues_[pivotIdx])
00492 upperIdx = pivotIdx;
00493 else
00494 segmentIdx = pivotIdx;
00495 }
00496
00497 assert(xValues_[segmentIdx] <= x);
00498 assert(x <= xValues_[segmentIdx + 1]);
00499 return segmentIdx;
00500 }
00501 }
00502
00503 template <class Evaluation>
00504 Evaluation evalDerivative_(const Evaluation& x OPM_UNUSED, size_t segIdx) const
00505 {
00506 Scalar x0 = xValues_[segIdx];
00507 Scalar x1 = xValues_[segIdx + 1];
00508
00509 Scalar y0 = yValues_[segIdx];
00510 Scalar y1 = yValues_[segIdx + 1];
00511
00512 return (y1 - y0)/(x1 - x0);
00513 }
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 int updateMonotonicity_(size_t i, int& r) const
00524 {
00525 if (yValues_[i] < yValues_[i + 1]) {
00526
00527 if (r == 3 || r == 1)
00528 r = 1;
00529 else
00530 r = 0;
00531 return 1;
00532 }
00533 else if (yValues_[i] > yValues_[i + 1]) {
00534
00535 if (r == 3 || r == -1)
00536 r = -1;
00537 else
00538 r = 0;
00539 return -1;
00540 }
00541
00542 return 3;
00543 }
00544
00548 struct ComparatorX_
00549 {
00550 ComparatorX_(const std::vector<Scalar>& x)
00551 : x_(x)
00552 {}
00553
00554 bool operator ()(size_t idxA, size_t idxB) const
00555 { return x_.at(idxA) < x_.at(idxB); }
00556
00557 const std::vector<Scalar>& x_;
00558 };
00559
00563 void sortInput_()
00564 {
00565 size_t n = numSamples();
00566
00567
00568 std::vector<unsigned> idxVector(n);
00569 for (unsigned i = 0; i < n; ++i)
00570 idxVector[i] = i;
00571
00572
00573
00574 ComparatorX_ cmp(xValues_);
00575 std::sort(idxVector.begin(), idxVector.end(), cmp);
00576
00577
00578 std::vector<Scalar> tmpX(n), tmpY(n);
00579 for (size_t i = 0; i < idxVector.size(); ++ i) {
00580 tmpX[i] = xValues_[idxVector[i]];
00581 tmpY[i] = yValues_[idxVector[i]];
00582 }
00583 xValues_ = tmpX;
00584 yValues_ = tmpY;
00585 }
00586
00591 void reverseSamplingPoints_()
00592 {
00593
00594 size_t n = numSamples();
00595 for (size_t i = 0; i <= (n - 1)/2; ++i) {
00596 std::swap(xValues_[i], xValues_[n - i - 1]);
00597 std::swap(yValues_[i], yValues_[n - i - 1]);
00598 }
00599 }
00600
00604 void resizeArrays_(size_t nSamples)
00605 {
00606 xValues_.resize(nSamples);
00607 yValues_.resize(nSamples);
00608 }
00609
00610 std::vector<Scalar> xValues_;
00611 std::vector<Scalar> yValues_;
00612 };
00613 }
00614
00615 #endif