27#ifndef OPM_TABULATED_1D_FUNCTION_HPP
28#define OPM_TABULATED_1D_FUNCTION_HPP
30#include <opm/common/OpmLog/OpmLog.hpp>
46template <
class Scalar>
65 template <
class ScalarArrayX,
class ScalarArrayY>
67 const ScalarArrayX& x,
68 const ScalarArrayY& y,
69 bool sortInputs =
true)
80 template <
class ScalarContainer>
82 const ScalarContainer& y,
83 bool sortInputs =
true)
92 template <
class Po
intContainer>
94 bool sortInputs =
true)
102 template <
class ScalarArrayX,
class ScalarArrayY>
104 const ScalarArrayX& x,
105 const ScalarArrayY& y,
106 bool sortInputs =
true)
108 assert(nSamples > 1);
110 resizeArrays_(nSamples);
111 for (
size_t i = 0; i < nSamples; ++i) {
118 else if (xValues_[0] > xValues_[
numSamples() - 1])
119 reverseSamplingPoints_();
127 template <
class ScalarContainerX,
class ScalarContainerY>
129 const ScalarContainerY& y,
130 bool sortInputs =
true)
132 assert(x.size() == y.size());
134 resizeArrays_(x.size());
136 std::copy(x.begin(), x.end(), xValues_.begin());
137 std::copy(y.begin(), y.end(), yValues_.begin());
141 else if (xValues_[0] > xValues_[
numSamples() - 1])
142 reverseSamplingPoints_();
149 template <
class Po
intArray>
151 const PointArray& points,
152 bool sortInputs =
true)
156 assert(nSamples > 1);
158 resizeArrays_(nSamples);
159 for (
size_t i = 0; i < nSamples; ++i) {
160 xValues_[i] = points[i][0];
161 yValues_[i] = points[i][1];
166 else if (xValues_[0] > xValues_[
numSamples() - 1])
167 reverseSamplingPoints_();
184 template <
class XYContainer>
186 bool sortInputs =
true)
190 assert(points.size() > 1);
192 resizeArrays_(points.size());
193 typename XYContainer::const_iterator it = points.begin();
194 typename XYContainer::const_iterator endIt = points.end();
195 for (
unsigned i = 0; it != endIt; ++i, ++it) {
196 xValues_[i] = std::get<0>(*it);
197 yValues_[i] = std::get<1>(*it);
202 else if (xValues_[0] > xValues_[
numSamples() - 1])
203 reverseSamplingPoints_();
210 {
return xValues_.size(); }
216 {
return xValues_[0]; }
227 Scalar
xAt(
size_t i)
const
228 {
return xValues_[i]; }
240 {
return yValues_[i]; }
245 template <
class Evaluation>
247 {
return xValues_[0] <= x && x <= xValues_[
numSamples() - 1]; }
258 template <
class Evaluation>
259 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const
261 size_t segIdx = findSegmentIndex_(x, extrapolate);
263 Scalar x0 = xValues_[segIdx];
264 Scalar x1 = xValues_[segIdx + 1];
266 Scalar y0 = yValues_[segIdx];
267 Scalar y1 = yValues_[segIdx + 1];
269 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
283 template <
class Evaluation>
286 unsigned segIdx = findSegmentIndex_(x, extrapolate);
287 return evalDerivative_(x, segIdx);
304 template <
class Evaluation>
322 template <
class Evaluation>
335 [[maybe_unused]]
bool extrapolate =
false)
const
352 size_t i = findSegmentIndex_(x0, extrapolate);
353 if (xValues_[i + 1] >= x1) {
356 updateMonotonicity_(i, r);
362 updateMonotonicity_(i, r);
367 size_t iEnd = findSegmentIndex_(x1, extrapolate);
368 for (; i < iEnd - 1; ++i) {
369 updateMonotonicity_(i, r);
382 return (r < 0 || r==3)?-1:0;
384 return (r > 0 || r==3)?1:0;
390 updateMonotonicity_(iEnd, r);
418 void printCSV(Scalar xi0, Scalar xi1,
unsigned k, std::ostream& os = std::cout)
const
423 for (
unsigned i = 0; i <= k; ++i) {
424 double x = i*(x1 - x0)/k + x0;
428 if (x < xValues_[0]) {
430 y = (x - xValues_[0])*dy_dx + yValues_[0];
432 else if (x > xValues_[n]) {
434 y = (x - xValues_[n])*dy_dx + yValues_[n];
437 throw std::runtime_error(
"The sampling points given to a function must be sorted by their x value!");
445 os << x <<
" " << y <<
" " << dy_dx <<
"\n";
450 return xValues_ == data.xValues_ &&
451 yValues_ == data.yValues_;
455 template <
class Evaluation>
456 size_t findSegmentIndex_(
const Evaluation& x,
bool extrapolate =
false)
const
459 std::ostringstream sstream;
460 sstream <<
"We can not search for extrapolation/interpolation segment in an 1D table for non-finite value " <<
getValue(x) <<
" .";
461 throw std::runtime_error(sstream.str());
464 if (!extrapolate && !
applies(x))
465 throw std::logic_error(
"Trying to evaluate a tabulated function outside of its range");
469 std::ostringstream sstream;
470 sstream <<
"We need at least two sampling points to do interpolation/extrapolation,"
471 "and the table only contains {} sampling points" <<
numSamples();
472 throw std::logic_error(sstream.str());
475 if (x <= xValues_[1])
477 else if (x >= xValues_[xValues_.size() - 2])
478 return xValues_.size() - 2;
482 size_t upperIdx = xValues_.size() - 2;
483 while (lowerIdx + 1 < upperIdx) {
484 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
485 if (x < xValues_[pivotIdx])
491 if (xValues_[lowerIdx] > x || x > xValues_[lowerIdx + 1]) {
492 std::ostringstream sstream;
493 sstream <<
"Problematic interpolation/extrapolation segment is found for the input value " <<
Opm::getValue(x)
494 <<
"\nthe lower index of the found segment is " << lowerIdx <<
", the size of the table is " <<
numSamples()
495 <<
",\nand the end values of the found segment are " << xValues_[lowerIdx] <<
" and " << xValues_[lowerIdx + 1]
496 <<
", respectively.";
497 std::ostringstream sstream2;
498 sstream2 <<
" Outputting the problematic table for more information(with *** marking the found segment):";
504 sstream2 <<
" " << xValues_[i];
505 if (i == lowerIdx + 1)
509 OpmLog::debug(sstream.str() +
"\n" + sstream2.str());
510 throw std::runtime_error(sstream.str());
516 template <
class Evaluation>
517 Evaluation evalDerivative_(
const Evaluation& x,
size_t segIdx)
const
519 Scalar x0 = xValues_[segIdx];
520 Scalar x1 = xValues_[segIdx + 1];
522 Scalar y0 = yValues_[segIdx];
523 Scalar y1 = yValues_[segIdx + 1];
525 Evaluation ret =
blank(x);
526 ret = (y1 - y0)/(x1 - x0);
538 int updateMonotonicity_(
size_t i,
int& r)
const
540 if (yValues_[i] < yValues_[i + 1]) {
542 if (r == 3 || r == 1)
548 else if (yValues_[i] > yValues_[i + 1]) {
550 if (r == 3 || r == -1)
565 ComparatorX_(
const std::vector<Scalar>& x)
569 bool operator ()(
size_t idxA,
size_t idxB)
const
570 {
return x_.at(idxA) < x_.at(idxB); }
572 const std::vector<Scalar>& x_;
583 std::vector<unsigned> idxVector(n);
584 for (
unsigned i = 0; i < n; ++i)
589 ComparatorX_ cmp(xValues_);
590 std::sort(idxVector.begin(), idxVector.end(), cmp);
593 std::vector<Scalar> tmpX(n), tmpY(n);
594 for (
size_t i = 0; i < idxVector.size(); ++ i) {
595 tmpX[i] = xValues_[idxVector[i]];
596 tmpY[i] = yValues_[idxVector[i]];
606 void reverseSamplingPoints_()
610 for (
size_t i = 0; i <= (n - 1)/2; ++i) {
611 std::swap(xValues_[i], xValues_[n - i - 1]);
612 std::swap(yValues_[i], yValues_[n - i - 1]);
619 void resizeArrays_(
size_t nSamples)
621 xValues_.resize(nSamples);
622 yValues_.resize(nSamples);
625 std::vector<Scalar> xValues_;
626 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...
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:48
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:259
const std::vector< Scalar > & yValues() const
Definition: Tabulated1DFunction.hpp:233
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:209
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Tabulated1DFunction.hpp:284
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:150
bool operator==(const Tabulated1DFunction< Scalar > &data) const
Definition: Tabulated1DFunction.hpp:449
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:221
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:81
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:215
Evaluation evalThirdDerivative(const Evaluation &, bool=false) const
Evaluate the function's third derivative at a given position.
Definition: Tabulated1DFunction.hpp:323
Evaluation evalSecondDerivative(const Evaluation &, bool=false) const
Evaluate the function's second derivative at a given position.
Definition: Tabulated1DFunction.hpp:305
Tabulated1DFunction(size_t nSamples, const ScalarArrayX &x, const ScalarArrayY &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:66
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:418
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:185
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:227
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:93
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:55
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval.
Definition: Tabulated1DFunction.hpp:399
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:103
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:128
int monotonic(Scalar x0, Scalar x1, bool extrapolate=false) const
Returns 1 if the function is monotonically increasing, -1 if the function is mononously decreasing an...
Definition: Tabulated1DFunction.hpp:334
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:246
const std::vector< Scalar > & xValues() const
Definition: Tabulated1DFunction.hpp:230
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:239
Definition: Air_Mesitylene.hpp:34
bool isfinite(const Evaluation &value)
Definition: MathToolbox.hpp:420
Evaluation blank(const Evaluation &x)
Definition: MathToolbox.hpp:297
ReturnEval_< Evaluation1, Evaluation2 >::type min(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:346
ReturnEval_< Evaluation1, Evaluation2 >::type max(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:341
auto getValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::value(val))
Definition: MathToolbox.hpp:330