OGS 6.1.0-1721-g6382411ad
PiecewiseLinearMonotonicCurve.cpp
Go to the documentation of this file.
1 
14 
15 #include <algorithm>
16 #include <cmath>
17 #include <limits>
18 
19 #include "BaseLib/Error.h"
20 
21 namespace MathLib
22 {
24  std::vector<double>&& x, std::vector<double>&& y)
25  : PiecewiseLinearInterpolation(std::move(x), std::move(y), false)
26 {
27  if (!isStrongMonotonic())
28  {
29  OGS_FATAL("The given curve is not strong monotonic.");
30  }
31 }
32 
34 {
35  const double gradient0 = getDerivative(_supp_pnts[0]);
36 
37  if (std::fabs(gradient0) < std::numeric_limits<double>::min())
38  return false;
39 
40  return std::none_of(_supp_pnts.begin(), _supp_pnts.end(),
41  [&](const double p) {
42  return this->getDerivative(p) * gradient0 <= 0.;
43  });
44 }
45 
47 {
48  std::size_t interval_idx = 0;
49  if (_values_at_supp_pnts.front() < _values_at_supp_pnts.back())
50  {
51  if (y <= _values_at_supp_pnts.front())
52  {
53  return _supp_pnts[0];
54  }
55  if (y >= _values_at_supp_pnts.back())
56  {
57  return _supp_pnts[_supp_pnts.size() - 1];
58  }
59 
60  // search interval that has the point inside
61  auto const& it(std::lower_bound(_values_at_supp_pnts.begin(),
62  _values_at_supp_pnts.end(), y));
63  interval_idx = std::distance(_values_at_supp_pnts.begin(), it) - 1;
64  }
65  else
66  {
67  if (y >= _values_at_supp_pnts.front())
68  {
69  return _supp_pnts[0];
70  }
71  if (y <= _values_at_supp_pnts.back())
72  {
73  return _supp_pnts[_supp_pnts.size() - 1];
74  }
75 
76  // search interval in the reverse direction for the point inside
77  auto const& it(std::lower_bound(_values_at_supp_pnts.rbegin(),
78  _values_at_supp_pnts.rend(), y));
79  interval_idx = _values_at_supp_pnts.size() -
80  (std::distance(_values_at_supp_pnts.rbegin(), it)) - 1;
81  }
82 
83  const double xi_1 = _supp_pnts[interval_idx + 1];
84  const double xi = _supp_pnts[interval_idx];
85  const double yi_1 = _values_at_supp_pnts[interval_idx + 1];
86  const double yi = _values_at_supp_pnts[interval_idx];
87 
88  // compute gradient: m = (x_{i+1} - x_i) / (y_{i+1} - y_i)
89  const double m = (xi_1 - xi) / (yi_1 - yi);
90 
91  // compute the variable by linear interpolation: x = m * (y - y_i) + x_i,
92  // and then return the result.
93  return m * (y - yi) + xi;
94 }
95 
96 } // namespace
double getDerivative(double const pnt_to_interpolate) const
Calculates derivative using quadratic interpolation and central difference quotient.
static const double p
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
PiecewiseLinearMonotonicCurve(std::vector< double > &&x, std::vector< double > &&y)