My Project
IntervalTabulated2DFunction.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
29 #define OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
30 
35 
36 #include <vector>
37 #include <limits>
38 #include <sstream>
39 #include <cassert>
40 #include <algorithm>
41 
42 namespace Opm {
43 
50 template <class Scalar>
52 {
53 public:
55  { }
56 
57  template <class DataContainer>
58  IntervalTabulated2DFunction(const std::vector<Scalar>& xPos,
59  const std::vector<Scalar>& yPos,
60  const DataContainer& data,
61  const bool xExtrapolate = false,
62  const bool yExtrapolate = false)
63  : xPos_(xPos)
64  , yPos_(yPos)
65  , samples_(data)
66  , xExtrapolate_(xExtrapolate)
67  , yExtrapolate_(yExtrapolate)
68  {
69 #ifndef NDEBUG
70  // in debug mode, ensure that the x and y positions arrays are strictly
71  // mononically increasing.
72  for (unsigned i = 0; i < xPos.size() - 1; ++ i) {
73  if (xPos[i + 1] <= xPos[i])
74  throw std::runtime_error("The array for the x-positions is not strictly increasing!");
75  }
76 
77  for (unsigned i = 0; i < yPos.size() - 1; ++ i) {
78  if (yPos[i + 1] <= yPos[i])
79  throw std::runtime_error("The array for the y-positions is not strictly increasing!");
80  }
81 #endif
82 
83  // make sure the size is correct
84  if (numX() != samples_.size())
85  throw std::runtime_error("numX() is not equal to the number of rows of the sampling points");
86 
87  for (unsigned xIdx = 0; xIdx < numX(); ++xIdx) {
88  if (samples_[xIdx].size() != numY()) {
89  std::ostringstream oss;
90  oss << "The " << xIdx << "-th row of the sampling points has different size than numY() ";
91  throw std::runtime_error(oss.str());
92  }
93  }
94  }
95 
99  size_t numX() const
100  { return xPos_.size(); }
101 
105  size_t numY() const
106  { return yPos_.size(); }
107 
111  Scalar xMin() const
112  { return xPos_.front(); }
113 
117  Scalar xMax() const
118  { return xPos_.back(); }
119 
123  Scalar yMin() const
124  { return yPos_.front(); }
125 
129  Scalar yMax() const
130  { return yPos_.back(); }
131 
132  const std::vector<Scalar>& xPos() const
133  { return xPos_; }
134 
135  const std::vector<Scalar>& yPos() const
136  { return yPos_; }
137 
138  const std::vector<std::vector<Scalar>>& samples() const
139  { return samples_; }
140 
141  bool xExtrapolate() const
142  { return xExtrapolate_; }
143 
144  bool yExtrapolate() const
145  { return yExtrapolate_; }
146 
147  bool operator==(const IntervalTabulated2DFunction<Scalar>& data) const {
148  return this->xPos() == data.xPos() &&
149  this->yPos() == data.yPos() &&
150  this->samples() == data.samples() &&
151  this->xExtrapolate() == data.xExtrapolate() &&
152  this->yExtrapolate() == data.yExtrapolate();
153  }
154 
158  Scalar valueAt(size_t i, size_t j) const
159  { return samples_[i][j]; }
160 
164  template <class Evaluation>
165  bool applies(const Evaluation& x, const Evaluation& y) const
166  { return appliesX(x) && appliesY(y); }
167 
171  template <class Evaluation>
172  bool appliesX(const Evaluation& x) const
173  { return xMin() <= x && x <= xMax(); }
174 
178  template <class Evaluation>
179  bool appliesY(const Evaluation& y) const
180  { return yMin() <= y && y <= yMax(); }
181 
182 
190  template <typename Evaluation>
191  Evaluation eval(const Evaluation& x, const Evaluation& y) const
192  {
193  if ((!xExtrapolate_ && !appliesX(x)) || (!yExtrapolate_ && !appliesY(y))) {
194  std::ostringstream oss;
195  oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
196  throw NumericalIssue(oss.str());
197  };
198 
199  // bi-linear interpolation: first, calculate the x and y indices in the lookup
200  // table ...
201  const unsigned i = xSegmentIndex_(x);
202  const unsigned j = ySegmentIndex_(y);
203 
204  // bi-linear interpolation / extrapolation
205  const Evaluation alpha = xToAlpha(x, i);
206  const Evaluation beta = yToBeta(y, j);
207 
208  const Evaluation s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta;
209  const Evaluation s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta;
210 
211  Valgrind::CheckDefined(s1);
212  Valgrind::CheckDefined(s2);
213 
214  // ... and combine them using the x position
215  return s1*(1.0 - alpha) + s2*alpha;
216  }
217 
218 private:
219  // the sampling points in the x-drection
220  std::vector<Scalar> xPos_;
221  // the sampling points in the y-drection
222  std::vector<Scalar> yPos_;
223  // data at the sampling points
224  std::vector<std::vector<Scalar> > samples_;
225 
226  bool xExtrapolate_ = false;
227  bool yExtrapolate_ = false;
228 
232  template <class Evaluation>
233  unsigned xSegmentIndex_(const Evaluation& x) const
234  {
235  assert(xExtrapolate_ || appliesX(x) );
236 
237  return segmentIndex_(x, xPos_);
238  }
239 
243  template <class Evaluation>
244  unsigned ySegmentIndex_(const Evaluation& y) const
245  {
246  assert(yExtrapolate_ || appliesY(y) );
247 
248  return segmentIndex_(y, yPos_);
249  }
250 
251 
252  template <class Evaluation>
253  static unsigned segmentIndex_(const Evaluation& v, const std::vector<Scalar>& vPos)
254  {
255  const unsigned n = vPos.size();
256  assert(n >= 2);
257 
258  if (v <= vPos.front() || n == 2)
259  return 0;
260  else if (v >= vPos.back())
261  return n - 2;
262 
263  assert(n > 2 && v > vPos.front() && v < vPos.back());
264 
265  // bisection. this assumes that the vPos array is strictly mononically
266  // increasing.
267  size_t lowerIdx = 0;
268  size_t upperIdx = vPos.size() - 1;
269  while (lowerIdx + 1 < upperIdx) {
270  size_t pivotIdx = (lowerIdx + upperIdx) / 2;
271  if (v < vPos[pivotIdx])
272  upperIdx = pivotIdx;
273  else
274  lowerIdx = pivotIdx;
275  }
276 
277  assert(vPos[lowerIdx] <= v);
278  assert(v <= vPos[lowerIdx + 1]);
279  return lowerIdx;
280  }
281 
288  template <class Evaluation>
289  Evaluation xToAlpha(const Evaluation& x, unsigned xSegmentIdx) const
290  {
291  Scalar x1 = xPos_[xSegmentIdx];
292  Scalar x2 = xPos_[xSegmentIdx + 1];
293  return (x - x1)/(x2 - x1);
294  }
295 
302  template <class Evaluation>
303  Evaluation yToBeta(const Evaluation& y, unsigned ySegmentIdx) const
304  {
305  Scalar y1 = yPos_[ySegmentIdx];
306  Scalar y2 = yPos_[ySegmentIdx + 1];
307  return (y - y1)/(y2 - y1);
308  }
309 
310 };
311 } // namespace Opm
312 
313 #endif
Provides the opm-material specific exception classes.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Provides the OPM_UNUSED macro.
Some templates to wrap the valgrind client request macros.
Implements a function that depends on two variables.
Definition: IntervalTabulated2DFunction.hpp:52
size_t numY() const
Returns the number of sampling points in Y direction.
Definition: IntervalTabulated2DFunction.hpp:105
Evaluation eval(const Evaluation &x, const Evaluation &y) const
Evaluate the function at a given (x,y) position.
Definition: IntervalTabulated2DFunction.hpp:191
bool appliesX(const Evaluation &x) const
Returns true if a coordinate lies in the tabulated range on the x direction.
Definition: IntervalTabulated2DFunction.hpp:172
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:111
Scalar yMin() const
Returns the minimum of the Y coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:123
size_t numX() const
Returns the number of sampling points in X direction.
Definition: IntervalTabulated2DFunction.hpp:99
bool appliesY(const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range on the y direction.
Definition: IntervalTabulated2DFunction.hpp:179
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range.
Definition: IntervalTabulated2DFunction.hpp:165
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:117
Scalar yMax() const
Returns the maximum of the Y coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:129
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition: IntervalTabulated2DFunction.hpp:158
Definition: Exceptions.hpp:46