25 #ifndef OPM_TABULATED_1D_FUNCTION_HPP
26 #define OPM_TABULATED_1D_FUNCTION_HPP
28 #include <opm/common/ErrorMacros.hpp>
29 #include <opm/common/Exceptions.hpp>
43 template <
class Scalar>
62 template <
class ScalarArrayX,
class ScalarArrayY>
64 const ScalarArrayX &x,
65 const ScalarArrayY &y,
66 bool sortInputs =
false)
77 template <
class ScalarContainer>
79 const ScalarContainer &y,
80 bool sortInputs =
false)
89 template <
class Po
intContainer>
91 bool sortInputs =
false)
99 template <
class ScalarArrayX,
class ScalarArrayY>
101 const ScalarArrayX &x,
102 const ScalarArrayY &y,
103 bool sortInputs =
false)
105 assert(nSamples > 1);
107 resizeArrays_(nSamples);
108 for (
size_t i = 0; i < nSamples; ++i) {
115 else if (xValues_[0] > xValues_[
numSamples() - 1])
116 reverseSamplingPoints_();
124 template <
class ScalarContainerX,
class ScalarContainerY>
126 const ScalarContainerY &y,
127 bool sortInputs =
false)
129 assert(x.size() == y.size());
130 assert(x.size() > 1);
132 resizeArrays_(x.size());
133 std::copy(x.begin(), x.end(), xValues_.begin());
134 std::copy(y.begin(), y.end(), yValues_.begin());
138 else if (xValues_[0] > xValues_[
numSamples() - 1])
139 reverseSamplingPoints_();
145 template <
class Po
intArray>
147 const PointArray &points,
148 bool sortInputs =
false)
152 assert(nSamples > 1);
154 resizeArrays_(nSamples);
155 for (
size_t i = 0; i < nSamples; ++i) {
156 xValues_[i] = points[i][0];
157 yValues_[i] = points[i][1];
162 else if (xValues_[0] > xValues_[
numSamples() - 1])
163 reverseSamplingPoints_();
180 template <
class XYContainer>
182 bool sortInputs =
false)
186 assert(points.size() > 1);
188 resizeArrays_(points.size());
189 typename XYContainer::const_iterator it = points.begin();
190 typename XYContainer::const_iterator endIt = points.end();
191 for (
int i = 0; it != endIt; ++i, ++it) {
192 xValues_[i] = std::get<0>(*it);
193 yValues_[i] = std::get<1>(*it);
198 else if (xValues_[0] > xValues_[
numSamples() - 1])
199 reverseSamplingPoints_();
206 {
return xValues_.size(); }
212 {
return xValues_[0]; }
223 Scalar
xAt(
size_t i)
const
224 {
return xValues_[i]; }
230 {
return yValues_[i]; }
236 {
return xValues_[0] <= x && x <= xValues_[
numSamples() - 1]; }
247 Scalar
eval(Scalar x,
bool extrapolate=
false)
const
250 if (extrapolate && x < xValues_.front())
252 else if (extrapolate && x > xValues_.back())
255 assert(xValues_.front() <= x && x <= xValues_.back());
256 segIdx = findSegmentIndex_(x);
259 Scalar x0 = xValues_[segIdx];
260 Scalar x1 = xValues_[segIdx + 1];
262 Scalar y0 = yValues_[segIdx];
263 Scalar y1 = yValues_[segIdx + 1];
265 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
277 template <
class Evaluation>
278 Evaluation
eval(
const Evaluation& x,
bool extrapolate=
false)
const
281 if (extrapolate && x.value < xValues_.front())
283 else if (extrapolate && x.value > xValues_.back())
286 assert(xValues_.front() <= x.value && x.value <= xValues_.back());
287 segIdx = findSegmentIndex_(x.value);
290 Scalar x0 = xValues_[segIdx];
291 Scalar x1 = xValues_[segIdx + 1];
293 Scalar y0 = yValues_[segIdx];
294 Scalar y1 = yValues_[segIdx + 1];
296 Scalar m = (y1 - y0)/(x1 - x0);
299 result.value = y0 + m*(x.value - x0);
300 for (
unsigned varIdx = 0; varIdx < result.derivatives.size(); ++varIdx)
301 result.derivatives[varIdx] = m*x.derivatives[varIdx];
319 int segIdx = findSegmentIndex_(x);
340 assert(extrapolate ||
applies(x));
360 assert(extrapolate ||
applies(x));
389 size_t i = findSegmentIndex_(x0);
390 if (xValues_[i + 1] >= x1) {
393 updateMonotonicity_(i, r);
399 updateMonotonicity_(i, r);
404 size_t iEnd = findSegmentIndex_(x1);
405 for (; i < iEnd - 1; ++i) {
406 updateMonotonicity_(i, r);
419 return (r < 0 || r==3)?-1:0;
421 return (r > 0 || r==3)?1:0;
427 updateMonotonicity_(iEnd, r);
455 void printCSV(Scalar xi0, Scalar xi1,
int k, std::ostream &os = std::cout)
const
460 for (
int i = 0; i <= k; ++i) {
461 double x = i*(x1 - x0)/k + x0;
465 if (x < xValues_[0]) {
467 y = (x - xValues_[0])*dy_dx + yValues_[0];
469 else if (x > xValues_[n]) {
471 y = (x - xValues_[n])*dy_dx + yValues_[n];
474 OPM_THROW(std::runtime_error,
475 "The sampling points given to a function must be sorted by their x value!");
483 os << x <<
" " << y <<
" " << dy_dx <<
"\n";
488 size_t findSegmentIndex_(Scalar x)
const
491 assert(xValues_.size() >= 2);
493 if (x <= xValues_[1])
495 else if (x >= xValues_[xValues_.size() - 2])
496 return xValues_.size() - 2;
499 size_t segmentIdx = 1;
500 size_t upperIdx = xValues_.size() - 2;
501 while (segmentIdx + 1 < upperIdx) {
502 size_t pivotIdx = (segmentIdx + upperIdx) / 2;
503 if (x < xValues_[pivotIdx])
506 segmentIdx = pivotIdx;
509 assert(xValues_[segmentIdx] <= x);
510 assert(x <= xValues_[segmentIdx + 1]);
515 Scalar evalDerivative_(
OPM_UNUSED Scalar x,
int segIdx)
const
517 Scalar x0 = xValues_[segIdx];
518 Scalar x1 = xValues_[segIdx + 1];
520 Scalar y0 = yValues_[segIdx];
521 Scalar y1 = yValues_[segIdx + 1];
523 return (y1 - y0)/(x1 - x0);
534 int updateMonotonicity_(
int i,
int &r)
const
536 if (yValues_[i] < yValues_[i + 1]) {
538 if (r == 3 || r == 1)
544 else if (yValues_[i] > yValues_[i + 1]) {
546 if (r == 3 || r == -1)
561 ComparatorX_(
const std::vector<Scalar> &x)
565 bool operator ()(
size_t idxA,
size_t idxB)
const
566 {
return x_.at(idxA) < x_.at(idxB); }
568 const std::vector<Scalar> &x_;
579 std::vector<unsigned> idxVector(n);
580 for (
unsigned i = 0; i < n; ++i)
585 ComparatorX_ cmp(xValues_);
586 std::sort(idxVector.begin(), idxVector.end(), cmp);
589 std::vector<Scalar> tmpX(n), tmpY(n);
590 for (
size_t i = 0; i < idxVector.size(); ++ i) {
591 tmpX[i] = xValues_[idxVector[i]];
592 tmpY[i] = yValues_[idxVector[i]];
602 void reverseSamplingPoints_()
606 for (
size_t i = 0; i <= (n - 1)/2; ++i) {
607 std::swap(xValues_[i], xValues_[n - i - 1]);
608 std::swap(yValues_[i], yValues_[n - i - 1]);
615 void resizeArrays_(
size_t nSamples)
617 xValues_.resize(nSamples);
618 yValues_.resize(nSamples);
621 std::vector<Scalar> xValues_;
622 std::vector<Scalar> yValues_;
void printCSV(Scalar xi0, Scalar xi1, int k, std::ostream &os=std::cout) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition: Tabulated1DFunction.hpp:455
void setContainerOfTuples(const XYContainer &points, bool sortInputs=false)
Set the sampling points of the piecewise linear function using a STL-compatible container of tuple-li...
Definition: Tabulated1DFunction.hpp:181
Definition: Air_Mesitylene.hpp:31
Scalar evalDerivative(Scalar x, bool=false) const
Evaluate the spline's derivative at a given position.
Definition: Tabulated1DFunction.hpp:317
Scalar eval(Scalar x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:247
bool applies(Scalar x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:235
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=false)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:100
Evaluation< Scalar, VarSetTag, numVars > max(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:114
int monotonic(Scalar x0, Scalar x1, bool extrapolate OPM_OPTIM_UNUSED=false) const
Returns 1 if the function is monotonically increasing, -1 if the function is mononously decreasing an...
Definition: Tabulated1DFunction.hpp:372
Scalar evalThirdDerivative(OPM_OPTIM_UNUSED Scalar x, OPM_OPTIM_UNUSED bool extrapolate=false) const
Evaluate the function's third derivative at a given position.
Definition: Tabulated1DFunction.hpp:358
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:211
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=false)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:125
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval. ...
Definition: Tabulated1DFunction.hpp:436
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:278
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:223
Evaluation< Scalar, VarSetTag, numVars > min(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:61
#define OPM_UNUSED
Definition: Unused.hpp:32
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:217
Scalar evalSecondDerivative(OPM_OPTIM_UNUSED Scalar x, OPM_OPTIM_UNUSED bool extrapolate=false) const
Evaluate the function's second derivative at a given position.
Definition: Tabulated1DFunction.hpp:338
#define OPM_OPTIM_UNUSED
Definition: Unused.hpp:42
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:205
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:229
Tabulated1DFunction(const PointContainer &points, bool sortInputs=false)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:90
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=false)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:146
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=false)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:63
Provides the OPM_UNUSED macro.
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:44
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=false)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:78
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:52