My Project
NonuniformTableLinear.hpp
1 /*
2  Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
3  Copyright 2009, 2010, 2011, 2012 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_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
22 #define OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
23 
24 #include <cmath>
25 #include <exception>
26 #include <vector>
27 #include <utility>
28 
29 #include <opm/common/ErrorMacros.hpp>
30 #include <opm/common/utility/numeric/linearInterpolation.hpp>
31 
32 
33 namespace Opm
34 {
35 
36 
42  template<typename T>
44  {
45  public:
48 
52  template<class XContainer, class YContainer>
53  NonuniformTableLinear(const XContainer& x_values,
54  const YContainer& y_values);
55 
58  std::pair<double, double> domain();
59 
62  void rescaleDomain(std::pair<double, double> new_domain);
63 
67  double operator()(const double x) const;
68 
72  double derivative(const double x) const;
73 
77  double inverse(const double y) const;
78 
82  bool operator==(const NonuniformTableLinear& other) const;
83 
84  protected:
85  std::vector<double> x_values_;
86  std::vector<T> y_values_;
87  mutable std::vector<T> x_values_reversed_;
88  mutable std::vector<T> y_values_reversed_;
89  };
90 
91 
92  // A utility function
99  template <typename FI>
100  bool isNondecreasing(const FI beg, const FI end)
101  {
102  if (beg == end) return true;
103  FI it = beg;
104  ++it;
105  FI prev = beg;
106  for (; it != end; ++it, ++prev) {
107  if (*it < *prev) {
108  return false;
109  }
110  }
111  return true;
112  }
113 
114 
115 
116  // Member implementations.
117 
118  template<typename T>
119  inline
122  {
123  }
124 
125 
126  template<typename T>
127  template<class XContainer, class YContainer>
128  inline
130  ::NonuniformTableLinear(const XContainer& x_column,
131  const YContainer& y_column)
132  : x_values_( x_column.begin() , x_column.end()),
133  y_values_( y_column.begin() , y_column.end())
134  {
135  assert(isNondecreasing(x_values_.begin(), x_values_.end()));
136  }
137 
138 
139  template<typename T>
140  inline std::pair<double, double>
143  {
144  return std::make_pair(x_values_[0], x_values_.back());
145  }
146 
147  template<typename T>
148  inline void
150  ::rescaleDomain(std::pair<double, double> new_domain)
151  {
152  const double a = x_values_[0];
153  const double b = x_values_.back();
154  const double c = new_domain.first;
155  const double d = new_domain.second;
156  // x in [a, b] -> x in [c, d]
157  for (int i = 0; i < int(x_values_.size()); ++i) {
158  x_values_[i] = (x_values_[i] - a)*(d - c)/(b - a) + c;
159  }
160  }
161 
162  template<typename T>
163  inline double
165  ::operator()(const double x) const
166  {
167  return Opm::linearInterpolation(x_values_, y_values_, x);
168  }
169 
170  template<typename T>
171  inline double
173  ::derivative(const double x) const
174  {
175  return Opm::linearInterpolationDerivative(x_values_, y_values_, x);
176  }
177 
178  template<typename T>
179  inline double
181  ::inverse(const double y) const
182  {
183  if (y_values_.front() < y_values_.back()) {
184  return Opm::linearInterpolation(y_values_, x_values_, y);
185  } else {
186  if (y_values_reversed_.empty()) {
187  y_values_reversed_ = y_values_;
188  std::reverse(y_values_reversed_.begin(), y_values_reversed_.end());
189  assert(isNondecreasing(y_values_reversed_.begin(), y_values_reversed_.end()));
190  x_values_reversed_ = x_values_;
191  std::reverse(x_values_reversed_.begin(), x_values_reversed_.end());
192  }
193  return Opm::linearInterpolation(y_values_reversed_, x_values_reversed_, y);
194  }
195  }
196 
197  template<typename T>
198  inline bool
201  {
202  return x_values_ == other.x_values_
203  && y_values_ == other.y_values_;
204  }
205 
206 } // namespace Opm
207 
208 #endif // OPM_NONUNIFORMTABLELINEAR_HEADER_INCLUDED
This class uses linear interpolation to compute the value (and its derivative) of a function f sample...
Definition: NonuniformTableLinear.hpp:44
double derivative(const double x) const
Evaluate the derivative at x.
Definition: NonuniformTableLinear.hpp:173
double operator()(const double x) const
Evaluate the value at x.
Definition: NonuniformTableLinear.hpp:165
void rescaleDomain(std::pair< double, double > new_domain)
Rescale the domain.
Definition: NonuniformTableLinear.hpp:150
double inverse(const double y) const
Evaluate the inverse at y.
Definition: NonuniformTableLinear.hpp:181
bool operator==(const NonuniformTableLinear &other) const
Equality operator.
Definition: NonuniformTableLinear.hpp:200
NonuniformTableLinear()
Default constructor.
Definition: NonuniformTableLinear.hpp:121
std::pair< double, double > domain()
Get the domain.
Definition: NonuniformTableLinear.hpp:142
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition: Exceptions.hpp:29
bool isNondecreasing(const FI beg, const FI end)
Detect if a sequence is nondecreasing.
Definition: NonuniformTableLinear.hpp:100