OGS 6.2.1-97-g73d1aeda3
PiecewiseLinearInterpolation.cpp
Go to the documentation of this file.
1 
14 #include <cmath>
15 #include <limits>
16 #include <utility>
17 
18 #include "BaseLib/Error.h"
19 #include "BaseLib/quicksort.h"
20 
22 
23 namespace MathLib
24 {
26  std::vector<double>&& supporting_points,
27  std::vector<double>&& values_at_supp_pnts,
28  bool supp_pnts_sorted)
29  : _supp_pnts(std::move(supporting_points)),
30  _values_at_supp_pnts(std::move(values_at_supp_pnts))
31 {
32  if (!supp_pnts_sorted)
33  {
34  BaseLib::quicksort<double, double>(
35  _supp_pnts, static_cast<std::size_t>(0), _supp_pnts.size(),
37  }
38 
39  const auto it = std::adjacent_find(_supp_pnts.begin(), _supp_pnts.end());
40  if (it != _supp_pnts.end())
41  {
42  const std::size_t i = std::distance(_supp_pnts.begin(), it);
43  OGS_FATAL(
44  "Variable %d and variable %d are the same. "
45  "Piecewise linear interpolation is not possible\n",
46  i, i + 1);
47  }
48 }
49 
50 double PiecewiseLinearInterpolation::getValue(double pnt_to_interpolate) const
51 {
52  // search interval that has the point inside
53  if (pnt_to_interpolate <= _supp_pnts.front())
54  {
55  return _values_at_supp_pnts[0];
56  }
57 
58  if (_supp_pnts.back() <= pnt_to_interpolate)
59  {
60  return _values_at_supp_pnts[_supp_pnts.size() - 1];
61  }
62 
63  auto const& it(std::lower_bound(_supp_pnts.begin(), _supp_pnts.end(),
64  pnt_to_interpolate));
65  std::size_t const interval_idx = std::distance(_supp_pnts.begin(), it) - 1;
66 
67  // support points.
68  double const x = _supp_pnts[interval_idx];
69  double const x_r = _supp_pnts[interval_idx + 1];
70 
71  // values.
72  double const f = _values_at_supp_pnts[interval_idx];
73  double const f_r = _values_at_supp_pnts[interval_idx + 1];
74 
75  // compute linear interpolation polynom: y = m * (x - support[i]) + value[i]
76  const double m = (f_r - f) / (x_r - x);
77 
78  return m * (pnt_to_interpolate - x) + f;
79 }
80 
82  double const pnt_to_interpolate) const
83 {
84  if (pnt_to_interpolate < _supp_pnts.front() ||
85  _supp_pnts.back() < pnt_to_interpolate)
86  {
87  return 0;
88  }
89 
90  auto const& it(std::lower_bound(_supp_pnts.begin(), _supp_pnts.end(),
91  pnt_to_interpolate));
92  std::size_t interval_idx = std::distance(_supp_pnts.begin(), it);
93 
94  if (pnt_to_interpolate == _supp_pnts.front())
95  {
96  interval_idx = 1;
97  }
98 
99  if (interval_idx > 1 && interval_idx < _supp_pnts.size() - 2)
100  {
101  // left and right support points.
102  double const x_ll = _supp_pnts[interval_idx - 2];
103  double const x_l = _supp_pnts[interval_idx - 1];
104  double const x = _supp_pnts[interval_idx];
105  double const x_r = _supp_pnts[interval_idx + 1];
106 
107  // left and right values.
108  double const f_ll = _values_at_supp_pnts[interval_idx - 2];
109  double const f_l = _values_at_supp_pnts[interval_idx - 1];
110  double const f = _values_at_supp_pnts[interval_idx];
111  double const f_r = _values_at_supp_pnts[interval_idx + 1];
112 
113  double const tangent_right = (f_l - f_r) / (x_l - x_r);
114  double const tangent_left = (f_ll - f) / (x_ll - x);
115  double const w = (pnt_to_interpolate - x) / (x_l - x);
116  return (1. - w) * tangent_right + w * tangent_left;
117  }
118 
119  return (_values_at_supp_pnts[interval_idx] -
120  _values_at_supp_pnts[interval_idx - 1]) /
121  (_supp_pnts[interval_idx] - _supp_pnts[interval_idx - 1]);
122 }
123 
125 {
126  assert(!_supp_pnts.empty());
127  return _supp_pnts.back();
128 }
130 {
131  assert(!_supp_pnts.empty());
132  return _supp_pnts.front();
133 }
134 } // namespace MathLib
double getDerivative(double const pnt_to_interpolate) const
Calculates derivative using quadratic interpolation and central difference quotient.
Definition of the PiecewiseLinearInterpolation class.
PiecewiseLinearInterpolation(std::vector< double > &&supporting_points, std::vector< double > &&values_at_supp_pnts, bool supp_pnts_sorted=false)
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
double getValue(double pnt_to_interpolate) const
Calculates the interpolation value.