27 #ifndef OPM_TABULATED_1D_FUNCTION_HPP
28 #define OPM_TABULATED_1D_FUNCTION_HPP
45 template <
class Scalar>
64 template <
class ScalarArrayX,
class ScalarArrayY>
66 const ScalarArrayX& x,
67 const ScalarArrayY& y,
68 bool sortInputs =
true)
79 template <
class ScalarContainer>
81 const ScalarContainer& y,
82 bool sortInputs =
true)
91 template <
class Po
intContainer>
93 bool sortInputs =
true)
101 template <
class ScalarArrayX,
class ScalarArrayY>
103 const ScalarArrayX& x,
104 const ScalarArrayY& y,
105 bool sortInputs =
true)
107 assert(nSamples > 1);
109 resizeArrays_(nSamples);
110 for (
size_t i = 0; i < nSamples; ++i) {
117 else if (xValues_[0] > xValues_[
numSamples() - 1])
118 reverseSamplingPoints_();
126 template <
class ScalarContainerX,
class ScalarContainerY>
128 const ScalarContainerY& y,
129 bool sortInputs =
true)
131 assert(x.size() == y.size());
133 resizeArrays_(x.size());
135 std::copy(x.begin(), x.end(), xValues_.begin());
136 std::copy(y.begin(), y.end(), yValues_.begin());
140 else if (xValues_[0] > xValues_[
numSamples() - 1])
141 reverseSamplingPoints_();
148 template <
class Po
intArray>
150 const PointArray& points,
151 bool sortInputs =
true)
155 assert(nSamples > 1);
157 resizeArrays_(nSamples);
158 for (
size_t i = 0; i < nSamples; ++i) {
159 xValues_[i] = points[i][0];
160 yValues_[i] = points[i][1];
165 else if (xValues_[0] > xValues_[
numSamples() - 1])
166 reverseSamplingPoints_();
183 template <
class XYContainer>
185 bool sortInputs =
true)
189 assert(points.size() > 1);
191 resizeArrays_(points.size());
192 typename XYContainer::const_iterator it = points.begin();
193 typename XYContainer::const_iterator endIt = points.end();
194 for (
unsigned i = 0; it != endIt; ++i, ++it) {
195 xValues_[i] = std::get<0>(*it);
196 yValues_[i] = std::get<1>(*it);
201 else if (xValues_[0] > xValues_[
numSamples() - 1])
202 reverseSamplingPoints_();
209 {
return xValues_.size(); }
215 {
return xValues_[0]; }
226 Scalar
xAt(
size_t i)
const
227 {
return xValues_[i]; }
229 const std::vector<Scalar>& xValues()
const
232 const std::vector<Scalar>& yValues()
const
239 {
return yValues_[i]; }
244 template <
class Evaluation>
246 {
return xValues_[0] <= x && x <= xValues_[
numSamples() - 1]; }
257 template <
class Evaluation>
258 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const
260 size_t segIdx = findSegmentIndex_(x, extrapolate);
262 Scalar x0 = xValues_[segIdx];
263 Scalar x1 = xValues_[segIdx + 1];
265 Scalar y0 = yValues_[segIdx];
266 Scalar y1 = yValues_[segIdx + 1];
268 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
282 template <
class Evaluation>
285 unsigned segIdx = findSegmentIndex_(x, extrapolate);
286 return evalDerivative_(x, segIdx);
303 template <
class Evaluation>
321 template <
class Evaluation>
333 int monotonic(Scalar x0, Scalar x1,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const
350 size_t i = findSegmentIndex_(x0, extrapolate);
351 if (xValues_[i + 1] >= x1) {
354 updateMonotonicity_(i, r);
360 updateMonotonicity_(i, r);
365 size_t iEnd = findSegmentIndex_(x1, extrapolate);
366 for (; i < iEnd - 1; ++i) {
367 updateMonotonicity_(i, r);
380 return (r < 0 || r==3)?-1:0;
382 return (r > 0 || r==3)?1:0;
388 updateMonotonicity_(iEnd, r);
416 void printCSV(Scalar xi0, Scalar xi1,
unsigned k, std::ostream& os = std::cout)
const
418 Scalar x0 = std::min(xi0, xi1);
419 Scalar x1 = std::max(xi0, xi1);
421 for (
unsigned i = 0; i <= k; ++i) {
422 double x = i*(x1 - x0)/k + x0;
426 if (x < xValues_[0]) {
428 y = (x - xValues_[0])*dy_dx + yValues_[0];
430 else if (x > xValues_[n]) {
432 y = (x - xValues_[n])*dy_dx + yValues_[n];
435 throw std::runtime_error(
"The sampling points given to a function must be sorted by their x value!");
443 os << x <<
" " << y <<
" " << dy_dx <<
"\n";
448 return xValues_ == data.xValues_ &&
449 yValues_ == data.yValues_;
453 template <
class Evaluation>
454 size_t findSegmentIndex_(
const Evaluation& x,
bool extrapolate =
false)
const
456 if (!extrapolate && !
applies(x))
457 throw NumericalIssue(
"Tried to evaluate a tabulated function outside of its range");
460 assert(xValues_.size() >= 2);
462 if (x <= xValues_[1])
464 else if (x >= xValues_[xValues_.size() - 2])
465 return xValues_.size() - 2;
469 size_t upperIdx = xValues_.size() - 2;
470 while (lowerIdx + 1 < upperIdx) {
471 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
472 if (x < xValues_[pivotIdx])
478 assert(xValues_[lowerIdx] <= x);
479 assert(x <= xValues_[lowerIdx + 1]);
484 template <
class Evaluation>
485 Evaluation evalDerivative_(
const Evaluation& x,
size_t segIdx)
const
487 Scalar x0 = xValues_[segIdx];
488 Scalar x1 = xValues_[segIdx + 1];
490 Scalar y0 = yValues_[segIdx];
491 Scalar y1 = yValues_[segIdx + 1];
493 Evaluation ret = blank(x);
494 ret = (y1 - y0)/(x1 - x0);
506 int updateMonotonicity_(
size_t i,
int& r)
const
508 if (yValues_[i] < yValues_[i + 1]) {
510 if (r == 3 || r == 1)
516 else if (yValues_[i] > yValues_[i + 1]) {
518 if (r == 3 || r == -1)
533 ComparatorX_(
const std::vector<Scalar>& x)
537 bool operator ()(
size_t idxA,
size_t idxB)
const
538 {
return x_.at(idxA) < x_.at(idxB); }
540 const std::vector<Scalar>& x_;
551 std::vector<unsigned> idxVector(n);
552 for (
unsigned i = 0; i < n; ++i)
557 ComparatorX_ cmp(xValues_);
558 std::sort(idxVector.begin(), idxVector.end(), cmp);
561 std::vector<Scalar> tmpX(n), tmpY(n);
562 for (
size_t i = 0; i < idxVector.size(); ++ i) {
563 tmpX[i] = xValues_[idxVector[i]];
564 tmpY[i] = yValues_[idxVector[i]];
574 void reverseSamplingPoints_()
578 for (
size_t i = 0; i <= (n - 1)/2; ++i) {
579 std::swap(xValues_[i], xValues_[n - i - 1]);
580 std::swap(yValues_[i], yValues_[n - i - 1]);
587 void resizeArrays_(
size_t nSamples)
589 xValues_.resize(nSamples);
590 yValues_.resize(nSamples);
593 std::vector<Scalar> xValues_;
594 std::vector<Scalar> yValues_;
Provides the opm-material specific exception classes.
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Provides the OPM_UNUSED macro.
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:258
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:208
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Tabulated1DFunction.hpp:283
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:149
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:220
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:80
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:214
Evaluation evalThirdDerivative(const Evaluation &, bool=false) const
Evaluate the function's third derivative at a given position.
Definition: Tabulated1DFunction.hpp:322
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:333
Evaluation evalSecondDerivative(const Evaluation &, bool=false) const
Evaluate the function's second derivative at a given position.
Definition: Tabulated1DFunction.hpp:304
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:65
void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream &os=std::cout) const
Prints k tuples of the format (x, y, dx/dy, isMonotonic) to stdout.
Definition: Tabulated1DFunction.hpp:416
void setContainerOfTuples(const XYContainer &points, bool sortInputs=true)
Set the sampling points of the piecewise linear function using a STL-compatible container of tuple-li...
Definition: Tabulated1DFunction.hpp:184
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:226
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:92
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:54
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval.
Definition: Tabulated1DFunction.hpp:397
void setXYArrays(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:102
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:127
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:245
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:238