diff --git a/dumux/common/splinecommon_.hh b/dumux/common/splinecommon_.hh index 2a1529f971b2ddfbca1e0eed2154967a157cb75b..691bf1efa769f32aabaff17a1641256d23ce5ed6 100644 --- a/dumux/common/splinecommon_.hh +++ b/dumux/common/splinecommon_.hh @@ -145,12 +145,12 @@ public: // handle extrapolation if (extrapolate) { if (x < xMin()) { - Scalar m = evalDerivative(xMin(), /*segmentIdx=*/0); + Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0); Scalar y0 = y_(0); return y0 + m*(x - xMin()); } else if (x > xMax()) { - Scalar m = evalDerivative(xMax(), /*segmentIdx=*/numSamples_()-1); + Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples_()-2); Scalar y0 = y_(numSamples_()-1); return y0 + m*(x - xMax()); } @@ -179,7 +179,7 @@ public: if (Dune::FloatCmp::le(x, xMin())) return evalDerivative_(xMin(), 0); else if (Dune::FloatCmp::ge(x, xMax())) - return evalDerivative_(xMax(), numSamples_() - 1); + return evalDerivative_(xMax(), numSamples_() - 2); } return evalDerivative_(x, segmentIdx_(x)); diff --git a/test/common/spline/test_spline.cc b/test/common/spline/test_spline.cc index 34ae0fc6c5bb99cc6a91627668c1c8739f680ef5..b55adf99560806d8bfb835bf06828221d838da7e 100644 --- a/test/common/spline/test_spline.cc +++ b/test/common/spline/test_spline.cc @@ -131,6 +131,39 @@ void testNatural(const Spline &sp, << (d3 - d2)/eps << " ought to be 0"); } +/*! + * \brief Tests extrapolation, see issue 925 + * + * Note: tests only Spline<double, -1> since the implementation is common to + * to all spline classes. + */ +void testExtrapolation() +{ + std::vector<double> x = {0.0, 0.25, 0.5, 0.75, 1.0}; + std::vector<double> y = {1.0, 1.0, 1.0, 1.0, 1.0}; + + Dumux::Spline<double, -1> sp(x, y); + + double xi = sp.eval(0.6); + if (Dune::FloatCmp::ne(xi, 1.0)) + DUNE_THROW(Dune::InvalidStateException, + "Spline interpolation gives " + << xi << " ought to be 1"); + + xi = sp.eval(-0.1, /*interpolate=*/true); + if (Dune::FloatCmp::ne(xi, 1.0)) + DUNE_THROW(Dune::InvalidStateException, + "Spline extrapolation at -0.1 gives " + << xi << " ought to be 1"); + + xi = sp.eval(1.1, /*interpolate=*/true); + if (Dune::FloatCmp::ne(xi, 1.0)) + DUNE_THROW(Dune::InvalidStateException, + "Spline extrapolation at 1.1 gives " + << xi << " ought to be 1"); +} + + void testAll() { double x[] = { 0, 5, 7.5, 8.75, 9.375 }; @@ -238,6 +271,7 @@ void plot() int main(int argc, char** argv) { testAll(); + testExtrapolation(); plot(); return 0;