27 #ifndef OPM_SPLINE_HPP
28 #define OPM_SPLINE_HPP
90 template<
class Scalar>
94 typedef std::vector<Scalar> Vector;
128 Scalar y0, Scalar y1,
129 Scalar m0, Scalar m1)
130 {
set(x0, x1, y0, y1, m0, m1); }
140 template <
class ScalarArrayX,
class ScalarArrayY>
142 const ScalarArrayX& x,
143 const ScalarArrayY& y,
145 bool sortInputs =
true)
146 { this->
setXYArrays(nSamples, x, y, splineType, sortInputs); }
155 template <
class Po
intArray>
157 const PointArray& points,
159 bool sortInputs =
true)
169 template <
class ScalarContainer>
171 const ScalarContainer& y,
173 bool sortInputs =
true)
182 template <
class Po
intContainer>
185 bool sortInputs =
true)
198 template <
class ScalarArray>
200 const ScalarArray& x,
201 const ScalarArray& y,
204 bool sortInputs =
true)
205 { this->
setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
216 template <
class Po
intArray>
218 const PointArray& points,
221 bool sortInputs =
true)
233 template <
class ScalarContainerX,
class ScalarContainerY>
235 const ScalarContainerY& y,
238 bool sortInputs =
true)
249 template <
class Po
intContainer>
253 bool sortInputs =
true)
267 void set(Scalar x0, Scalar x1,
268 Scalar y0, Scalar y1,
269 Scalar m0, Scalar m1)
322 template <
class ScalarArrayX,
class ScalarArrayY>
324 const ScalarArrayX& x,
325 const ScalarArrayY& y,
326 Scalar m0, Scalar m1,
327 bool sortInputs =
true)
329 assert(nSamples > 1);
332 for (
size_t i = 0; i < nSamples; ++i) {
356 template <
class ScalarContainerX,
class ScalarContainerY>
358 const ScalarContainerY& y,
359 Scalar m0, Scalar m1,
360 bool sortInputs =
true)
362 assert(x.size() == y.size());
363 assert(x.size() > 1);
367 std::copy(x.begin(), x.end(), xPos_.begin());
368 std::copy(y.begin(), y.end(), yPos_.begin());
390 template <
class Po
intArray>
392 const PointArray& points,
395 bool sortInputs =
true)
399 assert(nSamples > 1);
402 for (
size_t i = 0; i < nSamples; ++i) {
403 xPos_[i] = points[i][0];
404 yPos_[i] = points[i][1];
427 template <
class XYContainer>
431 bool sortInputs =
true)
435 assert(points.size() > 1);
438 typename XYContainer::const_iterator it = points.begin();
439 typename XYContainer::const_iterator endIt = points.end();
440 for (
size_t i = 0; it != endIt; ++i, ++it) {
469 template <
class XYContainer>
473 bool sortInputs =
true)
477 typename XYContainer::const_iterator it = points.begin();
478 typename XYContainer::const_iterator endIt = points.end();
479 for (
unsigned i = 0; it != endIt; ++i, ++it) {
480 xPos_[i] = std::get<0>(*it);
481 yPos_[i] = std::get<1>(*it);
510 template <
class ScalarArrayX,
class ScalarArrayY>
512 const ScalarArrayX& x,
513 const ScalarArrayY& y,
515 bool sortInputs =
true)
517 assert(nSamples > 1);
520 for (
size_t i = 0; i < nSamples; ++i) {
530 if (splineType == Periodic)
532 else if (splineType == Natural)
534 else if (splineType == Monotonic)
537 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
551 template <
class ScalarContainerX,
class ScalarContainerY>
553 const ScalarContainerY& y,
555 bool sortInputs =
true)
557 assert(x.size() == y.size());
558 assert(x.size() > 1);
561 std::copy(x.begin(), x.end(), xPos_.begin());
562 std::copy(y.begin(), y.end(), yPos_.begin());
569 if (splineType == Periodic)
571 else if (splineType == Natural)
573 else if (splineType == Monotonic)
576 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
591 template <
class Po
intArray>
593 const PointArray& points,
595 bool sortInputs =
true)
599 assert(nSamples > 1);
602 for (
size_t i = 0; i < nSamples; ++i) {
603 xPos_[i] = points[i][0];
604 yPos_[i] = points[i][1];
612 if (splineType == Periodic)
614 else if (splineType == Natural)
616 else if (splineType == Monotonic)
619 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
633 template <
class XYContainer>
636 bool sortInputs =
true)
640 assert(points.size() > 1);
643 typename XYContainer::const_iterator it = points.begin();
644 typename XYContainer::const_iterator endIt = points.end();
645 for (
size_t i = 0; it != endIt; ++ i, ++it) {
655 if (splineType == Periodic)
657 else if (splineType == Natural)
659 else if (splineType == Monotonic)
662 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
679 template <
class XYContainer>
682 bool sortInputs =
true)
686 typename XYContainer::const_iterator it = points.begin();
687 typename XYContainer::const_iterator endIt = points.end();
688 for (
unsigned i = 0; it != endIt; ++i, ++it) {
689 xPos_[i] = std::get<0>(*it);
690 yPos_[i] = std::get<1>(*it);
698 if (splineType == Periodic)
700 else if (splineType == Natural)
702 else if (splineType == Monotonic)
705 throw std::runtime_error(
"Spline type "+std::to_string(
int(splineType))+
" not supported at this place");
711 template <
class Evaluation>
719 {
return xPos_.size(); }
724 Scalar
xAt(
size_t sampleIdx)
const
725 {
return x_(sampleIdx); }
731 {
return y_(sampleIdx); }
750 void printCSV(Scalar xi0, Scalar xi1,
size_t k, std::ostream& os = std::cout)
const
752 Scalar x0 = std::min(xi0, xi1);
753 Scalar x1 = std::max(xi0, xi1);
755 for (
size_t i = 0; i <= k; ++i) {
756 double x = i*(x1 - x0)/k + x0;
757 double x_p1 = x + (x1 - x0)/k;
764 y = (x -
x_(0))*dy_dx +
y_(0);
765 mono = (dy_dx>0)?1:-1;
767 else if (x >
x_(n)) {
769 y = (x -
x_(n))*dy_dx +
y_(n);
770 mono = (dy_dx>0)?1:-1;
773 throw std::runtime_error(
"The sampling points given to a spline must be sorted by their x value!");
782 os << x <<
" " << y <<
" " << dy_dx <<
" " << mono <<
"\n";
797 template <
class Evaluation>
798 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const
800 if (!extrapolate && !
applies(x))
801 throw NumericalIssue(
"Tried to evaluate a spline outside of its range");
806 Scalar m = evalDerivative_(
xAt(0), 0);
808 return y0 + m*(x -
xAt(0));
810 else if (x >
xAt(
static_cast<size_t>(
static_cast<long int>(
numSamples()) - 1))) {
811 Scalar m = evalDerivative_(
xAt(
static_cast<size_t>(
numSamples() - 1)),
814 return y0 + m*(x -
xAt(
static_cast<size_t>(
numSamples() - 1)));
818 return eval_(x, segmentIdx_(scalarValue(x)));
833 template <
class Evaluation>
836 if (!extrapolate && !
applies(x))
837 throw NumericalIssue(
"Tried to evaluate the derivative of a spline outside of its range");
842 return evalDerivative_(
xAt(0), 0);
847 return evalDerivative_(x, segmentIdx_(scalarValue(x)));
862 template <
class Evaluation>
865 if (!extrapolate && !
applies(x))
866 throw NumericalIssue(
"Tried to evaluate the second derivative of a spline outside of its range");
867 else if (extrapolate)
870 return evalDerivative2_(x, segmentIdx_(scalarValue(x)));
885 template <
class Evaluation>
888 if (!extrapolate && !
applies(x))
889 throw NumericalIssue(
"Tried to evaluate the third derivative of a spline outside of its range");
890 else if (extrapolate)
893 return evalDerivative3_(x, segmentIdx_(scalarValue(x)));
902 template <
class Evaluation>
906 const Evaluation& d)
const
917 template <
class Evaluation>
922 const Evaluation& d)
const
926 Evaluation tmpSol[3], sol = 0;
928 size_t iFirst = segmentIdx_(x0);
929 size_t iLast = segmentIdx_(x1);
930 for (
size_t i = iFirst; i <= iLast; ++i)
938 throw std::runtime_error(
"Spline has more than one intersection");
943 throw std::runtime_error(
"Spline has no intersection");
956 int monotonic(Scalar x0, Scalar x1,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const
958 assert(std::abs(x0 - x1) > 1e-30);
969 Scalar m = evalDerivative_(
xAt(0), 0);
970 if (std::abs(m) < 1e-20)
975 size_t i = segmentIdx_(x0);
976 if (
x_(i + 1) >= x1) {
979 monotonic_(i, x0, x1, r);
985 monotonic_(i, x0,
x_(i+1), r);
990 size_t iEnd = segmentIdx_(x1);
991 for (; i < iEnd - 1; ++i) {
992 monotonic_(i,
x_(i),
x_(i + 1), r);
1001 assert(extrapolate);
1005 return (r < 0 || r==3)?-1:0;
1007 return (r > 0 || r==3)?1:0;
1013 monotonic_(iEnd,
x_(iEnd), x1, r);
1035 bool operator ()(
unsigned idxA,
unsigned idxB)
const
1036 {
return x_.at(idxA) < x_.at(idxB); }
1038 const std::vector<Scalar>& x_;
1049 std::vector<unsigned> idxVector(n);
1050 for (
unsigned i = 0; i < n; ++i)
1056 std::sort(idxVector.begin(), idxVector.end(), cmp);
1059 std::vector<Scalar> tmpX(n), tmpY(n);
1060 for (
size_t i = 0; i < idxVector.size(); ++ i) {
1061 tmpX[i] = xPos_[idxVector[i]];
1062 tmpY[i] = yPos_[idxVector[i]];
1076 for (
unsigned i = 0; i <= (n - 1)/2; ++i) {
1077 std::swap(xPos_[i], xPos_[n - i - 1]);
1078 std::swap(yPos_[i], yPos_[n - i - 1]);
1087 xPos_.resize(nSamples);
1088 yPos_.resize(nSamples);
1089 slopeVec_.resize(nSamples);
1107 M.
solve(moments, d);
1128 M.
solve(moments, d);
1152 M.
solve(moments, d);
1155 for (
int i =
static_cast<int>(
numSamples()) - 2; i >= 0; --i) {
1156 unsigned ui =
static_cast<unsigned>(i);
1157 moments[ui+1] = moments[ui];
1171 template <
class DestVector,
class SourceVector>
1174 const SourceVector& srcX,
1175 const SourceVector& srcY,
1178 assert(nSamples >= 2);
1182 for (
unsigned i = 0; i < nSamples; ++i) {
1184 if (srcX[0] > srcX[nSamples - 1])
1185 idx = nSamples - i - 1;
1186 destX[i] = srcX[idx];
1187 destY[i] = srcY[idx];
1191 template <
class DestVector,
class ListIterator>
1192 void assignFromArrayList_(DestVector& destX,
1194 const ListIterator& srcBegin,
1195 const ListIterator& srcEnd,
1198 assert(nSamples >= 2);
1201 ListIterator it = srcBegin;
1203 bool reverse =
false;
1204 if ((*srcBegin)[0] > (*it)[0])
1209 for (
unsigned i = 0; it != srcEnd; ++i, ++it) {
1212 idx = nSamples - i - 1;
1213 destX[i] = (*it)[0];
1214 destY[i] = (*it)[1];
1225 template <
class DestVector,
class ListIterator>
1228 ListIterator srcBegin,
1229 ListIterator srcEnd,
1232 assert(nSamples >= 2);
1238 ListIterator it = srcBegin;
1240 bool reverse =
false;
1241 if (std::get<0>(*srcBegin) > std::get<0>(*it))
1246 for (
unsigned i = 0; it != srcEnd; ++i, ++it) {
1249 idx = nSamples - i - 1;
1250 destX[i] = std::get<0>(*it);
1251 destY[i] = std::get<1>(*it);
1259 template <
class Vector,
class Matrix>
1267 d[0] = 6/
h_(1) * ( (
y_(1) -
y_(0))/
h_(1) - m0);
1276 (m1 - (
y_(n) -
y_(n - 1))/
h_(n));
1283 template <
class Vector,
class Matrix>
1293 for (
size_t i = 1; i < n; ++i) {
1294 Scalar lambda_i =
h_(i + 1) / (
h_(i) +
h_(i + 1));
1295 Scalar mu_i = 1 - lambda_i;
1297 6 / (
h_(i) +
h_(i + 1))
1299 ( (
y_(i + 1) -
y_(i))/
h_(i + 1) - (
y_(i) -
y_(i - 1))/
h_(i));
1303 M[i][i + 1] = lambda_i;
1308 Scalar lambda_0 = 0;
1329 template <
class Matrix,
class Vector>
1338 assert(M.
rows() == n);
1341 for (
size_t i = 2; i < n; ++i) {
1342 Scalar lambda_i =
h_(i + 1) / (
h_(i) +
h_(i + 1));
1343 Scalar mu_i = 1 - lambda_i;
1345 6 / (
h_(i) +
h_(i + 1))
1347 ( (
y_(i + 1) -
y_(i))/
h_(i + 1) - (
y_(i) -
y_(i - 1))/
h_(i));
1351 M[i-1][i] = lambda_i;
1355 Scalar lambda_n =
h_(1) / (
h_(n) +
h_(1));
1356 Scalar lambda_1 =
h_(2) / (
h_(1) +
h_(2));;
1357 Scalar mu_1 = 1 - lambda_1;
1358 Scalar mu_n = 1 - lambda_n;
1377 M[n-1][0] = lambda_n;
1391 template <
class Vector>
1397 std::vector<Scalar> delta(n);
1398 for (
size_t k = 0; k < n - 1; ++k)
1399 delta[k] = (
y_(k + 1) -
y_(k))/(
x_(k + 1) -
x_(k));
1402 for (
size_t k = 1; k < n - 1; ++k)
1403 slopes[k] = (delta[k - 1] + delta[k])/2;
1404 slopes[0] = delta[0];
1405 slopes[n - 1] = delta[n - 2];
1408 for (
size_t k = 0; k < n - 1; ++k) {
1409 if (std::abs(delta[k]) < 1e-50) {
1417 Scalar alpha = slopes[k] / delta[k];
1418 Scalar beta = slopes[k + 1] / delta[k];
1420 if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
1424 else if (alpha*alpha + beta*beta > 3*3) {
1425 Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
1426 slopes[k] = tau*alpha*delta[k];
1427 slopes[k + 1] = tau*beta*delta[k];
1440 template <
class MomentsVector,
class SlopeVector>
1451 Scalar h = this->
h_(n - 1);
1456 (
y_(n - 1) -
y_(n - 2))/h
1458 h/6*(moments[n-1] - moments[n - 2]);
1463 moments[n - 1] * x*x / (2 * h)
1469 for (
size_t i = 0; i < n - 1; ++ i) {
1472 Scalar h_i = this->
h_(i + 1);
1477 (
y_(i+1) -
y_(i))/h_i
1479 h_i/6*(moments[i+1] - moments[i]);
1482 - moments[i] * x_i1*x_i1 / (2 * h_i)
1489 slopes[n - 1] = mRight;
1495 template <
class Evaluation>
1496 Evaluation eval_(
const Evaluation& x,
size_t i)
const
1499 Scalar delta =
h_(i + 1);
1500 Evaluation t = (x -
x_(i))/delta;
1504 + h10_(t) *
slope_(i)*delta
1505 + h01_(t) *
y_(i + 1)
1506 + h11_(t) *
slope_(i + 1)*delta;
1511 template <
class Evaluation>
1512 Evaluation evalDerivative_(
const Evaluation& x,
size_t i)
const
1515 Scalar delta =
h_(i + 1);
1516 Evaluation t = (x -
x_(i))/delta;
1517 Evaluation alpha = 1 / delta;
1521 (h00_prime_(t) *
y_(i)
1522 + h10_prime_(t) *
slope_(i)*delta
1523 + h01_prime_(t) *
y_(i + 1)
1524 + h11_prime_(t) *
slope_(i + 1)*delta);
1529 template <
class Evaluation>
1530 Evaluation evalDerivative2_(
const Evaluation& x,
size_t i)
const
1533 Scalar delta =
h_(i + 1);
1534 Evaluation t = (x -
x_(i))/delta;
1535 Evaluation alpha = 1 / delta;
1539 *(h00_prime2_(t) *
y_(i)
1540 + h10_prime2_(t) *
slope_(i)*delta
1541 + h01_prime2_(t) *
y_(i + 1)
1542 + h11_prime2_(t) *
slope_(i + 1)*delta);
1547 template <
class Evaluation>
1548 Evaluation evalDerivative3_(
const Evaluation& x,
size_t i)
const
1551 Scalar delta =
h_(i + 1);
1552 Evaluation t = (x -
x_(i))/delta;
1553 Evaluation alpha = 1 / delta;
1557 *(h00_prime3_(t)*
y_(i)
1558 + h10_prime3_(t)*
slope_(i)*delta
1559 + h01_prime3_(t)*
y_(i + 1)
1560 + h11_prime3_(t)*
slope_(i + 1)*delta);
1564 template <
class Evaluation>
1565 Evaluation h00_(
const Evaluation& t)
const
1566 {
return (2*t - 3)*t*t + 1; }
1568 template <
class Evaluation>
1569 Evaluation h10_(
const Evaluation& t)
const
1570 {
return ((t - 2)*t + 1)*t; }
1572 template <
class Evaluation>
1573 Evaluation h01_(
const Evaluation& t)
const
1574 {
return (-2*t + 3)*t*t; }
1576 template <
class Evaluation>
1577 Evaluation h11_(
const Evaluation& t)
const
1578 {
return (t - 1)*t*t; }
1581 template <
class Evaluation>
1582 Evaluation h00_prime_(
const Evaluation& t)
const
1583 {
return (3*2*t - 2*3)*t; }
1585 template <
class Evaluation>
1586 Evaluation h10_prime_(
const Evaluation& t)
const
1587 {
return (3*t - 2*2)*t + 1; }
1589 template <
class Evaluation>
1590 Evaluation h01_prime_(
const Evaluation& t)
const
1591 {
return (-3*2*t + 2*3)*t; }
1593 template <
class Evaluation>
1594 Evaluation h11_prime_(
const Evaluation& t)
const
1595 {
return (3*t - 2)*t; }
1598 template <
class Evaluation>
1599 Evaluation h00_prime2_(
const Evaluation& t)
const
1600 {
return 2*3*2*t - 2*3; }
1602 template <
class Evaluation>
1603 Evaluation h10_prime2_(
const Evaluation& t)
const
1604 {
return 2*3*t - 2*2; }
1606 template <
class Evaluation>
1607 Evaluation h01_prime2_(
const Evaluation& t)
const
1608 {
return -2*3*2*t + 2*3; }
1610 template <
class Evaluation>
1611 Evaluation h11_prime2_(
const Evaluation& t)
const
1612 {
return 2*3*t - 2; }
1615 template <
class Evaluation>
1616 Scalar h00_prime3_(
const Evaluation&)
const
1619 template <
class Evaluation>
1620 Scalar h10_prime3_(
const Evaluation&)
const
1623 template <
class Evaluation>
1624 Scalar h01_prime3_(
const Evaluation&)
const
1627 template <
class Evaluation>
1628 Scalar h11_prime3_(
const Evaluation&)
const
1639 int monotonic_(
size_t i, Scalar x0, Scalar x1,
int& r)
const
1646 if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
1649 Scalar disc = b*b - 4*a*c;
1653 if (x0*(x0*a + b) + c > 0) {
1654 r = (r==3 || r == 1)?1:0;
1658 r = (r==3 || r == -1)?-1:0;
1662 disc = std::sqrt(disc);
1663 Scalar xE1 = (-b + disc)/(2*a);
1664 Scalar xE2 = (-b - disc)/(2*a);
1666 if (std::abs(disc) < 1e-30) {
1668 if (std::abs(xE1 - x0) < 1e-30)
1673 if (x0*(x0*a + b) + c > 0) {
1674 r = (r==3 || r == 1)?1:0;
1678 r = (r==3 || r == -1)?-1:0;
1682 if ((x0 < xE1 && xE1 < x1) ||
1683 (x0 < xE2 && xE2 < x1))
1692 if (x0*(x0*a + b) + c > 0) {
1693 r = (r==3 || r == 1)?1:0;
1697 r = (r==3 || r == -1)?-1:0;
1706 template <
class Evaluation>
1709 const Evaluation& a,
1710 const Evaluation& b,
1711 const Evaluation& c,
1712 const Evaluation& d,
1713 Scalar x0 = -1e30, Scalar x1 = 1e30)
const
1721 x0 = std::max(
x_(segIdx), x0);
1722 x1 = std::min(
x_(segIdx+1), x1);
1726 for (
unsigned j = 0; j < n; ++j) {
1727 if (x0 <= sol[j] && sol[j] <= x1) {
1736 size_t segmentIdx_(Scalar x)
const
1742 while (iLow + 1 < iHigh) {
1743 size_t i = (iLow + iHigh) / 2;
1755 Scalar
h_(
size_t i)
const
1757 assert(
x_(i) >
x_(i-1));
1759 return x_(i) -
x_(i - 1);
1765 Scalar
x_(
size_t i)
const
1766 {
return xPos_[i]; }
1771 Scalar
y_(
size_t i)
const
1772 {
return yPos_[i]; }
1779 {
return slopeVec_[i]; }
1783 Scalar a_(
size_t i)
const
1784 {
return evalDerivative3_(Scalar(0.0), i)/6.0; }
1788 Scalar b_(
size_t i)
const
1789 {
return evalDerivative2_(Scalar(0.0), i)/2.0; }
1793 Scalar c_(
size_t i)
const
1794 {
return evalDerivative_(Scalar(0.0), i); }
1798 Scalar d_(
size_t i)
const
1799 {
return eval_(Scalar(0.0), i); }
Provides the opm-material specific exception classes.
Provides free functions to invert polynomials of degree 1, 2 and 3.
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition: PolynomialUtils.hpp:148
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Provides the OPM_UNUSED macro.
Definition: Exceptions.hpp:46
Class implementing cubic splines.
Definition: Spline.hpp:92
size_t numSamples() const
Return the number of (x, y) values.
Definition: Spline.hpp:718
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:1260
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:592
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:323
Scalar valueAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:730
void makeNaturalSpline_()
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1118
void sortInput_()
Sort the sample points in ascending order of their x value.
Definition: Spline.hpp:1044
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:141
Spline(const PointContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:250
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
Convert the moments at the sample points to slopes.
Definition: Spline.hpp:1441
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:470
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:1330
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Spline.hpp:712
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:170
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:127
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,...
Definition: Spline.hpp:903
Scalar y_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1771
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Spline.hpp:834
SplineType
The type of the spline to be created.
Definition: Spline.hpp:102
Spline(const PointContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition: Spline.hpp:183
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:1284
Scalar h_(size_t i) const
Returns x[i] - x[i - 1].
Definition: Spline.hpp:1755
void set(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)
Set the sampling points and the boundary slopes of the spline with two sampling points.
Definition: Spline.hpp:267
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:956
Spline(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:234
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:634
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's third derivative at a given position.
Definition: Spline.hpp:886
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:750
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:199
void assignSamplingPoints_(DestVector &destX, DestVector &destY, const SourceVector &srcX, const SourceVector &srcY, unsigned nSamples)
Set the sampling point vectors.
Definition: Spline.hpp:1172
Scalar xAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition: Spline.hpp:724
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the spline as interval.
Definition: Spline.hpp:1022
void assignFromTupleList_(DestVector &destX, DestVector &destY, ListIterator srcBegin, ListIterator srcEnd, unsigned nSamples)
Set the sampling points.
Definition: Spline.hpp:1226
Spline(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition: Spline.hpp:217
Scalar slope_(size_t i) const
Returns the slope (i.e.
Definition: Spline.hpp:1778
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Spline.hpp:798
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:1707
Scalar x_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition: Spline.hpp:1765
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:391
void makeFullSpline_(Scalar m0, Scalar m1)
Create a natural spline from the already set sampling points.
Definition: Spline.hpp:1097
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's second derivative at a given position.
Definition: Spline.hpp:863
void makeMonotonicSpline_(Vector &slopes)
Create a monotonic spline from the already set sampling points.
Definition: Spline.hpp:1392
void setNumSamples_(size_t nSamples)
Resizes the internal vectors to store the sample points.
Definition: Spline.hpp:1085
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:511
Spline()
Default constructor for a spline.
Definition: Spline.hpp:114
void reverseSamplingPoints_()
Reverse order of the elements in the arrays which contain the sampling points.
Definition: Spline.hpp:1072
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:680
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:552
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:357
void makePeriodicSpline_()
Create a periodic spline from the already set sampling points.
Definition: Spline.hpp:1139
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:428
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:918
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:156
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition: TridiagonalMatrix.hpp:51
size_t rows() const
Return the number of rows of the matrix.
Definition: TridiagonalMatrix.hpp:140
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition: TridiagonalMatrix.hpp:714
Helper class needed to sort the input sampling points.
Definition: Spline.hpp:1030