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

Finished unification of linear interpolation.

Browse Source 
The functions of linInt.hpp are now used everywhere, but:
 - linInt.hpp -> linearInterpolation.hpp (better name)
 - linearInterpolationExtrap() -> linearInterpolation() (extrapolate by default)
This commit is contained in:
Atgeirr Flø Rasmussen
2013-03-22 15:33:07 +01:00
parent a8097317a5
commit c1657b427a
6 changed files with 53 additions and 196 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
opm/core/io/eclipse/EclipseGridParserHelpers.hpp
View File

@@ -41,7 +41,7 @@
#include <istream>
#include <vector>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <boost/date_time/gregorian/gregorian.hpp>
namespace Opm

2
opm/core/props/pvt/BlackoilPvtProperties.cpp
View File

@@ -29,7 +29,7 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
namespace Opm

2
opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@@ -30,7 +30,7 @@
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>

2
opm/core/props/pvt/SinglePvtLiveOil.cpp
View File

@@ -19,7 +19,7 @@
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>

106
opm/core/utility/linInt.hpp
View File

@@ -1,106 +0,0 @@
//===========================================================================
//
// File: linInt.hpp
//
// Created: Tue Feb 16 14:44:10 2010
//
// Author(s): Bj<42>rn Spjelkavik <bsp@sintef.no>
// Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM 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 3 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_LININT_HEADER
#define OPM_LININT_HEADER
#include <vector>
#include <algorithm>
namespace Opm
{
inline int tableIndex(const std::vector<double>& table, double x)
{
// Returns an index in an ordered table such that x is between
// table[j] and table[j+1]. If x is out of range, first or last
// interval is returned; Binary search.
int n = table.size() - 1;
if (n < 2) {
return 0;
}
int jl = 0;
int ju = n;
bool ascend = (table[n] > table[0]);
while (ju - jl > 1) {
int jm = (ju + jl)/2; // Compute a midpoint
if ( (x >= table[jm]) == ascend) {
jl = jm; // Replace lower limit
} else {
ju = jm; // Replace upper limit
}
}
return jl;
}
inline double linearInterpolationDerivative(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
// Extrapolates if x is outside xv
int ix1 = tableIndex(xv, x);
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1]);
}
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
// Extrapolates if x is outside xv
int ix1 = tableIndex(xv, x);
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)
{
// Extrapolates if x is outside xv
ix1 = tableIndex(xv, x);
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
} // namespace Opm
#endif // OPM_LININT_HEADER

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

@@ -1,19 +1,5 @@
//===========================================================================
//
// File: linearInterpolation.hpp
//
// Created: Tue Sep 9 12:49:39 2008
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2013 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
@@ -32,97 +18,74 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_LINEARINTERPOLATION_HEADER
#define OPM_LINEARINTERPOLATION_HEADER
#ifndef OPM_LINEARINTERPOLATION_HEADER_INCLUDED
#define OPM_LINEARINTERPOLATION_HEADER_INCLUDED
#include <opm/core/utility/linInt.hpp>
#if 0
#include <vector>
#include <algorithm>
namespace Opm
{
/** Linear interpolation.
* Given an increasing vector xv of parameter values and
* a vector yv of point values of the same size,
* the function returns ...
*/
template <typename T>
T linearInterpolation(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
inline int tableIndex(const std::vector<double>& table, double x)
{
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
int lb_ix = lb - xv.begin();
if (lb_ix == 0) {
return yv[0];
} else if (lb_ix == int(xv.size())) {
return yv.back();
} else {
double w = (x - xv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
return (1.0 - w)*yv[lb_ix - 1] + w*yv[lb_ix];
// Returns an index in an ordered table such that x is between
// table[j] and table[j+1]. If x is out of range, first or last
// interval is returned; Binary search.
int n = table.size() - 1;
if (n < 2) {
return 0;
}
int jl = 0;
int ju = n;
bool ascend = (table[n] > table[0]);
while (ju - jl > 1) {
int jm = (ju + jl)/2; // Compute a midpoint
if ( (x >= table[jm]) == ascend) {
jl = jm; // Replace lower limit
} else {
ju = jm; // Replace upper limit
}
}
return jl;
}
/// @brief
/// @todo Doc me!
/// @tparam
/// @param
/// @return
template <typename T>
T linearInterpolationDerivative(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
inline double linearInterpolationDerivative(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
int lb_ix = lb - xv.begin();
if (lb_ix == 0) {
return 0.;
} else if (lb_ix == int(xv.size())) {
return 0.;
} else {
return (yv[lb_ix] - yv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
}
// Extrapolates if x is outside xv
int ix1 = tableIndex(xv, x);
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1]);
}
template <typename T>
T linearInterpolationDerivativeExtrap(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv, double x)
{
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
int lb_ix = lb - xv.begin();
int nend = int(xv.size());
if (lb_ix == 0) {
return (yv[1]-yv[0])/(xv[1]-xv[0]);
} else if (lb_ix == int(xv.size())) {
return (yv[nend-1]-yv[nend-2])/(xv[nend-1]-xv[nend-2]);
} else {
return (yv[lb_ix] - yv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
}
// Extrapolates if x is outside xv
int ix1 = tableIndex(xv, x);
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
template <typename T>
T linearInterpolation(const std::vector<double>& xv,
const std::vector<T>& yv,
double x)
inline double linearInterpolation(const std::vector<double>& xv,
const std::vector<double>& yv,
double x, int& ix1)
{
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
int lb_ix = lb - xv.begin();
if (lb_ix == 0) {
lb_ix=1;
} else if (lb_ix == int(xv.size())){
lb_ix = int(xv.size())-1;
}
double w = (x - xv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
return (1.0 - w)*yv[lb_ix - 1] + w*yv[lb_ix];
// Extrapolates if x is outside xv
ix1 = tableIndex(xv, x);
int ix2 = ix1 + 1;
return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
}
} // namespace Opm
#endif
#endif // OPM_LINEARINTERPOLATION_HEADER
#endif // OPM_LINEARINTERPOLATION_HEADER_INCLUDED