Commit cc39be0a authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'fix/issue-925' into 'master'

[common] fix segment index issue in spline extrapolation

Closes #925

See merge request !2235
parents 6d430f90 dfae9ed1
......@@ -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));
......
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment