28 #ifndef OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
29 #define OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
50 template <
class Scalar>
57 template <
class DataContainer>
59 const std::vector<Scalar>& yPos,
60 const DataContainer& data,
61 const bool xExtrapolate =
false,
62 const bool yExtrapolate =
false)
66 , xExtrapolate_(xExtrapolate)
67 , yExtrapolate_(yExtrapolate)
72 for (
unsigned i = 0; i < xPos.size() - 1; ++ i) {
73 if (xPos[i + 1] <= xPos[i])
74 throw std::runtime_error(
"The array for the x-positions is not strictly increasing!");
77 for (
unsigned i = 0; i < yPos.size() - 1; ++ i) {
78 if (yPos[i + 1] <= yPos[i])
79 throw std::runtime_error(
"The array for the y-positions is not strictly increasing!");
84 if (
numX() != samples_.size())
85 throw std::runtime_error(
"numX() is not equal to the number of rows of the sampling points");
87 for (
unsigned xIdx = 0; xIdx <
numX(); ++xIdx) {
88 if (samples_[xIdx].size() !=
numY()) {
89 std::ostringstream oss;
90 oss <<
"The " << xIdx <<
"-th row of the sampling points has different size than numY() ";
91 throw std::runtime_error(oss.str());
100 {
return xPos_.size(); }
106 {
return yPos_.size(); }
112 {
return xPos_.front(); }
118 {
return xPos_.back(); }
124 {
return yPos_.front(); }
130 {
return yPos_.back(); }
132 const std::vector<Scalar>& xPos()
const
135 const std::vector<Scalar>& yPos()
const
138 const std::vector<std::vector<Scalar>>& samples()
const
141 bool xExtrapolate()
const
142 {
return xExtrapolate_; }
144 bool yExtrapolate()
const
145 {
return yExtrapolate_; }
147 bool operator==(
const IntervalTabulated2DFunction<Scalar>& data)
const {
148 return this->xPos() == data.xPos() &&
149 this->yPos() == data.yPos() &&
150 this->samples() == data.samples() &&
151 this->xExtrapolate() == data.xExtrapolate() &&
152 this->yExtrapolate() == data.yExtrapolate();
159 {
return samples_[i][j]; }
164 template <
class Evaluation>
165 bool applies(
const Evaluation& x,
const Evaluation& y)
const
171 template <
class Evaluation>
173 {
return xMin() <= x && x <=
xMax(); }
178 template <
class Evaluation>
180 {
return yMin() <= y && y <=
yMax(); }
190 template <
typename Evaluation>
191 Evaluation
eval(
const Evaluation& x,
const Evaluation& y)
const
194 std::ostringstream oss;
195 oss <<
"Attempt to get undefined table value (" << x <<
", " << y <<
")";
201 const unsigned i = xSegmentIndex_(x);
202 const unsigned j = ySegmentIndex_(y);
205 const Evaluation alpha = xToAlpha(x, i);
206 const Evaluation beta = yToBeta(y, j);
208 const Evaluation s1 =
valueAt(i, j) * (1.0 - beta) +
valueAt(i, j + 1) * beta;
209 const Evaluation s2 =
valueAt(i + 1, j) * (1.0 - beta) +
valueAt(i + 1, j + 1) * beta;
211 Valgrind::CheckDefined(s1);
212 Valgrind::CheckDefined(s2);
215 return s1*(1.0 - alpha) + s2*alpha;
220 std::vector<Scalar> xPos_;
222 std::vector<Scalar> yPos_;
224 std::vector<std::vector<Scalar> > samples_;
226 bool xExtrapolate_ =
false;
227 bool yExtrapolate_ =
false;
232 template <
class Evaluation>
233 unsigned xSegmentIndex_(
const Evaluation& x)
const
235 assert(xExtrapolate_ ||
appliesX(x) );
237 return segmentIndex_(x, xPos_);
243 template <
class Evaluation>
244 unsigned ySegmentIndex_(
const Evaluation& y)
const
246 assert(yExtrapolate_ ||
appliesY(y) );
248 return segmentIndex_(y, yPos_);
252 template <
class Evaluation>
253 static unsigned segmentIndex_(
const Evaluation& v,
const std::vector<Scalar>& vPos)
255 const unsigned n = vPos.size();
258 if (v <= vPos.front() || n == 2)
260 else if (v >= vPos.back())
263 assert(n > 2 && v > vPos.front() && v < vPos.back());
268 size_t upperIdx = vPos.size() - 1;
269 while (lowerIdx + 1 < upperIdx) {
270 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
271 if (v < vPos[pivotIdx])
277 assert(vPos[lowerIdx] <= v);
278 assert(v <= vPos[lowerIdx + 1]);
288 template <
class Evaluation>
289 Evaluation xToAlpha(
const Evaluation& x,
unsigned xSegmentIdx)
const
291 Scalar x1 = xPos_[xSegmentIdx];
292 Scalar x2 = xPos_[xSegmentIdx + 1];
293 return (x - x1)/(x2 - x1);
302 template <
class Evaluation>
303 Evaluation yToBeta(
const Evaluation& y,
unsigned ySegmentIdx)
const
305 Scalar y1 = yPos_[ySegmentIdx];
306 Scalar y2 = yPos_[ySegmentIdx + 1];
307 return (y - y1)/(y2 - y1);
Provides the opm-material specific exception classes.
Provides the OPM_UNUSED macro.
Some templates to wrap the valgrind client request macros.
Implements a function that depends on two variables.
Definition: IntervalTabulated2DFunction.hpp:52
size_t numY() const
Returns the number of sampling points in Y direction.
Definition: IntervalTabulated2DFunction.hpp:105
Evaluation eval(const Evaluation &x, const Evaluation &y) const
Evaluate the function at a given (x,y) position.
Definition: IntervalTabulated2DFunction.hpp:191
bool appliesX(const Evaluation &x) const
Returns true if a coordinate lies in the tabulated range on the x direction.
Definition: IntervalTabulated2DFunction.hpp:172
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:111
Scalar yMin() const
Returns the minimum of the Y coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:123
size_t numX() const
Returns the number of sampling points in X direction.
Definition: IntervalTabulated2DFunction.hpp:99
bool appliesY(const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range on the y direction.
Definition: IntervalTabulated2DFunction.hpp:179
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range.
Definition: IntervalTabulated2DFunction.hpp:165
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:117
Scalar yMax() const
Returns the maximum of the Y coordinate of the sampling points.
Definition: IntervalTabulated2DFunction.hpp:129
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition: IntervalTabulated2DFunction.hpp:158
Definition: Exceptions.hpp:46