diff --git a/opm/material/common/Spline.hpp b/opm/material/common/Spline.hpp index 8b06ecd66..8816afad4 100644 --- a/opm/material/common/Spline.hpp +++ b/opm/material/common/Spline.hpp @@ -251,12 +251,6 @@ public: bool sortInputs = false) { this->setContainerOfPoints(points, m0, m1, sortInputs); } - /*! - * \brief Returns the number of sampling points. - */ - size_t numSamples() const - { return xPos_.size(); } - /*! * \brief Set the sampling points and the boundary slopes of the * spline with two sampling points. @@ -717,29 +711,24 @@ public: return x_(0) <= x && x <= x_(numSamples() - 1); } - /*! - * \brief Return the x value of the leftmost sampling point. - */ - Scalar xMin() const - { return x_(0); } /*! - * \brief Return the x value of the rightmost sampling point. + * \brief Return the number of (x, y) values. */ - Scalar xMax() const - { return x_(numSamples() - 1); } + size_t numSamples() const + { return xPos_.size(); } /*! - * \brief Return the y value of the leftmost sampling point. + * \brief Return the x value of a given sampling point. */ - Scalar yFirst() const - { return y_(0); } + Scalar xAt(int sampleIdx) const + { return x_(sampleIdx); } /*! - * \brief Return the y value of the rightmost sampling point. + * \brief Return the x value of a given sampling point. */ - Scalar yLast() const - { return y_(numSamples() - 1); } + Scalar valueAt(int sampleIdx) const + { return y_(sampleIdx); } /*! * \brief Prints k tuples of the format (x, y, dx/dy, isMonotonic) @@ -812,15 +801,15 @@ public: // handle extrapolation if (extrapolate) { - if (x < xMin()) { - Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0); + if (x < xAt(0)) { + Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0); Scalar y0 = y_(0); - return y0 + m*(x - xMin()); + return y0 + m*(x - xAt(0)); } - else if (x > xMax()) { - Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples()-2); + else if (x > xAt(numSamples() - 1)) { + Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples()-2); Scalar y0 = y_(numSamples() - 1); - return y0 + m*(x - xMax()); + return y0 + m*(x - xAt(numSamples() - 1)); } } @@ -845,15 +834,15 @@ public: // handle extrapolation if (extrapolate) { - if (x.value < xMin()) { - Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0); + if (x.value < xAt(0)) { + Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0); Scalar y0 = y_(0); - return Evaluation::createConstant(y0 + m*(x.value - xMin())); + return Evaluation::createConstant(y0 + m*(x.value - xAt(0))); } - else if (x > xMax()) { - Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples()-2); + else if (x > xAt(numSamples() - 1)) { + Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples()-2); Scalar y0 = y_(numSamples() - 1); - return Evaluation::createConstant(y0 + m*(x.value - xMax())); + return Evaluation::createConstant(y0 + m*(x.value - xAt(numSamples() - 1))); } } @@ -876,10 +865,10 @@ public: { assert(extrapolate || applies(x)); if (extrapolate) { - if (x < xMin()) - return evalDerivative_(xMin(), /*segmentIdx=*/0); - else if (x > xMax()) - return evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2); + if (x < xAt(0)) + return evalDerivative_(xAt(0), /*segmentIdx=*/0); + else if (x > xAt(numSamples() - 1)) + return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2); } return evalDerivative_(x, segmentIdx_(x)); @@ -939,7 +928,7 @@ public: const Evaluation& c, const Evaluation& d) const { - return intersectInterval(xMin(), xMax(), a, b, c, d); + return intersectInterval(xAt(0), xAt(numSamples() - 1), a, b, c, d); } /*! @@ -1000,12 +989,12 @@ public: assert(x0 < x1); int r = 3; - if (x0 < xMin()) { + if (x0 < xAt(0)) { assert(extrapolate); - Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0); + Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0); if (std::abs(m) < 1e-20) r = (m < 0)?-1:1; - x0 = xMin(); + x0 = xAt(0); }; size_t i = segmentIdx_(x0); @@ -1033,10 +1022,10 @@ public: // if the user asked for a part of the spline which is // extrapolated, we need to check the slope at the spline's // endpoint - if (x1 > xMax()) { + if (x1 > xAt(numSamples() - 1)) { assert(extrapolate); - Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2); + Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2); if (m < 0) return (r < 0 || r==3)?-1:0; else if (m > 0) @@ -1056,7 +1045,7 @@ public: * spline as interval. */ int monotonic() const - { return monotonic(xMin(), xMax()); } + { return monotonic(xAt(0), xAt(numSamples() - 1)); } protected: /*! diff --git a/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp index bb1af9575..5ecd5020b 100644 --- a/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp @@ -139,12 +139,13 @@ public: static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) { // this assumes that the capillary pressure is monotonically decreasing - if (Sw <= params.pcnwSpline().xMin()) - return Evaluation(params.pcnwSpline().yFirst()); - if (Sw >= params.pcnwSpline().xMax()) - return Evaluation(params.pcnwSpline().yLast()); + const auto& pcnwSpline = params.pcnwSpline(); + if (Sw <= pcnwSpline.xAt(0)) + return Evaluation(pcnwSpline.valueAt(0)); + if (Sw >= pcnwSpline.xAt(pcnwSpline.numSamples() - 1)) + return Evaluation(pcnwSpline.valueAt(pcnwSpline.numSamples() - 1)); - return params.pcnwSpline().eval(Sw); + return pcnwSpline.eval(Sw); } template @@ -153,14 +154,15 @@ public: static const Evaluation nil(0.0); // this assumes that the capillary pressure is monotonically decreasing - if (pcnw >= params.pcnwSpline().yFirst()) - return Evaluation(params.pcnwSpline().xMin()); - if (pcnw <= params.pcnwSpline().yLast()) - return Evaluation(params.pcnwSpline().xMax()); + const auto& pcnwSpline = params.pcnwSpline(); + if (pcnw >= pcnwSpline.valueAt(0)) + return Evaluation(pcnwSpline.xAt(0)); + if (pcnw <= pcnwSpline.y(pcnwSpline.numSamples() - 1)) + return Evaluation(pcnwSpline.xAt(pcnwSpline.numSamples() - 1)); // the intersect() method of splines is a bit slow, but this code path is not too // time critical... - return params.pcnwSpline().intersect(/*a=*/nil, /*b=*/nil, /*c=*/nil, /*d=*/pcnw); + return pcnwSpline.intersect(/*a=*/nil, /*b=*/nil, /*c=*/nil, /*d=*/pcnw); } /*! @@ -204,12 +206,13 @@ public: template static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { - if (Sw <= params.krnSpline().xMin()) - return Evaluation(params.krwSpline().yFirst()); - if (Sw >= params.krnSpline().xMax()) - return Evaluation(params.krwSpline().yLast()); + const auto& krwSpline = params.krwSpline(); + if (Sw <= krwSpline.xAt(0)) + return Evaluation(krwSpline.valueAt(0)); + if (Sw >= krwSpline.xAt(krwSpline.numSamples() - 1)) + return Evaluation(krwSpline.valueAt(krwSpline.numSamples() - 1)); - return params.krwSpline().eval(Sw); + return krwSpline.eval(Sw); } template @@ -217,12 +220,13 @@ public: { static const Evaluation nil(0.0); - if (krw <= params.krwSpline().yFirst()) - return Evaluation(params.krwSpline().xMin()); - if (krw >= params.krwSpline().yLast()) - return Evaluation(params.krwSpline().xMax()); + const auto& krwSpline = params.krwSpline(); + if (krw <= krwSpline.valueAt(0)) + return Evaluation(krwSpline.xAt(0)); + if (krw >= krwSpline.valueAt(krwSpline.numSamples() - 1)) + return Evaluation(krwSpline.xAt(krwSpline.numSamples() - 1)); - return params.krwSpline().intersect(/*a=*/nil, /*b=*/nil, /*c=*/nil, /*d=*/krw); + return krwSpline.intersect(/*a=*/nil, /*b=*/nil, /*c=*/nil, /*d=*/krw); } /*! @@ -243,12 +247,13 @@ public: template static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) { - if (Sw <= params.krnSpline().xMin()) - return Evaluation(params.krnSpline().yFirst()); - if (Sw >= params.krnSpline().xMax()) - return Evaluation(params.krnSpline().yLast()); + const auto& krnSpline = params.krnSpline(); + if (Sw <= krnSpline.xAt(0)) + return Evaluation(krnSpline.valueAt(0)); + if (Sw >= krnSpline.xAt(krnSpline.numSamples() - 1)) + return Evaluation(krnSpline.valueAt(krnSpline.numSamples() - 1)); - return params.krnSpline().eval(Sw); + return krnSpline.eval(Sw); } template @@ -256,12 +261,13 @@ public: { static const Evaluation nil(0.0); - if (krn >= params.krnSpline().yFirst()) - return Evaluation(params.krnSpline().xMin()); - if (krn <= params.krnSpline().yLast()) - return Evaluation(params.krnSpline().xMax()); + const auto& krnSpline = params.krnSpline(); + if (krn >= krnSpline.valueAt(0)) + return Evaluation(krnSpline.xAt(0)); + if (krn <= krnSpline.valueAt(krnSpline.numSamples() - 1)) + return Evaluation(krnSpline.xAt(krnSpline.numSamples() - 1)); - return params.krnSpline().intersect(/*a=*/nil, /*b=*/nil, /*c=*/nil, /*d=*/krn); + return krnSpline.intersect(/*a=*/nil, /*b=*/nil, /*c=*/nil, /*d=*/krn); } }; } // namespace Opm diff --git a/tests/test_spline.cpp b/tests/test_spline.cpp index 2aa24a15c..c9b11625b 100644 --- a/tests/test_spline.cpp +++ b/tests/test_spline.cpp @@ -82,7 +82,9 @@ void testCommon(const Spline &sp, // make sure the derivatives are consistent with the curve size_t np = 3*n; for (size_t i = 0; i < np; ++i) { - double xval = sp.xMin() + (sp.xMax() - sp.xMin())*i/np; + double xMin = sp.xAt(0); + double xMax = sp.xAt(sp.numSamples() - 1); + double xval = xMin + (xMax - xMin)*i/np; // first derivative double y1 = sp.eval(xval+epsFD);