Merged in UniformTableLinear.hpp and friends.
This commit is contained in:
Atgeirr Flø Rasmussen
commit
3914dcb8e9
1
.hgtags
Normal file
1
.hgtags
Normal file
@ -0,0 +1 @@
|
||||
b3cc85c2b5ab3f5b283244624d06a1afab9f8d67 Version-InitialImport-svn935
|
||||
260
dune/porsol/common/UniformTableLinear.hpp
Normal file
260
dune/porsol/common/UniformTableLinear.hpp
Normal file
@ -0,0 +1,260 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
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_UNIFORMTABLELINEAR_HEADER_INCLUDED
|
||||
#define OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED
|
||||
|
||||
#include <cmath>
|
||||
#include <exception>
|
||||
#include <vector>
|
||||
#include <utility>
|
||||
#include <iostream>
|
||||
|
||||
#include <dune/common/ErrorMacros.hpp>
|
||||
#include <dune/porsol/common/linearInterpolation.hpp>
|
||||
|
||||
namespace Dune {
|
||||
namespace utils {
|
||||
|
||||
|
||||
/// @brief This class uses linear interpolation to compute the value
|
||||
/// (and its derivative) of a function f sampled at uniform points.
|
||||
/// @tparam T the range type of the function (should be an algebraic ring type)
|
||||
template<typename T>
|
||||
class UniformTableLinear
|
||||
{
|
||||
public:
|
||||
/// @brief Default constructor.
|
||||
UniformTableLinear();
|
||||
|
||||
/// @brief Construct from vector of y-values.
|
||||
/// @param xmin the x value corresponding to the first y value.
|
||||
/// @param xmax the x value corresponding to the last y value.
|
||||
/// @param y_values vector of range values.
|
||||
UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const std::vector<T>& y_values);
|
||||
|
||||
/// @brief Construct from array of y-values.
|
||||
/// @param xmin the x value corresponding to the first y value.
|
||||
/// @param xmax the x value corresponding to the last y value.
|
||||
/// @param y_values array of range values.
|
||||
/// @param num_y_values the number of values in y_values.
|
||||
UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const T* y_values,
|
||||
int num_y_values);
|
||||
|
||||
/// @brief Get the domain.
|
||||
/// @return the domain as a pair of doubles.
|
||||
std::pair<double, double> domain();
|
||||
|
||||
/// @brief Rescale the domain.
|
||||
/// @param new_domain the new domain as a pair of doubles.
|
||||
void rescaleDomain(std::pair<double, double> new_domain);
|
||||
|
||||
/// @brief Evaluate the value at x.
|
||||
/// @param x a domain value
|
||||
/// @return f(x)
|
||||
double operator()(const double x) const;
|
||||
|
||||
/// @brief Evaluate the derivative at x.
|
||||
/// @param x a domain value
|
||||
/// @return f'(x)
|
||||
double derivative(const double x) const;
|
||||
|
||||
/// @brief Equality operator.
|
||||
/// @param other another UniformTableLinear.
|
||||
/// @return true if they are represented exactly alike.
|
||||
bool operator==(const UniformTableLinear& other) const;
|
||||
|
||||
/// @brief Policies for how to behave when trying to evaluate outside the domain.
|
||||
enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2};
|
||||
|
||||
/// @brief Sets the behavioural policy for evaluation to the left of the domain.
|
||||
/// @param rp the policy
|
||||
void setLeftPolicy(RangePolicy rp);
|
||||
|
||||
/// @brief Sets the behavioural policy for evaluation to the right of the domain.
|
||||
/// @param rp the policy
|
||||
void setRightPolicy(RangePolicy rp);
|
||||
|
||||
protected:
|
||||
double xmin_;
|
||||
double xmax_;
|
||||
double xdelta_;
|
||||
std::vector<T> y_values_;
|
||||
RangePolicy left_;
|
||||
RangePolicy right_;
|
||||
template <typename U>
|
||||
friend std::ostream& operator<<(std::ostream& os, const UniformTableLinear<U>& t);
|
||||
};
|
||||
|
||||
|
||||
// Member implementations.
|
||||
|
||||
template<typename T>
|
||||
inline
|
||||
UniformTableLinear<T>
|
||||
::UniformTableLinear()
|
||||
: left_(ClosestValue), right_(ClosestValue)
|
||||
{
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline
|
||||
UniformTableLinear<T>
|
||||
::UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const std::vector<T>& y_values)
|
||||
: xmin_(xmin), xmax_(xmax), y_values_(y_values),
|
||||
left_(ClosestValue), right_(ClosestValue)
|
||||
{
|
||||
ASSERT(xmax > xmin);
|
||||
ASSERT(y_values.size() > 1);
|
||||
xdelta_ = (xmax - xmin)/(y_values.size() - 1);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline
|
||||
UniformTableLinear<T>
|
||||
::UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const T* y_values,
|
||||
int num_y_values)
|
||||
: xmin_(xmin), xmax_(xmax),
|
||||
y_values_(y_values, y_values + num_y_values),
|
||||
left_(ClosestValue), right_(ClosestValue)
|
||||
{
|
||||
ASSERT(xmax > xmin);
|
||||
ASSERT(y_values_.size() > 1);
|
||||
xdelta_ = (xmax - xmin)/(y_values_.size() - 1);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline std::pair<double, double>
|
||||
UniformTableLinear<T>
|
||||
::domain()
|
||||
{
|
||||
return std::make_pair(xmin_, xmax_);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline void
|
||||
UniformTableLinear<T>
|
||||
::rescaleDomain(std::pair<double, double> new_domain)
|
||||
{
|
||||
xmin_ = new_domain.first;
|
||||
xmax_ = new_domain.second;
|
||||
xdelta_ = (xmax_ - xmin_)/(y_values_.size() - 1);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline double
|
||||
UniformTableLinear<T>
|
||||
::operator()(const double xparam) const
|
||||
{
|
||||
// Implements ClosestValue policy.
|
||||
double x = std::min(xparam, xmax_);
|
||||
x = std::max(x, xmin_);
|
||||
|
||||
// Lookup is easy since we are uniform in x.
|
||||
double pos = (x - xmin_)/xdelta_;
|
||||
double posi = std::floor(pos);
|
||||
int left = int(posi);
|
||||
if (left == int(y_values_.size()) - 1) {
|
||||
// We are at xmax_
|
||||
return y_values_.back();
|
||||
}
|
||||
double w = pos - posi;
|
||||
return (1.0 - w)*y_values_[left] + w*y_values_[left + 1];
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline double
|
||||
UniformTableLinear<T>
|
||||
::derivative(const double xparam) const
|
||||
{
|
||||
// Implements ClosestValue policy.
|
||||
double x = std::min(xparam, xmax_);
|
||||
x = std::max(x, xmin_);
|
||||
|
||||
// Lookup is easy since we are uniform in x.
|
||||
double pos = (x - xmin_)/xdelta_;
|
||||
double posi = std::floor(pos);
|
||||
int left = int(posi);
|
||||
if (left == int(y_values_.size()) - 1) {
|
||||
// We are at xmax_
|
||||
--left;
|
||||
}
|
||||
return (y_values_[left + 1] - y_values_[left])/xdelta_;
|
||||
}
|
||||
|
||||
|
||||
template<typename T>
|
||||
inline bool
|
||||
UniformTableLinear<T>
|
||||
::operator==(const UniformTableLinear<T>& other) const
|
||||
{
|
||||
return xmin_ == other.xmin_
|
||||
&& xdelta_ == other.xdelta_
|
||||
&& y_values_ == other.y_values_
|
||||
&& left_ == other.left_
|
||||
&& right_ == other.right_;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline void
|
||||
UniformTableLinear<T>
|
||||
::setLeftPolicy(RangePolicy rp)
|
||||
{
|
||||
if (rp != ClosestValue) {
|
||||
THROW("Only ClosestValue RangePolicy implemented.");
|
||||
}
|
||||
left_ = rp;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline void
|
||||
UniformTableLinear<T>
|
||||
::setRightPolicy(RangePolicy rp)
|
||||
{
|
||||
if (rp != ClosestValue) {
|
||||
THROW("Only ClosestValue RangePolicy implemented.");
|
||||
}
|
||||
right_ = rp;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
inline std::ostream& operator<<(std::ostream& os, const UniformTableLinear<T>& t)
|
||||
{
|
||||
int n = t.y_values_.size();
|
||||
for (int i = 0; i < n; ++i) {
|
||||
double f = double(i)/double(n - 1);
|
||||
os << (1.0 - f)*t.xmin_ + f*t.xmax_
|
||||
<< " " << t.y_values_[i] << '\n';
|
||||
}
|
||||
return os;
|
||||
}
|
||||
|
||||
} // namespace utils
|
||||
} // namespace Dune
|
||||
|
||||
#endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED
|
||||
52
dune/porsol/common/buildUniformMonotoneTable.hpp
Normal file
52
dune/porsol/common/buildUniformMonotoneTable.hpp
Normal file
@ -0,0 +1,52 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
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_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED
|
||||
#define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED
|
||||
|
||||
#include <dune/common/MonotCubicInterpolator.hpp>
|
||||
#include <dune/porsol/common/UniformTableLinear.hpp>
|
||||
|
||||
namespace Dune {
|
||||
namespace utils {
|
||||
|
||||
template <typename T>
|
||||
void buildUniformMonotoneTable(const std::vector<double>& xv,
|
||||
const std::vector<T>& yv,
|
||||
const int samples,
|
||||
UniformTableLinear<T>& table)
|
||||
{
|
||||
MonotCubicInterpolator interp(xv, yv);
|
||||
std::vector<T> uniform_yv(samples);
|
||||
double xmin = xv[0];
|
||||
double xmax = xv.back();
|
||||
for (int i = 0; i < samples; ++i) {
|
||||
double w = double(i)/double(samples - 1);
|
||||
double x = (1.0 - w)*xmin + w*xmax;
|
||||
uniform_yv[i] = interp(x);
|
||||
}
|
||||
table = UniformTableLinear<T>(xmin, xmax, uniform_yv);
|
||||
}
|
||||
|
||||
} // namespace utils
|
||||
} // namespace Dune
|
||||
|
||||
|
||||
|
||||
#endif // OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED
|
||||
90
dune/porsol/common/linearInterpolation.hpp
Normal file
90
dune/porsol/common/linearInterpolation.hpp
Normal file
@ -0,0 +1,90 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// 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 Statoil ASA.
|
||||
|
||||
This file is part of The Open Reservoir Simulator Project (OpenRS).
|
||||
|
||||
OpenRS 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.
|
||||
|
||||
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPENRS_LINEARINTERPOLATION_HEADER
|
||||
#define OPENRS_LINEARINTERPOLATION_HEADER
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
|
||||
namespace Dune
|
||||
{
|
||||
|
||||
/** 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)
|
||||
{
|
||||
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];
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief
|
||||
/// @todo Doc me!
|
||||
/// @tparam
|
||||
/// @param
|
||||
/// @return
|
||||
template <typename T>
|
||||
T linearInterpolationDerivative(const std::vector<double>& xv,
|
||||
const std::vector<T>& yv,
|
||||
double x)
|
||||
{
|
||||
double epsilon = 1e-4; // @@ Ad hoc, should choose based on xv.
|
||||
double x_low = std::max(xv[0], x - epsilon);
|
||||
double x_high = std::min(xv.back(), x + epsilon);
|
||||
T low = linearInterpolation(xv, yv, x_low);
|
||||
T high = linearInterpolation(xv, yv, x_high);
|
||||
return (high - low)/(x_high - x_low);
|
||||
}
|
||||
|
||||
|
||||
} // namespace Dune
|
||||
|
||||
|
||||
|
||||
#endif // OPENRS_LINEARINTERPOLATION_HEADER
|
||||
Loading…
Reference in New Issue
Block a user