OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Remove trailing whitespaces

Browse Source
This commit is contained in:
Júlio Hoffimann
2013-07-28 08:34:13 -03:00
parent 78bed975ed
commit f4244e9382
3 changed files with 228 additions and 228 deletions
Download Patch File Download Diff File Expand all files Collapse all files

236
opm/core/utility/MonotCubicInterpolator.cpp
View File

@@ -1,17 +1,17 @@
/*
MonotCubicInterpolator
Copyright (C) 2006 Statoil ASA
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
@@ -24,7 +24,7 @@
Class to represent a one-dimensional function with single-valued
argument. Cubic interpolation, preserving the monotonicity of the
discrete known function values
@see MonotCubicInterpolator.h for further documentation.
*/
@@ -50,8 +50,8 @@ using namespace std;
Internal data structure of points and values:
vector(s):
- Easier coding
vector(s):
- Easier coding
- Faster vector operations when setting up derivatives.
- sorting slightly more complex.
- insertion of further values bad.
@@ -61,11 +61,11 @@ using namespace std;
- code complexity almost as for map.
- insertion of additional values bad
vector over struct datapoint { x, f, d}
vector over struct datapoint { x, f, d}
- nice code
- not as sortable, insertion is cumbersome.
** This is used currently: **
** This is used currently: **
map<double, double> one for (x,f) and one for (x,d)
- Naturally sorted on x-values (done by the map-construction)
- Slower to set up, awkward loop coding (?)
@@ -79,7 +79,7 @@ using namespace std;
- even more unreadable(??) code?
- higher skills needed by programmer.
MONOTONE CUBIC INTERPOLATION:
It is a local algorithm. Adding one point only incur recomputation
@@ -87,7 +87,7 @@ using namespace std;
divided).
NOTE: We do not currently make use of this local fact. Upon
insertion of a new data pair, everything is recomputed. Revisit
insertion of a new data pair, everything is recomputed. Revisit
this when needed.
*/
@@ -102,7 +102,7 @@ MonotCubicInterpolator(const vector<double> & x, const vector<double> & f) {
if (x.size() != f.size()) {
throw("Unable to constuct MonotCubicInterpolator from vectors.") ;
}
// Add the contents of the input vectors to our map of values.
vector<double>::const_iterator posx, posf;
for (posx = x.begin(), posf = f.begin(); posx != x.end(); ++posx, ++posf) {
@@ -114,9 +114,9 @@ MonotCubicInterpolator(const vector<double> & x, const vector<double> & f) {
bool
bool
MonotCubicInterpolator::
read(const std::string & datafilename, int xColumn, int fColumn)
read(const std::string & datafilename, int xColumn, int fColumn)
{
data.clear() ;
ddata.clear() ;
@@ -125,7 +125,7 @@ read(const std::string & datafilename, int xColumn, int fColumn)
if (!datafile_fs) {
return false ;
}
string linestring;
while (!datafile_fs.eof()) {
getline(datafile_fs, linestring);
@@ -140,7 +140,7 @@ read(const std::string & datafilename, int xColumn, int fColumn)
stringstream strs(linestring);
int columnindex = 0;
std::vector<double> value;
std::vector<double> value;
if (linestring.size() > 0 && linestring.at(0) != '#') {
while (!(strs.rdstate() & std::ios::failbit)) {
double tmp;
@@ -152,34 +152,34 @@ read(const std::string & datafilename, int xColumn, int fColumn)
if (columnindex >= (max(xColumn, fColumn))) {
data[value[xColumn-1]] = value[fColumn-1] ;
}
}
}
datafile_fs.close();
if (data.size() == 0) {
return false ;
}
computeInternalFunctionData();
return true ;
}
void
void
MonotCubicInterpolator::
addPair(double newx, double newf) throw(const char*) {
if (std::isnan(newx) || std::isinf(newx) || std::isnan(newf) || std::isinf(newf)) {
throw("MonotCubicInterpolator: addPair() received inf/nan input.");
}
data[newx] = newf ;
// In a critical application, we should only update the
// internal function data for the offended interval,
// internal function data for the offended interval,
// not for all function values over again.
computeInternalFunctionData();
}
double
double
MonotCubicInterpolator::
evaluate(double x) const throw(const char*){
@@ -189,7 +189,7 @@ evaluate(double x) const throw(const char*){
// xf becomes the first (xdata,fdata) pair where xdata >= x
map<double,double>::const_iterator xf_iterator = data.lower_bound(x);
// First check if we must extrapolate:
if (xf_iterator == data.begin()) {
if (data.begin()->first == x) { // Just on the interval limit
@@ -204,36 +204,36 @@ evaluate(double x) const throw(const char*){
// Constant extrapolation (!!)
return data.rbegin()->second;
}
// Ok, we have x_min < x < x_max
pair<double,double> xf2 = *xf_iterator;
pair<double,double> xf1 = *(--xf_iterator);
// we now have: xf2.first > x
// we now have: xf2.first > x
// Linear interpolation if derivative data is not available:
if (ddata.size() != data.size()) {
double finterp = xf1.second +
(xf2.second - xf1.second) / (xf2.first - xf1.first)
double finterp = xf1.second +
(xf2.second - xf1.second) / (xf2.first - xf1.first)
* (x - xf1.first);
return finterp;
}
else { // Do Cubic Hermite spline
double t = (x - xf1.first)/(xf2.first - xf1.first); // t \in [0,1]
double h = xf2.first - xf1.first;
double finterp
= xf1.second * H00(t)
double finterp
= xf1.second * H00(t)
+ ddata[xf1.first] * H10(t) * h
+ xf2.second * H01(t)
+ xf2.second * H01(t)
+ ddata[xf2.first] * H11(t) * h ;
return finterp;
}
}
// double
// double
// MonotCubicInterpolator::
// evaluate(double x, double& errorestimate_output) {
// cout << "Error: errorestimate not implemented" << endl;
@@ -241,12 +241,12 @@ evaluate(double x) const throw(const char*){
// return x;
// }
vector<double>
vector<double>
MonotCubicInterpolator::
get_xVector() const
get_xVector() const
{
map<double,double>::const_iterator xf_iterator = data.begin();
vector<double> outputvector;
outputvector.reserve(data.size());
for (xf_iterator = data.begin(); xf_iterator != data.end(); ++xf_iterator) {
@@ -258,11 +258,11 @@ get_xVector() const
vector<double>
MonotCubicInterpolator::
get_fVector() const
get_fVector() const
{
map<double,double>::const_iterator xf_iterator = data.begin();
vector<double> outputvector;
outputvector.reserve(data.size());
for (xf_iterator = data.begin(); xf_iterator != data.end(); ++xf_iterator) {
@@ -275,7 +275,7 @@ get_fVector() const
string
MonotCubicInterpolator::
toString() const
toString() const
{
const int precision = 20;
std::string dataString;
@@ -295,9 +295,9 @@ toString() const
dataStringStream << setprecision(precision) << it->second;
dataStringStream << '\n';
}
return dataStringStream.str();
}
@@ -306,31 +306,31 @@ MonotCubicInterpolator::
getMissingX() const throw(const char*)
{
if( data.size() < 2) {
throw("MonotCubicInterpolator::getMissingX() only one datapoint.");
throw("MonotCubicInterpolator::getMissingX() only one datapoint.");
}
// Search for biggest difference value in function-datavalues:
map<double,double>::const_iterator maxfDiffPair_iterator = data.begin();;
double maxfDiffValue = 0;
map<double,double>::const_iterator xf_iterator;
map<double,double>::const_iterator xf_next_iterator;
for (xf_iterator = data.begin(), xf_next_iterator = ++(data.begin());
xf_next_iterator != data.end();
xf_next_iterator != data.end();
++xf_iterator, ++xf_next_iterator) {
double absfDiff = fabs((double)(*xf_next_iterator).second
double absfDiff = fabs((double)(*xf_next_iterator).second
- (double)(*xf_iterator).second);
if (absfDiff > maxfDiffValue) {
maxfDiffPair_iterator = xf_iterator;
maxfDiffValue = absfDiff;
}
}
double newXvalue = ((*maxfDiffPair_iterator).first + ((*(++maxfDiffPair_iterator)).first))/2;
return make_pair(newXvalue, maxfDiffValue);
}
@@ -382,33 +382,33 @@ getMinimumF() const throw(const char*) {
}
void
void
MonotCubicInterpolator::
computeInternalFunctionData() const {
computeInternalFunctionData() const {
/* The contents of this function is meaningless if there is only one datapoint */
if (data.size() <= 1) {
return;
}
/* We do not check the caching flag if we are instructed
to do this computation */
/* We compute monotoneness and directions by assuming
monotoneness, and setting to false if the function is not for
some value */
map<double,double>::const_iterator xf_iterator;
map<double,double>::const_iterator xf_next_iterator;
strictlyMonotone = true; // We assume this is true, and will set to false if not
monotone = true;
strictlyDecreasing = true;
decreasing = true;
strictlyIncreasing = true;
increasing = true;
// Increasing or decreasing??
xf_iterator = data.begin();
xf_next_iterator = ++(data.begin());
@@ -419,12 +419,12 @@ computeInternalFunctionData() const {
strictlyMonotone = false;
strictlyIncreasing = false;
strictlyDecreasing = false;
++xf_iterator;
++xf_next_iterator;
}
if (xf_next_iterator != data.end()) {
if ( xf_iterator->second > xf_next_iterator->second) {
@@ -465,27 +465,27 @@ computeInternalFunctionData() const {
}
}
}
else {
// first two values must be equal if we get
else {
// first two values must be equal if we get
// here, but that should have been taken care
// of by the while loop above.
throw("Programming logic error.") ;
}
}
}
computeSimpleDerivatives();
// If our input data is monotone, we can do monotone cubic
// interpolation, so adjust the derivatives if so.
//
// If input data is not monotone, we should not touch
// If input data is not monotone, we should not touch
// the derivatives, as this code should reduce to a
// standard cubic interpolation algorithm.
if (monotone) {
adjustDerivativesForMonotoneness();
}
strictlyMonotoneCached = true;
monotoneCached = true;
}
@@ -495,7 +495,7 @@ computeInternalFunctionData() const {
//
// The notion of "flat" is determined by the input parameter "epsilon"
// Values whose difference are less than epsilon are regarded as equal.
//
//
// This is implemented to be able to obtain a strictly monotone
// curve from a data set that is strictly monotone except at the
// endpoints.
@@ -505,7 +505,7 @@ computeInternalFunctionData() const {
// (1,3), (2,3), (3,4), (4,5), (5,5), (6,5)
// will become
// (2,3), (3,4), (4,5)
//
//
// Assumes at least 3 datapoints. If less than three, this function is a noop.
void
MonotCubicInterpolator::
@@ -514,26 +514,26 @@ chopFlatEndpoints(const double epsilon) {
if (getSize() < 3) {
return;
}
map<double,double>::iterator xf_iterator;
map<double,double>::iterator xf_next_iterator;
// Clear flags:
strictlyMonotoneCached = false;
monotoneCached = false;
// Chop left end:
xf_iterator = data.begin();
xf_next_iterator = ++(data.begin());
// Erase data points that are similar to its right value from the left end.
while ((xf_next_iterator != data.end()) &&
while ((xf_next_iterator != data.end()) &&
(fabs(xf_iterator->second - xf_next_iterator->second) < epsilon )) {
xf_next_iterator++;
data.erase(xf_iterator);
xf_iterator++;
}
xf_iterator = data.end();
xf_iterator = data.end();
xf_iterator--; // (data.end() points beyond the last element)
xf_next_iterator = xf_iterator;
xf_next_iterator--;
@@ -544,15 +544,15 @@ chopFlatEndpoints(const double epsilon) {
data.erase(xf_iterator);
xf_iterator--;
}
// Finished chopping, so recompute function data:
computeInternalFunctionData();
}
//
//
// If function is monotone, but not strictly monotone,
// this function will remove datapoints from intervals
// this function will remove datapoints from intervals
// with zero derivative so that the curves become
// strictly monotone.
//
@@ -561,44 +561,44 @@ chopFlatEndpoints(const double epsilon) {
// (1,2), (2,3), (3,4), (4,4), (5,5), (6,6)
// will become
// (1,2), (2,3), (3,4), (5,5), (6,6)
//
//
// Assumes at least two datapoints, if one or zero datapoint, this is a noop.
//
//
//
void
MonotCubicInterpolator::
shrinkFlatAreas(const double epsilon) {
if (getSize() < 2) {
return;
}
map<double,double>::iterator xf_iterator;
map<double,double>::iterator xf_next_iterator;
// Nothing to do if we already are strictly monotone
if (isStrictlyMonotone()) {
return;
}
// Refuse to change a curve that is not monotone.
if (!isMonotone()) {
return;
}
// Clear flags, they are not to be trusted after we modify the
// data
strictlyMonotoneCached = false;
monotoneCached = false;
// Iterate through data values, if two data pairs
// have equal values, delete one of the data pair.
// Do not trust the source code on which data point is being
// removed (x-values of equal y-points might be averaged in the future)
xf_iterator = data.begin();
xf_next_iterator = ++(data.begin());
while (xf_next_iterator != data.end()) {
//cout << xf_iterator->first << "," << xf_iterator->second << " " << xf_next_iterator->first << "," << xf_next_iterator->second << "\n";
if (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon ) {
@@ -612,71 +612,71 @@ shrinkFlatAreas(const double epsilon) {
xf_next_iterator++;
}
}
}
void
void
MonotCubicInterpolator::
computeSimpleDerivatives() const {
ddata.clear();
// Do endpoints first:
map<double,double>::const_iterator xf_prev_iterator;
map<double,double>::const_iterator xf_iterator;
map<double,double>::const_iterator xf_next_iterator;
double diff;
// Leftmost interval:
xf_iterator = data.begin();
xf_next_iterator = ++(data.begin());
diff =
diff =
(xf_next_iterator->second - xf_iterator->second) /
(xf_next_iterator->first - xf_iterator->first);
ddata[xf_iterator->first] = diff ;
// Rightmost interval:
xf_iterator = --(--(data.end()));
xf_next_iterator = --(data.end());
diff =
diff =
(xf_next_iterator->second - xf_iterator->second) /
(xf_next_iterator->first - xf_iterator->first);
ddata[xf_next_iterator->first] = diff ;
// If we have more than two intervals, loop over internal points:
if (data.size() > 2) {
map<double,double>::const_iterator intpoint;
for (intpoint = ++data.begin(); intpoint != --data.end(); ++intpoint) {
/*
diff = (f2 - f1)/(x2-x1)/w + (f3-f1)/(x3-x2)/2
diff = (f2 - f1)/(x2-x1)/w + (f3-f1)/(x3-x2)/2
average of the forward and backward difference.
Weights are equal, should we weigh with h_i?
*/
map<double,double>::const_iterator lastpoint = intpoint; --lastpoint;
map<double,double>::const_iterator nextpoint = intpoint; ++nextpoint;
diff = (nextpoint->second - intpoint->second)/
(2*(nextpoint->first - intpoint->first))
(2*(nextpoint->first - intpoint->first))
+
(intpoint->second - lastpoint->second) /
(2*(intpoint->first - lastpoint->first));
ddata[intpoint->first] = diff ;
ddata[intpoint->first] = diff ;
}
}
}
void
void
MonotCubicInterpolator::
adjustDerivativesForMonotoneness() const {
map<double,double>::const_iterator point, dpoint;
/* Loop over all intervals, ie. loop over all points and look
at the interval to the right of the point */
for (point = data.begin(), dpoint = ddata.begin();
@@ -685,8 +685,8 @@ adjustDerivativesForMonotoneness() const {
map<double,double>::const_iterator nextpoint, nextdpoint;
nextpoint = point; ++nextpoint;
nextdpoint = dpoint; ++nextdpoint;
double delta =
double delta =
(nextpoint->second - point->second) /
(nextpoint->first - point->first);
if (fabs(delta) < 1e-14) {
@@ -695,25 +695,25 @@ adjustDerivativesForMonotoneness() const {
} else {
double alpha = ddata[point->first] / delta;
double beta = ddata[nextpoint->first] / delta;
if (! isMonotoneCoeff(alpha, beta)) {
double tau = 3/sqrt(alpha*alpha + beta*beta);
ddata[point->first] = tau*alpha*delta;
ddata[nextpoint->first] = tau*beta*delta;
}
}
}
}
void
void
MonotCubicInterpolator::
scaleData(double factor) {
map<double,double>::iterator it , itd ;
@@ -725,8 +725,8 @@ scaleData(double factor) {
} else {
for (it = data.begin() ; it != data.end() ; ++it ) {
it->second *= factor ;
}
}
}
}
}

218
opm/core/utility/MonotCubicInterpolator.hpp
View File

@@ -10,17 +10,17 @@
/*
MonotCubicInterpolator
Copyright (C) 2006 Statoil ASA
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
@@ -32,21 +32,21 @@ namespace Opm
/**
Class to represent a one-dimensional function f with single-valued
argument x. The function is represented by a table of function
values. Interpolation between table values is cubic and monotonicity
values. Interpolation between table values is cubic and monotonicity
preserving if input values are monotonous.
Outside x_min and x_max, the class will extrapolate using the
constant f(x_min) or f(x_max).
Extra functionality:
- Can return (x_1+x_2)/2 where x_1 and x_2 are such that
abs(f(x_1) - f(x_2)) is maximized. This is used to determine where
- Can return (x_1+x_2)/2 where x_1 and x_2 are such that
abs(f(x_1) - f(x_2)) is maximized. This is used to determine where
one should calculate a new value for increased accuracy in the
current function
Monotonicity preserving cubic interpolation algorithm is taken
from Fritsch and Carlson, "Monotone piecewise cubic interpolation",
SIAM J. Numer. Anal. 17, 238--246, no. 2,
Monotonicity preserving cubic interpolation algorithm is taken
from Fritsch and Carlson, "Monotone piecewise cubic interpolation",
SIAM J. Numer. Anal. 17, 238--246, no. 2,
$Id$
@@ -61,9 +61,9 @@ namespace Opm
class MonotCubicInterpolator {
public:
/**
@param datafilename A datafile with the x values and the corresponding f(x) values
@param datafilename A datafile with the x values and the corresponding f(x) values
Accepts a filename as input and parses this file for
two-column floating point data, interpreting the data as
@@ -71,36 +71,36 @@ class MonotCubicInterpolator {
Ignores all lines not conforming to \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
*/
MonotCubicInterpolator(const std::string & datafilename) throw (const char*)
MonotCubicInterpolator(const std::string & datafilename) throw (const char*)
{
if (!read(datafilename)) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
if (!read(datafilename)) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
} ;
}
/**
@param datafilename A datafile with the x values and the corresponding f(x) values
@param datafilename A datafile with the x values and the corresponding f(x) values
Accepts a filename as input and parses this file for
two-column floating point data, interpreting the data as
representing function values x and f(x).
Ignores all lines not conforming to \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
All commas in the file will be treated as spaces when parsing.
*/
MonotCubicInterpolator(const char* datafilename) throw (const char*)
MonotCubicInterpolator(const char* datafilename) throw (const char*)
{
if (!read(std::string(datafilename))) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
if (!read(std::string(datafilename))) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
} ;
}
/**
/**
@param datafilename data file
@param XColumn x values
@param fColumn f values
@@ -108,14 +108,14 @@ class MonotCubicInterpolator {
Accepts a filename as input, and parses the chosen columns in
that file.
*/
MonotCubicInterpolator(const char* datafilename, int xColumn, int fColumn) throw (const char*)
MonotCubicInterpolator(const char* datafilename, int xColumn, int fColumn) throw (const char*)
{
if (!read(std::string(datafilename),xColumn,fColumn)) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
if (!read(std::string(datafilename),xColumn,fColumn)) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
} ;
}
/**
/**
@param datafilename data file
@param XColumn x values
@param fColumn f values
@@ -123,13 +123,13 @@ class MonotCubicInterpolator {
Accepts a filename as input, and parses the chosen columns in
that file.
*/
MonotCubicInterpolator(const std::string & datafilename, int xColumn, int fColumn) throw (const char*)
MonotCubicInterpolator(const std::string & datafilename, int xColumn, int fColumn) throw (const char*)
{
if (!read(datafilename,xColumn,fColumn)) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
if (!read(datafilename,xColumn,fColumn)) {
throw("Unable to constuct MonotCubicInterpolator from file.") ;
} ;
}
/**
@param x vector of x values
@param f vector of corresponding f values
@@ -138,12 +138,12 @@ class MonotCubicInterpolator {
the interpolation object. First vector is the x-values, the
second vector is the function values
*/
MonotCubicInterpolator(const std::vector<double> & x ,
MonotCubicInterpolator(const std::vector<double> & x ,
const std::vector<double> & f);
/**
No input, an empty function object is created.
This object must be treated with care until
populated.
*/
@@ -152,7 +152,7 @@ class MonotCubicInterpolator {
/**
@param datafilename A datafile with the x values and the corresponding f(x) values
@param datafilename A datafile with the x values and the corresponding f(x) values
Accepts a filename as input and parses this file for
two-column floating point data, interpreting the data as
@@ -167,8 +167,8 @@ class MonotCubicInterpolator {
bool read(const std::string & datafilename) {
return read(datafilename,1,2) ;
}
/**
/**
@param datafilename data file
@param XColumn x values
@param fColumn f values
@@ -176,60 +176,60 @@ class MonotCubicInterpolator {
Accepts a filename as input, and parses the chosen columns in
that file.
*/
bool read(const std::string & datafilename, int xColumn, int fColumn) ;
bool read(const std::string & datafilename, int xColumn, int fColumn) ;
/**
@param x x value
Returns f(x) for given x (input). Interpolates (monotone cubic
Returns f(x) for given x (input). Interpolates (monotone cubic
or linearly) if necessary.
Extrapolates using the constants f(x_min) or f(x_max) if
input x is outside (x_min, x_max)
@return f(x) for a given x
*/
double operator () (double x) const { return evaluate(x) ; }
/**
@param x x value
Returns f(x) for given x (input). Interpolates (monotone cubic
Returns f(x) for given x (input). Interpolates (monotone cubic
or linearly) if necessary.
Extrapolates using the constants f(x_min) or f(x_max) if
input x is outside (x_min, x_max)
@return f(x) for a given x
*/
double evaluate(double x) const throw(const char*);
/**
@param x x value
@param errorestimate_output
@param errorestimate_output
Returns f(x) and an error estimate for given x (input).
Interpolates (linearly) if necessary.
Throws an exception if extrapolation would be necessary for
evaluation. We do not want to do extrapolation (yet).
The error estimate for x1 < x < x2 is
(x2 - x1)^2/8 * f''(x) where f''(x) is evaluated using
the stencil (1 -2 1) using either (x0, x1, x2) or (x1, x2, x3);
Throws an exception if the table contains only two x-values.
NOT IMPLEMENTED YET!
*/
double evaluate(double x, double & errorestimate_output ) const ;
double evaluate(double x, double & errorestimate_output ) const ;
/**
Minimum x-value, returns both x and f in a pair.
@return minimum x value
@return f(minimum x value)
*/
@@ -237,10 +237,10 @@ class MonotCubicInterpolator {
// Easy since the data is sorted on x:
return *data.begin();
}
/**
Maximum x-value, returns both x and f in a pair.
@return maximum x value
@return f(maximum x value)
*/
@@ -248,10 +248,10 @@ class MonotCubicInterpolator {
// Easy since the data is sorted on x:
return *data.rbegin();
}
/**
Maximum f-value, returns both x and f in a pair.
@return x value corresponding to maximum f value
@return maximum f value
*/
@@ -259,39 +259,39 @@ class MonotCubicInterpolator {
/**
Minimum f-value, returns both x and f in a pair
@return x value corresponding to minimal f value
@return minimum f value
*/
std::pair<double,double> getMinimumF() const throw(const char*) ;
/**
Provide a copy of the x-data as a vector
Unspecified order, but corresponds to get_fVector.
@return x values as a vector
*/
std::vector<double> get_xVector() const ;
/**
Provide a copy of tghe function data as a vector
Unspecified order, but corresponds to get_xVector
@return f values as a vector
*/
std::vector<double> get_fVector() const ;
/**
@param factor Scaling constant
Scale all the function value data by a constant
*/
void scaleData(double factor);
/**
Determines if the current function-value-data is strictly
monotone. This is a utility function for outsiders if they want
@@ -300,7 +300,7 @@ class MonotCubicInterpolator {
@return True if f(x) is strictly monotone, else False
*/
bool isStrictlyMonotone() {
/* Use cached value if it can be trusted */
if (strictlyMonotoneCached) {
return strictlyMonotone;
@@ -313,7 +313,7 @@ class MonotCubicInterpolator {
/**
Determines if the current function-value-data is monotone.
@return True if f(x) is monotone, else False
*/
bool isMonotone() const {
@@ -333,7 +333,7 @@ class MonotCubicInterpolator {
@return True if f(x) is strictly increasing, else False
*/
bool isStrictlyIncreasing() {
/* Use cached value if it can be trusted */
if (strictlyMonotoneCached) {
return (strictlyMonotone && strictlyIncreasing);
@@ -346,7 +346,7 @@ class MonotCubicInterpolator {
/**
Determines if the current function-value-data is monotone and increasing.
@return True if f(x) is monotone and increasing, else False
*/
bool isMonotoneIncreasing() const {
@@ -366,7 +366,7 @@ class MonotCubicInterpolator {
@return True if f(x) is strictly decreasing, else False
*/
bool isStrictlyDecreasing() {
/* Use cached value if it can be trusted */
if (strictlyMonotoneCached) {
return (strictlyMonotone && strictlyDecreasing);
@@ -379,7 +379,7 @@ class MonotCubicInterpolator {
/**
Determines if the current function-value-data is monotone and decreasing
@return True if f(x) is monotone and decreasing, else False
*/
bool isMonotoneDecreasing() const {
@@ -391,9 +391,9 @@ class MonotCubicInterpolator {
return (monotone && decreasing);
}
}
/**
@param newx New x point
@param newf New f(x) point
@@ -408,12 +408,12 @@ class MonotCubicInterpolator {
*/
void addPair(double newx, double newf) throw(const char*);
/**
Returns an x-value that is believed to yield the best
improvement in global accuracy for the interpolation if
computed.
Searches for the largest jump in f-values, and returns a x
value being the average of the two x-values representing the
f-value-jump.
@@ -422,14 +422,14 @@ class MonotCubicInterpolator {
@return Maximal difference
*/
std::pair<double,double> getMissingX() const throw(const char*) ;
/**
Constructs a string containing the data in a table
Constructs a string containing the data in a table
@return a string containing the data in a table
*/
std::string toString() const;
/**
@return Number of datapoint pairs in this object
*/
@@ -440,7 +440,7 @@ class MonotCubicInterpolator {
/**
Checks if the function curve is flat at the endpoints, chop off
endpoint data points if that is the case.
The notion of "flat" is determined by the input parameter "epsilon"
Values whose difference are less than epsilon are regarded as equal.
@@ -460,7 +460,7 @@ class MonotCubicInterpolator {
/**
Wrapper function for chopFlatEndpoints(const double)
providing a default epsilon parameter
providing a default epsilon parameter
*/
void chopFlatEndpoints() {
chopFlatEndpoints(1e-14);
@@ -468,7 +468,7 @@ class MonotCubicInterpolator {
/**
If function is monotone, but not strictly monotone,
this function will remove datapoints from intervals
this function will remove datapoints from intervals
with zero derivative so that the curve become
strictly monotone.
@@ -477,14 +477,14 @@ class MonotCubicInterpolator {
(1,2), (2,3), (3,4), (4,4), (5,5), (6,6)
will become
(1,2), (2,3), (3,4), (5,5), (6,6)
Assumes at least two datapoints, if one or zero datapoint, this is a noop.
*/
void shrinkFlatAreas(const double);
/**
Wrapper function for shrinkFlatAreas(const double)
providing a default epsilon parameter
providing a default epsilon parameter
*/
void shrinkFlatAreas() {
shrinkFlatAreas(1e-14);
@@ -493,17 +493,17 @@ class MonotCubicInterpolator {
private:
// Data structure to store x- and f-values
std::map<double, double> data;
// Data structure to store x- and d-values
mutable std::map<double, double> ddata;
mutable std::map<double, double> ddata;
// Storage containers for precomputed interpolation data
// std::vector<double> dvalues; // derivatives in Hermite interpolation.
// Flag to determine whether the boolean strictlyMonotone can be
// trusted.
mutable bool strictlyMonotoneCached;
@@ -517,13 +517,13 @@ private:
mutable bool strictlyIncreasing;
mutable bool decreasing;
mutable bool increasing;
/* Hermite basis functions, t \in [0,1] ,
notation from:
http://en.wikipedia.org/w/index.php?title=Cubic_Hermite_spline&oldid=84495502
*/
double H00(double t) const {
return 2*t*t*t - 3*t*t + 1;
}
@@ -536,35 +536,35 @@ private:
double H11(double t) const {
return t*t*t - t*t;
}
void computeInternalFunctionData() const ;
/**
/**
Computes initial derivative values using centered (second order) difference
for internal datapoints, and one-sided derivative for endpoints
The internal datastructure map<double,double> ddata is populated by this method.
*/
void computeSimpleDerivatives() const ;
/**
Adjusts the derivative values (ddata) so that we can guarantee that
the resulting piecewise Hermite polymial is monotone. This is
done according to the algorithm of Fritsch and Carlsson 1980,
done according to the algorithm of Fritsch and Carlsson 1980,
see Section 4, especially the two last lines.
*/
void adjustDerivativesForMonotoneness() const ;
/**
Checks if the coefficient alpha and beta is in
the region that guarantees monotoneness of the
derivative values they represent
See Fritsch and Carlson 1980, Lemma 2,
alternatively Step 5 in Wikipedia's article
/**
Checks if the coefficient alpha and beta is in
the region that guarantees monotoneness of the
derivative values they represent
See Fritsch and Carlson 1980, Lemma 2,
alternatively Step 5 in Wikipedia's article
on Monotone cubic interpolation.
*/
bool isMonotoneCoeff(double alpha, double beta) const {
@@ -574,8 +574,8 @@ private:
return false;
}
}
};
};
} // namespace Opm

2
opm/core/utility/linearInterpolation.hpp
View File

@@ -69,7 +69,7 @@ namespace Opm
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv,
double x, int& ix1)