OGS
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::abs(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 
48 {
49  std::size_t interval_idx = 0;
50  if (values_at_supp_pnts_.front() < values_at_supp_pnts_.back())
51  {
52  if (y <= values_at_supp_pnts_.front())
53  {
54  return supp_pnts_[0];
55  }
56  if (y >= values_at_supp_pnts_.back())
57  {
58  return supp_pnts_[supp_pnts_.size() - 1];
59  }
60 
61  // search interval that has the point inside
62  auto const& it(std::lower_bound(values_at_supp_pnts_.begin(),
63  values_at_supp_pnts_.end(), y));
64  interval_idx = std::distance(values_at_supp_pnts_.begin(), it) - 1;
65  }
66  else
67  {
68  if (y >= values_at_supp_pnts_.front())
69  {
70  return supp_pnts_[0];
71  }
72  if (y <= values_at_supp_pnts_.back())
73  {
74  return supp_pnts_[supp_pnts_.size() - 1];
75  }
76 
77  // search interval in the reverse direction for the point inside
78  auto const& it(std::lower_bound(values_at_supp_pnts_.rbegin(),
79  values_at_supp_pnts_.rend(), y));
80  interval_idx = values_at_supp_pnts_.size() -
81  (std::distance(values_at_supp_pnts_.rbegin(), it)) - 1;
82  }
83 
84  const double xi_1 = supp_pnts_[interval_idx + 1];
85  const double xi = supp_pnts_[interval_idx];
86  const double yi_1 = values_at_supp_pnts_[interval_idx + 1];
87  const double yi = values_at_supp_pnts_[interval_idx];
88 
89  // compute gradient: m = (x_{i+1} - x_i) / (y_{i+1} - y_i)
90  const double m = (xi_1 - xi) / (yi_1 - yi);
91 
92  // compute the variable by linear interpolation: x = m * (y - y_i) + x_i,
93  // and then return the result.
94  return m * (y - yi) + xi;
95 }
96 
97 } // namespace MathLib
#define OGS_FATAL(...)
Definition: Error.h:26
double getDerivative(double const pnt_to_interpolate) const
Calculates derivative using quadratic interpolation and central difference quotient.
PiecewiseLinearMonotonicCurve(std::vector< double > x, std::vector< double > y)
static const double p