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_SPLINE_HPP
00028 #define OPM_SPLINE_HPP
00029
00030 #include <opm/material/common/TridiagonalMatrix.hpp>
00031 #include <opm/material/common/PolynomialUtils.hpp>
00032 #include <opm/common/ErrorMacros.hpp>
00033 #include <opm/common/Exceptions.hpp>
00034 #include <opm/common/Unused.hpp>
00035
00036 #include <ostream>
00037 #include <vector>
00038 #include <tuple>
00039
00040 namespace Opm
00041 {
00091 template<class Scalar>
00092 class Spline
00093 {
00094 typedef Opm::TridiagonalMatrix<Scalar> Matrix;
00095 typedef std::vector<Scalar> Vector;
00096
00097 public:
00103 enum SplineType {
00104 Full,
00105 Natural,
00106 Periodic,
00107 Monotonic
00108 };
00109
00115 Spline()
00116 { }
00117
00128 Spline(Scalar x0, Scalar x1,
00129 Scalar y0, Scalar y1,
00130 Scalar m0, Scalar m1)
00131 { set(x0, x1, y0, y1, m0, m1); }
00132
00141 template <class ScalarArrayX, class ScalarArrayY>
00142 Spline(size_t nSamples,
00143 const ScalarArrayX& x,
00144 const ScalarArrayY& y,
00145 SplineType splineType = Natural,
00146 bool sortInputs = true)
00147 { this->setXYArrays(nSamples, x, y, splineType, sortInputs); }
00148
00156 template <class PointArray>
00157 Spline(size_t nSamples,
00158 const PointArray& points,
00159 SplineType splineType = Natural,
00160 bool sortInputs = true)
00161 { this->setArrayOfPoints(nSamples, points, splineType, sortInputs); }
00162
00170 template <class ScalarContainer>
00171 Spline(const ScalarContainer& x,
00172 const ScalarContainer& y,
00173 SplineType splineType = Natural,
00174 bool sortInputs = true)
00175 { this->setXYContainers(x, y, splineType, sortInputs); }
00176
00183 template <class PointContainer>
00184 Spline(const PointContainer& points,
00185 SplineType splineType = Natural,
00186 bool sortInputs = true)
00187 { this->setContainerOfPoints(points, splineType, sortInputs); }
00188
00199 template <class ScalarArray>
00200 Spline(size_t nSamples,
00201 const ScalarArray& x,
00202 const ScalarArray& y,
00203 Scalar m0,
00204 Scalar m1,
00205 bool sortInputs = true)
00206 { this->setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
00207
00217 template <class PointArray>
00218 Spline(size_t nSamples,
00219 const PointArray& points,
00220 Scalar m0,
00221 Scalar m1,
00222 bool sortInputs = true)
00223 { this->setArrayOfPoints(nSamples, points, m0, m1, sortInputs); }
00224
00234 template <class ScalarContainerX, class ScalarContainerY>
00235 Spline(const ScalarContainerX& x,
00236 const ScalarContainerY& y,
00237 Scalar m0,
00238 Scalar m1,
00239 bool sortInputs = true)
00240 { this->setXYContainers(x, y, m0, m1, sortInputs); }
00241
00250 template <class PointContainer>
00251 Spline(const PointContainer& points,
00252 Scalar m0,
00253 Scalar m1,
00254 bool sortInputs = true)
00255 { this->setContainerOfPoints(points, m0, m1, sortInputs); }
00256
00268 void set(Scalar x0, Scalar x1,
00269 Scalar y0, Scalar y1,
00270 Scalar m0, Scalar m1)
00271 {
00272 slopeVec_.resize(2);
00273 xPos_.resize(2);
00274 yPos_.resize(2);
00275
00276 if (x0 > x1) {
00277 xPos_[0] = x1;
00278 xPos_[1] = x0;
00279 yPos_[0] = y1;
00280 yPos_[1] = y0;
00281 }
00282 else {
00283 xPos_[0] = x0;
00284 xPos_[1] = x1;
00285 yPos_[0] = y0;
00286 yPos_[1] = y1;
00287 }
00288
00289 slopeVec_[0] = m0;
00290 slopeVec_[1] = m1;
00291
00292 Matrix M(numSamples());
00293 Vector d(numSamples());
00294 Vector moments(numSamples());
00295 this->makeFullSystem_(M, d, m0, m1);
00296
00297
00298 M.solve(moments, d);
00299
00300 this->setSlopesFromMoments_(slopeVec_, moments);
00301 }
00302
00303
00307
00311
00323 template <class ScalarArrayX, class ScalarArrayY>
00324 void setXYArrays(size_t nSamples,
00325 const ScalarArrayX& x,
00326 const ScalarArrayY& y,
00327 Scalar m0, Scalar m1,
00328 bool sortInputs = true)
00329 {
00330 assert(nSamples > 1);
00331
00332 setNumSamples_(nSamples);
00333 for (size_t i = 0; i < nSamples; ++i) {
00334 xPos_[i] = x[i];
00335 yPos_[i] = y[i];
00336 }
00337
00338 if (sortInputs)
00339 sortInput_();
00340 else if (xPos_[0] > xPos_[numSamples() - 1])
00341 reverseSamplingPoints_();
00342
00343 makeFullSpline_(m0, m1);
00344 }
00345
00357 template <class ScalarContainerX, class ScalarContainerY>
00358 void setXYContainers(const ScalarContainerX& x,
00359 const ScalarContainerY& y,
00360 Scalar m0, Scalar m1,
00361 bool sortInputs = true)
00362 {
00363 assert(x.size() == y.size());
00364 assert(x.size() > 1);
00365
00366 setNumSamples_(x.size());
00367
00368 std::copy(x.begin(), x.end(), xPos_.begin());
00369 std::copy(y.begin(), y.end(), yPos_.begin());
00370
00371 if (sortInputs)
00372 sortInput_();
00373 else if (xPos_[0] > xPos_[numSamples() - 1])
00374 reverseSamplingPoints_();
00375
00376 makeFullSpline_(m0, m1);
00377 }
00378
00391 template <class PointArray>
00392 void setArrayOfPoints(size_t nSamples,
00393 const PointArray& points,
00394 Scalar m0,
00395 Scalar m1,
00396 bool sortInputs = true)
00397 {
00398
00399
00400 assert(nSamples > 1);
00401
00402 setNumSamples_(nSamples);
00403 for (size_t i = 0; i < nSamples; ++i) {
00404 xPos_[i] = points[i][0];
00405 yPos_[i] = points[i][1];
00406 }
00407
00408 if (sortInputs)
00409 sortInput_();
00410 else if (xPos_[0] > xPos_[numSamples() - 1])
00411 reverseSamplingPoints_();
00412
00413 makeFullSpline_(m0, m1);
00414 }
00415
00428 template <class XYContainer>
00429 void setContainerOfPoints(const XYContainer& points,
00430 Scalar m0,
00431 Scalar m1,
00432 bool sortInputs = true)
00433 {
00434
00435
00436 assert(points.size() > 1);
00437
00438 setNumSamples_(points.size());
00439 typename XYContainer::const_iterator it = points.begin();
00440 typename XYContainer::const_iterator endIt = points.end();
00441 for (size_t i = 0; it != endIt; ++i, ++it) {
00442 xPos_[i] = (*it)[0];
00443 yPos_[i] = (*it)[1];
00444 }
00445
00446 if (sortInputs)
00447 sortInput_();
00448 else if (xPos_[0] > xPos_[numSamples() - 1])
00449 reverseSamplingPoints_();
00450
00451
00452 makeFullSpline_(m0, m1);
00453 }
00454
00470 template <class XYContainer>
00471 void setContainerOfTuples(const XYContainer& points,
00472 Scalar m0,
00473 Scalar m1,
00474 bool sortInputs = true)
00475 {
00476
00477 setNumSamples_(points.size());
00478 typename XYContainer::const_iterator it = points.begin();
00479 typename XYContainer::const_iterator endIt = points.end();
00480 for (unsigned i = 0; it != endIt; ++i, ++it) {
00481 xPos_[i] = std::get<0>(*it);
00482 yPos_[i] = std::get<1>(*it);
00483 }
00484
00485 if (sortInputs)
00486 sortInput_();
00487 else if (xPos_[0] > xPos_[numSamples() - 1])
00488 reverseSamplingPoints_();
00489
00490
00491 makeFullSpline_(m0, m1);
00492 }
00493
00497
00501
00511 template <class ScalarArrayX, class ScalarArrayY>
00512 void setXYArrays(size_t nSamples,
00513 const ScalarArrayX& x,
00514 const ScalarArrayY& y,
00515 SplineType splineType = Natural,
00516 bool sortInputs = true)
00517 {
00518 assert(nSamples > 1);
00519
00520 setNumSamples_(nSamples);
00521 for (size_t i = 0; i < nSamples; ++i) {
00522 xPos_[i] = x[i];
00523 yPos_[i] = y[i];
00524 }
00525
00526 if (sortInputs)
00527 sortInput_();
00528 else if (xPos_[0] > xPos_[numSamples() - 1])
00529 reverseSamplingPoints_();
00530
00531 if (splineType == Periodic)
00532 makePeriodicSpline_();
00533 else if (splineType == Natural)
00534 makeNaturalSpline_();
00535 else if (splineType == Monotonic)
00536 this->makeMonotonicSpline_(slopeVec_);
00537 else
00538 OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
00539 }
00540
00552 template <class ScalarContainerX, class ScalarContainerY>
00553 void setXYContainers(const ScalarContainerX& x,
00554 const ScalarContainerY& y,
00555 SplineType splineType = Natural,
00556 bool sortInputs = true)
00557 {
00558 assert(x.size() == y.size());
00559 assert(x.size() > 1);
00560
00561 setNumSamples_(x.size());
00562 std::copy(x.begin(), x.end(), xPos_.begin());
00563 std::copy(y.begin(), y.end(), yPos_.begin());
00564
00565 if (sortInputs)
00566 sortInput_();
00567 else if (xPos_[0] > xPos_[numSamples() - 1])
00568 reverseSamplingPoints_();
00569
00570 if (splineType == Periodic)
00571 makePeriodicSpline_();
00572 else if (splineType == Natural)
00573 makeNaturalSpline_();
00574 else if (splineType == Monotonic)
00575 this->makeMonotonicSpline_(slopeVec_);
00576 else
00577 OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
00578 }
00579
00592 template <class PointArray>
00593 void setArrayOfPoints(size_t nSamples,
00594 const PointArray& points,
00595 SplineType splineType = Natural,
00596 bool sortInputs = true)
00597 {
00598
00599
00600 assert(nSamples > 1);
00601
00602 setNumSamples_(nSamples);
00603 for (size_t i = 0; i < nSamples; ++i) {
00604 xPos_[i] = points[i][0];
00605 yPos_[i] = points[i][1];
00606 }
00607
00608 if (sortInputs)
00609 sortInput_();
00610 else if (xPos_[0] > xPos_[numSamples() - 1])
00611 reverseSamplingPoints_();
00612
00613 if (splineType == Periodic)
00614 makePeriodicSpline_();
00615 else if (splineType == Natural)
00616 makeNaturalSpline_();
00617 else if (splineType == Monotonic)
00618 this->makeMonotonicSpline_(slopeVec_);
00619 else
00620 OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
00621 }
00622
00634 template <class XYContainer>
00635 void setContainerOfPoints(const XYContainer& points,
00636 SplineType splineType = Natural,
00637 bool sortInputs = true)
00638 {
00639
00640
00641 assert(points.size() > 1);
00642
00643 setNumSamples_(points.size());
00644 typename XYContainer::const_iterator it = points.begin();
00645 typename XYContainer::const_iterator endIt = points.end();
00646 for (size_t i = 0; it != endIt; ++ i, ++it) {
00647 xPos_[i] = (*it)[0];
00648 yPos_[i] = (*it)[1];
00649 }
00650
00651 if (sortInputs)
00652 sortInput_();
00653 else if (xPos_[0] > xPos_[numSamples() - 1])
00654 reverseSamplingPoints_();
00655
00656 if (splineType == Periodic)
00657 makePeriodicSpline_();
00658 else if (splineType == Natural)
00659 makeNaturalSpline_();
00660 else if (splineType == Monotonic)
00661 this->makeMonotonicSpline_(slopeVec_);
00662 else
00663 OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
00664 }
00665
00680 template <class XYContainer>
00681 void setContainerOfTuples(const XYContainer& points,
00682 SplineType splineType = Natural,
00683 bool sortInputs = true)
00684 {
00685
00686 setNumSamples_(points.size());
00687 typename XYContainer::const_iterator it = points.begin();
00688 typename XYContainer::const_iterator endIt = points.end();
00689 for (unsigned i = 0; it != endIt; ++i, ++it) {
00690 xPos_[i] = std::get<0>(*it);
00691 yPos_[i] = std::get<1>(*it);
00692 }
00693
00694 if (sortInputs)
00695 sortInput_();
00696 else if (xPos_[0] > xPos_[numSamples() - 1])
00697 reverseSamplingPoints_();
00698
00699 if (splineType == Periodic)
00700 makePeriodicSpline_();
00701 else if (splineType == Natural)
00702 makeNaturalSpline_();
00703 else if (splineType == Monotonic)
00704 this->makeMonotonicSpline_(slopeVec_);
00705 else
00706 OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
00707 }
00708
00712 template <class Evaluation>
00713 bool applies(const Evaluation& x) const
00714 { return x_(0) <= x && x <= x_(numSamples() - 1); }
00715
00719 size_t numSamples() const
00720 { return xPos_.size(); }
00721
00725 Scalar xAt(size_t sampleIdx) const
00726 { return x_(sampleIdx); }
00727
00731 Scalar valueAt(size_t sampleIdx) const
00732 { return y_(sampleIdx); }
00733
00751 void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream& os = std::cout) const
00752 {
00753 Scalar x0 = std::min(xi0, xi1);
00754 Scalar x1 = std::max(xi0, xi1);
00755 const size_t n = numSamples() - 1;
00756 for (size_t i = 0; i <= k; ++i) {
00757 double x = i*(x1 - x0)/k + x0;
00758 double x_p1 = x + (x1 - x0)/k;
00759 double y;
00760 double dy_dx;
00761 double mono = 1;
00762 if (!applies(x)) {
00763 if (x < x_(0)) {
00764 dy_dx = evalDerivative(x_(0));
00765 y = (x - x_(0))*dy_dx + y_(0);
00766 mono = (dy_dx>0)?1:-1;
00767 }
00768 else if (x > x_(n)) {
00769 dy_dx = evalDerivative(x_(n));
00770 y = (x - x_(n))*dy_dx + y_(n);
00771 mono = (dy_dx>0)?1:-1;
00772 }
00773 else {
00774 OPM_THROW(std::runtime_error,
00775 "The sampling points given to a spline must be sorted by their x value!");
00776 }
00777 }
00778 else {
00779 y = eval(x);
00780 dy_dx = evalDerivative(x);
00781 mono = monotonic(x, x_p1, true);
00782 }
00783
00784 os << x << " " << y << " " << dy_dx << " " << mono << "\n";
00785 }
00786 }
00787
00799 template <class Evaluation>
00800 Evaluation eval(const Evaluation& x, bool extrapolate = false) const
00801 {
00802 if (!extrapolate && !applies(x))
00803 OPM_THROW(Opm::NumericalProblem,
00804 "Tried to evaluate a spline outside of its range");
00805
00806
00807 if (extrapolate) {
00808 if (x < xAt(0)) {
00809 Scalar m = evalDerivative_(xAt(0), 0);
00810 Scalar y0 = y_(0);
00811 return y0 + m*(x - xAt(0));
00812 }
00813 else if (x > xAt(static_cast<size_t>(static_cast<long int>(numSamples()) - 1))) {
00814 Scalar m = evalDerivative_(xAt(static_cast<size_t>(numSamples() - 1)),
00815 static_cast<size_t>(numSamples()-2));
00816 Scalar y0 = y_(static_cast<size_t>(numSamples() - 1));
00817 return y0 + m*(x - xAt(static_cast<size_t>(numSamples() - 1)));
00818 }
00819 }
00820
00821 return eval_(x, segmentIdx_(Opm::scalarValue(x)));
00822 }
00823
00836 template <class Evaluation>
00837 Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
00838 {
00839 if (!extrapolate && !applies(x))
00840 OPM_THROW(Opm::NumericalProblem,
00841 "Tried to evaluate the derivative of a spline outside of its range");
00842
00843
00844 if (extrapolate) {
00845 if (x < xAt(0))
00846 return evalDerivative_(xAt(0), 0);
00847 else if (x > xAt(numSamples() - 1))
00848 return evalDerivative_(xAt(numSamples() - 1), numSamples() - 2);
00849 }
00850
00851 return evalDerivative_(x, segmentIdx_(Opm::scalarValue(x)));
00852 }
00853
00866 template <class Evaluation>
00867 Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
00868 {
00869 if (!extrapolate && !applies(x))
00870 OPM_THROW(Opm::NumericalProblem,
00871 "Tried to evaluate the second derivative of a spline outside of its range");
00872 else if (extrapolate)
00873 return 0.0;
00874
00875 return evalDerivative2_(x, segmentIdx_(Opm::scalarValue(x)));
00876 }
00877
00890 template <class Evaluation>
00891 Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
00892 {
00893 if (!extrapolate && !applies(x))
00894 OPM_THROW(Opm::NumericalProblem,
00895 "Tried to evaluate the third derivative of a spline outside of its range");
00896 else if (extrapolate)
00897 return 0.0;
00898
00899 return evalDerivative3_(x, segmentIdx_(Opm::scalarValue(x)));
00900 }
00901
00908 template <class Evaluation>
00909 Evaluation intersect(const Evaluation& a,
00910 const Evaluation& b,
00911 const Evaluation& c,
00912 const Evaluation& d) const
00913 {
00914 return intersectInterval(xAt(0), xAt(numSamples() - 1), a, b, c, d);
00915 }
00916
00923 template <class Evaluation>
00924 Evaluation intersectInterval(Scalar x0, Scalar x1,
00925 const Evaluation& a,
00926 const Evaluation& b,
00927 const Evaluation& c,
00928 const Evaluation& d) const
00929 {
00930 assert(applies(x0) && applies(x1));
00931
00932 Evaluation tmpSol[3], sol = 0;
00933 size_t nSol = 0;
00934 size_t iFirst = segmentIdx_(x0);
00935 size_t iLast = segmentIdx_(x1);
00936 for (size_t i = iFirst; i <= iLast; ++i)
00937 {
00938 size_t nCur = intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
00939 if (nCur == 1)
00940 sol = tmpSol[0];
00941
00942 nSol += nCur;
00943 if (nSol > 1) {
00944 OPM_THROW(std::runtime_error,
00945 "Spline has more than one intersection");
00946 }
00947 }
00948
00949 if (nSol != 1)
00950 OPM_THROW(std::runtime_error,
00951 "Spline has no intersection");
00952
00953 return sol;
00954 }
00955
00964 int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED = false) const
00965 {
00966 assert(std::abs(x0 - x1) > 1e-30);
00967
00968
00969 if (x0 > x1)
00970 std::swap(x0, x1);
00971
00972 assert(x0 < x1);
00973
00974 int r = 3;
00975 if (x0 < xAt(0)) {
00976 assert(extrapolate);
00977 Scalar m = evalDerivative_(xAt(0), 0);
00978 if (std::abs(m) < 1e-20)
00979 r = (m < 0)?-1:1;
00980 x0 = xAt(0);
00981 };
00982
00983 size_t i = segmentIdx_(x0);
00984 if (x_(i + 1) >= x1) {
00985
00986
00987 monotonic_(i, x0, x1, r);
00988 return r;
00989 }
00990
00991
00992
00993 monotonic_(i, x0, x_(i+1), r);
00994 ++ i;
00995
00996
00997
00998 size_t iEnd = segmentIdx_(x1);
00999 for (; i < iEnd - 1; ++i) {
01000 monotonic_(i, x_(i), x_(i + 1), r);
01001 if (!r)
01002 return 0;
01003 }
01004
01005
01006
01007
01008 if (x1 > xAt(numSamples() - 1)) {
01009 assert(extrapolate);
01010
01011 Scalar m = evalDerivative_(xAt(numSamples() - 1), numSamples() - 2);
01012 if (m < 0)
01013 return (r < 0 || r==3)?-1:0;
01014 else if (m > 0)
01015 return (r > 0 || r==3)?1:0;
01016
01017 return r;
01018 }
01019
01020
01021 monotonic_(iEnd, x_(iEnd), x1, r);
01022
01023 return r;
01024 }
01025
01030 int monotonic() const
01031 { return monotonic(xAt(0), xAt(numSamples() - 1)); }
01032
01033 protected:
01037 struct ComparatorX_
01038 {
01039 ComparatorX_(const std::vector<Scalar>& x)
01040 : x_(x)
01041 {}
01042
01043 bool operator ()(unsigned idxA, unsigned idxB) const
01044 { return x_.at(idxA) < x_.at(idxB); }
01045
01046 const std::vector<Scalar>& x_;
01047 };
01048
01052 void sortInput_()
01053 {
01054 size_t n = numSamples();
01055
01056
01057 std::vector<unsigned> idxVector(n);
01058 for (unsigned i = 0; i < n; ++i)
01059 idxVector[i] = i;
01060
01061
01062
01063 ComparatorX_ cmp(xPos_);
01064 std::sort(idxVector.begin(), idxVector.end(), cmp);
01065
01066
01067 std::vector<Scalar> tmpX(n), tmpY(n);
01068 for (size_t i = 0; i < idxVector.size(); ++ i) {
01069 tmpX[i] = xPos_[idxVector[i]];
01070 tmpY[i] = yPos_[idxVector[i]];
01071 }
01072 xPos_ = tmpX;
01073 yPos_ = tmpY;
01074 }
01075
01080 void reverseSamplingPoints_()
01081 {
01082
01083 size_t n = numSamples();
01084 for (unsigned i = 0; i <= (n - 1)/2; ++i) {
01085 std::swap(xPos_[i], xPos_[n - i - 1]);
01086 std::swap(yPos_[i], yPos_[n - i - 1]);
01087 }
01088 }
01089
01093 void setNumSamples_(size_t nSamples)
01094 {
01095 xPos_.resize(nSamples);
01096 yPos_.resize(nSamples);
01097 slopeVec_.resize(nSamples);
01098 }
01099
01105 void makeFullSpline_(Scalar m0, Scalar m1)
01106 {
01107 Matrix M(numSamples());
01108 std::vector<Scalar> d(numSamples());
01109 std::vector<Scalar> moments(numSamples());
01110
01111
01112 this->makeFullSystem_(M, d, m0, m1);
01113
01114
01115 M.solve(moments, d);
01116
01117
01118 this->setSlopesFromMoments_(slopeVec_, moments);
01119 }
01120
01126 void makeNaturalSpline_()
01127 {
01128 Matrix M(numSamples(), numSamples());
01129 Vector d(numSamples());
01130 Vector moments(numSamples());
01131
01132
01133 this->makeNaturalSystem_(M, d);
01134
01135
01136 M.solve(moments, d);
01137
01138
01139 this->setSlopesFromMoments_(slopeVec_, moments);
01140 }
01141
01147 void makePeriodicSpline_()
01148 {
01149 Matrix M(numSamples() - 1);
01150 Vector d(numSamples() - 1);
01151 Vector moments(numSamples() - 1);
01152
01153
01154
01155
01156
01157 this->makePeriodicSystem_(M, d);
01158
01159
01160 M.solve(moments, d);
01161
01162 moments.resize(numSamples());
01163 for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i) {
01164 unsigned ui = static_cast<unsigned>(i);
01165 moments[ui+1] = moments[ui];
01166 }
01167 moments[0] = moments[numSamples() - 1];
01168
01169
01170 this->setSlopesFromMoments_(slopeVec_, moments);
01171 }
01172
01179 template <class DestVector, class SourceVector>
01180 void assignSamplingPoints_(DestVector& destX,
01181 DestVector& destY,
01182 const SourceVector& srcX,
01183 const SourceVector& srcY,
01184 unsigned nSamples)
01185 {
01186 assert(nSamples >= 2);
01187
01188
01189
01190 for (unsigned i = 0; i < nSamples; ++i) {
01191 unsigned idx = i;
01192 if (srcX[0] > srcX[nSamples - 1])
01193 idx = nSamples - i - 1;
01194 destX[i] = srcX[idx];
01195 destY[i] = srcY[idx];
01196 }
01197 }
01198
01199 template <class DestVector, class ListIterator>
01200 void assignFromArrayList_(DestVector& destX,
01201 DestVector& destY,
01202 const ListIterator& srcBegin,
01203 const ListIterator& srcEnd,
01204 unsigned nSamples)
01205 {
01206 assert(nSamples >= 2);
01207
01208
01209 ListIterator it = srcBegin;
01210 ++it;
01211 bool reverse = false;
01212 if ((*srcBegin)[0] > (*it)[0])
01213 reverse = true;
01214 --it;
01215
01216
01217 for (unsigned i = 0; it != srcEnd; ++i, ++it) {
01218 unsigned idx = i;
01219 if (reverse)
01220 idx = nSamples - i - 1;
01221 destX[i] = (*it)[0];
01222 destY[i] = (*it)[1];
01223 }
01224 }
01225
01233 template <class DestVector, class ListIterator>
01234 void assignFromTupleList_(DestVector& destX,
01235 DestVector& destY,
01236 ListIterator srcBegin,
01237 ListIterator srcEnd,
01238 unsigned nSamples)
01239 {
01240 assert(nSamples >= 2);
01241
01242
01243
01244
01245
01246 ListIterator it = srcBegin;
01247 ++it;
01248 bool reverse = false;
01249 if (std::get<0>(*srcBegin) > std::get<0>(*it))
01250 reverse = true;
01251 --it;
01252
01253
01254 for (unsigned i = 0; it != srcEnd; ++i, ++it) {
01255 unsigned idx = i;
01256 if (reverse)
01257 idx = nSamples - i - 1;
01258 destX[i] = std::get<0>(*it);
01259 destY[i] = std::get<1>(*it);
01260 }
01261 }
01262
01267 template <class Vector, class Matrix>
01268 void makeFullSystem_(Matrix& M, Vector& d, Scalar m0, Scalar m1)
01269 {
01270 makeNaturalSystem_(M, d);
01271
01272 size_t n = numSamples() - 1;
01273
01274 M[0][1] = 1;
01275 d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0);
01276
01277
01278 M[n][n - 1] = 1;
01279
01280
01281 d[n] =
01282 6/h_(n)
01283 *
01284 (m1 - (y_(n) - y_(n - 1))/h_(n));
01285 }
01286
01291 template <class Vector, class Matrix>
01292 void makeNaturalSystem_(Matrix& M, Vector& d)
01293 {
01294 M = 0.0;
01295
01296
01297
01298 size_t n = numSamples() - 1;
01299
01300
01301 for (size_t i = 1; i < n; ++i) {
01302 Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
01303 Scalar mu_i = 1 - lambda_i;
01304 Scalar d_i =
01305 6 / (h_(i) + h_(i + 1))
01306 *
01307 ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
01308
01309 M[i][i-1] = mu_i;
01310 M[i][i] = 2;
01311 M[i][i + 1] = lambda_i;
01312 d[i] = d_i;
01313 };
01314
01315
01316 Scalar lambda_0 = 0;
01317 Scalar d_0 = 0;
01318
01319 Scalar mu_n = 0;
01320 Scalar d_n = 0;
01321
01322
01323 M[0][0] = 2;
01324 M[0][1] = lambda_0;
01325 d[0] = d_0;
01326
01327
01328 M[n][n-1] = mu_n;
01329 M[n][n] = 2;
01330 d[n] = d_n;
01331 }
01332
01337 template <class Matrix, class Vector>
01338 void makePeriodicSystem_(Matrix& M, Vector& d)
01339 {
01340 M = 0.0;
01341
01342
01343
01344 size_t n = numSamples() - 1;
01345
01346 assert(M.rows() == n);
01347
01348
01349 for (size_t i = 2; i < n; ++i) {
01350 Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
01351 Scalar mu_i = 1 - lambda_i;
01352 Scalar d_i =
01353 6 / (h_(i) + h_(i + 1))
01354 *
01355 ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
01356
01357 M[i-1][i-2] = mu_i;
01358 M[i-1][i-1] = 2;
01359 M[i-1][i] = lambda_i;
01360 d[i-1] = d_i;
01361 };
01362
01363 Scalar lambda_n = h_(1) / (h_(n) + h_(1));
01364 Scalar lambda_1 = h_(2) / (h_(1) + h_(2));;
01365 Scalar mu_1 = 1 - lambda_1;
01366 Scalar mu_n = 1 - lambda_n;
01367
01368 Scalar d_1 =
01369 6 / (h_(1) + h_(2))
01370 *
01371 ( (y_(2) - y_(1))/h_(2) - (y_(1) - y_(0))/h_(1));
01372 Scalar d_n =
01373 6 / (h_(n) + h_(1))
01374 *
01375 ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n));
01376
01377
01378
01379 M[0][0] = 2;
01380 M[0][1] = lambda_1;
01381 M[0][n-1] = mu_1;
01382 d[0] = d_1;
01383
01384
01385 M[n-1][0] = lambda_n;
01386 M[n-1][n-2] = mu_n;
01387 M[n-1][n-1] = 2;
01388 d[n-1] = d_n;
01389 }
01390
01399 template <class Vector>
01400 void makeMonotonicSpline_(Vector& slopes)
01401 {
01402 auto n = numSamples();
01403
01404
01405 std::vector<Scalar> delta(n);
01406 for (size_t k = 0; k < n - 1; ++k)
01407 delta[k] = (y_(k + 1) - y_(k))/(x_(k + 1) - x_(k));
01408
01409
01410 for (size_t k = 1; k < n - 1; ++k)
01411 slopes[k] = (delta[k - 1] + delta[k])/2;
01412 slopes[0] = delta[0];
01413 slopes[n - 1] = delta[n - 2];
01414
01415
01416 for (size_t k = 0; k < n - 1; ++k) {
01417 if (std::abs(delta[k]) < 1e-50) {
01418
01419 slopes[k] = 0;
01420 slopes[k + 1] = 0;
01421 ++ k;
01422 continue;
01423 }
01424 else {
01425 Scalar alpha = slopes[k] / delta[k];
01426 Scalar beta = slopes[k + 1] / delta[k];
01427
01428 if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
01429 slopes[k] = 0;
01430 }
01431
01432 else if (alpha*alpha + beta*beta > 3*3) {
01433 Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
01434 slopes[k] = tau*alpha*delta[k];
01435 slopes[k + 1] = tau*beta*delta[k];
01436 }
01437 }
01438 }
01439 }
01440
01448 template <class MomentsVector, class SlopeVector>
01449 void setSlopesFromMoments_(SlopeVector& slopes, const MomentsVector& moments)
01450 {
01451 size_t n = numSamples();
01452
01453
01454
01455
01456 Scalar mRight;
01457
01458 {
01459 Scalar h = this->h_(n - 1);
01460 Scalar x = h;
01461
01462
01463 Scalar A =
01464 (y_(n - 1) - y_(n - 2))/h
01465 -
01466 h/6*(moments[n-1] - moments[n - 2]);
01467
01468 mRight =
01469
01470
01471 moments[n - 1] * x*x / (2 * h)
01472 +
01473 A;
01474 }
01475
01476
01477 for (size_t i = 0; i < n - 1; ++ i) {
01478
01479
01480 Scalar h_i = this->h_(i + 1);
01481
01482 Scalar x_i1 = h_i;
01483
01484 Scalar A_i =
01485 (y_(i+1) - y_(i))/h_i
01486 -
01487 h_i/6*(moments[i+1] - moments[i]);
01488
01489 slopes[i] =
01490 - moments[i] * x_i1*x_i1 / (2 * h_i)
01491 +
01492
01493
01494 A_i;
01495
01496 }
01497 slopes[n - 1] = mRight;
01498 }
01499
01500
01501
01502
01503 template <class Evaluation>
01504 Evaluation eval_(const Evaluation& x, size_t i) const
01505 {
01506
01507 Scalar delta = h_(i + 1);
01508 Evaluation t = (x - x_(i))/delta;
01509
01510 return
01511 h00_(t) * y_(i)
01512 + h10_(t) * slope_(i)*delta
01513 + h01_(t) * y_(i + 1)
01514 + h11_(t) * slope_(i + 1)*delta;
01515 }
01516
01517
01518
01519 template <class Evaluation>
01520 Evaluation evalDerivative_(const Evaluation& x, size_t i) const
01521 {
01522
01523 Scalar delta = h_(i + 1);
01524 Evaluation t = (x - x_(i))/delta;
01525 Evaluation alpha = 1 / delta;
01526
01527 return
01528 alpha *
01529 (h00_prime_(t) * y_(i)
01530 + h10_prime_(t) * slope_(i)*delta
01531 + h01_prime_(t) * y_(i + 1)
01532 + h11_prime_(t) * slope_(i + 1)*delta);
01533 }
01534
01535
01536
01537 template <class Evaluation>
01538 Evaluation evalDerivative2_(const Evaluation& x, size_t i) const
01539 {
01540
01541 Scalar delta = h_(i + 1);
01542 Evaluation t = (x - x_(i))/delta;
01543 Evaluation alpha = 1 / delta;
01544
01545 return
01546 alpha*alpha
01547 *(h00_prime2_(t) * y_(i)
01548 + h10_prime2_(t) * slope_(i)*delta
01549 + h01_prime2_(t) * y_(i + 1)
01550 + h11_prime2_(t) * slope_(i + 1)*delta);
01551 }
01552
01553
01554
01555 template <class Evaluation>
01556 Evaluation evalDerivative3_(const Evaluation& x, size_t i) const
01557 {
01558
01559 Scalar delta = h_(i + 1);
01560 Evaluation t = (x - x_(i))/delta;
01561 Evaluation alpha = 1 / delta;
01562
01563 return
01564 alpha*alpha*alpha
01565 *(h00_prime3_(t)*y_(i)
01566 + h10_prime3_(t)*slope_(i)*delta
01567 + h01_prime3_(t)*y_(i + 1)
01568 + h11_prime3_(t)*slope_(i + 1)*delta);
01569 }
01570
01571
01572 template <class Evaluation>
01573 Evaluation h00_(const Evaluation& t) const
01574 { return (2*t - 3)*t*t + 1; }
01575
01576 template <class Evaluation>
01577 Evaluation h10_(const Evaluation& t) const
01578 { return ((t - 2)*t + 1)*t; }
01579
01580 template <class Evaluation>
01581 Evaluation h01_(const Evaluation& t) const
01582 { return (-2*t + 3)*t*t; }
01583
01584 template <class Evaluation>
01585 Evaluation h11_(const Evaluation& t) const
01586 { return (t - 1)*t*t; }
01587
01588
01589 template <class Evaluation>
01590 Evaluation h00_prime_(const Evaluation& t) const
01591 { return (3*2*t - 2*3)*t; }
01592
01593 template <class Evaluation>
01594 Evaluation h10_prime_(const Evaluation& t) const
01595 { return (3*t - 2*2)*t + 1; }
01596
01597 template <class Evaluation>
01598 Evaluation h01_prime_(const Evaluation& t) const
01599 { return (-3*2*t + 2*3)*t; }
01600
01601 template <class Evaluation>
01602 Evaluation h11_prime_(const Evaluation& t) const
01603 { return (3*t - 2)*t; }
01604
01605
01606 template <class Evaluation>
01607 Evaluation h00_prime2_(const Evaluation& t) const
01608 { return 2*3*2*t - 2*3; }
01609
01610 template <class Evaluation>
01611 Evaluation h10_prime2_(const Evaluation& t) const
01612 { return 2*3*t - 2*2; }
01613
01614 template <class Evaluation>
01615 Evaluation h01_prime2_(const Evaluation& t) const
01616 { return -2*3*2*t + 2*3; }
01617
01618 template <class Evaluation>
01619 Evaluation h11_prime2_(const Evaluation& t) const
01620 { return 2*3*t - 2; }
01621
01622
01623 template <class Evaluation>
01624 Scalar h00_prime3_(const Evaluation& t OPM_UNUSED) const
01625 { return 2*3*2; }
01626
01627 template <class Evaluation>
01628 Scalar h10_prime3_(const Evaluation& t OPM_UNUSED) const
01629 { return 2*3; }
01630
01631 template <class Evaluation>
01632 Scalar h01_prime3_(const Evaluation& t OPM_UNUSED) const
01633 { return -2*3*2; }
01634
01635 template <class Evaluation>
01636 Scalar h11_prime3_(const Evaluation& t OPM_UNUSED) const
01637 { return 2*3; }
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647 int monotonic_(size_t i, Scalar x0, Scalar x1, int& r) const
01648 {
01649
01650 Scalar a = 3*a_(i);
01651 Scalar b = 2*b_(i);
01652 Scalar c = c_(i);
01653
01654 if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
01655 return 3;
01656
01657 Scalar disc = b*b - 4*a*c;
01658 if (disc < 0) {
01659
01660
01661 if (x0*(x0*a + b) + c > 0) {
01662 r = (r==3 || r == 1)?1:0;
01663 return 1;
01664 }
01665 else {
01666 r = (r==3 || r == -1)?-1:0;
01667 return -1;
01668 }
01669 }
01670 disc = std::sqrt(disc);
01671 Scalar xE1 = (-b + disc)/(2*a);
01672 Scalar xE2 = (-b - disc)/(2*a);
01673
01674 if (std::abs(disc) < 1e-30) {
01675
01676 if (std::abs(xE1 - x0) < 1e-30)
01677
01678
01679
01680 x0 = x1;
01681 if (x0*(x0*a + b) + c > 0) {
01682 r = (r==3 || r == 1)?1:0;
01683 return 1;
01684 }
01685 else {
01686 r = (r==3 || r == -1)?-1:0;
01687 return -1;
01688 }
01689 };
01690 if ((x0 < xE1 && xE1 < x1) ||
01691 (x0 < xE2 && xE2 < x1))
01692 {
01693
01694 r = 0;
01695 return 0;
01696 }
01697
01698 x0 = (x0 + x1)/2;
01699
01700 if (x0*(x0*a + b) + c > 0) {
01701 r = (r==3 || r == 1)?1:0;
01702 return 1;
01703 }
01704 else {
01705 r = (r==3 || r == -1)?-1:0;
01706 return -1;
01707 }
01708 }
01709
01714 template <class Evaluation>
01715 size_t intersectSegment_(Evaluation* sol,
01716 size_t segIdx,
01717 const Evaluation& a,
01718 const Evaluation& b,
01719 const Evaluation& c,
01720 const Evaluation& d,
01721 Scalar x0 = -1e30, Scalar x1 = 1e30) const
01722 {
01723 unsigned n =
01724 Opm::invertCubicPolynomial(sol,
01725 a_(segIdx) - a,
01726 b_(segIdx) - b,
01727 c_(segIdx) - c,
01728 d_(segIdx) - d);
01729 x0 = std::max(x_(segIdx), x0);
01730 x1 = std::min(x_(segIdx+1), x1);
01731
01732
01733 size_t k = 0;
01734 for (unsigned j = 0; j < n; ++j) {
01735 if (x0 <= sol[j] && sol[j] <= x1) {
01736 sol[k] = sol[j];
01737 ++k;
01738 }
01739 }
01740 return k;
01741 }
01742
01743
01744 size_t segmentIdx_(Scalar x) const
01745 {
01746
01747 size_t iLow = 0;
01748 size_t iHigh = numSamples() - 1;
01749
01750 while (iLow + 1 < iHigh) {
01751 size_t i = (iLow + iHigh) / 2;
01752 if (x_(i) > x)
01753 iHigh = i;
01754 else
01755 iLow = i;
01756 };
01757 return iLow;
01758 }
01759
01763 Scalar h_(size_t i) const
01764 {
01765 assert(x_(i) > x_(i-1));
01766
01767 return x_(i) - x_(i - 1);
01768 }
01769
01773 Scalar x_(size_t i) const
01774 { return xPos_[i]; }
01775
01779 Scalar y_(size_t i) const
01780 { return yPos_[i]; }
01781
01786 Scalar slope_(size_t i) const
01787 { return slopeVec_[i]; }
01788
01789
01790
01791 Scalar a_(size_t i) const
01792 { return evalDerivative3_(Scalar(0.0), i)/6.0; }
01793
01794
01795
01796 Scalar b_(size_t i) const
01797 { return evalDerivative2_(Scalar(0.0), i)/2.0; }
01798
01799
01800
01801 Scalar c_(size_t i) const
01802 { return evalDerivative_(Scalar(0.0), i); }
01803
01804
01805
01806 Scalar d_(size_t i) const
01807 { return eval_(Scalar(0.0), i); }
01808
01809 Vector xPos_;
01810 Vector yPos_;
01811 Vector slopeVec_;
01812 };
01813 }
01814
01815 #endif