My Project
Loading...
Searching...
No Matches
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
32
35
36#include <cassert>
37#include <cmath>
38#include <iosfwd>
39#include <tuple>
40#include <type_traits>
41#include <vector>
42
43namespace Opm {
52template <class Scalar>
54{
55public:
56 typedef std::tuple</*x=*/Scalar, /*y=*/Scalar, /*value=*/Scalar> SamplePoint;
57
69 LeftExtreme,
70 RightExtreme,
71 Vertical
72 };
73
74 explicit UniformXTabulated2DFunction(const InterpolationPolicy interpolationGuide = Vertical)
75 : interpolationGuide_(interpolationGuide)
76 { }
77
78 UniformXTabulated2DFunction(const std::vector<Scalar>& xPos,
79 const std::vector<Scalar>& yPos,
80 const std::vector<std::vector<SamplePoint>>& samples,
81 InterpolationPolicy interpolationGuide)
82 : samples_(samples)
83 , xPos_(xPos)
84 , yPos_(yPos)
85 , interpolationGuide_(interpolationGuide)
86 { }
87
91 Scalar xMin() const
92 { return xPos_.front(); }
93
97 Scalar xMax() const
98 { return xPos_.back(); }
99
103 Scalar xAt(size_t i) const
104 { return xPos_[i]; }
105
109 Scalar yAt(size_t i, size_t j) const
110 { return std::get<1>(samples_[i][j]); }
111
115 Scalar valueAt(size_t i, size_t j) const
116 { return std::get<2>(samples_[i][j]); }
117
121 size_t numX() const
122 { return xPos_.size(); }
123
127 Scalar yMin(unsigned i) const
128 { return std::get<1>(samples_.at(i).front()); }
129
133 Scalar yMax(unsigned i) const
134 { return std::get<1>(samples_.at(i).back()); }
135
139 size_t numY(unsigned i) const
140 { return samples_.at(i).size(); }
141
145 Scalar iToX(unsigned i) const
146 {
147 assert(i < numX());
148
149 return xPos_.at(i);
150 }
151
152 const std::vector<std::vector<SamplePoint>>& samples() const
153 {
154 return samples_;
155 }
156
157 const std::vector<Scalar>& xPos() const
158 {
159 return xPos_;
160 }
161
162 const std::vector<Scalar>& yPos() const
163 {
164 return yPos_;
165 }
166
167 InterpolationPolicy interpolationGuide() const
168 {
169 return interpolationGuide_;
170 }
171
175 Scalar jToY(unsigned i, unsigned j) const
176 {
177 assert(i < numX());
178 assert(size_t(j) < samples_[i].size());
179
180 return std::get<1>(samples_.at(i).at(j));
181 }
182
186 template <class Evaluation>
187 unsigned xSegmentIndex(const Evaluation& x,
188 [[maybe_unused]] bool extrapolate = 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,
236 [[maybe_unused]] bool extrapolate = false) const
237 {
238 assert(xSampleIdx < numX());
239 const auto& colSamplePoints = samples_.at(xSampleIdx);
240
241 assert(colSamplePoints.size() >= 2);
242 assert(extrapolate || (yMin(xSampleIdx) <= y && y <= yMax(xSampleIdx)));
243
244 if (y <= std::get<1>(colSamplePoints[1]))
245 return 0;
246 else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
247 return colSamplePoints.size() - 2;
248 else {
249 assert(colSamplePoints.size() >= 3);
250
251 // bisection
252 unsigned lowerIdx = 1;
253 unsigned upperIdx = colSamplePoints.size() - 2;
254 while (lowerIdx + 1 < upperIdx) {
255 unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
256 if (y < std::get<1>(colSamplePoints[pivotIdx]))
257 upperIdx = pivotIdx;
258 else
259 lowerIdx = pivotIdx;
260 }
261
262 return lowerIdx;
263 }
264 }
265
272 template <class Evaluation>
273 Evaluation yToBeta(const Evaluation& y, unsigned xSampleIdx, unsigned ySegmentIdx) const
274 {
275 assert(xSampleIdx < numX());
276 assert(ySegmentIdx < numY(xSampleIdx) - 1);
277
278 const auto& colSamplePoints = samples_.at(xSampleIdx);
279
280 Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
281 Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
282
283 return (y - y1)/(y2 - y1);
284 }
285
289 template <class Evaluation>
290 bool applies(const Evaluation& x, const Evaluation& y) const
291 {
292 if (x < xMin() || xMax() < x)
293 return false;
294
295 unsigned i = xSegmentIndex(x, /*extrapolate=*/false);
296 Scalar alpha = xToAlpha(decay<Scalar>(x), i);
297
298 const auto& col1SamplePoints = samples_.at(i);
299 const auto& col2SamplePoints = samples_.at(i + 1);
300
301 Scalar minY =
302 alpha*std::get<1>(col1SamplePoints.front()) +
303 (1 - alpha)*std::get<1>(col2SamplePoints.front());
304
305 Scalar maxY =
306 alpha*std::get<1>(col1SamplePoints.back()) +
307 (1 - alpha)*std::get<1>(col2SamplePoints.back());
308
309 return minY <= y && y <= maxY;
310 }
317 template <class Evaluation>
318 Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
319 {
320 Evaluation alpha, beta1, beta2;
321 unsigned i, j1, j2;
322 findPoints(i, j1, j2, alpha, beta1, beta2, x, y, extrapolate);
323 return eval(i, j1, j2, alpha, beta1, beta2);
324 }
325
326 template <class Evaluation>
327 void findPoints(unsigned& i,
328 unsigned& j1,
329 unsigned& j2,
330 Evaluation& alpha,
331 Evaluation& beta1,
332 Evaluation& beta2,
333 const Evaluation& x,
334 const Evaluation& y,
335 bool extrapolate) const
336 {
337#ifndef NDEBUG
338 if (!extrapolate && !applies(x, y)) {
339 if constexpr (std::is_floating_point_v<Evaluation>) {
340 throw NumericalProblem("Attempt to get undefined table value (" +
341 std::to_string(x) + ", " +
342 std::to_string(y) + ")");
343 } else {
344 throw NumericalProblem("Attempt to get undefined table value (" +
345 std::to_string(x.value()) + ", " +
346 std::to_string(y.value()) + ")");
347 }
348 };
349#endif
350
351 // bi-linear interpolation: first, calculate the x and y indices in the lookup
352 // table ...
353 i = xSegmentIndex(x, extrapolate);
354 alpha = xToAlpha(x, i);
355 // The 'shift' is used to shift the points used to interpolate within
356 // the (i) and (i+1) sets of sample points, so that when approaching
357 // the boundary of the domain given by the samples, one gets the same
358 // value as one would get by interpolating along the boundary curve
359 // itself.
360 Evaluation shift = 0.0;
361 if (interpolationGuide_ == InterpolationPolicy::Vertical) {
362 // Shift is zero, no need to reset it.
363 } else {
364 // find upper and lower y value
365 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
366 // The domain is above the boundary curve, up to y = infinity.
367 // The shift is therefore the same for all values of y.
368 shift = yPos_[i+1] - yPos_[i];
369 } else {
370 assert(interpolationGuide_ == InterpolationPolicy::RightExtreme);
371 // The domain is below the boundary curve, down to y = 0.
372 // The shift is therefore no longer the the same for all
373 // values of y, since at y = 0 the shift must be zero.
374 // The shift is computed by linear interpolation between
375 // the maximal value at the domain boundary curve, and zero.
376 shift = yPos_[i+1] - yPos_[i];
377 auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
378 if (yEnd > 0.) {
379 shift = shift * y / yEnd;
380 } else {
381 shift = 0.;
382 }
383 }
384 }
385 auto yLower = y - alpha*shift;
386 auto yUpper = y + (1-alpha)*shift;
387
388 j1 = ySegmentIndex(yLower, i, extrapolate);
389 j2 = ySegmentIndex(yUpper, i + 1, extrapolate);
390 beta1 = yToBeta(yLower, i, j1);
391 beta2 = yToBeta(yUpper, i + 1, j2);
392 }
393
394 template <class Evaluation>
395 Evaluation eval(const unsigned& i, const unsigned& j1, const unsigned& j2, const Evaluation& alpha,const Evaluation& beta1,const Evaluation& beta2) const
396 {
397 // evaluate the two function values for the same y value ...
398 const Evaluation& s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
399 const Evaluation& s2 = valueAt(i + 1, j2)*(1.0 - beta2) + valueAt(i + 1, j2 + 1)*beta2;
400
401 Valgrind::CheckDefined(s1);
402 Valgrind::CheckDefined(s2);
403
404 // ... and combine them using the x position
405 const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
406 Valgrind::CheckDefined(result);
407
408 return result;
409 }
410
416 size_t appendXPos(Scalar nextX)
417 {
418 if (xPos_.empty() || xPos_.back() < nextX) {
419 xPos_.push_back(nextX);
420 yPos_.push_back(-1e100);
421 samples_.push_back({});
422 return xPos_.size() - 1;
423 }
424 else if (xPos_.front() > nextX) {
425 // this is slow, but so what?
426 xPos_.insert(xPos_.begin(), nextX);
427 yPos_.insert(yPos_.begin(), -1e100);
428 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
429 return 0;
430 }
431 throw std::invalid_argument("Sampling points should be specified either monotonically "
432 "ascending or descending.");
433 }
434
440 size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
441 {
442 assert(i < numX());
443 Scalar x = iToX(i);
444 if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
445 samples_[i].push_back(SamplePoint(x, y, value));
446 if (interpolationGuide_ == InterpolationPolicy::RightExtreme) {
447 yPos_[i] = y;
448 }
449 return samples_[i].size() - 1;
450 }
451 else if (std::get<1>(samples_[i].front()) > y) {
452 // slow, but we still don't care...
453 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
454 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
455 yPos_[i] = y;
456 }
457 return 0;
458 }
459
460 throw std::invalid_argument("Sampling points must be specified in either monotonically "
461 "ascending or descending order.");
462 }
463
470 void print(std::ostream& os) const;
471
472 bool operator==(const UniformXTabulated2DFunction<Scalar>& data) const {
473 return this->xPos() == data.xPos() &&
474 this->yPos() == data.yPos() &&
475 this->samples() == data.samples() &&
476 this->interpolationGuide() == data.interpolationGuide();
477 }
478
479private:
480 // the vector which contains the values of the sample points
481 // f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
482 // instead!
483 std::vector<std::vector<SamplePoint> > samples_;
484
485 // the position of each vertical line on the x-axis
486 std::vector<Scalar> xPos_;
487 // the position on the y-axis of the guide point
488 std::vector<Scalar> yPos_;
489 InterpolationPolicy interpolationGuide_;
490};
491} // namespace Opm
492
493#endif
Provides the OPM 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.
Definition Exceptions.hpp:40
Implements a scalar function that depends on two variables and which is sampled uniformly in the X di...
Definition UniformXTabulated2DFunction.hpp:54
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:97
size_t appendSamplePoint(size_t i, Scalar y, Scalar value)
Append a sample point.
Definition UniformXTabulated2DFunction.hpp:440
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:273
Scalar xAt(size_t i) const
Returns the value of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:103
Scalar jToY(unsigned i, unsigned j) const
Return the position on the y-axis of the j-th interval.
Definition UniformXTabulated2DFunction.hpp:175
InterpolationPolicy
Indicates how interpolation will be performed.
Definition UniformXTabulated2DFunction.hpp:68
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition UniformXTabulated2DFunction.hpp:115
unsigned ySegmentIndex(const Evaluation &y, unsigned xSampleIdx, bool extrapolate=false) const
Return the interval index of a given position on the y-axis.
Definition UniformXTabulated2DFunction.hpp:235
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition UniformXTabulated2DFunction.hpp:91
void print(std::ostream &os) const
Print the table for debugging purposes.
Definition UniformXTabulated2DFunction.cpp:32
Scalar yMin(unsigned i) const
Returns the minimum of the Y coordinate of the sampling points for a given column.
Definition UniformXTabulated2DFunction.hpp:127
unsigned xSegmentIndex(const Evaluation &x, bool extrapolate=false) const
Return the interval index of a given position on the x-axis.
Definition UniformXTabulated2DFunction.hpp:187
size_t numY(unsigned i) const
Returns the number of sampling points in Y direction a given column.
Definition UniformXTabulated2DFunction.hpp:139
size_t numX() const
Returns the number of sampling points in X direction.
Definition UniformXTabulated2DFunction.hpp:121
Scalar yMax(unsigned i) const
Returns the maximum of the Y coordinate of the sampling points for a given column.
Definition UniformXTabulated2DFunction.hpp:133
Scalar iToX(unsigned i) const
Return the position on the x-axis of the i-th interval.
Definition UniformXTabulated2DFunction.hpp:145
Scalar yAt(size_t i, size_t j) const
Returns the value of the Y coordinate of a sampling point.
Definition UniformXTabulated2DFunction.hpp:109
size_t appendXPos(Scalar nextX)
Set the x-position of a vertical line.
Definition UniformXTabulated2DFunction.hpp:416
Evaluation eval(const Evaluation &x, const Evaluation &y, bool extrapolate=false) const
Evaluate the function at a given (x,y) position.
Definition UniformXTabulated2DFunction.hpp:318
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true iff a coordinate lies in the tabulated range.
Definition UniformXTabulated2DFunction.hpp:290
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30