28#ifndef OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
29#define OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
51template <
class Scalar>
58 template <
class DataContainer>
60 const std::vector<Scalar>& yPos,
61 const DataContainer& data,
62 const bool xExtrapolate =
false,
63 const bool yExtrapolate =
false)
67 , xExtrapolate_(xExtrapolate)
68 , yExtrapolate_(yExtrapolate)
73 for (
unsigned i = 0; i < xPos.size() - 1; ++ i) {
74 if (xPos[i + 1] <= xPos[i])
75 throw std::runtime_error(
"The array for the x-positions is not strictly increasing!");
78 for (
unsigned i = 0; i < yPos.size() - 1; ++ i) {
79 if (yPos[i + 1] <= yPos[i])
80 throw std::runtime_error(
"The array for the y-positions is not strictly increasing!");
85 if (
numX() != samples_.size())
86 throw std::runtime_error(
"numX() is not equal to the number of rows of the sampling points");
88 for (
unsigned xIdx = 0; xIdx <
numX(); ++xIdx) {
89 if (samples_[xIdx].size() !=
numY()) {
90 throw std::runtime_error(
"The " + std::to_string(xIdx) +
91 "-th row of the sampling points has "
92 "different size than numY() ");
101 {
return xPos_.size(); }
107 {
return yPos_.size(); }
113 {
return xPos_.front(); }
119 {
return xPos_.back(); }
125 {
return yPos_.front(); }
131 {
return yPos_.back(); }
133 const std::vector<Scalar>& xPos()
const
136 const std::vector<Scalar>& yPos()
const
139 const std::vector<std::vector<Scalar>>& samples()
const
142 bool xExtrapolate()
const
143 {
return xExtrapolate_; }
145 bool yExtrapolate()
const
146 {
return yExtrapolate_; }
148 bool operator==(
const IntervalTabulated2DFunction<Scalar>& data)
const {
149 return this->xPos() == data.xPos() &&
150 this->yPos() == data.yPos() &&
151 this->samples() == data.samples() &&
152 this->xExtrapolate() == data.xExtrapolate() &&
153 this->yExtrapolate() == data.yExtrapolate();
160 {
return samples_[i][j]; }
165 template <
class Evaluation>
166 bool applies(
const Evaluation& x,
const Evaluation& y)
const
172 template <
class Evaluation>
174 {
return xMin() <= x && x <=
xMax(); }
179 template <
class Evaluation>
181 {
return yMin() <= y && y <=
yMax(); }
191 template <
typename Evaluation>
192 Evaluation
eval(
const Evaluation& x,
const Evaluation& y)
const
195 if constexpr (std::is_floating_point_v<Evaluation>) {
197 std::to_string(x) +
", " +
198 std::to_string(y) +
")");
201 std::to_string(x.value()) +
", " +
202 std::to_string(y.value()) +
")");
208 const unsigned i = xSegmentIndex_(x);
209 const unsigned j = ySegmentIndex_(y);
212 const Evaluation alpha = xToAlpha(x, i);
213 const Evaluation beta = yToBeta(y, j);
215 const Evaluation s1 =
valueAt(i, j) * (1.0 - beta) +
valueAt(i, j + 1) * beta;
216 const Evaluation s2 =
valueAt(i + 1, j) * (1.0 - beta) +
valueAt(i + 1, j + 1) * beta;
218 Valgrind::CheckDefined(s1);
219 Valgrind::CheckDefined(s2);
222 return s1*(1.0 - alpha) + s2*alpha;
227 std::vector<Scalar> xPos_;
229 std::vector<Scalar> yPos_;
231 std::vector<std::vector<Scalar> > samples_;
233 bool xExtrapolate_ =
false;
234 bool yExtrapolate_ =
false;
239 template <
class Evaluation>
240 unsigned xSegmentIndex_(
const Evaluation& x)
const
242 assert(xExtrapolate_ ||
appliesX(x) );
244 return segmentIndex_(x, xPos_);
250 template <
class Evaluation>
251 unsigned ySegmentIndex_(
const Evaluation& y)
const
253 assert(yExtrapolate_ ||
appliesY(y) );
255 return segmentIndex_(y, yPos_);
259 template <
class Evaluation>
260 static unsigned segmentIndex_(
const Evaluation& v,
const std::vector<Scalar>& vPos)
262 const unsigned n = vPos.size();
265 if (v <= vPos.front() || n == 2)
267 else if (v >= vPos.back())
270 assert(n > 2 && v > vPos.front() && v < vPos.back());
275 size_t upperIdx = vPos.size() - 1;
276 while (lowerIdx + 1 < upperIdx) {
277 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
278 if (v < vPos[pivotIdx])
284 assert(vPos[lowerIdx] <= v);
285 assert(v <= vPos[lowerIdx + 1]);
295 template <
class Evaluation>
296 Evaluation xToAlpha(
const Evaluation& x,
unsigned xSegmentIdx)
const
298 Scalar x1 = xPos_[xSegmentIdx];
299 Scalar x2 = xPos_[xSegmentIdx + 1];
300 return (x - x1)/(x2 - x1);
309 template <
class Evaluation>
310 Evaluation yToBeta(
const Evaluation& y,
unsigned ySegmentIdx)
const
312 Scalar y1 = yPos_[ySegmentIdx];
313 Scalar y2 = yPos_[ySegmentIdx + 1];
314 return (y - y1)/(y2 - y1);
Provides the OPM specific exception classes.
Some templates to wrap the valgrind client request macros.
Implements a function that depends on two variables.
Definition IntervalTabulated2DFunction.hpp:53
size_t numY() const
Returns the number of sampling points in Y direction.
Definition IntervalTabulated2DFunction.hpp:106
Evaluation eval(const Evaluation &x, const Evaluation &y) const
Evaluate the function at a given (x,y) position.
Definition IntervalTabulated2DFunction.hpp:192
bool appliesX(const Evaluation &x) const
Returns true if a coordinate lies in the tabulated range on the x direction.
Definition IntervalTabulated2DFunction.hpp:173
Scalar xMin() const
Returns the minimum of the X coordinate of the sampling points.
Definition IntervalTabulated2DFunction.hpp:112
Scalar yMin() const
Returns the minimum of the Y coordinate of the sampling points.
Definition IntervalTabulated2DFunction.hpp:124
size_t numX() const
Returns the number of sampling points in X direction.
Definition IntervalTabulated2DFunction.hpp:100
bool appliesY(const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range on the y direction.
Definition IntervalTabulated2DFunction.hpp:180
bool applies(const Evaluation &x, const Evaluation &y) const
Returns true if a coordinate lies in the tabulated range.
Definition IntervalTabulated2DFunction.hpp:166
Scalar xMax() const
Returns the maximum of the X coordinate of the sampling points.
Definition IntervalTabulated2DFunction.hpp:118
Scalar yMax() const
Returns the maximum of the Y coordinate of the sampling points.
Definition IntervalTabulated2DFunction.hpp:130
Scalar valueAt(size_t i, size_t j) const
Returns the value of a sampling point.
Definition IntervalTabulated2DFunction.hpp:159
Definition Exceptions.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30