My Project
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 
35 
36 #include <iostream>
37 #include <vector>
38 #include <limits>
39 #include <tuple>
40 #include <sstream>
41 #include <cassert>
42 #include <cmath>
43 
44 namespace Opm {
53 template <class Scalar>
55 {
56 public:
57  typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
58 
70  LeftExtreme,
71  RightExtreme,
72  Vertical
73  };
74 
75  explicit UniformXTabulated2DFunction(const InterpolationPolicy interpolationGuide = Vertical)
76  : interpolationGuide_(interpolationGuide)
77  { }
78 
79  UniformXTabulated2DFunction(const std::vector<Scalar>& xPos,
80  const std::vector<Scalar>& yPos,
81  const std::vector<std::vector<SamplePoint>>& samples,
82  InterpolationPolicy interpolationGuide)
83  : samples_(samples)
84  , xPos_(xPos)
85  , yPos_(yPos)
86  , interpolationGuide_(interpolationGuide)
87  { }
88 
92  Scalar xMin() const
93  { return xPos_.front(); }
94 
98  Scalar xMax() const
99  { return xPos_.back(); }
100 
104  Scalar xAt(size_t i) const
105  { return xPos_[i]; }
106 
110  Scalar yAt(size_t i, size_t j) const
111  { return std::get<1>(samples_[i][j]); }
112 
116  Scalar valueAt(size_t i, size_t j) const
117  { return std::get<2>(samples_[i][j]); }
118 
122  size_t numX() const
123  { return xPos_.size(); }
124 
128  Scalar yMin(unsigned i) const
129  { return std::get<1>(samples_.at(i).front()); }
130 
134  Scalar yMax(unsigned i) const
135  { return std::get<1>(samples_.at(i).back()); }
136 
140  size_t numY(unsigned i) const
141  { return samples_.at(i).size(); }
142 
146  Scalar iToX(unsigned i) const
147  {
148  assert(i < numX());
149 
150  return xPos_.at(i);
151  }
152 
153  const std::vector<std::vector<SamplePoint>>& samples() const
154  {
155  return samples_;
156  }
157 
158  const std::vector<Scalar>& xPos() const
159  {
160  return xPos_;
161  }
162 
163  const std::vector<Scalar>& yPos() const
164  {
165  return yPos_;
166  }
167 
168  InterpolationPolicy interpolationGuide() const
169  {
170  return interpolationGuide_;
171  }
172 
176  Scalar jToY(unsigned i, unsigned j) const
177  {
178  assert(i < numX());
179  assert(size_t(j) < samples_[i].size());
180 
181  return std::get<1>(samples_.at(i).at(j));
182  }
183 
187  template <class Evaluation>
188  unsigned xSegmentIndex(const Evaluation& x, bool extrapolate OPM_OPTIM_UNUSED = false) const
189  {
190  assert(extrapolate || (xMin() <= x && x <= xMax()));
191 
192  // we need at least two sampling points!
193  assert(xPos_.size() >= 2);
194 
195  if (x <= xPos_[1])
196  return 0;
197  else if (x >= xPos_[xPos_.size() - 2])
198  return xPos_.size() - 2;
199  else {
200  assert(xPos_.size() >= 3);
201 
202  // bisection
203  unsigned lowerIdx = 1;
204  unsigned upperIdx = xPos_.size() - 2;
205  while (lowerIdx + 1 < upperIdx) {
206  unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
207  if (x < xPos_[pivotIdx])
208  upperIdx = pivotIdx;
209  else
210  lowerIdx = pivotIdx;
211  }
212 
213  return lowerIdx;
214  }
215  }
216 
223  template <class Evaluation>
224  Evaluation xToAlpha(const Evaluation& x, unsigned segmentIdx) const
225  {
226  Scalar x1 = xPos_[segmentIdx];
227  Scalar x2 = xPos_[segmentIdx + 1];
228  return (x - x1)/(x2 - x1);
229  }
230 
234  template <class Evaluation>
235  unsigned ySegmentIndex(const Evaluation& y, unsigned xSampleIdx, bool extrapolate OPM_OPTIM_UNUSED = false) const
236  {
237  assert(xSampleIdx < numX());
238  const auto& colSamplePoints = samples_.at(xSampleIdx);
239 
240  assert(colSamplePoints.size() >= 2);
241  assert(extrapolate || (yMin(xSampleIdx) <= y && y <= yMax(xSampleIdx)));
242 
243  if (y <= std::get<1>(colSamplePoints[1]))
244  return 0;
245  else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
246  return colSamplePoints.size() - 2;
247  else {
248  assert(colSamplePoints.size() >= 3);
249 
250  // bisection
251  unsigned lowerIdx = 1;
252  unsigned upperIdx = colSamplePoints.size() - 2;
253  while (lowerIdx + 1 < upperIdx) {
254  unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
255  if (y < std::get<1>(colSamplePoints[pivotIdx]))
256  upperIdx = pivotIdx;
257  else
258  lowerIdx = pivotIdx;
259  }
260 
261  return lowerIdx;
262  }
263  }
264 
271  template <class Evaluation>
272  Evaluation yToBeta(const Evaluation& y, unsigned xSampleIdx, unsigned ySegmentIdx) const
273  {
274  assert(xSampleIdx < numX());
275  assert(ySegmentIdx < numY(xSampleIdx) - 1);
276 
277  const auto& colSamplePoints = samples_.at(xSampleIdx);
278 
279  Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
280  Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
281 
282  return (y - y1)/(y2 - y1);
283  }
284 
288  template <class Evaluation>
289  bool applies(const Evaluation& x, const Evaluation& y) const
290  {
291  if (x < xMin() || xMax() < x)
292  return false;
293 
294  unsigned i = xSegmentIndex(x, /*extrapolate=*/false);
295  Scalar alpha = xToAlpha(decay<Scalar>(x), i);
296 
297  const auto& col1SamplePoints = samples_.at(i);
298  const auto& col2SamplePoints = samples_.at(i + 1);
299 
300  Scalar minY =
301  alpha*std::get<1>(col1SamplePoints.front()) +
302  (1 - alpha)*std::get<1>(col2SamplePoints.front());
303 
304  Scalar maxY =
305  alpha*std::get<1>(col1SamplePoints.back()) +
306  (1 - alpha)*std::get<1>(col2SamplePoints.back());
307 
308  return minY <= y && y <= maxY;
309  }
316  template <class Evaluation>
317  Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
318  {
319 #ifndef NDEBUG
320  if (!extrapolate && !applies(x, y)) {
321  std::ostringstream oss;
322  oss << "Attempt to get undefined table value (" << x << ", " << y << ")";
323  throw NumericalIssue(oss.str());
324  };
325 #endif
326 
327  // bi-linear interpolation: first, calculate the x and y indices in the lookup
328  // table ...
329  unsigned i = xSegmentIndex(x, extrapolate);
330  const Evaluation& alpha = xToAlpha(x, i);
331  // The 'shift' is used to shift the points used to interpolate within
332  // the (i) and (i+1) sets of sample points, so that when approaching
333  // the boundary of the domain given by the samples, one gets the same
334  // value as one would get by interpolating along the boundary curve
335  // itself.
336  Evaluation shift = 0.0;
337  if (interpolationGuide_ == InterpolationPolicy::Vertical) {
338  // Shift is zero, no need to reset it.
339  } else {
340  // find upper and lower y value
341  if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
342  // The domain is above the boundary curve, up to y = infinity.
343  // The shift is therefore the same for all values of y.
344  shift = yPos_[i+1] - yPos_[i];
345  } else {
346  assert(interpolationGuide_ == InterpolationPolicy::RightExtreme);
347  // The domain is below the boundary curve, down to y = 0.
348  // The shift is therefore no longer the the same for all
349  // values of y, since at y = 0 the shift must be zero.
350  // The shift is computed by linear interpolation between
351  // the maximal value at the domain boundary curve, and zero.
352  shift = yPos_[i+1] - yPos_[i];
353  auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
354  if (yEnd > 0.) {
355  shift = shift * y / yEnd;
356  } else {
357  shift = 0.;
358  }
359  }
360  }
361  auto yLower = y - alpha*shift;
362  auto yUpper = y + (1-alpha)*shift;
363 
364  unsigned j1 = ySegmentIndex(yLower, i, extrapolate);
365  unsigned j2 = ySegmentIndex(yUpper, i + 1, extrapolate);
366  const Evaluation& beta1 = yToBeta(yLower, i, j1);
367  const Evaluation& beta2 = yToBeta(yUpper, i + 1, j2);
368 
369  // evaluate the two function values for the same y value ...
370  const Evaluation& s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
371  const Evaluation& s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
372 
373  Valgrind::CheckDefined(s1);
374  Valgrind::CheckDefined(s2);
375 
376  // ... and combine them using the x position
377  const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
378  Valgrind::CheckDefined(result);
379 
380  return result;
381  }
382 
388  size_t appendXPos(Scalar nextX)
389  {
390  if (xPos_.empty() || xPos_.back() < nextX) {
391  xPos_.push_back(nextX);
392  yPos_.push_back(-1e100);
393  samples_.push_back({});
394  return xPos_.size() - 1;
395  }
396  else if (xPos_.front() > nextX) {
397  // this is slow, but so what?
398  xPos_.insert(xPos_.begin(), nextX);
399  yPos_.insert(yPos_.begin(), -1e100);
400  samples_.insert(samples_.begin(), std::vector<SamplePoint>());
401  return 0;
402  }
403  throw std::invalid_argument("Sampling points should be specified either monotonically "
404  "ascending or descending.");
405  }
406 
412  size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
413  {
414  assert(i < numX());
415  Scalar x = iToX(i);
416  if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
417  samples_[i].push_back(SamplePoint(x, y, value));
418  if (interpolationGuide_ == InterpolationPolicy::RightExtreme) {
419  yPos_[i] = y;
420  }
421  return samples_[i].size() - 1;
422  }
423  else if (std::get<1>(samples_[i].front()) > y) {
424  // slow, but we still don't care...
425  samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
426  if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
427  yPos_[i] = y;
428  }
429  return 0;
430  }
431 
432  throw std::invalid_argument("Sampling points must be specified in either monotonically "
433  "ascending or descending order.");
434  }
435 
442  void print(std::ostream& os = std::cout) const
443  {
444  Scalar x0 = xMin();
445  Scalar x1 = xMax();
446  int m = numX();
447 
448  Scalar y0 = 1e30;
449  Scalar y1 = -1e30;
450  size_t n = 0;
451  for (int i = 0; i < m; ++ i) {
452  y0 = std::min(y0, yMin(i));
453  y1 = std::max(y1, yMax(i));
454  n = std::max(n, numY(i));
455  }
456 
457  m *= 3;
458  n *= 3;
459  for (int i = 0; i <= m; ++i) {
460  Scalar x = x0 + (x1 - x0)*i/m;
461  for (size_t j = 0; j <= n; ++j) {
462  Scalar y = y0 + (y1 - y0)*j/n;
463  os << x << " " << y << " " << eval(x, y) << "\n";
464  }
465  os << "\n";
466  }
467  }
468 
469  bool operator==(const UniformXTabulated2DFunction<Scalar>& data) const {
470  return this->xPos() == data.xPos() &&
471  this->yPos() == data.yPos() &&
472  this->samples() == data.samples() &&
473  this->interpolationGuide() == data.interpolationGuide();
474  }
475 
476 private:
477  // the vector which contains the values of the sample points
478  // f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
479  // instead!
480  std::vector<std::vector<SamplePoint> > samples_;
481 
482  // the position of each vertical line on the x-axis
483  std::vector<Scalar> xPos_;
484  // the position on the y-axis of the guide point
485  std::vector<Scalar> yPos_;
486  InterpolationPolicy interpolationGuide_;
487 };
488 } // namespace Opm
489 
490 #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.
Definition: Exceptions.hpp:46
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition: UniformXTabulated2DFunction.hpp:55
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition: UniformXTabulated2DFunction.hpp:98
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
Append a sample point.
Definition: UniformXTabulated2DFunction.hpp:412
unsigned ySegmentIndex(const Evaluation &y, unsigned xSampleIdx, bool extrapolate OPM_OPTIM_UNUSED=false) const
Return the interval index of a given position on the y-axis.
Definition: UniformXTabulated2DFunction.hpp:235
Evaluation xToAlpha(const Evaluation &x, unsigned segmentIdx) const
Return the relative position of an x value in an intervall.
Definition: UniformXTabulated2DFunction.hpp:224
Evaluation yToBeta(const Evaluation &y, unsigned xSampleIdx, unsigned ySegmentIdx) const
Return the relative position of an y value in an interval.
Definition: UniformXTabulated2DFunction.hpp:272
Scalar xAt(size_t i) const
Returns the value of the X coordinate of the sampling points.
Definition: UniformXTabulated2DFunction.hpp:104
void print(std::ostream &os=std::cout) const
Print the table for debugging purposes.
Definition: UniformXTabulated2DFunction.hpp:442
Scalar jToY(unsigned i, unsigned j) const
Return the position on the y-axis of the j-th interval.
Definition: UniformXTabulated2DFunction.hpp:176
InterpolationPolicy
Indicates how interpolation will be performed.
Definition: UniformXTabulated2DFunction.hpp:69
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition: UniformXTabulated2DFunction.hpp:116
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition: UniformXTabulated2DFunction.hpp:92
Scalar yMin(unsigned i) const
Returns the minimum of the Y coordinate of the sampling points for a given column.
Definition: UniformXTabulated2DFunction.hpp:128
size_t numY(unsigned i) const
Returns the number of sampling points in Y direction a given column.
Definition: UniformXTabulated2DFunction.hpp:140
size_t numX() const
Returns the number of sampling points in X direction.
Definition: UniformXTabulated2DFunction.hpp:122
Scalar yMax(unsigned i) const
Returns the maximum of the Y coordinate of the sampling points for a given column.
Definition: UniformXTabulated2DFunction.hpp:134
Scalar iToX(unsigned i) const
Return the position on the x-axis of the i-th interval.
Definition: UniformXTabulated2DFunction.hpp:146
Scalar yAt(size_t i, size_t j) const
Returns the value of the Y coordinate of a sampling point.
Definition: UniformXTabulated2DFunction.hpp:110
size_t appendXPos(Scalar nextX)
Set the x-position of a vertical line.
Definition: UniformXTabulated2DFunction.hpp:388
Evaluation eval(const Evaluation &x, const Evaluation &y, bool extrapolate=false) const
Evaluate the function at a given (x,y) position.
Definition: UniformXTabulated2DFunction.hpp:317
unsigned xSegmentIndex(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:188
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true iff a coordinate lies in the tabulated range.
Definition: UniformXTabulated2DFunction.hpp:289