My Project
Loading...
Searching...
No Matches
Spline.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef OPM_SPLINE_HPP
28#define OPM_SPLINE_HPP
29
31
34
35#include <iosfwd>
36#include <vector>
37
38namespace Opm
39{
89template<class Scalar>
90class Spline
91{
93 typedef std::vector<Scalar> Vector;
94
95public:
102 Full,
103 Natural,
104 Periodic,
105 Monotonic
106 };
107
114 { }
115
126 Spline(Scalar x0, Scalar x1,
127 Scalar y0, Scalar y1,
128 Scalar m0, Scalar m1)
129 { set(x0, x1, y0, y1, m0, m1); }
130
139 template <class ScalarArrayX, class ScalarArrayY>
140 Spline(size_t nSamples,
141 const ScalarArrayX& x,
142 const ScalarArrayY& y,
143 SplineType splineType = Natural,
144 bool sortInputs = true)
145 { this->setXYArrays(nSamples, x, y, splineType, sortInputs); }
146
154 template <class PointArray>
155 Spline(size_t nSamples,
156 const PointArray& points,
157 SplineType splineType = Natural,
158 bool sortInputs = true)
159 { this->setArrayOfPoints(nSamples, points, splineType, sortInputs); }
160
168 template <class ScalarContainer>
169 Spline(const ScalarContainer& x,
170 const ScalarContainer& y,
171 SplineType splineType = Natural,
172 bool sortInputs = true)
173 { this->setXYContainers(x, y, splineType, sortInputs); }
174
181 template <class PointContainer>
182 Spline(const PointContainer& points,
183 SplineType splineType = Natural,
184 bool sortInputs = true)
185 { this->setContainerOfPoints(points, splineType, sortInputs); }
186
197 template <class ScalarArray>
198 Spline(size_t nSamples,
199 const ScalarArray& x,
200 const ScalarArray& y,
201 Scalar m0,
202 Scalar m1,
203 bool sortInputs = true)
204 { this->setXYArrays(nSamples, x, y, m0, m1, sortInputs); }
205
215 template <class PointArray>
216 Spline(size_t nSamples,
217 const PointArray& points,
218 Scalar m0,
219 Scalar m1,
220 bool sortInputs = true)
221 { this->setArrayOfPoints(nSamples, points, m0, m1, sortInputs); }
222
232 template <class ScalarContainerX, class ScalarContainerY>
233 Spline(const ScalarContainerX& x,
234 const ScalarContainerY& y,
235 Scalar m0,
236 Scalar m1,
237 bool sortInputs = true)
238 { this->setXYContainers(x, y, m0, m1, sortInputs); }
239
248 template <class PointContainer>
249 Spline(const PointContainer& points,
250 Scalar m0,
251 Scalar m1,
252 bool sortInputs = true)
253 { this->setContainerOfPoints(points, m0, m1, sortInputs); }
254
266 void set(Scalar x0, Scalar x1,
267 Scalar y0, Scalar y1,
268 Scalar m0, Scalar m1)
269 {
270 slopeVec_.resize(2);
271 xPos_.resize(2);
272 yPos_.resize(2);
273
274 if (x0 > x1) {
275 xPos_[0] = x1;
276 xPos_[1] = x0;
277 yPos_[0] = y1;
278 yPos_[1] = y0;
279 }
280 else {
281 xPos_[0] = x0;
282 xPos_[1] = x1;
283 yPos_[0] = y0;
284 yPos_[1] = y1;
285 }
286
287 slopeVec_[0] = m0;
288 slopeVec_[1] = m1;
289
290 Matrix M(numSamples());
291 Vector d(numSamples());
292 Vector moments(numSamples());
293 this->makeFullSystem_(M, d, m0, m1);
294
295 // solve for the moments
296 M.solve(moments, d);
297
298 this->setSlopesFromMoments_(slopeVec_, moments);
299 }
300
301
305 // Full splines //
309
321 template <class ScalarArrayX, class ScalarArrayY>
322 void setXYArrays(size_t nSamples,
323 const ScalarArrayX& x,
324 const ScalarArrayY& y,
325 Scalar m0, Scalar m1,
326 bool sortInputs = true)
327 {
328 assert(nSamples > 1);
329
330 setNumSamples_(nSamples);
331 for (size_t i = 0; i < nSamples; ++i) {
332 xPos_[i] = x[i];
333 yPos_[i] = y[i];
334 }
335
336 if (sortInputs)
337 sortInput_();
338 else if (xPos_[0] > xPos_[numSamples() - 1])
340
341 makeFullSpline_(m0, m1);
342 }
343
355 template <class ScalarContainerX, class ScalarContainerY>
356 void setXYContainers(const ScalarContainerX& x,
357 const ScalarContainerY& y,
358 Scalar m0, Scalar m1,
359 bool sortInputs = true)
360 {
361 assert(x.size() == y.size());
362 assert(x.size() > 1);
363
364 setNumSamples_(x.size());
365
366 std::copy(x.begin(), x.end(), xPos_.begin());
367 std::copy(y.begin(), y.end(), yPos_.begin());
368
369 if (sortInputs)
370 sortInput_();
371 else if (xPos_[0] > xPos_[numSamples() - 1])
373
374 makeFullSpline_(m0, m1);
375 }
376
389 template <class PointArray>
390 void setArrayOfPoints(size_t nSamples,
391 const PointArray& points,
392 Scalar m0,
393 Scalar m1,
394 bool sortInputs = true)
395 {
396 // a spline with no or just one sampling points? what an
397 // incredible bad idea!
398 assert(nSamples > 1);
399
400 setNumSamples_(nSamples);
401 for (size_t i = 0; i < nSamples; ++i) {
402 xPos_[i] = points[i][0];
403 yPos_[i] = points[i][1];
404 }
405
406 if (sortInputs)
407 sortInput_();
408 else if (xPos_[0] > xPos_[numSamples() - 1])
410
411 makeFullSpline_(m0, m1);
412 }
413
426 template <class XYContainer>
427 void setContainerOfPoints(const XYContainer& points,
428 Scalar m0,
429 Scalar m1,
430 bool sortInputs = true)
431 {
432 // a spline with no or just one sampling points? what an
433 // incredible bad idea!
434 assert(points.size() > 1);
435
436 setNumSamples_(points.size());
437 typename XYContainer::const_iterator it = points.begin();
438 typename XYContainer::const_iterator endIt = points.end();
439 for (size_t i = 0; it != endIt; ++i, ++it) {
440 xPos_[i] = (*it)[0];
441 yPos_[i] = (*it)[1];
442 }
443
444 if (sortInputs)
445 sortInput_();
446 else if (xPos_[0] > xPos_[numSamples() - 1])
448
449 // make a full spline
450 makeFullSpline_(m0, m1);
451 }
452
468 template <class XYContainer>
469 void setContainerOfTuples(const XYContainer& points,
470 Scalar m0,
471 Scalar m1,
472 bool sortInputs = true)
473 {
474 // resize internal arrays
475 setNumSamples_(points.size());
476 typename XYContainer::const_iterator it = points.begin();
477 typename XYContainer::const_iterator endIt = points.end();
478 for (unsigned i = 0; it != endIt; ++i, ++it) {
479 xPos_[i] = std::get<0>(*it);
480 yPos_[i] = std::get<1>(*it);
481 }
482
483 if (sortInputs)
484 sortInput_();
485 else if (xPos_[0] > xPos_[numSamples() - 1])
487
488 // make a full spline
489 makeFullSpline_(m0, m1);
490 }
491
495 // Natural/Periodic splines //
499
509 template <class ScalarArrayX, class ScalarArrayY>
510 void setXYArrays(size_t nSamples,
511 const ScalarArrayX& x,
512 const ScalarArrayY& y,
513 SplineType splineType = Natural,
514 bool sortInputs = true)
515 {
516 assert(nSamples > 1);
517
518 setNumSamples_(nSamples);
519 for (size_t i = 0; i < nSamples; ++i) {
520 xPos_[i] = x[i];
521 yPos_[i] = y[i];
522 }
523
524 if (sortInputs)
525 sortInput_();
526 else if (xPos_[0] > xPos_[numSamples() - 1])
528
529 if (splineType == Periodic)
531 else if (splineType == Natural)
533 else if (splineType == Monotonic)
534 this->makeMonotonicSpline_(slopeVec_);
535 else
536 throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
537 }
538
550 template <class ScalarContainerX, class ScalarContainerY>
551 void setXYContainers(const ScalarContainerX& x,
552 const ScalarContainerY& y,
553 SplineType splineType = Natural,
554 bool sortInputs = true)
555 {
556 assert(x.size() == y.size());
557 assert(x.size() > 1);
558
559 setNumSamples_(x.size());
560 std::copy(x.begin(), x.end(), xPos_.begin());
561 std::copy(y.begin(), y.end(), yPos_.begin());
562
563 if (sortInputs)
564 sortInput_();
565 else if (xPos_[0] > xPos_[numSamples() - 1])
567
568 if (splineType == Periodic)
570 else if (splineType == Natural)
572 else if (splineType == Monotonic)
573 this->makeMonotonicSpline_(slopeVec_);
574 else
575 throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
576 }
577
590 template <class PointArray>
591 void setArrayOfPoints(size_t nSamples,
592 const PointArray& points,
593 SplineType splineType = Natural,
594 bool sortInputs = true)
595 {
596 // a spline with no or just one sampling points? what an
597 // incredible bad idea!
598 assert(nSamples > 1);
599
600 setNumSamples_(nSamples);
601 for (size_t i = 0; i < nSamples; ++i) {
602 xPos_[i] = points[i][0];
603 yPos_[i] = points[i][1];
604 }
605
606 if (sortInputs)
607 sortInput_();
608 else if (xPos_[0] > xPos_[numSamples() - 1])
610
611 if (splineType == Periodic)
613 else if (splineType == Natural)
615 else if (splineType == Monotonic)
616 this->makeMonotonicSpline_(slopeVec_);
617 else
618 throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
619 }
620
632 template <class XYContainer>
633 void setContainerOfPoints(const XYContainer& points,
634 SplineType splineType = Natural,
635 bool sortInputs = true)
636 {
637 // a spline with no or just one sampling points? what an
638 // incredible bad idea!
639 assert(points.size() > 1);
640
641 setNumSamples_(points.size());
642 typename XYContainer::const_iterator it = points.begin();
643 typename XYContainer::const_iterator endIt = points.end();
644 for (size_t i = 0; it != endIt; ++ i, ++it) {
645 xPos_[i] = (*it)[0];
646 yPos_[i] = (*it)[1];
647 }
648
649 if (sortInputs)
650 sortInput_();
651 else if (xPos_[0] > xPos_[numSamples() - 1])
653
654 if (splineType == Periodic)
656 else if (splineType == Natural)
658 else if (splineType == Monotonic)
659 this->makeMonotonicSpline_(slopeVec_);
660 else
661 throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
662 }
663
678 template <class XYContainer>
679 void setContainerOfTuples(const XYContainer& points,
680 SplineType splineType = Natural,
681 bool sortInputs = true)
682 {
683 // resize internal arrays
684 setNumSamples_(points.size());
685 typename XYContainer::const_iterator it = points.begin();
686 typename XYContainer::const_iterator endIt = points.end();
687 for (unsigned i = 0; it != endIt; ++i, ++it) {
688 xPos_[i] = std::get<0>(*it);
689 yPos_[i] = std::get<1>(*it);
690 }
691
692 if (sortInputs)
693 sortInput_();
694 else if (xPos_[0] > xPos_[numSamples() - 1])
696
697 if (splineType == Periodic)
699 else if (splineType == Natural)
701 else if (splineType == Monotonic)
702 this->makeMonotonicSpline_(slopeVec_);
703 else
704 throw std::runtime_error("Spline type "+std::to_string(int(splineType))+" not supported at this place");
705 }
706
710 template <class Evaluation>
711 bool applies(const Evaluation& x) const
712 { return x_(0) <= x && x <= x_(numSamples() - 1); }
713
717 size_t numSamples() const
718 { return xPos_.size(); }
719
723 Scalar xAt(size_t sampleIdx) const
724 { return x_(sampleIdx); }
725
729 Scalar valueAt(size_t sampleIdx) const
730 { return y_(sampleIdx); }
731
749 void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream& os) const;
750
762 template <class Evaluation>
763 Evaluation eval(const Evaluation& x, bool extrapolate = false) const
764 {
765 if (!extrapolate && !applies(x))
766 throw NumericalProblem("Tried to evaluate a spline outside of its range");
767
768 // handle extrapolation
769 if (extrapolate) {
770 if (x < xAt(0)) {
771 Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
772 Scalar y0 = y_(0);
773 return y0 + m*(x - xAt(0));
774 }
775 else if (x > xAt(static_cast<size_t>(static_cast<long int>(numSamples()) - 1))) {
776 Scalar m = evalDerivative_(xAt(static_cast<size_t>(numSamples() - 1)),
777 /*segmentIdx=*/static_cast<size_t>(numSamples()-2));
778 Scalar y0 = y_(static_cast<size_t>(numSamples() - 1));
779 return y0 + m*(x - xAt(static_cast<size_t>(numSamples() - 1)));
780 }
781 }
782
783 return eval_(x, segmentIdx_(scalarValue(x)));
784 }
785
798 template <class Evaluation>
799 Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
800 {
801 if (!extrapolate && !applies(x))
802 throw NumericalProblem("Tried to evaluate the derivative of a spline outside of its range");
803
804 // handle extrapolation
805 if (extrapolate) {
806 if (x < xAt(0))
807 return evalDerivative_(xAt(0), /*segmentIdx=*/0);
808 else if (x > xAt(numSamples() - 1))
809 return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
810 }
811
812 return evalDerivative_(x, segmentIdx_(scalarValue(x)));
813 }
814
827 template <class Evaluation>
828 Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
829 {
830 if (!extrapolate && !applies(x))
831 throw NumericalProblem("Tried to evaluate the second derivative of a spline outside of its range");
832 else if (extrapolate)
833 return 0.0;
834
835 return evalDerivative2_(x, segmentIdx_(scalarValue(x)));
836 }
837
850 template <class Evaluation>
851 Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
852 {
853 if (!extrapolate && !applies(x))
854 throw NumericalProblem("Tried to evaluate the third derivative of a spline outside of its range");
855 else if (extrapolate)
856 return 0.0;
857
858 return evalDerivative3_(x, segmentIdx_(scalarValue(x)));
859 }
860
867 template <class Evaluation>
868 Evaluation intersect(const Evaluation& a,
869 const Evaluation& b,
870 const Evaluation& c,
871 const Evaluation& d) const
872 {
873 return intersectInterval(xAt(0), xAt(numSamples() - 1), a, b, c, d);
874 }
875
882 template <class Evaluation>
883 Evaluation intersectInterval(Scalar x0, Scalar x1,
884 const Evaluation& a,
885 const Evaluation& b,
886 const Evaluation& c,
887 const Evaluation& d) const
888 {
889 assert(applies(x0) && applies(x1));
890
891 Evaluation tmpSol[3], sol = 0;
892 size_t nSol = 0;
893 size_t iFirst = segmentIdx_(x0);
894 size_t iLast = segmentIdx_(x1);
895 for (size_t i = iFirst; i <= iLast; ++i)
896 {
897 size_t nCur = intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
898 if (nCur == 1)
899 sol = tmpSol[0];
900
901 nSol += nCur;
902 if (nSol > 1) {
903 throw std::runtime_error("Spline has more than one intersection"); //<<a<<"x^3 + "<<b<"x^2 + "<<c<"x + "<<d);
904 }
905 }
906
907 if (nSol != 1)
908 throw std::runtime_error("Spline has no intersection"); //<<a<"x^3 + " <<b<"x^2 + "<<c<"x + "<<d<<"!");
909
910 return sol;
911 }
912
921 int monotonic(Scalar x0, Scalar x1,
922 [[maybe_unused]] bool extrapolate = false) const
923 {
924 assert(std::abs(x0 - x1) > 1e-30);
925
926 // make sure that x0 is smaller than x1
927 if (x0 > x1)
928 std::swap(x0, x1);
929
930 assert(x0 < x1);
931
932 int r = 3;
933 if (x0 < xAt(0)) {
934 assert(extrapolate);
935 Scalar m = evalDerivative_(xAt(0), /*segmentIdx=*/0);
936 if (std::abs(m) < 1e-20)
937 r = (m < 0)?-1:1;
938 x0 = xAt(0);
939 };
940
941 size_t i = segmentIdx_(x0);
942 if (x_(i + 1) >= x1) {
943 // interval is fully contained within a single spline
944 // segment
945 monotonic_(i, x0, x1, r);
946 return r;
947 }
948
949 // the first segment overlaps with the specified interval
950 // partially
951 monotonic_(i, x0, x_(i+1), r);
952 ++ i;
953
954 // make sure that the segments which are completly in the
955 // interval [x0, x1] all exhibit the same monotonicity.
956 size_t iEnd = segmentIdx_(x1);
957 for (; i < iEnd - 1; ++i) {
958 monotonic_(i, x_(i), x_(i + 1), r);
959 if (!r)
960 return 0;
961 }
962
963 // if the user asked for a part of the spline which is
964 // extrapolated, we need to check the slope at the spline's
965 // endpoint
966 if (x1 > xAt(numSamples() - 1)) {
967 assert(extrapolate);
968
969 Scalar m = evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
970 if (m < 0)
971 return (r < 0 || r==3)?-1:0;
972 else if (m > 0)
973 return (r > 0 || r==3)?1:0;
974
975 return r;
976 }
977
978 // check for the last segment
979 monotonic_(iEnd, x_(iEnd), x1, r);
980
981 return r;
982 }
983
988 int monotonic() const
989 { return monotonic(xAt(0), xAt(numSamples() - 1)); }
990
991protected:
996 {
997 ComparatorX_(const std::vector<Scalar>& x)
998 : x_(x)
999 {}
1000
1001 bool operator ()(unsigned idxA, unsigned idxB) const
1002 { return x_.at(idxA) < x_.at(idxB); }
1003
1004 const std::vector<Scalar>& x_;
1005 };
1006
1011 {
1012 size_t n = numSamples();
1013
1014 // create a vector containing 0...n-1
1015 std::vector<unsigned> idxVector(n);
1016 for (unsigned i = 0; i < n; ++i)
1017 idxVector[i] = i;
1018
1019 // sort the indices according to the x values of the sample
1020 // points
1021 ComparatorX_ cmp(xPos_);
1022 std::sort(idxVector.begin(), idxVector.end(), cmp);
1023
1024 // reorder the sample points
1025 std::vector<Scalar> tmpX(n), tmpY(n);
1026 for (size_t i = 0; i < idxVector.size(); ++ i) {
1027 tmpX[i] = xPos_[idxVector[i]];
1028 tmpY[i] = yPos_[idxVector[i]];
1029 }
1030 xPos_ = tmpX;
1031 yPos_ = tmpY;
1032 }
1033
1039 {
1040 // reverse the arrays
1041 size_t n = numSamples();
1042 for (unsigned i = 0; i <= (n - 1)/2; ++i) {
1043 std::swap(xPos_[i], xPos_[n - i - 1]);
1044 std::swap(yPos_[i], yPos_[n - i - 1]);
1045 }
1046 }
1047
1051 void setNumSamples_(size_t nSamples)
1052 {
1053 xPos_.resize(nSamples);
1054 yPos_.resize(nSamples);
1055 slopeVec_.resize(nSamples);
1056 }
1057
1063 void makeFullSpline_(Scalar m0, Scalar m1)
1064 {
1065 Matrix M(numSamples());
1066 std::vector<Scalar> d(numSamples());
1067 std::vector<Scalar> moments(numSamples());
1068
1069 // create linear system of equations
1070 this->makeFullSystem_(M, d, m0, m1);
1071
1072 // solve for the moments (-> second derivatives)
1073 M.solve(moments, d);
1074
1075 // convert the moments to slopes at the sample points
1076 this->setSlopesFromMoments_(slopeVec_, moments);
1077 }
1078
1085 {
1087 Vector d(numSamples());
1088 Vector moments(numSamples());
1089
1090 // create linear system of equations
1091 this->makeNaturalSystem_(M, d);
1092
1093 // solve for the moments (-> second derivatives)
1094 M.solve(moments, d);
1095
1096 // convert the moments to slopes at the sample points
1097 this->setSlopesFromMoments_(slopeVec_, moments);
1098 }
1099
1106 {
1107 Matrix M(numSamples() - 1);
1108 Vector d(numSamples() - 1);
1109 Vector moments(numSamples() - 1);
1110
1111 // create linear system of equations. This is a bit hacky,
1112 // because it assumes that std::vector internally stores its
1113 // data as a big C-style array, but it saves us from yet
1114 // another copy operation
1115 this->makePeriodicSystem_(M, d);
1116
1117 // solve for the moments (-> second derivatives)
1118 M.solve(moments, d);
1119
1120 moments.resize(numSamples());
1121 for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i) {
1122 unsigned ui = static_cast<unsigned>(i);
1123 moments[ui+1] = moments[ui];
1124 }
1125 moments[0] = moments[numSamples() - 1];
1126
1127 // convert the moments to slopes at the sample points
1128 this->setSlopesFromMoments_(slopeVec_, moments);
1129 }
1130
1137 template <class DestVector, class SourceVector>
1138 void assignSamplingPoints_(DestVector& destX,
1139 DestVector& destY,
1140 const SourceVector& srcX,
1141 const SourceVector& srcY,
1142 unsigned nSamples)
1143 {
1144 assert(nSamples >= 2);
1145
1146 // copy sample points, make sure that the first x value is
1147 // smaller than the last one
1148 for (unsigned i = 0; i < nSamples; ++i) {
1149 unsigned idx = i;
1150 if (srcX[0] > srcX[nSamples - 1])
1151 idx = nSamples - i - 1;
1152 destX[i] = srcX[idx];
1153 destY[i] = srcY[idx];
1154 }
1155 }
1156
1157 template <class DestVector, class ListIterator>
1158 void assignFromArrayList_(DestVector& destX,
1159 DestVector& destY,
1160 const ListIterator& srcBegin,
1161 const ListIterator& srcEnd,
1162 unsigned nSamples)
1163 {
1164 assert(nSamples >= 2);
1165
1166 // find out wether the x values are in reverse order
1167 ListIterator it = srcBegin;
1168 ++it;
1169 bool reverse = false;
1170 if ((*srcBegin)[0] > (*it)[0])
1171 reverse = true;
1172 --it;
1173
1174 // loop over all sampling points
1175 for (unsigned i = 0; it != srcEnd; ++i, ++it) {
1176 unsigned idx = i;
1177 if (reverse)
1178 idx = nSamples - i - 1;
1179 destX[idx] = (*it)[0];
1180 destY[idx] = (*it)[1];
1181 }
1182 }
1183
1191 template <class DestVector, class ListIterator>
1192 void assignFromTupleList_(DestVector& destX,
1193 DestVector& destY,
1194 ListIterator srcBegin,
1195 ListIterator srcEnd,
1196 unsigned nSamples)
1197 {
1198 assert(nSamples >= 2);
1199
1200 // copy sample points, make sure that the first x value is
1201 // smaller than the last one
1202
1203 // find out wether the x values are in reverse order
1204 ListIterator it = srcBegin;
1205 ++it;
1206 bool reverse = false;
1207 if (std::get<0>(*srcBegin) > std::get<0>(*it))
1208 reverse = true;
1209 --it;
1210
1211 // loop over all sampling points
1212 for (unsigned i = 0; it != srcEnd; ++i, ++it) {
1213 unsigned idx = i;
1214 if (reverse)
1215 idx = nSamples - i - 1;
1216 destX[idx] = std::get<0>(*it);
1217 destY[idx] = std::get<1>(*it);
1218 }
1219 }
1220
1225 template <class Vector, class Matrix>
1226 void makeFullSystem_(Matrix& M, Vector& d, Scalar m0, Scalar m1)
1227 {
1228 makeNaturalSystem_(M, d);
1229
1230 size_t n = numSamples() - 1;
1231 // first row
1232 M[0][1] = 1;
1233 d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0);
1234
1235 // last row
1236 M[n][n - 1] = 1;
1237
1238 // right hand side
1239 d[n] =
1240 6/h_(n)
1241 *
1242 (m1 - (y_(n) - y_(n - 1))/h_(n));
1243 }
1244
1249 template <class Vector, class Matrix>
1250 void makeNaturalSystem_(Matrix& M, Vector& d)
1251 {
1252 M = 0.0;
1253
1254 // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1255 // Springer, 2005, p. 111
1256 size_t n = numSamples() - 1;
1257
1258 // second to next to last rows
1259 for (size_t i = 1; i < n; ++i) {
1260 Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
1261 Scalar mu_i = 1 - lambda_i;
1262 Scalar d_i =
1263 6 / (h_(i) + h_(i + 1))
1264 *
1265 ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
1266
1267 M[i][i-1] = mu_i;
1268 M[i][i] = 2;
1269 M[i][i + 1] = lambda_i;
1270 d[i] = d_i;
1271 };
1272
1273 // See Stroer, equation (2.5.2.7)
1274 Scalar lambda_0 = 0;
1275 Scalar d_0 = 0;
1276
1277 Scalar mu_n = 0;
1278 Scalar d_n = 0;
1279
1280 // first row
1281 M[0][0] = 2;
1282 M[0][1] = lambda_0;
1283 d[0] = d_0;
1284
1285 // last row
1286 M[n][n-1] = mu_n;
1287 M[n][n] = 2;
1288 d[n] = d_n;
1289 }
1290
1295 template <class Matrix, class Vector>
1296 void makePeriodicSystem_(Matrix& M, Vector& d)
1297 {
1298 M = 0.0;
1299
1300 // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1301 // Springer, 2005, p. 111
1302 size_t n = numSamples() - 1;
1303
1304 assert(M.rows() == n);
1305
1306 // second to next to last rows
1307 for (size_t i = 2; i < n; ++i) {
1308 Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1));
1309 Scalar mu_i = 1 - lambda_i;
1310 Scalar d_i =
1311 6 / (h_(i) + h_(i + 1))
1312 *
1313 ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i));
1314
1315 M[i-1][i-2] = mu_i;
1316 M[i-1][i-1] = 2;
1317 M[i-1][i] = lambda_i;
1318 d[i-1] = d_i;
1319 };
1320
1321 Scalar lambda_n = h_(1) / (h_(n) + h_(1));
1322 Scalar lambda_1 = h_(2) / (h_(1) + h_(2));;
1323 Scalar mu_1 = 1 - lambda_1;
1324 Scalar mu_n = 1 - lambda_n;
1325
1326 Scalar d_1 =
1327 6 / (h_(1) + h_(2))
1328 *
1329 ( (y_(2) - y_(1))/h_(2) - (y_(1) - y_(0))/h_(1));
1330 Scalar d_n =
1331 6 / (h_(n) + h_(1))
1332 *
1333 ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n));
1334
1335
1336 // first row
1337 M[0][0] = 2;
1338 M[0][1] = lambda_1;
1339 M[0][n-1] = mu_1;
1340 d[0] = d_1;
1341
1342 // last row
1343 M[n-1][0] = lambda_n;
1344 M[n-1][n-2] = mu_n;
1345 M[n-1][n-1] = 2;
1346 d[n-1] = d_n;
1347 }
1348
1357 template <class Vector>
1358 void makeMonotonicSpline_(Vector& slopes)
1359 {
1360 auto n = numSamples();
1361
1362 // calculate the slopes of the secant lines
1363 std::vector<Scalar> delta(n);
1364 for (size_t k = 0; k < n - 1; ++k)
1365 delta[k] = (y_(k + 1) - y_(k))/(x_(k + 1) - x_(k));
1366
1367 // calculate the "raw" slopes at the sample points
1368 for (size_t k = 1; k < n - 1; ++k)
1369 slopes[k] = (delta[k - 1] + delta[k])/2;
1370 slopes[0] = delta[0];
1371 slopes[n - 1] = delta[n - 2];
1372
1373 // post-process the "raw" slopes at the sample points
1374 for (size_t k = 0; k < n - 1; ++k) {
1375 if (std::abs(delta[k]) < 1e-50) {
1376 // make the spline flat if the inputs are equal
1377 slopes[k] = 0;
1378 slopes[k + 1] = 0;
1379 ++ k;
1380 continue;
1381 }
1382 else {
1383 Scalar alpha = slopes[k] / delta[k];
1384 Scalar beta = slopes[k + 1] / delta[k];
1385
1386 if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
1387 slopes[k] = 0;
1388 }
1389 // limit (alpha, beta) to a circle of radius 3
1390 else if (alpha*alpha + beta*beta > 3*3) {
1391 Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta);
1392 slopes[k] = tau*alpha*delta[k];
1393 slopes[k + 1] = tau*beta*delta[k];
1394 }
1395 }
1396 }
1397 }
1398
1406 template <class MomentsVector, class SlopeVector>
1407 void setSlopesFromMoments_(SlopeVector& slopes, const MomentsVector& moments)
1408 {
1409 size_t n = numSamples();
1410
1411 // evaluate slope at the rightmost point.
1412 // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1413 // Springer, 2005, p. 109
1414 Scalar mRight;
1415
1416 {
1417 Scalar h = this->h_(n - 1);
1418 Scalar x = h;
1419 //Scalar x_1 = 0;
1420
1421 Scalar A =
1422 (y_(n - 1) - y_(n - 2))/h
1423 -
1424 h/6*(moments[n-1] - moments[n - 2]);
1425
1426 mRight =
1427 //- moments[n - 2] * x_1*x_1 / (2 * h)
1428 //+
1429 moments[n - 1] * x*x / (2 * h)
1430 +
1431 A;
1432 }
1433
1434 // evaluate the slope for the first n-1 sample points
1435 for (size_t i = 0; i < n - 1; ++ i) {
1436 // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
1437 // Springer, 2005, p. 109
1438 Scalar h_i = this->h_(i + 1);
1439 //Scalar x_i = 0;
1440 Scalar x_i1 = h_i;
1441
1442 Scalar A_i =
1443 (y_(i+1) - y_(i))/h_i
1444 -
1445 h_i/6*(moments[i+1] - moments[i]);
1446
1447 slopes[i] =
1448 - moments[i] * x_i1*x_i1 / (2 * h_i)
1449 +
1450 //moments[i + 1] * x_i*x_i / (2 * h_i)
1451 //+
1452 A_i;
1453
1454 }
1455 slopes[n - 1] = mRight;
1456 }
1457
1458
1459 // evaluate the spline at a given the position and given the
1460 // segment index
1461 template <class Evaluation>
1462 Evaluation eval_(const Evaluation& x, size_t i) const
1463 {
1464 // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1465 Scalar delta = h_(i + 1);
1466 Evaluation t = (x - x_(i))/delta;
1467
1468 return
1469 h00_(t) * y_(i)
1470 + h10_(t) * slope_(i)*delta
1471 + h01_(t) * y_(i + 1)
1472 + h11_(t) * slope_(i + 1)*delta;
1473 }
1474
1475 // evaluate the derivative of a spline given the actual position
1476 // and the segment index
1477 template <class Evaluation>
1478 Evaluation evalDerivative_(const Evaluation& x, size_t i) const
1479 {
1480 // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1481 Scalar delta = h_(i + 1);
1482 Evaluation t = (x - x_(i))/delta;
1483 Evaluation alpha = 1 / delta;
1484
1485 return
1486 alpha *
1487 (h00_prime_(t) * y_(i)
1488 + h10_prime_(t) * slope_(i)*delta
1489 + h01_prime_(t) * y_(i + 1)
1490 + h11_prime_(t) * slope_(i + 1)*delta);
1491 }
1492
1493 // evaluate the second derivative of a spline given the actual
1494 // position and the segment index
1495 template <class Evaluation>
1496 Evaluation evalDerivative2_(const Evaluation& x, size_t i) const
1497 {
1498 // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1499 Scalar delta = h_(i + 1);
1500 Evaluation t = (x - x_(i))/delta;
1501 Evaluation alpha = 1 / delta;
1502
1503 return
1504 alpha*alpha
1505 *(h00_prime2_(t) * y_(i)
1506 + h10_prime2_(t) * slope_(i)*delta
1507 + h01_prime2_(t) * y_(i + 1)
1508 + h11_prime2_(t) * slope_(i + 1)*delta);
1509 }
1510
1511 // evaluate the third derivative of a spline given the actual
1512 // position and the segment index
1513 template <class Evaluation>
1514 Evaluation evalDerivative3_(const Evaluation& x, size_t i) const
1515 {
1516 // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline
1517 Scalar delta = h_(i + 1);
1518 Evaluation t = (x - x_(i))/delta;
1519 Evaluation alpha = 1 / delta;
1520
1521 return
1522 alpha*alpha*alpha
1523 *(h00_prime3_(t)*y_(i)
1524 + h10_prime3_(t)*slope_(i)*delta
1525 + h01_prime3_(t)*y_(i + 1)
1526 + h11_prime3_(t)*slope_(i + 1)*delta);
1527 }
1528
1529 // hermite basis functions
1530 template <class Evaluation>
1531 Evaluation h00_(const Evaluation& t) const
1532 { return (2*t - 3)*t*t + 1; }
1533
1534 template <class Evaluation>
1535 Evaluation h10_(const Evaluation& t) const
1536 { return ((t - 2)*t + 1)*t; }
1537
1538 template <class Evaluation>
1539 Evaluation h01_(const Evaluation& t) const
1540 { return (-2*t + 3)*t*t; }
1541
1542 template <class Evaluation>
1543 Evaluation h11_(const Evaluation& t) const
1544 { return (t - 1)*t*t; }
1545
1546 // first derivative of the hermite basis functions
1547 template <class Evaluation>
1548 Evaluation h00_prime_(const Evaluation& t) const
1549 { return (3*2*t - 2*3)*t; }
1550
1551 template <class Evaluation>
1552 Evaluation h10_prime_(const Evaluation& t) const
1553 { return (3*t - 2*2)*t + 1; }
1554
1555 template <class Evaluation>
1556 Evaluation h01_prime_(const Evaluation& t) const
1557 { return (-3*2*t + 2*3)*t; }
1558
1559 template <class Evaluation>
1560 Evaluation h11_prime_(const Evaluation& t) const
1561 { return (3*t - 2)*t; }
1562
1563 // second derivative of the hermite basis functions
1564 template <class Evaluation>
1565 Evaluation h00_prime2_(const Evaluation& t) const
1566 { return 2*3*2*t - 2*3; }
1567
1568 template <class Evaluation>
1569 Evaluation h10_prime2_(const Evaluation& t) const
1570 { return 2*3*t - 2*2; }
1571
1572 template <class Evaluation>
1573 Evaluation h01_prime2_(const Evaluation& t) const
1574 { return -2*3*2*t + 2*3; }
1575
1576 template <class Evaluation>
1577 Evaluation h11_prime2_(const Evaluation& t) const
1578 { return 2*3*t - 2; }
1579
1580 // third derivative of the hermite basis functions
1581 template <class Evaluation>
1582 Scalar h00_prime3_(const Evaluation&) const
1583 { return 2*3*2; }
1584
1585 template <class Evaluation>
1586 Scalar h10_prime3_(const Evaluation&) const
1587 { return 2*3; }
1588
1589 template <class Evaluation>
1590 Scalar h01_prime3_(const Evaluation&) const
1591 { return -2*3*2; }
1592
1593 template <class Evaluation>
1594 Scalar h11_prime3_(const Evaluation&) const
1595 { return 2*3; }
1596
1597 // returns the monotonicality of an interval of a spline segment
1598 //
1599 // The return value have the following meaning:
1600 //
1601 // 3: spline is constant within interval [x0, x1]
1602 // 1: spline is monotonously increasing in the specified interval
1603 // 0: spline is not monotonic (or constant) in the specified interval
1604 // -1: spline is monotonously decreasing in the specified interval
1605 int monotonic_(size_t i, Scalar x0, Scalar x1, int& r) const
1606 {
1607 // coefficients of derivative in monomial basis
1608 Scalar a = 3*a_(i);
1609 Scalar b = 2*b_(i);
1610 Scalar c = c_(i);
1611
1612 if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20)
1613 return 3; // constant in interval, r stays unchanged!
1614
1615 Scalar disc = b*b - 4*a*c;
1616 if (disc < 0) {
1617 // discriminant of derivative is smaller than 0, i.e. the
1618 // segment's derivative does not exhibit any extrema.
1619 if (x0*(x0*a + b) + c > 0) {
1620 r = (r==3 || r == 1)?1:0;
1621 return 1;
1622 }
1623 else {
1624 r = (r==3 || r == -1)?-1:0;
1625 return -1;
1626 }
1627 }
1628 disc = std::sqrt(disc);
1629 Scalar xE1 = (-b + disc)/(2*a);
1630 Scalar xE2 = (-b - disc)/(2*a);
1631
1632 if (std::abs(disc) < 1e-30) {
1633 // saddle point -> no extrema
1634 if (std::abs(xE1 - x0) < 1e-30)
1635 // make sure that we're not picking the saddle point
1636 // to determine whether we're monotonically increasing
1637 // or decreasing
1638 x0 = x1;
1639 if (x0*(x0*a + b) + c > 0) {
1640 r = (r==3 || r == 1)?1:0;
1641 return 1;
1642 }
1643 else {
1644 r = (r==3 || r == -1)?-1:0;
1645 return -1;
1646 }
1647 };
1648 if ((x0 < xE1 && xE1 < x1) ||
1649 (x0 < xE2 && xE2 < x1))
1650 {
1651 // there is an extremum in the range (x0, x1)
1652 r = 0;
1653 return 0;
1654 }
1655 // no extremum in range (x0, x1)
1656 x0 = (x0 + x1)/2; // pick point in the middle of the interval
1657 // to avoid extrema on the boundaries
1658 if (x0*(x0*a + b) + c > 0) {
1659 r = (r==3 || r == 1)?1:0;
1660 return 1;
1661 }
1662 else {
1663 r = (r==3 || r == -1)?-1:0;
1664 return -1;
1665 }
1666 }
1667
1672 template <class Evaluation>
1673 size_t intersectSegment_(Evaluation* sol,
1674 size_t segIdx,
1675 const Evaluation& a,
1676 const Evaluation& b,
1677 const Evaluation& c,
1678 const Evaluation& d,
1679 Scalar x0 = -1e30, Scalar x1 = 1e30) const
1680 {
1681 unsigned n =
1683 a_(segIdx) - a,
1684 b_(segIdx) - b,
1685 c_(segIdx) - c,
1686 d_(segIdx) - d);
1687 x0 = std::max(x_(segIdx), x0);
1688 x1 = std::min(x_(segIdx+1), x1);
1689
1690 // filter the intersections outside of the specified interval
1691 size_t k = 0;
1692 for (unsigned j = 0; j < n; ++j) {
1693 if (x0 <= sol[j] && sol[j] <= x1) {
1694 sol[k] = sol[j];
1695 ++k;
1696 }
1697 }
1698 return k;
1699 }
1700
1701 // find the segment index for a given x coordinate
1702 size_t segmentIdx_(Scalar x) const
1703 {
1704 // bisection
1705 size_t iLow = 0;
1706 size_t iHigh = numSamples() - 1;
1707
1708 while (iLow + 1 < iHigh) {
1709 size_t i = (iLow + iHigh) / 2;
1710 if (x_(i) > x)
1711 iHigh = i;
1712 else
1713 iLow = i;
1714 };
1715 return iLow;
1716 }
1717
1721 Scalar h_(size_t i) const
1722 {
1723 assert(x_(i) > x_(i-1)); // the sampling points must be given
1724 // in ascending order
1725 return x_(i) - x_(i - 1);
1726 }
1727
1731 Scalar x_(size_t i) const
1732 { return xPos_[i]; }
1733
1737 Scalar y_(size_t i) const
1738 { return yPos_[i]; }
1739
1744 Scalar slope_(size_t i) const
1745 { return slopeVec_[i]; }
1746
1747 // returns the coefficient in front of the x^3 term. In Stoer this
1748 // is delta.
1749 Scalar a_(size_t i) const
1750 { return evalDerivative3_(/*x=*/Scalar(0.0), i)/6.0; }
1751
1752 // returns the coefficient in front of the x^2 term In Stoer this
1753 // is gamma.
1754 Scalar b_(size_t i) const
1755 { return evalDerivative2_(/*x=*/Scalar(0.0), i)/2.0; }
1756
1757 // returns the coefficient in front of the x^1 term. In Stoer this
1758 // is beta.
1759 Scalar c_(size_t i) const
1760 { return evalDerivative_(/*x=*/Scalar(0.0), i); }
1761
1762 // returns the coefficient in front of the x^0 term. In Stoer this
1763 // is alpha.
1764 Scalar d_(size_t i) const
1765 { return eval_(/*x=*/Scalar(0.0), i); }
1766
1767 Vector xPos_;
1768 Vector yPos_;
1769 Vector slopeVec_;
1770};
1771}
1772
1773#endif
Provides the OPM specific exception classes.
Provides free functions to invert polynomials of degree 1, 2 and 3.
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition Exceptions.hpp:40
Class implementing cubic splines.
Definition Spline.hpp:91
size_t numSamples() const
Return the number of (x, y) values.
Definition Spline.hpp:717
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:1226
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:591
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:322
Scalar valueAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition Spline.hpp:729
void makeNaturalSpline_()
Create a natural spline from the already set sampling points.
Definition Spline.hpp:1084
void sortInput_()
Sort the sample points in ascending order of their x value.
Definition Spline.hpp:1010
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:140
void printCSV(Scalar xi0, Scalar xi1, size_t k, std::ostream &os) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition Spline.cpp:33
Spline(const PointContainer &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition Spline.hpp:249
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
Convert the moments at the sample points to slopes.
Definition Spline.hpp:1407
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:469
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:1296
int monotonic(Scalar x0, Scalar x1, bool extrapolate=false) const
Returns 1 if the spline is monotonically increasing, -1 if the spline is mononously decreasing and 0 ...
Definition Spline.hpp:921
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition Spline.hpp:711
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:169
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:126
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:868
Scalar y_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition Spline.hpp:1737
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition Spline.hpp:799
SplineType
The type of the spline to be created.
Definition Spline.hpp:101
Spline(const PointContainer &points, SplineType splineType=Natural, bool sortInputs=true)
Convenience constructor for a natural or a periodic spline.
Definition Spline.hpp:182
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:1250
Scalar h_(size_t i) const
Returns x[i] - x[i - 1].
Definition Spline.hpp:1721
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:266
Spline(const ScalarContainerX &x, const ScalarContainerY &y, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition Spline.hpp:233
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:633
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's third derivative at a given position.
Definition Spline.hpp:851
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:198
void assignSamplingPoints_(DestVector &destX, DestVector &destY, const SourceVector &srcX, const SourceVector &srcY, unsigned nSamples)
Set the sampling point vectors.
Definition Spline.hpp:1138
Scalar xAt(size_t sampleIdx) const
Return the x value of a given sampling point.
Definition Spline.hpp:723
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the spline as interval.
Definition Spline.hpp:988
void assignFromTupleList_(DestVector &destX, DestVector &destY, ListIterator srcBegin, ListIterator srcEnd, unsigned nSamples)
Set the sampling points.
Definition Spline.hpp:1192
Spline(size_t nSamples, const PointArray &points, Scalar m0, Scalar m1, bool sortInputs=true)
Convenience constructor for a full spline.
Definition Spline.hpp:216
Scalar slope_(size_t i) const
Returns the slope (i.e.
Definition Spline.hpp:1744
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition Spline.hpp:763
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:1673
Scalar x_(size_t i) const
Returns the y coordinate of the i-th sampling point.
Definition Spline.hpp:1731
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:390
void makeFullSpline_(Scalar m0, Scalar m1)
Create a natural spline from the already set sampling points.
Definition Spline.hpp:1063
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's second derivative at a given position.
Definition Spline.hpp:828
void makeMonotonicSpline_(Vector &slopes)
Create a monotonic spline from the already set sampling points.
Definition Spline.hpp:1358
void setNumSamples_(size_t nSamples)
Resizes the internal vectors to store the sample points.
Definition Spline.hpp:1051
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:510
Spline()
Default constructor for a spline.
Definition Spline.hpp:113
void reverseSamplingPoints_()
Reverse order of the elements in the arrays which contain the sampling points.
Definition Spline.hpp:1038
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:679
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:551
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:356
void makePeriodicSpline_()
Create a periodic spline from the already set sampling points.
Definition Spline.hpp:1105
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:427
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:883
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:155
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition TridiagonalMatrix.hpp:50
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
unsigned invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition PolynomialUtils.hpp:148
Helper class needed to sort the input sampling points.
Definition Spline.hpp:996