ECLTableInterpolation1D.hpp

Go to the documentation of this file.

/*
  Copyright 2017 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_ECLTABLEINTERPOLATION1D_HEADER_INCLUDED
#define OPM_ECLTABLEINTERPOLATION1D_HEADER_INCLUDED
 
#include <cassert>
#include <type_traits>
#include <vector>
 
 
namespace Opm { namespace Interp1D {
 
    enum class PointCategory {
        InRange,
 
        LeftOfRange,
 
        RightOfRange,
    };
 
    namespace PiecewisePolynomial {
 
        namespace ExtrapolationPolicy {
            class Constant
            {
            public:
                template <class TabulatedFunction, class Index>
                double left(const std::vector<double>& /* xi */,
                            const double               /* x */,
                            const Index                col ,
                            TabulatedFunction&&        f) const
                {
                    return f(0, col);
                }
 
                template <class TabulatedFunction, class Index>
                double right(const std::vector<double>& xi,
                             const double               /* x */,
                             const Index                col,
                             TabulatedFunction&&        f) const
                {
                    return f(xi.size() - 1, col);
                }
            };
 
            class Linearly
            {
            public:
                template <class TabulatedFunction, class Index>
                double left(const std::vector<double>& xi,
                            const double               x,
                            const Index                col,
                            TabulatedFunction&&        f) const
                {
                    // Derivative of f(i,col)
                    const auto f0   = f(0, col);
                    const auto f1   = f(1, col);
                    const auto dfdx = (f1 - f0) / (xi[1] - xi[0]);
 
                    // <= 0 if ascending range.
                    const auto dx   = x - xi.front();
 
                    return f0 + (dfdx * dx);
                }
 
                template <class TabulatedFunction, class Index>
                double right(const std::vector<double>& xi,
                             const double               x,
                             const Index                col,
                             TabulatedFunction&&        f) const
                {
                    const auto nRows = xi.size();
 
                    // Derivative of f(i,col)
                    const auto f0   = f(nRows - 2, col);
                    const auto f1   = f(nRows - 1, col);
                    const auto dfdx =
                        (f1 - f0) / (xi[nRows - 1] - xi[nRows - 2]);
 
                    // >= 0 if ascending range
                    const auto dx = x - xi.back();
 
                    return f1 + (dfdx * dx);
                }
            };
 
            class LinearlyWithDerivatives
            {
            public:
                explicit LinearlyWithDerivatives(const std::size_t nResCol)
                    : nResCol_(nResCol)
                {}
 
                template <class TabulatedFunction, class Index>
                double left(const std::vector<double>& xi,
                            const double               x,
                            const Index                col,
                            TabulatedFunction&&        f) const
                {
                    // Derivative of f(i,col) in f(i, nResCol_ + col)
                    const auto dfdx = f(0, this->nResCol_ + col);
                    const auto dx   = x - xi.front(); // <= 0 if ascending range
 
                    return f(0, col) + (dfdx * dx);
                }
 
                template <class TabulatedFunction, class Index>
                double right(const std::vector<double>& xi,
                             const double               x,
                             const Index                col,
                             TabulatedFunction&&        f) const
                {
                    const auto nRows = xi.size();
 
                    // Derivative of f(i,col) in f(i, nResCol_ + col)
                    const auto dfdx = f(nRows - 1, this->nResCol_ + col);
                    const auto dx   = x - xi.back(); // >= 0 if ascending range
 
                    return f(nRows - 1, col) + (dfdx * dx);
                }
 
            private:
                std::size_t nResCol_;
            };
        } // ExtrapolationPolicy
 
        struct LocalInterpPoint {
            PointCategory cat;
 
            std::vector<double>::size_type interval;
 
            double t;
 
            static LocalInterpPoint
            identify(const std::vector<double>& xi,
                     const double               x,
                     std::true_type is_ascending = std::true_type{});
 
            static LocalInterpPoint
            identify(const std::vector<double>& xi,
                     const double               x,
                     std::false_type is_ascending);
        };
 
    } // PiecewisePolynomial
 
}} // Opm::Interp1D
 
#endif // OPM_ECLTABLEINTERPOLATION1D_HEADER_INCLUDED