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