28 #ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
29 #define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
52 template <
class Scalar>
56 typedef std::tuple<Scalar, Scalar, Scalar> SamplePoint;
75 : interpolationGuide_(interpolationGuide)
78 UniformXTabulated2DFunction(
const std::vector<Scalar>& xPos,
79 const std::vector<Scalar>& yPos,
80 const std::vector<std::vector<SamplePoint>>& samples,
85 , interpolationGuide_(interpolationGuide)
92 {
return xPos_.front(); }
98 {
return xPos_.back(); }
103 Scalar
xAt(
size_t i)
const
109 Scalar
yAt(
size_t i,
size_t j)
const
110 {
return std::get<1>(samples_[i][j]); }
116 {
return std::get<2>(samples_[i][j]); }
122 {
return xPos_.size(); }
128 {
return std::get<1>(samples_.at(i).front()); }
134 {
return std::get<1>(samples_.at(i).back()); }
140 {
return samples_.at(i).size(); }
152 const std::vector<std::vector<SamplePoint>>& samples()
const
157 const std::vector<Scalar>& xPos()
const
162 const std::vector<Scalar>& yPos()
const
169 return interpolationGuide_;
175 Scalar
jToY(
unsigned i,
unsigned j)
const
178 assert(
size_t(j) < samples_[i].size());
180 return std::get<1>(samples_.at(i).at(j));
186 template <
class Evaluation>
188 [[maybe_unused]]
bool extrapolate =
false)
const
190 assert(extrapolate || (
xMin() <= x && x <=
xMax()));
193 assert(xPos_.size() >= 2);
197 else if (x >= xPos_[xPos_.size() - 2])
198 return xPos_.size() - 2;
200 assert(xPos_.size() >= 3);
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])
223 template <
class Evaluation>
224 Evaluation
xToAlpha(
const Evaluation& x,
unsigned segmentIdx)
const
226 Scalar x1 = xPos_[segmentIdx];
227 Scalar x2 = xPos_[segmentIdx + 1];
228 return (x - x1)/(x2 - x1);
234 template <
class Evaluation>
236 [[maybe_unused]]
bool extrapolate =
false)
const
238 assert(xSampleIdx <
numX());
239 const auto& colSamplePoints = samples_.at(xSampleIdx);
241 assert(colSamplePoints.size() >= 2);
242 assert(extrapolate || (
yMin(xSampleIdx) <= y && y <=
yMax(xSampleIdx)));
244 if (y <= std::get<1>(colSamplePoints[1]))
246 else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
247 return colSamplePoints.size() - 2;
249 assert(colSamplePoints.size() >= 3);
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]))
272 template <
class Evaluation>
273 Evaluation
yToBeta(
const Evaluation& y,
unsigned xSampleIdx,
unsigned ySegmentIdx)
const
275 assert(xSampleIdx <
numX());
276 assert(ySegmentIdx <
numY(xSampleIdx) - 1);
278 const auto& colSamplePoints = samples_.at(xSampleIdx);
280 Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
281 Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
283 return (y - y1)/(y2 - y1);
289 template <
class Evaluation>
290 bool applies(
const Evaluation& x,
const Evaluation& y)
const
296 Scalar alpha =
xToAlpha(decay<Scalar>(x), i);
298 const auto& col1SamplePoints = samples_.at(i);
299 const auto& col2SamplePoints = samples_.at(i + 1);
302 alpha*std::get<1>(col1SamplePoints.front()) +
303 (1 - alpha)*std::get<1>(col2SamplePoints.front());
306 alpha*std::get<1>(col1SamplePoints.back()) +
307 (1 - alpha)*std::get<1>(col2SamplePoints.back());
309 return minY <= y && y <= maxY;
317 template <
class Evaluation>
318 Evaluation
eval(
const Evaluation& x,
const Evaluation& y,
bool extrapolate=
false)
const
321 if (!extrapolate && !
applies(x, y)) {
322 std::ostringstream oss;
323 oss <<
"Attempt to get undefined table value (" << x <<
", " << y <<
")";
331 const Evaluation& alpha =
xToAlpha(x, i);
337 Evaluation shift = 0.0;
338 if (interpolationGuide_ == InterpolationPolicy::Vertical) {
342 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
345 shift = yPos_[i+1] - yPos_[i];
347 assert(interpolationGuide_ == InterpolationPolicy::RightExtreme);
353 shift = yPos_[i+1] - yPos_[i];
354 auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
356 shift = shift * y / yEnd;
362 auto yLower = y - alpha*shift;
363 auto yUpper = y + (1-alpha)*shift;
367 const Evaluation& beta1 =
yToBeta(yLower, i, j1);
368 const Evaluation& beta2 =
yToBeta(yUpper, i + 1, j2);
371 const Evaluation& s1 =
valueAt(i, j1)*(1.0 - beta1) +
valueAt(i, j1 + 1)*beta1;
372 const Evaluation& s2 =
valueAt(i + 1, j2)*(1.0 - beta2) +
valueAt(i + 1, j2 + 1)*beta2;
374 Valgrind::CheckDefined(s1);
375 Valgrind::CheckDefined(s2);
378 const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
379 Valgrind::CheckDefined(result);
391 if (xPos_.empty() || xPos_.back() < nextX) {
392 xPos_.push_back(nextX);
393 yPos_.push_back(-1e100);
394 samples_.push_back({});
395 return xPos_.size() - 1;
397 else if (xPos_.front() > nextX) {
399 xPos_.insert(xPos_.begin(), nextX);
400 yPos_.insert(yPos_.begin(), -1e100);
401 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
404 throw std::invalid_argument(
"Sampling points should be specified either monotonically "
405 "ascending or descending.");
417 if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
418 samples_[i].push_back(SamplePoint(x, y, value));
419 if (interpolationGuide_ == InterpolationPolicy::RightExtreme) {
422 return samples_[i].size() - 1;
424 else if (std::get<1>(samples_[i].front()) > y) {
426 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
427 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
433 throw std::invalid_argument(
"Sampling points must be specified in either monotonically "
434 "ascending or descending order.");
443 void print(std::ostream& os = std::cout)
const
452 for (
int i = 0; i < m; ++ i) {
453 y0 = std::min(y0,
yMin(i));
454 y1 = std::max(y1,
yMax(i));
455 n = std::max(n,
numY(i));
460 for (
int i = 0; i <= m; ++i) {
461 Scalar x = x0 + (x1 - x0)*i/m;
462 for (
size_t j = 0; j <= n; ++j) {
463 Scalar y = y0 + (y1 - y0)*j/n;
464 os << x <<
" " << y <<
" " <<
eval(x, y) <<
"\n";
471 return this->xPos() == data.xPos() &&
472 this->yPos() == data.yPos() &&
473 this->samples() == data.samples() &&
474 this->interpolationGuide() == data.interpolationGuide();
481 std::vector<std::vector<SamplePoint> > samples_;
484 std::vector<Scalar> xPos_;
486 std::vector<Scalar> yPos_;
Provides the opm-material specific exception classes.
Some templates to wrap the valgrind client request macros.
Definition: Exceptions.hpp:46