My Project
Loading...
Searching...
No Matches
IntervalTabulated2DFunction.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
29#define OPM_INTERVAL_TABULATED_2D_FUNCTION_HPP
30
32
35
36#include <algorithm>
37#include <cassert>
38#include <limits>
39#include <string>
40#include <type_traits>
41#include <vector>
42
43namespace Opm {
44
51template <class Scalar>
53{
54public:
56 { }
57
58 template <class DataContainer>
59 IntervalTabulated2DFunction(const std::vector<Scalar>& xPos,
60 const std::vector<Scalar>& yPos,
61 const DataContainer& data,
62 const bool xExtrapolate = false,
63 const bool yExtrapolate = false)
64 : xPos_(xPos)
65 , yPos_(yPos)
66 , samples_(data)
67 , xExtrapolate_(xExtrapolate)
68 , yExtrapolate_(yExtrapolate)
69 {
70#ifndef NDEBUG
71 // in debug mode, ensure that the x and y positions arrays are strictly
72 // mononically increasing.
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!");
76 }
77
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!");
81 }
82#endif
83
84 // make sure the size is correct
85 if (numX() != samples_.size())
86 throw std::runtime_error("numX() is not equal to the number of rows of the sampling points");
87
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() ");
93 }
94 }
95 }
96
100 size_t numX() const
101 { return xPos_.size(); }
102
106 size_t numY() const
107 { return yPos_.size(); }
108
112 Scalar xMin() const
113 { return xPos_.front(); }
114
118 Scalar xMax() const
119 { return xPos_.back(); }
120
124 Scalar yMin() const
125 { return yPos_.front(); }
126
130 Scalar yMax() const
131 { return yPos_.back(); }
132
133 const std::vector<Scalar>& xPos() const
134 { return xPos_; }
135
136 const std::vector<Scalar>& yPos() const
137 { return yPos_; }
138
139 const std::vector<std::vector<Scalar>>& samples() const
140 { return samples_; }
141
142 bool xExtrapolate() const
143 { return xExtrapolate_; }
144
145 bool yExtrapolate() const
146 { return yExtrapolate_; }
147
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();
154 }
155
159 Scalar valueAt(size_t i, size_t j) const
160 { return samples_[i][j]; }
161
165 template <class Evaluation>
166 bool applies(const Evaluation& x, const Evaluation& y) const
167 { return appliesX(x) && appliesY(y); }
168
172 template <class Evaluation>
173 bool appliesX(const Evaluation& x) const
174 { return xMin() <= x && x <= xMax(); }
175
179 template <class Evaluation>
180 bool appliesY(const Evaluation& y) const
181 { return yMin() <= y && y <= yMax(); }
182
183
191 template <typename Evaluation>
192 Evaluation eval(const Evaluation& x, const Evaluation& y) const
193 {
194 if ((!xExtrapolate_ && !appliesX(x)) || (!yExtrapolate_ && !appliesY(y))) {
195 if constexpr (std::is_floating_point_v<Evaluation>) {
196 throw NumericalProblem("Attempt to get undefined table value (" +
197 std::to_string(x) + ", " +
198 std::to_string(y) + ")");
199 } else {
200 throw NumericalProblem("Attempt to get undefined table value (" +
201 std::to_string(x.value()) + ", " +
202 std::to_string(y.value()) + ")");
203 }
204 };
205
206 // bi-linear interpolation: first, calculate the x and y indices in the lookup
207 // table ...
208 const unsigned i = xSegmentIndex_(x);
209 const unsigned j = ySegmentIndex_(y);
210
211 // bi-linear interpolation / extrapolation
212 const Evaluation alpha = xToAlpha(x, i);
213 const Evaluation beta = yToBeta(y, j);
214
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;
217
218 Valgrind::CheckDefined(s1);
219 Valgrind::CheckDefined(s2);
220
221 // ... and combine them using the x position
222 return s1*(1.0 - alpha) + s2*alpha;
223 }
224
225private:
226 // the sampling points in the x-drection
227 std::vector<Scalar> xPos_;
228 // the sampling points in the y-drection
229 std::vector<Scalar> yPos_;
230 // data at the sampling points
231 std::vector<std::vector<Scalar> > samples_;
232
233 bool xExtrapolate_ = false;
234 bool yExtrapolate_ = false;
235
239 template <class Evaluation>
240 unsigned xSegmentIndex_(const Evaluation& x) const
241 {
242 assert(xExtrapolate_ || appliesX(x) );
243
244 return segmentIndex_(x, xPos_);
245 }
246
250 template <class Evaluation>
251 unsigned ySegmentIndex_(const Evaluation& y) const
252 {
253 assert(yExtrapolate_ || appliesY(y) );
254
255 return segmentIndex_(y, yPos_);
256 }
257
258
259 template <class Evaluation>
260 static unsigned segmentIndex_(const Evaluation& v, const std::vector<Scalar>& vPos)
261 {
262 const unsigned n = vPos.size();
263 assert(n >= 2);
264
265 if (v <= vPos.front() || n == 2)
266 return 0;
267 else if (v >= vPos.back())
268 return n - 2;
269
270 assert(n > 2 && v > vPos.front() && v < vPos.back());
271
272 // bisection. this assumes that the vPos array is strictly mononically
273 // increasing.
274 size_t lowerIdx = 0;
275 size_t upperIdx = vPos.size() - 1;
276 while (lowerIdx + 1 < upperIdx) {
277 size_t pivotIdx = (lowerIdx + upperIdx) / 2;
278 if (v < vPos[pivotIdx])
279 upperIdx = pivotIdx;
280 else
281 lowerIdx = pivotIdx;
282 }
283
284 assert(vPos[lowerIdx] <= v);
285 assert(v <= vPos[lowerIdx + 1]);
286 return lowerIdx;
287 }
288
295 template <class Evaluation>
296 Evaluation xToAlpha(const Evaluation& x, unsigned xSegmentIdx) const
297 {
298 Scalar x1 = xPos_[xSegmentIdx];
299 Scalar x2 = xPos_[xSegmentIdx + 1];
300 return (x - x1)/(x2 - x1);
301 }
302
309 template <class Evaluation>
310 Evaluation yToBeta(const Evaluation& y, unsigned ySegmentIdx) const
311 {
312 Scalar y1 = yPos_[ySegmentIdx];
313 Scalar y2 = yPos_[ySegmentIdx + 1];
314 return (y - y1)/(y2 - y1);
315 }
316
317};
318} // namespace Opm
319
320#endif
Provides the OPM specific exception classes.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
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