My Project
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.
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
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 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.
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
size_t rows() const
Return the number of rows of the matrix.
Definition: TridiagonalMatrix.hpp:139
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition: TridiagonalMatrix.hpp:713
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