28 #ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
29 #define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP
53 template <
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>
188 unsigned xSegmentIndex(
const Evaluation& x,
bool extrapolate OPM_OPTIM_UNUSED =
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>
235 unsigned ySegmentIndex(
const Evaluation& y,
unsigned xSampleIdx,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const
237 assert(xSampleIdx <
numX());
238 const auto& colSamplePoints = samples_.at(xSampleIdx);
240 assert(colSamplePoints.size() >= 2);
241 assert(extrapolate || (
yMin(xSampleIdx) <= y && y <=
yMax(xSampleIdx)));
243 if (y <= std::get<1>(colSamplePoints[1]))
245 else if (y >= std::get<1>(colSamplePoints[colSamplePoints.size() - 2]))
246 return colSamplePoints.size() - 2;
248 assert(colSamplePoints.size() >= 3);
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]))
271 template <
class Evaluation>
272 Evaluation
yToBeta(
const Evaluation& y,
unsigned xSampleIdx,
unsigned ySegmentIdx)
const
274 assert(xSampleIdx <
numX());
275 assert(ySegmentIdx <
numY(xSampleIdx) - 1);
277 const auto& colSamplePoints = samples_.at(xSampleIdx);
279 Scalar y1 = std::get<1>(colSamplePoints[ySegmentIdx]);
280 Scalar y2 = std::get<1>(colSamplePoints[ySegmentIdx + 1]);
282 return (y - y1)/(y2 - y1);
288 template <
class Evaluation>
289 bool applies(
const Evaluation& x,
const Evaluation& y)
const
295 Scalar alpha =
xToAlpha(decay<Scalar>(x), i);
297 const auto& col1SamplePoints = samples_.at(i);
298 const auto& col2SamplePoints = samples_.at(i + 1);
301 alpha*std::get<1>(col1SamplePoints.front()) +
302 (1 - alpha)*std::get<1>(col2SamplePoints.front());
305 alpha*std::get<1>(col1SamplePoints.back()) +
306 (1 - alpha)*std::get<1>(col2SamplePoints.back());
308 return minY <= y && y <= maxY;
316 template <
class Evaluation>
317 Evaluation
eval(
const Evaluation& x,
const Evaluation& y,
bool extrapolate=
false)
const
320 if (!extrapolate && !
applies(x, y)) {
321 std::ostringstream oss;
322 oss <<
"Attempt to get undefined table value (" << x <<
", " << y <<
")";
330 const Evaluation& alpha =
xToAlpha(x, i);
336 Evaluation shift = 0.0;
337 if (interpolationGuide_ == InterpolationPolicy::Vertical) {
341 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
344 shift = yPos_[i+1] - yPos_[i];
346 assert(interpolationGuide_ == InterpolationPolicy::RightExtreme);
352 shift = yPos_[i+1] - yPos_[i];
353 auto yEnd = yPos_[i]*(1.0 - alpha) + yPos_[i+1]*alpha;
355 shift = shift * y / yEnd;
361 auto yLower = y - alpha*shift;
362 auto yUpper = y + (1-alpha)*shift;
366 const Evaluation& beta1 =
yToBeta(yLower, i, j1);
367 const Evaluation& beta2 =
yToBeta(yUpper, i + 1, j2);
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;
373 Valgrind::CheckDefined(s1);
374 Valgrind::CheckDefined(s2);
377 const Evaluation& result = s1*(1.0 - alpha) + s2*alpha;
378 Valgrind::CheckDefined(result);
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;
396 else if (xPos_.front() > nextX) {
398 xPos_.insert(xPos_.begin(), nextX);
399 yPos_.insert(yPos_.begin(), -1e100);
400 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
403 throw std::invalid_argument(
"Sampling points should be specified either monotonically "
404 "ascending or descending.");
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) {
421 return samples_[i].size() - 1;
423 else if (std::get<1>(samples_[i].front()) > y) {
425 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
426 if (interpolationGuide_ == InterpolationPolicy::LeftExtreme) {
432 throw std::invalid_argument(
"Sampling points must be specified in either monotonically "
433 "ascending or descending order.");
442 void print(std::ostream& os = std::cout)
const
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));
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";
470 return this->xPos() == data.xPos() &&
471 this->yPos() == data.yPos() &&
472 this->samples() == data.samples() &&
473 this->interpolationGuide() == data.interpolationGuide();
480 std::vector<std::vector<SamplePoint> > samples_;
483 std::vector<Scalar> xPos_;
485 std::vector<Scalar> yPos_;
Provides the opm-material specific exception classes.
Provides the OPM_UNUSED macro.
Some templates to wrap the valgrind client request macros.
Definition: Exceptions.hpp:46