diff --git a/dune/porsol/common/UniformTableLinear.hpp b/dune/porsol/common/UniformTableLinear.hpp new file mode 100644 index 000000000..cf3bef256 --- /dev/null +++ b/dune/porsol/common/UniformTableLinear.hpp @@ -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 . +*/ + +#ifndef OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED +#define OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED + +#include +#include +#include +#include +#include + +#include +#include + +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 + 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& 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 domain(); + + /// @brief Rescale the domain. + /// @param new_domain the new domain as a pair of doubles. + void rescaleDomain(std::pair 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 y_values_; + RangePolicy left_; + RangePolicy right_; + template + friend std::ostream& operator<<(std::ostream& os, const UniformTableLinear& t); + }; + + + // Member implementations. + + template + inline + UniformTableLinear + ::UniformTableLinear() + : left_(ClosestValue), right_(ClosestValue) + { + } + + template + inline + UniformTableLinear + ::UniformTableLinear(double xmin, + double xmax, + const std::vector& 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 + inline + UniformTableLinear + ::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 + inline std::pair + UniformTableLinear + ::domain() + { + return std::make_pair(xmin_, xmax_); + } + + template + inline void + UniformTableLinear + ::rescaleDomain(std::pair new_domain) + { + xmin_ = new_domain.first; + xmax_ = new_domain.second; + xdelta_ = (xmax_ - xmin_)/(y_values_.size() - 1); + } + + template + inline double + UniformTableLinear + ::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 + inline double + UniformTableLinear + ::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 + inline bool + UniformTableLinear + ::operator==(const UniformTableLinear& other) const + { + return xmin_ == other.xmin_ + && xdelta_ == other.xdelta_ + && y_values_ == other.y_values_ + && left_ == other.left_ + && right_ == other.right_; + } + + template + inline void + UniformTableLinear + ::setLeftPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + left_ = rp; + } + + template + inline void + UniformTableLinear + ::setRightPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + right_ = rp; + } + + + template + inline std::ostream& operator<<(std::ostream& os, const UniformTableLinear& 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 diff --git a/dune/porsol/common/buildUniformMonotoneTable.hpp b/dune/porsol/common/buildUniformMonotoneTable.hpp new file mode 100644 index 000000000..03c060bec --- /dev/null +++ b/dune/porsol/common/buildUniformMonotoneTable.hpp @@ -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 . +*/ + +#ifndef OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED +#define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED + +#include +#include + +namespace Dune { + namespace utils { + + template + void buildUniformMonotoneTable(const std::vector& xv, + const std::vector& yv, + const int samples, + UniformTableLinear& table) + { + MonotCubicInterpolator interp(xv, yv); + std::vector 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(xmin, xmax, uniform_yv); + } + + } // namespace utils +} // namespace Dune + + + +#endif // OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED