My Project
Loading...
Searching...
No Matches
linearInterpolation.hpp
1/*
2 Copyright 2009, 2010, 2013 SINTEF ICT, Applied Mathematics.
3 Copyright 2009, 2010 Statoil ASA.
4
5 This file is part of the Open Porous Media project (OPM).
6
7 OPM is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 OPM is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with OPM. If not, see <http://www.gnu.org/licenses/>.
19*/
20
21#ifndef OPM_LINEARINTERPOLATION_HEADER_INCLUDED
22#define OPM_LINEARINTERPOLATION_HEADER_INCLUDED
23
24
25#include <algorithm>
26#include <tuple>
27#include <vector>
28
29namespace Opm
30{
31
35inline int tableIndex(const std::vector<double>& table, double x)
36{
37 if (table.size() < 2)
38 return 0;
39
40 const auto lower = std::lower_bound(table.begin(), table.end(), x);
41
42 if (lower == table.end())
43 return table.size()-2;
44 else if (lower == table.begin())
45 return 0;
46 else
47 return std::distance(table.begin(), lower)-1;
48}
49
50inline std::pair<double, int>
51linearInterpolationSlope(const std::vector<double>& xv,
52 const std::vector<double>& yv,
53 const double x)
54{
55 const auto i = Opm::tableIndex(xv, x);
56 return { (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]), i };
57}
58
59
60inline double linearInterpolationDerivative(const std::vector<double>& xv,
61 const std::vector<double>& yv, double x)
62{
63 // Extrapolates if x is outside xv
64 return linearInterpolationSlope(xv, yv, x).first;
65}
66
67inline double linearInterpolation(const std::vector<double>& xv,
68 const std::vector<double>& yv, double x)
69{
70 // Extrapolates if x is outside xv
71 const auto& [t, i] = linearInterpolationSlope(xv, yv, x);
72 return t * (x - xv[i]) + yv[i];
73}
74
75inline double linearInterpolationNoExtrapolation(const std::vector<double>& xv,
76 const std::vector<double>& yv, double x)
77{
78 // Return end values if x is outside xv
79 if (x < xv.front()) {
80 return yv.front();
81 }
82 if (x > xv.back()) {
83 return yv.back();
84 }
85
86 return linearInterpolation(xv, yv, x);
87}
88
89inline double linearInterpolation(const std::vector<double>& xv,
90 const std::vector<double>& yv,
91 double x, int& ix1)
92{
93 // Extrapolates if x is outside xv
94 double t;
95 std::tie(t, ix1) = linearInterpolationSlope(xv, yv, x);
96 return t * (x - xv[ix1]) + yv[ix1];
97}
98
99} // namespace Opm
100
101#endif // OPM_LINEARINTERPOLATION_HEADER_INCLUDED
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
int tableIndex(const std::vector< double > &table, double x)
Returns an index in an ordered table such that x is between table[j] and table[j+1].
Definition linearInterpolation.hpp:35