ECLPiecewiseLinearInterpolant.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_ECLSIMPLE1DINTERPOLANT_HEADER_INCLUDED
#define OPM_ECLSIMPLE1DINTERPOLANT_HEADER_INCLUDED
 
#include <opm/utility/ECLTableInterpolation1D.hpp>
 
#include <cassert>
#include <cmath>
#include <exception>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace Opm { namespace Interp1D { namespace PiecewisePolynomial {
 
    template <class Extrapolation, bool IsAscendingRange = true>
    class Linear
    {
    public:
        explicit Linear(Extrapolation&& extrap)
            : extrap_(std::forward<Extrapolation>(extrap))
            , nCols_ (0)
        {}
 
        template <class ElmIterator, class ValueTransform>
        Linear(Extrapolation&&                    extrap,
               ElmIterator                        xBegin,
               ElmIterator                        xEnd,
               std::vector<ElmIterator>&          colIt,
               const ValueTransform&              xTransform,
               const std::vector<ValueTransform>& colTransform);
 
        LocalInterpPoint classifyPoint(const double x) const
        {
            return LocalInterpPoint::identify(this->x_, x, Ascending_P{});
        }
 
        double evaluate(const std::size_t       col,
                        const LocalInterpPoint& pt) const
        {
            if (pt.cat == PointCategory::InRange) {
                // Common case.  Input point is within range of x_.  Placed
                // first to enable early return.
                return this->interpolate(pt, col);
            }
 
            auto yval = [this](const std::size_t i,
                               const std::size_t j) -> double
            {
                return this->y(i, j);
            };
 
            if (pt.cat == PointCategory::LeftOfRange) {
                // Extrapolate function ot the left of input range.
                const auto xmin = this->x_.front();
 
                return this->extrap_.left(this->x_, xmin + pt.t, col, yval);
            }
 
            assert (pt.cat == PointCategory::RightOfRange);
 
            // Extrapolate function ot the right of input range.
            const auto xmax = this->x_.back();
            return this->extrap_.right(this->x_, xmax + pt.t, col, yval);
        }
 
        const std::vector<double>& independentVariable() const
        {
            return this->x_;
        }
 
        std::vector<double> resultVariable(const std::size_t col) const
        {
            auto result = std::vector<double>{};
 
            if (col >= this->nCols_) {
                throw std::domain_error {
                    "Result Column Identifier Ouf of Bounds"
                };
            }
 
            result.reserve(this->x_.size());
            for (auto n = this->x_.size(), i = 0*n; i < n; ++i) {
                result.push_back(this->y(i, col));
            }
 
            return result;
        }
 
    private:
        using Ascending_P =
            std::integral_constant<bool, IsAscendingRange>;
 
        Extrapolation extrap_;
 
        std::size_t nCols_;
 
        std::vector<double> x_;
 
        std::vector<double> y_;
 
        double interpolate(const LocalInterpPoint& pt,
                           const std::size_t       col) const
        {
            assert (pt.cat == PointCategory::InRange);
            assert (pt.interval + 1 < this->x_.size());
 
            const auto left = pt.interval + 0;
            const auto xl   = this->x_[left];
            const auto yl   = this->y (left, col);
 
            const auto right = pt.interval + 1;
            const auto xr    = this->x_[right];
            const auto yr    = this->y (right, col);
 
            const auto t = pt.t / (xr - xl);
 
            return t*yr + (1.0 - t)*yl;
        }
 
        double y(const std::size_t row,
                 const std::size_t col) const
        {
            // Recall: this->y_ stored with column index cycling the most
            // rapidly.
 
            assert (col < this->nCols_);
            assert (row*this->nCols_ + col < this->y_.size());
 
            return this->y_[row*this->nCols_ + col];
        }
    };
 
    template <class Extrapolation, bool IsAscendingRange>
    template <class ElmIterator, class ValueTransform>
    Linear<Extrapolation, IsAscendingRange>::
    Linear(Extrapolation&&                    extrap,
           ElmIterator                        xBegin,
           ElmIterator                        xEnd,
           std::vector<ElmIterator>&          colIt,
           const ValueTransform&              xTransform,
           const std::vector<ValueTransform>& colTransform)
        : extrap_(std::forward<Extrapolation>(extrap))
        , nCols_ (colIt.size())
    {
        // There must be at least one dependent variable/result variable.
        assert (colIt.size() >= 1);
 
        const auto nRows = std::distance(xBegin, xEnd);
 
        this->x_.reserve(nRows);
        this->y_.reserve(nRows * colIt.size());
 
        auto keyValid = [](const double xi)
        {
            // Indep. variable values <= -1.0e20 or >= 1.0e20 signal
            // "unused" table nodes (rows).  These nodes are in the table to
            // fill out the allocated size if one particular sub-table does
            // not use all nodes.  The magic value 1.0e20 is documented in
            // the Fileformats Reference Manual.
            return std::abs(xi) < 1.0e20;
        };
 
        while (xBegin != xEnd) {
            // Extract relevant portion of the table.  Preallocated rows
            // that are not actually part of the result set (i.e., those
            // that are set to a sentinel value) are discarded.
            if (keyValid(*xBegin)) {
                this->x_.push_back(xTransform(*xBegin));
 
                auto colID = 0*colTransform.size();
                for (auto ci : colIt) {
                    // Store 'y_' with column index cycling most rapidly.
                    this->y_.push_back(colTransform[colID++](*ci));
                }
            }
 
            // -------------------------------------------------------------
            // Advance iterators.
 
            // 1) Independent variable.
            ++xBegin;
 
            // 2) Dependent/result/columns.
            for (auto& ci : colIt) {
                ++ci;
            }
        }
 
        // Dispose of any excess capacity.
        if (this->x_.size() < static_cast<decltype(this->x_.size())>(nRows)) {
            this->x_.shrink_to_fit();
            this->y_.shrink_to_fit();
        }
 
        if (this->x_.size() < 2) {
            // Table has no interval that supports interpolation.  Either
            // just a single node or no nodes at all.  We can't do anything
            // useful here, so don't pretend that this is okay.
 
            throw std::invalid_argument {
                "No Interpolation Intervals of Non-Zero Size"
            };
        }
    }
 
}}} // Opm::Interp1D::PiecewisePolynomial
 
#endif // OPM_ECLSIMPLE1DINTERPOLANT_HEADER_INCLUDED