27 #ifndef OPM_TABULATED_1D_FUNCTION_HPP
28 #define OPM_TABULATED_1D_FUNCTION_HPP
31 #include <opm/common/ErrorMacros.hpp>
32 #include <opm/common/Exceptions.hpp>
33 #include <opm/common/Unused.hpp>
46 template <
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());
133 assert(x.size() > 1);
135 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_();
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]; }
233 {
return yValues_[i]; }
238 template <
class Evaluation>
240 {
return xValues_[0] <= x && x <= xValues_[
numSamples() - 1]; }
251 template <
class Evaluation>
252 Evaluation
eval(
const Evaluation& x,
bool extrapolate =
false)
const
254 if (!extrapolate && !
applies(x))
255 OPM_THROW(Opm::NumericalProblem,
256 "Tried to evaluate a tabulated function outside of its range");
260 if (extrapolate && x < xValues_.front())
262 else if (extrapolate && x > xValues_.back())
265 segIdx = findSegmentIndex_(Opm::scalarValue(x));
267 Scalar x0 = xValues_[segIdx];
268 Scalar x1 = xValues_[segIdx + 1];
270 Scalar y0 = yValues_[segIdx];
271 Scalar y1 = yValues_[segIdx + 1];
273 return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
287 template <
class Evaluation>
290 if (!extrapolate && !
applies(x)) {
291 OPM_THROW(Opm::NumericalProblem,
292 "Tried to evaluate a derivative of a tabulated"
293 " function outside of its range");
296 unsigned segIdx = findSegmentIndex_(Opm::scalarValue(x));
297 return evalDerivative_(x, segIdx);
314 template <
class Evaluation>
317 if (!extrapolate && !
applies(x)) {
318 OPM_THROW(Opm::NumericalProblem,
319 "Tried to evaluate a second derivative of a tabulated "
320 " function outside of its range");
340 template <
class Evaluation>
343 if (!extrapolate && !
applies(x)) {
344 OPM_THROW(Opm::NumericalProblem,
345 "Tried to evaluate a third derivative of a tabulated "
346 " function outside of its range");
360 int monotonic(Scalar x0, Scalar x1,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const
377 size_t i = findSegmentIndex_(x0);
378 if (xValues_[i + 1] >= x1) {
381 updateMonotonicity_(i, r);
387 updateMonotonicity_(i, r);
392 size_t iEnd = findSegmentIndex_(x1);
393 for (; i < iEnd - 1; ++i) {
394 updateMonotonicity_(i, r);
407 return (r < 0 || r==3)?-1:0;
409 return (r > 0 || r==3)?1:0;
415 updateMonotonicity_(iEnd, r);
443 void printCSV(Scalar xi0, Scalar xi1,
unsigned k, std::ostream& os = std::cout)
const
445 Scalar x0 = std::min(xi0, xi1);
446 Scalar x1 = std::max(xi0, xi1);
448 for (
int i = 0; i <= k; ++i) {
449 double x = i*(x1 - x0)/k + x0;
453 if (x < xValues_[0]) {
455 y = (x - xValues_[0])*dy_dx + yValues_[0];
457 else if (x > xValues_[n]) {
459 y = (x - xValues_[n])*dy_dx + yValues_[n];
462 OPM_THROW(std::runtime_error,
463 "The sampling points given to a function must be sorted by their x value!");
471 os << x <<
" " << y <<
" " << dy_dx <<
"\n";
476 size_t findSegmentIndex_(Scalar x)
const
479 assert(xValues_.size() >= 2);
481 if (x <= xValues_[1])
483 else if (x >= xValues_[xValues_.size() - 2])
484 return xValues_.size() - 2;
487 size_t segmentIdx = 1;
488 size_t upperIdx = xValues_.size() - 2;
489 while (segmentIdx + 1 < upperIdx) {
490 size_t pivotIdx = (segmentIdx + upperIdx) / 2;
491 if (x < xValues_[pivotIdx])
494 segmentIdx = pivotIdx;
497 assert(xValues_[segmentIdx] <= x);
498 assert(x <= xValues_[segmentIdx + 1]);
503 template <
class Evaluation>
504 Evaluation evalDerivative_(
const Evaluation& x OPM_UNUSED,
size_t segIdx)
const
506 Scalar x0 = xValues_[segIdx];
507 Scalar x1 = xValues_[segIdx + 1];
509 Scalar y0 = yValues_[segIdx];
510 Scalar y1 = yValues_[segIdx + 1];
512 return (y1 - y0)/(x1 - x0);
523 int updateMonotonicity_(
size_t i,
int& r)
const
525 if (yValues_[i] < yValues_[i + 1]) {
527 if (r == 3 || r == 1)
533 else if (yValues_[i] > yValues_[i + 1]) {
535 if (r == 3 || r == -1)
550 ComparatorX_(
const std::vector<Scalar>& x)
554 bool operator ()(
size_t idxA,
size_t idxB)
const
555 {
return x_.at(idxA) < x_.at(idxB); }
557 const std::vector<Scalar>& x_;
568 std::vector<unsigned> idxVector(n);
569 for (
unsigned i = 0; i < n; ++i)
574 ComparatorX_ cmp(xValues_);
575 std::sort(idxVector.begin(), idxVector.end(), cmp);
578 std::vector<Scalar> tmpX(n), tmpY(n);
579 for (
size_t i = 0; i < idxVector.size(); ++ i) {
580 tmpX[i] = xValues_[idxVector[i]];
581 tmpY[i] = yValues_[idxVector[i]];
591 void reverseSamplingPoints_()
595 for (
size_t i = 0; i <= (n - 1)/2; ++i) {
596 std::swap(xValues_[i], xValues_[n - i - 1]);
597 std::swap(yValues_[i], yValues_[n - i - 1]);
604 void resizeArrays_(
size_t nSamples)
606 xValues_.resize(nSamples);
607 yValues_.resize(nSamples);
610 std::vector<Scalar> xValues_;
611 std::vector<Scalar> yValues_;
int monotonic() const
Same as monotonic(x0, x1), but with the entire range of the function as interval. ...
Definition: Tabulated1DFunction.hpp:424
Evaluation evalDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline's derivative at a given position.
Definition: Tabulated1DFunction.hpp:288
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:360
void setXYContainers(const ScalarContainerX &x, const ScalarContainerY &y, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:128
Scalar xAt(size_t i) const
Return the x value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:226
Evaluation evalThirdDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the function's third derivative at a given position.
Definition: Tabulated1DFunction.hpp:341
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
Tabulated1DFunction()
Default constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:55
bool applies(const Evaluation &x) const
Return true iff the given x is in range [x1, xn].
Definition: Tabulated1DFunction.hpp:239
Scalar xMax() const
Return the x value of the rightmost sampling point.
Definition: Tabulated1DFunction.hpp:220
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Tabulated1DFunction(const PointContainer &points, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:93
Scalar xMin() const
Return the x value of the leftmost sampling point.
Definition: Tabulated1DFunction.hpp:214
Evaluation evalSecondDerivative(const Evaluation &x, bool extrapolate=false) const
Evaluate the function's second derivative at a given position.
Definition: Tabulated1DFunction.hpp:315
Evaluation eval(const Evaluation &x, bool extrapolate=false) const
Evaluate the spline at a given position.
Definition: Tabulated1DFunction.hpp:252
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 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
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:443
size_t numSamples() const
Returns the number of sampling points.
Definition: Tabulated1DFunction.hpp:208
Implements a linearly interpolated scalar function that depends on one variable.
Definition: Tabulated1DFunction.hpp:47
Scalar valueAt(size_t i) const
Return the value of the a sample point with a given index.
Definition: Tabulated1DFunction.hpp:232
void setArrayOfPoints(size_t nSamples, const PointArray &points, bool sortInputs=true)
Set the sampling points for the piecewise linear function.
Definition: Tabulated1DFunction.hpp:149
Tabulated1DFunction(const ScalarContainer &x, const ScalarContainer &y, bool sortInputs=true)
Convenience constructor for a piecewise linear function.
Definition: Tabulated1DFunction.hpp:81