OGS 6.2.1-97-g73d1aeda3
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  {
39  return false;
40  }
41 
42  return std::none_of(_supp_pnts.begin(), _supp_pnts.end(),
43  [&](const double p) {
44  return this->getDerivative(p) * gradient0 <= 0.;
45  });
46 }
47 
49 {
50  std::size_t interval_idx = 0;
51  if (_values_at_supp_pnts.front() < _values_at_supp_pnts.back())
52  {
53  if (y <= _values_at_supp_pnts.front())
54  {
55  return _supp_pnts[0];
56  }
57  if (y >= _values_at_supp_pnts.back())
58  {
59  return _supp_pnts[_supp_pnts.size() - 1];
60  }
61 
62  // search interval that has the point inside
63  auto const& it(std::lower_bound(_values_at_supp_pnts.begin(),
64  _values_at_supp_pnts.end(), y));
65  interval_idx = std::distance(_values_at_supp_pnts.begin(), it) - 1;
66  }
67  else
68  {
69  if (y >= _values_at_supp_pnts.front())
70  {
71  return _supp_pnts[0];
72  }
73  if (y <= _values_at_supp_pnts.back())
74  {
75  return _supp_pnts[_supp_pnts.size() - 1];
76  }
77 
78  // search interval in the reverse direction for the point inside
79  auto const& it(std::lower_bound(_values_at_supp_pnts.rbegin(),
80  _values_at_supp_pnts.rend(), y));
81  interval_idx = _values_at_supp_pnts.size() -
82  (std::distance(_values_at_supp_pnts.rbegin(), it)) - 1;
83  }
84 
85  const double xi_1 = _supp_pnts[interval_idx + 1];
86  const double xi = _supp_pnts[interval_idx];
87  const double yi_1 = _values_at_supp_pnts[interval_idx + 1];
88  const double yi = _values_at_supp_pnts[interval_idx];
89 
90  // compute gradient: m = (x_{i+1} - x_i) / (y_{i+1} - y_i)
91  const double m = (xi_1 - xi) / (yi_1 - yi);
92 
93  // compute the variable by linear interpolation: x = m * (y - y_i) + x_i,
94  // and then return the result.
95  return m * (y - yi) + xi;
96 }
97 
98 } // namespace MathLib
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:63
PiecewiseLinearMonotonicCurve(std::vector< double > &&x, std::vector< double > &&y)