27 #ifndef OPM_SPLINE_HPP 28 #define OPM_SPLINE_HPP 32 #include <opm/common/ErrorMacros.hpp> 33 #include <opm/common/Exceptions.hpp> 34 #include <opm/common/Unused.hpp> 91 template<
class Scalar>
95 typedef std::vector<Scalar> Vector;
129 Scalar y0, Scalar y1,
130 Scalar m0, Scalar m1)
131 {
set(x0, x1, y0, y1, m0, m1); }
141 template <
class ScalarArrayX,
class ScalarArrayY>
143 const ScalarArrayX& x,
144 const ScalarArrayY& y,
146 bool sortInputs =
true)
147 { this->
setXYArrays(nSamples, x, y, splineType, sortInputs); }
156 template <
class Po
intArray>
158 const PointArray& points,
160 bool sortInputs =
true)
170 template <
class ScalarContainer>
172 const ScalarContainer& y,
174 bool sortInputs =
true)
183 template <
class Po
intContainer>
186 bool sortInputs =
true)
199 template <
class ScalarArray>
201 const ScalarArray& x,
202 const ScalarArray& y,
205 bool sortInputs =
true)
206 { this->
setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
217 template <
class Po
intArray>
219 const PointArray& points,
222 bool sortInputs =
true)
234 template <
class ScalarContainerX,
class ScalarContainerY>
236 const ScalarContainerY& y,
239 bool sortInputs =
true)
250 template <
class Po
intContainer>
254 bool sortInputs =
true)
268 void set(Scalar x0, Scalar x1,
269 Scalar y0, Scalar y1,
270 Scalar m0, Scalar m1)
323 template <
class ScalarArrayX,
class ScalarArrayY>
325 const ScalarArrayX& x,
326 const ScalarArrayY& y,
327 Scalar m0, Scalar m1,
328 bool sortInputs =
true)
330 assert(nSamples > 1);
333 for (
size_t i = 0; i < nSamples; ++i) {
357 template <
class ScalarContainerX,
class ScalarContainerY>
359 const ScalarContainerY& y,
360 Scalar m0, Scalar m1,
361 bool sortInputs =
true)
363 assert(x.size() == y.size());
364 assert(x.size() > 1);
368 std::copy(x.begin(), x.end(), xPos_.begin());
369 std::copy(y.begin(), y.end(), yPos_.begin());
391 template <
class Po
intArray>
393 const PointArray& points,
396 bool sortInputs =
true)
400 assert(nSamples > 1);
403 for (
size_t i = 0; i < nSamples; ++i) {
404 xPos_[i] = points[i][0];
405 yPos_[i] = points[i][1];
428 template <
class XYContainer>
432 bool sortInputs =
true)
436 assert(points.size() > 1);
439 typename XYContainer::const_iterator it = points.begin();
440 typename XYContainer::const_iterator endIt = points.end();
441 for (
size_t i = 0; it != endIt; ++i, ++it) {
470 template <
class XYContainer>
474 bool sortInputs =
true)
478 typename XYContainer::const_iterator it = points.begin();
479 typename XYContainer::const_iterator endIt = points.end();
480 for (
unsigned i = 0; it != endIt; ++i, ++it) {
481 xPos_[i] = std::get<0>(*it);
482 yPos_[i] = std::get<1>(*it);
511 template <
class ScalarArrayX,
class ScalarArrayY>
513 const ScalarArrayX& x,
514 const ScalarArrayY& y,
516 bool sortInputs =
true)
518 assert(nSamples > 1);
521 for (
size_t i = 0; i < nSamples; ++i) {
531 if (splineType == Periodic)
533 else if (splineType == Natural)
535 else if (splineType == Monotonic)
538 OPM_THROW(std::runtime_error,
"Spline type " << splineType <<
" not supported at this place");
552 template <
class ScalarContainerX,
class ScalarContainerY>
554 const ScalarContainerY& y,
556 bool sortInputs =
true)
558 assert(x.size() == y.size());
559 assert(x.size() > 1);
562 std::copy(x.begin(), x.end(), xPos_.begin());
563 std::copy(y.begin(), y.end(), yPos_.begin());
570 if (splineType == Periodic)
572 else if (splineType == Natural)
574 else if (splineType == Monotonic)
577 OPM_THROW(std::runtime_error,
"Spline type " << splineType <<
" not supported at this place");
592 template <
class Po
intArray>
594 const PointArray& points,
596 bool sortInputs =
true)
600 assert(nSamples > 1);
603 for (
size_t i = 0; i < nSamples; ++i) {
604 xPos_[i] = points[i][0];
605 yPos_[i] = points[i][1];
613 if (splineType == Periodic)
615 else if (splineType == Natural)
617 else if (splineType == Monotonic)
620 OPM_THROW(std::runtime_error,
"Spline type " << splineType <<
" not supported at this place");
634 template <
class XYContainer>
637 bool sortInputs =
true)
641 assert(points.size() > 1);
644 typename XYContainer::const_iterator it = points.begin();
645 typename XYContainer::const_iterator endIt = points.end();
646 for (
size_t i = 0; it != endIt; ++ i, ++it) {
656 if (splineType == Periodic)
658 else if (splineType == Natural)
660 else if (splineType == Monotonic)
663 OPM_THROW(std::runtime_error,
"Spline type " << splineType <<
" not supported at this place");
680 template <
class XYContainer>
683 bool sortInputs =
true)
687 typename XYContainer::const_iterator it = points.begin();
688 typename XYContainer::const_iterator endIt = points.end();
689 for (
unsigned i = 0; it != endIt; ++i, ++it) {
690 xPos_[i] = std::get<0>(*it);
691 yPos_[i] = std::get<1>(*it);
699 if (splineType == Periodic)
701 else if (splineType == Natural)
703 else if (splineType == Monotonic)
706 OPM_THROW(std::runtime_error,
"Spline type " << splineType <<
" not supported at this place");
712 template <
class Evaluation>
720 {
return xPos_.size(); }
725 Scalar
xAt(
size_t sampleIdx)
const 726 {
return x_(sampleIdx); }
732 {
return y_(sampleIdx); }
751 void printCSV(Scalar xi0, Scalar xi1,
size_t k, std::ostream& os = std::cout)
const 753 Scalar x0 = std::min(xi0, xi1);
754 Scalar x1 = std::max(xi0, xi1);
756 for (
size_t i = 0; i <= k; ++i) {
757 double x = i*(x1 - x0)/k + x0;
758 double x_p1 = x + (x1 - x0)/k;
765 y = (x -
x_(0))*dy_dx +
y_(0);
766 mono = (dy_dx>0)?1:-1;
768 else if (x >
x_(n)) {
770 y = (x -
x_(n))*dy_dx +
y_(n);
771 mono = (dy_dx>0)?1:-1;
774 OPM_THROW(std::runtime_error,
775 "The sampling points given to a spline must be sorted by their x value!");
784 os << x <<
" " << y <<
" " << dy_dx <<
" " << mono <<
"\n";
799 template <
class Evaluation>
800 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const 802 if (!extrapolate && !
applies(x))
803 OPM_THROW(Opm::NumericalProblem,
804 "Tried to evaluate a spline outside of its range");
809 Scalar m = evalDerivative_(
xAt(0), 0);
811 return y0 + m*(x -
xAt(0));
813 else if (x >
xAt(static_cast<size_t>(static_cast<long int>(
numSamples()) - 1))) {
814 Scalar m = evalDerivative_(
xAt(static_cast<size_t>(
numSamples() - 1)),
817 return y0 + m*(x -
xAt(static_cast<size_t>(
numSamples() - 1)));
821 return eval_(x, segmentIdx_(Opm::scalarValue(x)));
836 template <
class Evaluation>
839 if (!extrapolate && !
applies(x))
840 OPM_THROW(Opm::NumericalProblem,
841 "Tried to evaluate the derivative of a spline outside of its range");
846 return evalDerivative_(
xAt(0), 0);
851 return evalDerivative_(x, segmentIdx_(Opm::scalarValue(x)));
866 template <
class Evaluation>
869 if (!extrapolate && !
applies(x))
870 OPM_THROW(Opm::NumericalProblem,
871 "Tried to evaluate the second derivative of a spline outside of its range");
872 else if (extrapolate)
875 return evalDerivative2_(x, segmentIdx_(Opm::scalarValue(x)));
890 template <
class Evaluation>
893 if (!extrapolate && !
applies(x))
894 OPM_THROW(Opm::NumericalProblem,
895 "Tried to evaluate the third derivative of a spline outside of its range");
896 else if (extrapolate)
899 return evalDerivative3_(x, segmentIdx_(Opm::scalarValue(x)));
908 template <
class Evaluation>
912 const Evaluation& d)
const 923 template <
class Evaluation>
928 const Evaluation& d)
const 932 Evaluation tmpSol[3], sol = 0;
934 size_t iFirst = segmentIdx_(x0);
935 size_t iLast = segmentIdx_(x1);
936 for (
size_t i = iFirst; i <= iLast; ++i)
944 OPM_THROW(std::runtime_error,
945 "Spline has more than one intersection");
950 OPM_THROW(std::runtime_error,
951 "Spline has no intersection");
964 int monotonic(Scalar x0, Scalar x1,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const 966 assert(std::abs(x0 - x1) > 1e-30);
977 Scalar m = evalDerivative_(
xAt(0), 0);
978 if (std::abs(m) < 1e-20)
983 size_t i = segmentIdx_(x0);
984 if (
x_(i + 1) >= x1) {
987 monotonic_(i, x0, x1, r);
993 monotonic_(i, x0,
x_(i+1), r);
998 size_t iEnd = segmentIdx_(x1);
999 for (; i < iEnd - 1; ++i) {
1000 monotonic_(i,
x_(i),
x_(i + 1), r);
1009 assert(extrapolate);
1013 return (r < 0 || r==3)?-1:0;
1015 return (r > 0 || r==3)?1:0;
1021 monotonic_(iEnd,
x_(iEnd), x1, r);
1043 bool operator ()(
unsigned idxA,
unsigned idxB)
const 1044 {
return x_.at(idxA) < x_.at(idxB); }
1046 const std::vector<Scalar>& x_;
1057 std::vector<unsigned> idxVector(n);
1058 for (
unsigned i = 0; i < n; ++i)
1064 std::sort(idxVector.begin(), idxVector.end(), cmp);
1067 std::vector<Scalar> tmpX(n), tmpY(n);
1068 for (
size_t i = 0; i < idxVector.size(); ++ i) {
1069 tmpX[i] = xPos_[idxVector[i]];
1070 tmpY[i] = yPos_[idxVector[i]];
1084 for (
unsigned i = 0; i <= (n - 1)/2; ++i) {
1085 std::swap(xPos_[i], xPos_[n - i - 1]);
1086 std::swap(yPos_[i], yPos_[n - i - 1]);
1095 xPos_.resize(nSamples);
1096 yPos_.resize(nSamples);
1097 slopeVec_.resize(nSamples);
1115 M.
solve(moments, d);
1136 M.
solve(moments, d);
1160 M.
solve(moments, d);
1163 for (
int i = static_cast<int>(
numSamples()) - 2; i >= 0; --i) {
1164 unsigned ui =
static_cast<unsigned>(i);
1165 moments[ui+1] = moments[ui];
1179 template <
class DestVector,
class SourceVector>
1182 const SourceVector& srcX,
1183 const SourceVector& srcY,
1186 assert(nSamples >= 2);
1190 for (
unsigned i = 0; i < nSamples; ++i) {
1192 if (srcX[0] > srcX[nSamples - 1])
1193 idx = nSamples - i - 1;
1194 destX[i] = srcX[idx];
1195 destY[i] = srcY[idx];
1199 template <
class DestVector,
class ListIterator>
1200 void assignFromArrayList_(DestVector& destX,
1202 const ListIterator& srcBegin,
1203 const ListIterator& srcEnd,
1206 assert(nSamples >= 2);
1209 ListIterator it = srcBegin;
1211 bool reverse =
false;
1212 if ((*srcBegin)[0] > (*it)[0])
1217 for (
unsigned i = 0; it != srcEnd; ++i, ++it) {
1220 idx = nSamples - i - 1;
1221 destX[i] = (*it)[0];
1222 destY[i] = (*it)[1];
1233 template <
class DestVector,
class ListIterator>
1236 ListIterator srcBegin,
1237 ListIterator srcEnd,
1240 assert(nSamples >= 2);
1246 ListIterator it = srcBegin;
1248 bool reverse =
false;
1249 if (std::get<0>(*srcBegin) > std::get<0>(*it))
1254 for (
unsigned i = 0; it != srcEnd; ++i, ++it) {
1257 idx = nSamples - i - 1;
1258 destX[i] = std::get<0>(*it);
1259 destY[i] = std::get<1>(*it);
1267 template <
class Vector,
class Matrix>
1275 d[0] = 6/
h_(1) * ( (
y_(1) -
y_(0))/
h_(1) - m0);
1284 (m1 - (
y_(n) -
y_(n - 1))/
h_(n));
1291 template <
class Vector,
class Matrix>
1301 for (
size_t i = 1; i < n; ++i) {
1302 Scalar lambda_i =
h_(i + 1) / (
h_(i) +
h_(i + 1));
1303 Scalar mu_i = 1 - lambda_i;
1305 6 / (
h_(i) +
h_(i + 1))
1307 ( (
y_(i + 1) -
y_(i))/
h_(i + 1) - (
y_(i) -
y_(i - 1))/
h_(i));
1311 M[i][i + 1] = lambda_i;
1316 Scalar lambda_0 = 0;
1337 template <
class Matrix,
class Vector>
1346 assert(M.
rows() == n);
1349 for (
size_t i = 2; i < n; ++i) {
1350 Scalar lambda_i =
h_(i + 1) / (
h_(i) +
h_(i + 1));
1351 Scalar mu_i = 1 - lambda_i;
1353 6 / (
h_(i) +
h_(i + 1))
1355 ( (
y_(i + 1) -
y_(i))/
h_(i + 1) - (
y_(i) -
y_(i - 1))/
h_(i));
1359 M[i-1][i] = lambda_i;
1363 Scalar lambda_n =
h_(1) / (
h_(n) +
h_(1));
1364 Scalar lambda_1 =
h_(2) / (
h_(1) +
h_(2));;
1365 Scalar mu_1 = 1 - lambda_1;
1366 Scalar mu_n = 1 - lambda_n;
1385 M[n-1][0] = lambda_n;
1399 template <
class Vector>
1405 std::vector<Scalar> delta(n);
1406 for (
size_t k = 0; k < n - 1; ++k)
1407 delta[k] = (
y_(k + 1) -
y_(k))/(
x_(k + 1) -
x_(k));
1410 for (
size_t k = 1; k < n - 1; ++k)
1411 slopes[k] = (delta[k - 1] + delta[k])/2;
1412 slopes[0] = delta[0];
1413 slopes[n - 1] = delta[n - 2];
1416 for (
size_t k = 0; k < n - 1; ++k) {
1417 if (std::abs(delta[k]) < 1e-50) {
1425 Scalar alpha = slopes[k] / delta[k];
1426 Scalar beta = slopes[k + 1] / delta[k];
1428 if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
1432 else if (alpha*alpha + beta*beta > 3*3) {
1433 Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
1434 slopes[k] = tau*alpha*delta[k];
1435 slopes[k + 1] = tau*beta*delta[k];
1448 template <
class MomentsVector,
class SlopeVector>
1459 Scalar h = this->
h_(n - 1);
1464 (
y_(n - 1) -
y_(n - 2))/h
1466 h/6*(moments[n-1] - moments[n - 2]);
1471 moments[n - 1] * x*x / (2 * h)
1477 for (
size_t i = 0; i < n - 1; ++ i) {
1480 Scalar h_i = this->
h_(i + 1);
1485 (
y_(i+1) -
y_(i))/h_i
1487 h_i/6*(moments[i+1] - moments[i]);
1490 - moments[i] * x_i1*x_i1 / (2 * h_i)
1497 slopes[n - 1] = mRight;
1503 template <
class Evaluation>
1504 Evaluation eval_(
const Evaluation& x,
size_t i)
const 1507 Scalar delta =
h_(i + 1);
1508 Evaluation t = (x -
x_(i))/delta;
1512 + h10_(t) *
slope_(i)*delta
1513 + h01_(t) *
y_(i + 1)
1514 + h11_(t) *
slope_(i + 1)*delta;
1519 template <
class Evaluation>
1520 Evaluation evalDerivative_(
const Evaluation& x,
size_t i)
const 1523 Scalar delta =
h_(i + 1);
1524 Evaluation t = (x -
x_(i))/delta;
1525 Evaluation alpha = 1 / delta;
1529 (h00_prime_(t) *
y_(i)
1530 + h10_prime_(t) *
slope_(i)*delta
1531 + h01_prime_(t) *
y_(i + 1)
1532 + h11_prime_(t) *
slope_(i + 1)*delta);
1537 template <
class Evaluation>
1538 Evaluation evalDerivative2_(
const Evaluation& x,
size_t i)
const 1541 Scalar delta =
h_(i + 1);
1542 Evaluation t = (x -
x_(i))/delta;
1543 Evaluation alpha = 1 / delta;
1547 *(h00_prime2_(t) *
y_(i)
1548 + h10_prime2_(t) *
slope_(i)*delta
1549 + h01_prime2_(t) *
y_(i + 1)
1550 + h11_prime2_(t) *
slope_(i + 1)*delta);
1555 template <
class Evaluation>
1556 Evaluation evalDerivative3_(
const Evaluation& x,
size_t i)
const 1559 Scalar delta =
h_(i + 1);
1560 Evaluation t = (x -
x_(i))/delta;
1561 Evaluation alpha = 1 / delta;
1565 *(h00_prime3_(t)*
y_(i)
1566 + h10_prime3_(t)*
slope_(i)*delta
1567 + h01_prime3_(t)*
y_(i + 1)
1568 + h11_prime3_(t)*
slope_(i + 1)*delta);
1572 template <
class Evaluation>
1573 Evaluation h00_(
const Evaluation& t)
const 1574 {
return (2*t - 3)*t*t + 1; }
1576 template <
class Evaluation>
1577 Evaluation h10_(
const Evaluation& t)
const 1578 {
return ((t - 2)*t + 1)*t; }
1580 template <
class Evaluation>
1581 Evaluation h01_(
const Evaluation& t)
const 1582 {
return (-2*t + 3)*t*t; }
1584 template <
class Evaluation>
1585 Evaluation h11_(
const Evaluation& t)
const 1586 {
return (t - 1)*t*t; }
1589 template <
class Evaluation>
1590 Evaluation h00_prime_(
const Evaluation& t)
const 1591 {
return (3*2*t - 2*3)*t; }
1593 template <
class Evaluation>
1594 Evaluation h10_prime_(
const Evaluation& t)
const 1595 {
return (3*t - 2*2)*t + 1; }
1597 template <
class Evaluation>
1598 Evaluation h01_prime_(
const Evaluation& t)
const 1599 {
return (-3*2*t + 2*3)*t; }
1601 template <
class Evaluation>
1602 Evaluation h11_prime_(
const Evaluation& t)
const 1603 {
return (3*t - 2)*t; }
1606 template <
class Evaluation>
1607 Evaluation h00_prime2_(
const Evaluation& t)
const 1608 {
return 2*3*2*t - 2*3; }
1610 template <
class Evaluation>
1611 Evaluation h10_prime2_(
const Evaluation& t)
const 1612 {
return 2*3*t - 2*2; }
1614 template <
class Evaluation>
1615 Evaluation h01_prime2_(
const Evaluation& t)
const 1616 {
return -2*3*2*t + 2*3; }
1618 template <
class Evaluation>
1619 Evaluation h11_prime2_(
const Evaluation& t)
const 1620 {
return 2*3*t - 2; }
1623 template <
class Evaluation>
1624 Scalar h00_prime3_(
const Evaluation& t OPM_UNUSED)
const 1627 template <
class Evaluation>
1628 Scalar h10_prime3_(
const Evaluation& t OPM_UNUSED)
const 1631 template <
class Evaluation>
1632 Scalar h01_prime3_(
const Evaluation& t OPM_UNUSED)
const 1635 template <
class Evaluation>
1636 Scalar h11_prime3_(
const Evaluation& t OPM_UNUSED)
const 1647 int monotonic_(
size_t i, Scalar x0, Scalar x1,
int& r)
const 1654 if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
1657 Scalar disc = b*b - 4*a*c;
1661 if (x0*(x0*a + b) + c > 0) {
1662 r = (r==3 || r == 1)?1:0;
1666 r = (r==3 || r == -1)?-1:0;
1670 disc = std::sqrt(disc);
1671 Scalar xE1 = (-b + disc)/(2*a);
1672 Scalar xE2 = (-b - disc)/(2*a);
1674 if (std::abs(disc) < 1e-30) {
1676 if (std::abs(xE1 - x0) < 1e-30)
1681 if (x0*(x0*a + b) + c > 0) {
1682 r = (r==3 || r == 1)?1:0;
1686 r = (r==3 || r == -1)?-1:0;
1690 if ((x0 < xE1 && xE1 < x1) ||
1691 (x0 < xE2 && xE2 < x1))
1700 if (x0*(x0*a + b) + c > 0) {
1701 r = (r==3 || r == 1)?1:0;
1705 r = (r==3 || r == -1)?-1:0;
1714 template <
class Evaluation>
1717 const Evaluation& a,
1718 const Evaluation& b,
1719 const Evaluation& c,
1720 const Evaluation& d,
1721 Scalar x0 = -1e30, Scalar x1 = 1e30)
const 1729 x0 = std::max(
x_(segIdx), x0);
1730 x1 = std::min(
x_(segIdx+1), x1);
1734 for (
unsigned j = 0; j < n; ++j) {
1735 if (x0 <= sol[j] && sol[j] <= x1) {
1744 size_t segmentIdx_(Scalar x)
const 1750 while (iLow + 1 < iHigh) {
1751 size_t i = (iLow + iHigh) / 2;
1763 Scalar
h_(
size_t i)
const 1765 assert(
x_(i) >
x_(i-1));
1767 return x_(i) -
x_(i - 1);
1773 Scalar
x_(
size_t i)
const 1774 {
return xPos_[i]; }
1779 Scalar
y_(
size_t i)
const 1780 {
return yPos_[i]; }
1787 {
return slopeVec_[i]; }
1791 Scalar a_(
size_t i)
const 1792 {
return evalDerivative3_(Scalar(0.0), i)/6.0; }
1796 Scalar b_(
size_t i)
const 1797 {
return evalDerivative2_(Scalar(0.0), i)/2.0; }
1801 Scalar c_(
size_t i)
const 1802 {
return evalDerivative_(Scalar(0.0), i); }
1806 Scalar d_(
size_t i)
const 1807 {
return eval_(Scalar(0.0), i); }
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using C-style arrays.
Definition: Spline.hpp:324
Spline(size_t nSamples, const ScalarArray &x, const ScalarArray &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:200
void assignSamplingPoints_(DestVector &destX, DestVector &destY, const SourceVector &srcX, const SourceVector &srcY, unsigned nSamples)
Set the sampling point vectors.
Definition: Spline.hpp:1180
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition: TridiagonalMatrix.hpp:714
void setNumSamples_(size_t nSamples)
Resizes the internal vectors to store the sample points.
Definition: Spline.hpp:1093
void makeMonotonicSpline_(Vector &slopes)
Create a monotonic spline from the already set sampling points.
Definition: Spline.hpp:1400
Scalar valueAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:731
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Spline.hpp:837
void sortInput_()
Sort the sample points in ascending order of their x value.
Definition: Spline.hpp:1052
SplineType
The type of the spline to be created.
Definition: Spline.hpp:103
Helper class needed to sort the input sampling points.
Definition: Spline.hpp:1037
void setArrayOfPoints(size_t nSamples, const PointArray &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a C-style array.
Definition: Spline.hpp:593
void makeFullSpline_(Scalar m0, Scalar m1)
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1105
Scalar y_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1779
Spline()
Default constructor for a spline.
Definition: Spline.hpp:115
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's third derivative at a given position.
Definition: Spline.hpp:891
Provides free functions to invert polynomials of degree 1, 2 and 3.
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using STL-compatible containers...
Definition: Spline.hpp:358
Scalar slope_(size_t i) const
Returns the slope (i.e.
Definition: Spline.hpp:1786
Scalar h_(size_t i) const
Returns x[i] - x[i - 1].
Definition: Spline.hpp:1763
Definition: Air_Mesitylene.hpp:33
size_t numSamples() const
Return the number of (x, y) values.
Definition: Spline.hpp:719
Scalar xAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:725
Scalar x_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1773
void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream &os=std::cout) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition: Spline.hpp:751
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the spline as interval.
Definition: Spline.hpp:1030
Spline(const PointContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:184
int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED=false) const
Returns 1 if the spline is monotonically increasing, -1 if the spline is mononously decreasing and 0 ...
Definition: Spline.hpp:964
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points natural spline using C-style arrays.
Definition: Spline.hpp:512
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Spline.hpp:713
void reverseSamplingPoints_()
Reverse order of the elements in the arrays which contain the sampling points.
Definition: Spline.hpp:1080
size_t rows() const
Return the number of rows of the matrix.
Definition: TridiagonalMatrix.hpp:140
Spline(const ScalarContainer &x, const ScalarContainer &y, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:171
Evaluation intersectInterval(Scalar x0, Scalar x1, const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d) const
Find the intersections of the spline with a cubic polynomial in a sub-interval of the spline...
Definition: Spline.hpp:924
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's second derivative at a given position.
Definition: Spline.hpp:867
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using STL-compatible containers.
Definition: Spline.hpp:553
size_t intersectSegment_(Evaluation *sol, size_t segIdx, const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d, Scalar x0=-1e30, Scalar x1=1e30) const
Find all the intersections of a segment of the spline with a cubic polynomial within a specified inte...
Definition: Spline.hpp:1715
Spline(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:142
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left...
void setContainerOfPoints(const XYContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a STL-compatible container of array-like objects...
Definition: Spline.hpp:635
void makePeriodicSpline_()
Create a periodic spline from the already set sampling points.
Definition: Spline.hpp:1147
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Spline.hpp:800
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left...
Definition: TridiagonalMatrix.hpp:50
Spline(const PointContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:251
Spline(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:235
Spline(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)
Convenience constructor for a full spline with just two sampling points.
Definition: Spline.hpp:128
void makeFullSystem_(Matrix &M, Vector &d, Scalar m0, Scalar m1)
Make the linear system of equations Mx = d which results in the moments of the full spline...
Definition: Spline.hpp:1268
void setContainerOfTuples(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:471
void setContainerOfTuples(const XYContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Set the sampling points of a natural spline using a STL-compatible container of tuple-like objects...
Definition: Spline.hpp:681
Evaluation intersect(const Evaluation &a, const Evaluation &b, const Evaluation &c, const Evaluation &d) const
Find the intersections of the spline with a cubic polynomial in the whole interval, throws Opm::MathError exception if there is more or less than one solution.
Definition: Spline.hpp:909
void assignFromTupleList_(DestVector &destX, DestVector &destY, ListIterator srcBegin, ListIterator srcEnd, unsigned nSamples)
Set the sampling points.
Definition: Spline.hpp:1234
Spline(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:218
void setArrayOfPoints(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a C-style array.
Definition: Spline.hpp:392
Spline(size_t nSamples, const PointArray &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:157
void makePeriodicSystem_(Matrix &M, Vector &d)
Make the linear system of equations Mx = d which results in the moments of the periodic spline...
Definition: Spline.hpp:1338
void makeNaturalSpline_()
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1126
Class implementing cubic splines.
Definition: Spline.hpp:92
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
Convert the moments at the sample points to slopes.
Definition: Spline.hpp:1449
void setContainerOfPoints(const XYContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Set the sampling points and the boundary slopes of a full spline using a STL-compatible container of ...
Definition: Spline.hpp:429
void makeNaturalSystem_(Matrix &M, Vector &d)
Make the linear system of equations Mx = d which results in the moments of the natural spline...
Definition: Spline.hpp:1292