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 
34 
35 #include <vector>
36 #include <limits>
37 #include <sstream>
38 #include <cassert>
39 #include <algorithm>
40 
41 namespace Opm {
42 
49 template <class Scalar>
51 {
52 public:
54  { }
55 
56  template <class DataContainer>
57  IntervalTabulated2DFunction(const std::vector<Scalar>& xPos,
58  const std::vector<Scalar>& yPos,
59  const DataContainer& data,
60  const bool xExtrapolate = false,
61  const bool yExtrapolate = false)
62  : xPos_(xPos)
63  , yPos_(yPos)
64  , samples_(data)
65  , xExtrapolate_(xExtrapolate)
66  , yExtrapolate_(yExtrapolate)
67  {
68 #ifndef NDEBUG
69  // in debug mode, ensure that the x and y positions arrays are strictly
70  // mononically increasing.
71  for (unsigned i = 0; i < xPos.size() - 1; ++ i) {
72  if (xPos[i + 1] <= xPos[i])
73  throw std::runtime_error("The array for the x-positions is not strictly increasing!");
74  }
75 
76  for (unsigned i = 0; i < yPos.size() - 1; ++ i) {
77  if (yPos[i + 1] <= yPos[i])
78  throw std::runtime_error("The array for the y-positions is not strictly increasing!");
79  }
80 #endif
81 
82  // make sure the size is correct
83  if (numX() != samples_.size())
84  throw std::runtime_error("numX() is not equal to the number of rows of the sampling points");
85 
86  for (unsigned xIdx = 0; xIdx < numX(); ++xIdx) {
87  if (samples_[xIdx].size() != numY()) {
88  std::ostringstream oss;
89  oss << "The " << xIdx << "-th row of the sampling points has different size than numY() ";
90  throw std::runtime_error(oss.str());
91  }
92  }
93  }
94 
98  size_t numX() const
99  { return xPos_.size(); }
100 
104  size_t numY() const
105  { return yPos_.size(); }
106 
110  Scalar xMin() const
111  { return xPos_.front(); }
112 
116  Scalar xMax() const
117  { return xPos_.back(); }
118 
122  Scalar yMin() const
123  { return yPos_.front(); }
124 
128  Scalar yMax() const
129  { return yPos_.back(); }
130 
131  const std::vector<Scalar>& xPos() const
132  { return xPos_; }
133 
134  const std::vector<Scalar>& yPos() const
135  { return yPos_; }
136 
137  const std::vector<std::vector<Scalar>>& samples() const
138  { return samples_; }
139 
140  bool xExtrapolate() const
141  { return xExtrapolate_; }
142 
143  bool yExtrapolate() const
144  { return yExtrapolate_; }
145 
146  bool operator==(const IntervalTabulated2DFunction<Scalar>& data) const {
147  return this->xPos() == data.xPos() &&
148  this->yPos() == data.yPos() &&
149  this->samples() == data.samples() &&
150  this->xExtrapolate() == data.xExtrapolate() &&
151  this->yExtrapolate() == data.yExtrapolate();
152  }
153 
157  Scalar valueAt(size_t i, size_t j) const
158  { return samples_[i][j]; }
159 
163  template <class Evaluation>
164  bool applies(const Evaluation& x, const Evaluation& y) const
165  { return appliesX(x) && appliesY(y); }
166 
170  template <class Evaluation>
171  bool appliesX(const Evaluation& x) const
172  { return xMin() <= x && x <= xMax(); }
173 
177  template <class Evaluation>
178  bool appliesY(const Evaluation& y) const
179  { return yMin() <= y && y <= yMax(); }
180 
181 
189  template <typename Evaluation>
190  Evaluation eval(const Evaluation& x, const Evaluation& y) const
191  {
192  if ((!xExtrapolate_ && !appliesX(x)) || (!yExtrapolate_ && !appliesY(y))) {
193  std::ostringstream oss;
194  oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
195  throw NumericalIssue(oss.str());
196  };
197 
198  // bi-linear interpolation: first, calculate the x and y indices in the lookup
199  // table ...
200  const unsigned i = xSegmentIndex_(x);
201  const unsigned j = ySegmentIndex_(y);
202 
203  // bi-linear interpolation / extrapolation
204  const Evaluation alpha = xToAlpha(x, i);
205  const Evaluation beta = yToBeta(y, j);
206 
207  const Evaluation s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta;
208  const Evaluation s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta;
209 
210  Valgrind::CheckDefined(s1);
211  Valgrind::CheckDefined(s2);
212 
213  // ... and combine them using the x position
214  return s1*(1.0 - alpha) + s2*alpha;
215  }
216 
217 private:
218  // the sampling points in the x-drection
219  std::vector<Scalar> xPos_;
220  // the sampling points in the y-drection
221  std::vector<Scalar> yPos_;
222  // data at the sampling points
223  std::vector<std::vector<Scalar> > samples_;
224 
225  bool xExtrapolate_ = false;
226  bool yExtrapolate_ = false;
227 
231  template <class Evaluation>
232  unsigned xSegmentIndex_(const Evaluation& x) const
233  {
234  assert(xExtrapolate_ || appliesX(x) );
235 
236  return segmentIndex_(x, xPos_);
237  }
238 
242  template <class Evaluation>
243  unsigned ySegmentIndex_(const Evaluation& y) const
244  {
245  assert(yExtrapolate_ || appliesY(y) );
246 
247  return segmentIndex_(y, yPos_);
248  }
249 
250 
251  template <class Evaluation>
252  static unsigned segmentIndex_(const Evaluation& v, const std::vector<Scalar>& vPos)
253  {
254  const unsigned n = vPos.size();
255  assert(n >= 2);
256 
257  if (v <= vPos.front() || n == 2)
258  return 0;
259  else if (v >= vPos.back())
260  return n - 2;
261 
262  assert(n > 2 && v > vPos.front() && v < vPos.back());
263 
264  // bisection. this assumes that the vPos array is strictly mononically
265  // increasing.
266  size_t lowerIdx = 0;
267  size_t upperIdx = vPos.size() - 1;
268  while (lowerIdx + 1 < upperIdx) {
269  size_t pivotIdx = (lowerIdx + upperIdx) / 2;
270  if (v < vPos[pivotIdx])
271  upperIdx = pivotIdx;
272  else
273  lowerIdx = pivotIdx;
274  }
275 
276  assert(vPos[lowerIdx] <= v);
277  assert(v <= vPos[lowerIdx + 1]);
278  return lowerIdx;
279  }
280 
287  template <class Evaluation>
288  Evaluation xToAlpha(const Evaluation& x, unsigned xSegmentIdx) const
289  {
290  Scalar x1 = xPos_[xSegmentIdx];
291  Scalar x2 = xPos_[xSegmentIdx + 1];
292  return (x - x1)/(x2 - x1);
293  }
294 
301  template <class Evaluation>
302  Evaluation yToBeta(const Evaluation& y, unsigned ySegmentIdx) const
303  {
304  Scalar y1 = yPos_[ySegmentIdx];
305  Scalar y2 = yPos_[ySegmentIdx + 1];
306  return (y - y1)/(y2 - y1);
307  }
308 
309 };
310 } // namespace Opm
311 
312 #endif
Provides the opm-material specific exception classes.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Some templates to wrap the valgrind client request macros.
Implements a function that depends on two variables.
Definition: IntervalTabulated2DFunction.hpp:51
size_t numY() const
Returns the number of sampling points in Y direction.
Definition: IntervalTabulated2DFunction.hpp:104
Evaluation eval(const Evaluation &x, const Evaluation &y) const
Evaluate the function at a given (x,y) position.
Definition: IntervalTabulated2DFunction.hpp:190
bool appliesX(const Evaluation &x) const
Returns true if a coordinate lies in the tabulated range on the x direction.
Definition: IntervalTabulated2DFunction.hpp:171
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:110
Scalar yMin() const
Returns the minimum of the Y coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:122
size_t numX() const
Returns the number of sampling points in X direction.
Definition: IntervalTabulated2DFunction.hpp:98
bool appliesY(const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range on the y direction.
Definition: IntervalTabulated2DFunction.hpp:178
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range.
Definition: IntervalTabulated2DFunction.hpp:164
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:116
Scalar yMax() const
Returns the maximum of the Y coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:128
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition: IntervalTabulated2DFunction.hpp:157
Definition: Exceptions.hpp:46