28#ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
29#define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
53template <
class Scalar>
57 typedef std::tuple<Scalar, Scalar, Scalar> SamplePoint;
76 : interpolationGuide_(interpolationGuide)
79 UniformXTabulated2DFunction(
const std::vector<Scalar>& xPos,
80 const std::vector<Scalar>& yPos,
81 const std::vector<std::vector<SamplePoint>>& samples,
86 , interpolationGuide_(interpolationGuide)
93 {
return xPos_.front(); }
99 {
return xPos_.back(); }
104 Scalar
xAt(
size_t i)
const
110 Scalar
yAt(
size_t i,
size_t j)
const
111 {
return std::get<1>(samples_[i][j]); }
117 {
return std::get<2>(samples_[i][j]); }
123 {
return xPos_.size(); }
129 {
return std::get<1>(samples_.at(i).front()); }
135 {
return std::get<1>(samples_.at(i).back()); }
141 {
return samples_.at(i).size(); }
153 const std::vector<std::vector<SamplePoint>>& samples()
const
158 const std::vector<Scalar>& xPos()
const
163 const std::vector<Scalar>& yPos()
const
170 return interpolationGuide_;
176 Scalar
jToY(
unsigned i,
unsigned j)
const
179 assert(
size_t(j) < samples_[i].size());
181 return std::get<1>(samples_.at(i).at(j));
187 template <
class Evaluation>
189 [[maybe_unused]]
bool extrapolate =
false)
const
191 assert(extrapolate || (
xMin() <= x && x <=
xMax()));
194 assert(xPos_.size() >= 2);
198 else if (x >= xPos_[xPos_.size() - 2])
199 return xPos_.size() - 2;
201 assert(xPos_.size() >= 3);
204 unsigned lowerIdx = 1;
205 unsigned upperIdx = xPos_.size() - 2;
206 while (lowerIdx + 1 < upperIdx) {
207 unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
208 if (x < xPos_[pivotIdx])
224 template <
class Evaluation>
225 Evaluation
xToAlpha(
const Evaluation& x,
unsigned segmentIdx)
const
227 Scalar x1 = xPos_[segmentIdx];
228 Scalar x2 = xPos_[segmentIdx + 1];
229 return (x - x1)/(x2 - x1);
235 template <
class Evaluation>
237 [[maybe_unused]]
bool extrapolate =
false)
const
239 assert(xSampleIdx <
numX());
240 const auto& colSamplePoints = samples_.at(xSampleIdx);
242 assert(colSamplePoints.size() >= 2);
243 assert(extrapolate || (
yMin(xSampleIdx) <= y && y <=
yMax(xSampleIdx)));
245 if (y <= std::get<1>(colSamplePoints[1]))
247 else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
248 return colSamplePoints.size() - 2;
250 assert(colSamplePoints.size() >= 3);
253 unsigned lowerIdx = 1;
254 unsigned upperIdx = colSamplePoints.size() - 2;
255 while (lowerIdx + 1 < upperIdx) {
256 unsigned pivotIdx = (lowerIdx + upperIdx) / 2;
257 if (y < std::get<1>(colSamplePoints[pivotIdx]))
273 template <
class Evaluation>
274 Evaluation
yToBeta(
const Evaluation& y,
unsigned xSampleIdx,
unsigned ySegmentIdx)
const
276 assert(xSampleIdx <
numX());
277 assert(ySegmentIdx <
numY(xSampleIdx) - 1);
279 const auto& colSamplePoints = samples_.at(xSampleIdx);
281 Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
282 Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
284 return (y - y1)/(y2 - y1);
290 template <
class Evaluation>
291 bool applies(
const Evaluation& x,
const Evaluation& y)
const
297 Scalar alpha =
xToAlpha(decay<Scalar>(x), i);
299 const auto& col1SamplePoints = samples_.at(i);
300 const auto& col2SamplePoints = samples_.at(i + 1);
303 alpha*std::get<1>(col1SamplePoints.front()) +
304 (1 - alpha)*std::get<1>(col2SamplePoints.front());
307 alpha*std::get<1>(col1SamplePoints.back()) +
308 (1 - alpha)*std::get<1>(col2SamplePoints.back());
310 return minY <= y && y <= maxY;
318 template <
class Evaluation>
319 Evaluation
eval(
const Evaluation& x,
const Evaluation& y,
bool extrapolate=
false)
const
321 Evaluation alpha, beta1, beta2;
323 findPoints(i, j1, j2, alpha, beta1, beta2, x, y, extrapolate);
324 return eval(i, j1, j2, alpha, beta1, beta2);
327 template <
class Evaluation>
328 void findPoints(
unsigned& i,
336 bool extrapolate)
const
339 if (!extrapolate && !
applies(x, y)) {
340 if constexpr (std::is_floating_point_v<Evaluation>) {
342 std::to_string(x) +
", " +
343 std::to_string(y) +
")");
345 throw NumericalProblem(
"Attempt to get undefined table value (" +
346 std::to_string(x.value()) +
", " +
347 std::to_string(y.value()) +
")");
361 Evaluation shift = 0.0;
362 if (interpolationGuide_ == InterpolationPolicy::Vertical) {
366 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
369 shift = yPos_[i+1] - yPos_[i];
371 assert(interpolationGuide_ == InterpolationPolicy::RightExtreme);
377 shift = yPos_[i+1] - yPos_[i];
378 auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
380 shift = shift * y / yEnd;
386 auto yLower = y - alpha*shift;
387 auto yUpper = y + (1-alpha)*shift;
391 beta1 =
yToBeta(yLower, i, j1);
392 beta2 =
yToBeta(yUpper, i + 1, j2);
395 template <
class Evaluation>
396 Evaluation
eval(
const unsigned& i,
const unsigned& j1,
const unsigned& j2,
const Evaluation& alpha,
const Evaluation& beta1,
const Evaluation& beta2)
const
399 const Evaluation& s1 =
valueAt(i, j1)*(1.0 - beta1) +
valueAt(i, j1 + 1)*beta1;
400 const Evaluation& s2 =
valueAt(i + 1, j2)*(1.0 - beta2) +
valueAt(i + 1, j2 + 1)*beta2;
402 Valgrind::CheckDefined(s1);
403 Valgrind::CheckDefined(s2);
406 const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
407 Valgrind::CheckDefined(result);
419 if (xPos_.empty() || xPos_.back() < nextX) {
420 xPos_.push_back(nextX);
421 yPos_.push_back(std::numeric_limits<Scalar>::lowest() / 2);
422 samples_.push_back({});
423 return xPos_.size() - 1;
425 else if (xPos_.front() > nextX) {
427 xPos_.insert(xPos_.begin(), nextX);
428 yPos_.insert(yPos_.begin(), std::numeric_limits<Scalar>::lowest() / 2);
429 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
432 throw std::invalid_argument(
"Sampling points should be specified either monotonically "
433 "ascending or descending.");
445 if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
446 samples_[i].push_back(SamplePoint(x, y, value));
447 if (interpolationGuide_ == InterpolationPolicy::RightExtreme) {
450 return samples_[i].size() - 1;
452 else if (std::get<1>(samples_[i].front()) > y) {
454 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
455 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
461 throw std::invalid_argument(
"Sampling points must be specified in either monotonically "
462 "ascending or descending order.");
471 void print(std::ostream& os)
const;
474 return this->xPos() == data.xPos() &&
475 this->yPos() == data.yPos() &&
476 this->samples() == data.samples() &&
477 this->interpolationGuide() == data.interpolationGuide();
484 std::vector<std::vector<SamplePoint> > samples_;
487 std::vector<Scalar> xPos_;
489 std::vector<Scalar> yPos_;
Provides the OPM specific exception classes.
Some templates to wrap the valgrind client request macros.
Definition Exceptions.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30