UniformXTabulated2DFunction.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_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
29 #define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
30 
31 #include <opm/common/Valgrind.hpp>
32 #include <opm/common/Exceptions.hpp>
33 #include <opm/common/ErrorMacros.hpp>
34 #include <opm/common/Unused.hpp>
36 
37 #include <iostream>
38 #include <vector>
39 #include <limits>
40 #include <tuple>
41 
42 #include <assert.h>
43 
44 namespace Opm {
45 
54 template <class Scalar>
56 {
57  typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
58 
59 public:
61  { }
62 
66  Scalar xMin() const
67  { return xPos_.front(); }
68 
72  Scalar xMax() const
73  { return xPos_.back(); }
74 
78  Scalar xAt(size_t i) const
79  { return xPos_[i]; }
80 
84  Scalar yAt(size_t i, size_t j) const
85  { return std::get<1>(samples_[i][j]); }
86 
90  Scalar valueAt(size_t i, size_t j) const
91  { return std::get<2>(samples_[i][j]); }
92 
96  size_t numX() const
97  { return xPos_.size(); }
98 
102  Scalar yMin(size_t i) const
103  { return std::get<1>(samples_.at(i).front()); }
104 
108  Scalar yMax(size_t i) const
109  { return std::get<1>(samples_.at(i).back()); }
110 
114  size_t numY(size_t i) const
115  { return samples_.at(i).size(); }
116 
120  Scalar iToX(size_t i) const
121  {
122  assert(0 <= i && i < numX());
123 
124  return xPos_.at(i);
125  }
126 
130  Scalar jToY(size_t i, size_t j) const
131  {
132  assert(0 <= i && i < numX());
133  assert(0 <= j && size_t(j) < samples_[i].size());
134 
135  return std::get<1>(samples_.at(i).at(j));
136  }
137 
146  template <class Evaluation>
147  Evaluation xToI(const Evaluation& x, bool extrapolate OPM_OPTIM_UNUSED = false) const
148  {
149  assert(extrapolate || (xMin() <= x && x <= xMax()));
150 
151  // we need at least two sampling points!
152  assert(xPos_.size() >= 2);
153 
154  size_t segmentIdx;
155  if (x <= xPos_[1])
156  segmentIdx = 0;
157  else if (x >= xPos_[xPos_.size() - 2])
158  segmentIdx = xPos_.size() - 2;
159  else {
160  // bisection
161  segmentIdx = 1;
162  size_t upperIdx = xPos_.size() - 2;
163  while (segmentIdx + 1 < upperIdx) {
164  size_t pivotIdx = (segmentIdx + upperIdx) / 2;
165  if (x < xPos_[pivotIdx])
166  upperIdx = pivotIdx;
167  else
168  segmentIdx = pivotIdx;
169  }
170 
171  assert(xPos_[segmentIdx] <= x);
172  assert(x <= xPos_[segmentIdx + 1]);
173  }
174 
175  Scalar x1 = xPos_[segmentIdx];
176  Scalar x2 = xPos_[segmentIdx + 1];
177  return Scalar(segmentIdx) + (x - x1)/(x2 - x1);
178  }
179 
188  template <class Evaluation>
189  Evaluation yToJ(size_t i, const Evaluation& y, bool extrapolate OPM_OPTIM_UNUSED = false) const
190  {
191  assert(0 <= i && i < numX());
192  const auto& colSamplePoints = samples_.at(i);
193 
194  assert(extrapolate || (yMin(i) <= y && y <= yMax(i)));
195 
196  Scalar y1;
197  Scalar y2;
198 
199  // interval halving
200  size_t lowerIdx = 0;
201  size_t upperIdx = colSamplePoints.size() - 1;
202  size_t pivotIdx = (lowerIdx + upperIdx) / 2;
203  while (lowerIdx + 1 < upperIdx) {
204  if (y < std::get<1>(colSamplePoints[pivotIdx]))
205  upperIdx = pivotIdx;
206  else
207  lowerIdx = pivotIdx;
208  pivotIdx = (lowerIdx + upperIdx) / 2;
209  }
210 
211  y1 = std::get<1>(colSamplePoints[lowerIdx]);
212  y2 = std::get<1>(colSamplePoints[lowerIdx + 1]);
213 
214  assert(y1 <= y || (extrapolate && lowerIdx == 0));
215  assert(y <= y2 || (extrapolate && lowerIdx == colSamplePoints.size() - 2));
216 
217  return Scalar(lowerIdx) + (y - y1)/(y2 - y1);
218  }
219 
223  bool applies(Scalar x, Scalar y) const
224  {
225  if (x < xMin() || xMax() < x)
226  return false;
227 
228  Scalar i = xToI(x, /*extrapolate=*/false);
229  const auto& col1SamplePoints = samples_.at(unsigned(i));
230  const auto& col2SamplePoints = samples_.at(unsigned(i));
231  Scalar alpha = i - static_cast<int>(i);
232 
233  Scalar minY =
234  alpha*std::get<1>(col1SamplePoints.front()) +
235  (1 - alpha)*std::get<1>(col2SamplePoints.front());
236 
237  Scalar maxY =
238  alpha*std::get<1>(col1SamplePoints.back()) +
239  (1 - alpha)*std::get<1>(col2SamplePoints.back());
240 
241  return minY <= y && y <= maxY;
242  }
249  template <class Evaluation>
250  Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
251  {
252 #ifndef NDEBUG
253  if (!extrapolate && !applies(Opm::scalarValue(x), Opm::scalarValue(y))) {
254  OPM_THROW(NumericalProblem,
255  "Attempt to get undefined table value (" << x << ", " << y << ")");
256  };
257 #endif
258 
259  // bi-linear interpolation: first, calculate the x and y indices in the lookup
260  // table ...
261  Evaluation alpha = xToI(x, extrapolate);
262  size_t i =
263  static_cast<size_t>(std::max(0, std::min(static_cast<int>(numX()) - 2,
264  static_cast<int>(Opm::scalarValue(alpha)))));
265  alpha -= i;
266 
267  Evaluation beta1;
268  Evaluation beta2;
269 
270  beta1 = yToJ(i, y, extrapolate);
271  beta2 = yToJ(i + 1, y, extrapolate);
272 
273  size_t j1 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i)) - 2,
274  static_cast<int>(Opm::scalarValue(beta1)))));
275  size_t j2 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i + 1)) - 2,
276  static_cast<int>(Opm::scalarValue(beta2)))));
277 
278  beta1 -= j1;
279  beta2 -= j2;
280 
281  // evaluate the two function values for the same y value ...
282  Evaluation s1, s2;
283  s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
284  s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
285 
286  Valgrind::CheckDefined(s1);
287  Valgrind::CheckDefined(s2);
288 
289  // ... and finally combine them using x the position
290  Evaluation result;
291  result = s1*(1.0 - alpha) + s2*alpha;
292  Valgrind::CheckDefined(result);
293 
294  return result;
295  }
296 
302  size_t appendXPos(Scalar nextX)
303  {
304  if (xPos_.empty() || xPos_.back() < nextX) {
305  xPos_.push_back(nextX);
306  samples_.resize(xPos_.size());
307  return xPos_.size() - 1;
308  }
309  else if (xPos_.front() > nextX) {
310  // this is slow, but so what?
311  xPos_.insert(xPos_.begin(), nextX);
312  samples_.insert(samples_.begin(), std::vector<SamplePoint>());
313  return 0;
314  }
315  OPM_THROW(std::invalid_argument,
316  "Sampling points should be specified either monotonically "
317  "ascending or descending.");
318  }
319 
325  size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
326  {
327  assert(0 <= i && i < numX());
328 
329  Scalar x = iToX(i);
330  if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
331  samples_[i].push_back(SamplePoint(x, y, value));
332  return samples_[i].size() - 1;
333  }
334  else if (std::get<1>(samples_[i].front()) > y) {
335  // slow, but we still don't care...
336  samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
337  return 0;
338  }
339 
340  OPM_THROW(std::invalid_argument,
341  "Sampling points must be specified in either monotonically "
342  "ascending or descending order.");
343  }
344 
351  void print(std::ostream& os = std::cout) const
352  {
353  Scalar x0 = xMin();
354  Scalar x1 = xMax();
355  int m = numX();
356 
357  Scalar y0 = 1e30;
358  Scalar y1 = -1e30;
359  int n = 0;
360  for (int i = 0; i < m; ++ i) {
361  y0 = std::min(y0, yMin(i));
362  y1 = std::max(y1, yMax(i));
363  n = std::max(n, numY(i));
364  }
365 
366  m *= 3;
367  n *= 3;
368  for (int i = 0; i <= m; ++i) {
369  Scalar x = x0 + (x1 - x0)*i/m;
370  for (int j = 0; j <= n; ++j) {
371  Scalar y = y0 + (y1 - y0)*j/n;
372  os << x << " " << y << " " << eval(x, y) << "\n";
373  }
374  os << "\n";
375  }
376  }
377 
378 private:
379  // the vector which contains the values of the sample points
380  // f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
381  // instead!
382  std::vector<std::vector<SamplePoint> > samples_;
383 
384  // the position of each vertical line on the x-axis
385  std::vector<Scalar> xPos_;
386 };
387 } // namespace Opm
388 
389 #endif
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition: UniformXTabulated2DFunction.hpp:72
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
void print(std::ostream &os=std::cout) const
Print the table for debugging purposes.
Definition: UniformXTabulated2DFunction.hpp:351
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition: UniformXTabulated2DFunction.hpp:90
Scalar yAt(size_t i, size_t j) const
Returns the value of the Y coordinate of a sampling point.
Definition: UniformXTabulated2DFunction.hpp:84
Scalar yMin(size_t i) const
Returns the minimum of the Y coordinate of the sampling points for a given column.
Definition: UniformXTabulated2DFunction.hpp:102
Evaluation eval(const Evaluation &x, const Evaluation &y, bool extrapolate=false) const
Evaluate the function at a given (x,y) position.
Definition: UniformXTabulated2DFunction.hpp:250
Scalar xAt(size_t i) const
Returns the value of the X coordinate of the sampling points.
Definition: UniformXTabulated2DFunction.hpp:78
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition: UniformXTabulated2DFunction.hpp:66
size_t numY(size_t i) const
Returns the number of sampling points in Y direction a given column.
Definition: UniformXTabulated2DFunction.hpp:114
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition: UniformXTabulated2DFunction.hpp:55
Definition: Air_Mesitylene.hpp:33
size_t numX() const
Returns the number of sampling points in X direction.
Definition: UniformXTabulated2DFunction.hpp:96
Scalar iToX(size_t i) const
Return the position on the x-axis of the i-th interval.
Definition: UniformXTabulated2DFunction.hpp:120
Scalar jToY(size_t i, size_t j) const
Return the position on the y-axis of the j-th interval.
Definition: UniformXTabulated2DFunction.hpp:130
Scalar yMax(size_t i) const
Returns the maximum of the Y coordinate of the sampling points for a given column.
Definition: UniformXTabulated2DFunction.hpp:108
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
Append a sample point.
Definition: UniformXTabulated2DFunction.hpp:325
bool applies(Scalar x, Scalar y) const
Returns true iff a coordinate lies in the tabulated range.
Definition: UniformXTabulated2DFunction.hpp:223
Evaluation xToI(const Evaluation &x, bool extrapolate OPM_OPTIM_UNUSED=false) const
Return the interval index of a given position on the x-axis.
Definition: UniformXTabulated2DFunction.hpp:147
size_t appendXPos(Scalar nextX)
Set the x-position of a vertical line.
Definition: UniformXTabulated2DFunction.hpp:302
Evaluation yToJ(size_t i, const Evaluation &y, bool extrapolate OPM_OPTIM_UNUSED=false) const
Return the interval index of a given position on the y-axis.
Definition: UniformXTabulated2DFunction.hpp:189