28 #ifndef OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP 29 #define OPM_UNIFORM_X_TABULATED_2D_FUNCTION_HPP 31 #include <opm/common/Valgrind.hpp> 32 #include <opm/common/Exceptions.hpp> 33 #include <opm/common/ErrorMacros.hpp> 34 #include <opm/common/Unused.hpp> 54 template <
class Scalar>
57 typedef std::tuple<Scalar, Scalar, Scalar> SamplePoint;
67 {
return xPos_.front(); }
73 {
return xPos_.back(); }
78 Scalar
xAt(
size_t i)
const 84 Scalar
yAt(
size_t i,
size_t j)
const 85 {
return std::get<1>(samples_[i][j]); }
91 {
return std::get<2>(samples_[i][j]); }
97 {
return xPos_.size(); }
103 {
return std::get<1>(samples_.at(i).front()); }
109 {
return std::get<1>(samples_.at(i).back()); }
115 {
return samples_.at(i).size(); }
122 assert(0 <= i && i <
numX());
130 Scalar
jToY(
size_t i,
size_t j)
const 132 assert(0 <= i && i <
numX());
133 assert(0 <= j &&
size_t(j) < samples_[i].size());
135 return std::get<1>(samples_.at(i).at(j));
146 template <
class Evaluation>
147 Evaluation
xToI(
const Evaluation& x,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const 149 assert(extrapolate || (
xMin() <= x && x <=
xMax()));
152 assert(xPos_.size() >= 2);
157 else if (x >= xPos_[xPos_.size() - 2])
158 segmentIdx = xPos_.size() - 2;
162 size_t upperIdx = xPos_.size() - 2;
163 while (segmentIdx + 1 < upperIdx) {
164 size_t pivotIdx = (segmentIdx + upperIdx) / 2;
165 if (x < xPos_[pivotIdx])
168 segmentIdx = pivotIdx;
171 assert(xPos_[segmentIdx] <= x);
172 assert(x <= xPos_[segmentIdx + 1]);
175 Scalar x1 = xPos_[segmentIdx];
176 Scalar x2 = xPos_[segmentIdx + 1];
177 return Scalar(segmentIdx) + (x - x1)/(x2 - x1);
188 template <
class Evaluation>
189 Evaluation
yToJ(
size_t i,
const Evaluation& y,
bool extrapolate OPM_OPTIM_UNUSED =
false)
const 191 assert(0 <= i && i <
numX());
192 const auto& colSamplePoints = samples_.at(i);
194 assert(extrapolate || (
yMin(i) <= y && y <=
yMax(i)));
201 size_t upperIdx = colSamplePoints.size() - 1;
202 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
203 while (lowerIdx + 1 < upperIdx) {
204 if (y < std::get<1>(colSamplePoints[pivotIdx]))
208 pivotIdx = (lowerIdx + upperIdx) / 2;
211 y1 = std::get<1>(colSamplePoints[lowerIdx]);
212 y2 = std::get<1>(colSamplePoints[lowerIdx + 1]);
214 assert(y1 <= y || (extrapolate && lowerIdx == 0));
215 assert(y <= y2 || (extrapolate && lowerIdx == colSamplePoints.size() - 2));
217 return Scalar(lowerIdx) + (y - y1)/(y2 - y1);
228 Scalar i =
xToI(x,
false);
229 const auto& col1SamplePoints = samples_.at(
unsigned(i));
230 const auto& col2SamplePoints = samples_.at(
unsigned(i));
231 Scalar alpha = i -
static_cast<int>(i);
234 alpha*std::get<1>(col1SamplePoints.front()) +
235 (1 - alpha)*std::get<1>(col2SamplePoints.front());
238 alpha*std::get<1>(col1SamplePoints.back()) +
239 (1 - alpha)*std::get<1>(col2SamplePoints.back());
241 return minY <= y && y <= maxY;
249 template <
class Evaluation>
250 Evaluation
eval(
const Evaluation& x,
const Evaluation& y,
bool extrapolate=
false)
const 253 if (!extrapolate && !
applies(Opm::scalarValue(x), Opm::scalarValue(y))) {
254 OPM_THROW(NumericalProblem,
255 "Attempt to get undefined table value (" << x <<
", " << y <<
")");
261 Evaluation alpha =
xToI(x, extrapolate);
263 static_cast<size_t>(std::max(0, std::min(static_cast<int>(
numX()) - 2,
264 static_cast<int>(Opm::scalarValue(alpha)))));
270 beta1 =
yToJ(i, y, extrapolate);
271 beta2 =
yToJ(i + 1, y, extrapolate);
273 size_t j1 =
static_cast<size_t>(std::max(0, std::min(static_cast<int>(
numY(i)) - 2,
274 static_cast<int>(Opm::scalarValue(beta1)))));
275 size_t j2 =
static_cast<size_t>(std::max(0, std::min(static_cast<int>(
numY(i + 1)) - 2,
276 static_cast<int>(Opm::scalarValue(beta2)))));
284 s2 =
valueAt(i + 1, j2)*(1.0 - beta2) +
valueAt(i + 1, j2 + 1)*beta2;
286 Valgrind::CheckDefined(s1);
287 Valgrind::CheckDefined(s2);
291 result = s1*(1.0 - alpha) + s2*alpha;
292 Valgrind::CheckDefined(result);
304 if (xPos_.empty() || xPos_.back() < nextX) {
305 xPos_.push_back(nextX);
306 samples_.resize(xPos_.size());
307 return xPos_.size() - 1;
309 else if (xPos_.front() > nextX) {
311 xPos_.insert(xPos_.begin(), nextX);
312 samples_.insert(samples_.begin(), std::vector<SamplePoint>());
315 OPM_THROW(std::invalid_argument,
316 "Sampling points should be specified either monotonically " 317 "ascending or descending.");
327 assert(0 <= i && i <
numX());
330 if (samples_[i].empty() || std::get<1>(samples_[i].back()) < y) {
331 samples_[i].push_back(SamplePoint(x, y, value));
332 return samples_[i].size() - 1;
334 else if (std::get<1>(samples_[i].front()) > y) {
336 samples_[i].insert(samples_[i].begin(), SamplePoint(x, y, value));
340 OPM_THROW(std::invalid_argument,
341 "Sampling points must be specified in either monotonically " 342 "ascending or descending order.");
351 void print(std::ostream& os = std::cout)
const 360 for (
int i = 0; i < m; ++ i) {
361 y0 = std::min(y0,
yMin(i));
362 y1 = std::max(y1,
yMax(i));
363 n = std::max(n,
numY(i));
368 for (
int i = 0; i <= m; ++i) {
369 Scalar x = x0 + (x1 - x0)*i/m;
370 for (
int j = 0; j <= n; ++j) {
371 Scalar y = y0 + (y1 - y0)*j/n;
372 os << x <<
" " << y <<
" " <<
eval(x, y) <<
"\n";
382 std::vector<std::vector<SamplePoint> > samples_;
385 std::vector<Scalar> xPos_;
Definition: Air_Mesitylene.hpp:33