UniformTabulated2DFunction.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright (C) 2013 by Andreas Lauser

  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 2 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_UNIFORM_TABULATED_2D_FUNCTION_HPP
#define OPM_UNIFORM_TABULATED_2D_FUNCTION_HPP

#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/material/common/MathToolbox.hpp>


#include <vector>

#include <assert.h>

namespace Opm {

template <class Scalar>
class UniformTabulated2DFunction
{
public:
    UniformTabulated2DFunction()
    { }

    UniformTabulated2DFunction(Scalar xMin, Scalar xMax, unsigned m,
                               Scalar yMin, Scalar yMax, unsigned n)
    {
        resize(xMin, xMax, m, yMin, yMax, n);
    }

    void resize(Scalar xMin, Scalar xMax, unsigned m,
                Scalar yMin, Scalar yMax, unsigned n)
    {
        samples_.resize(m*n);

        m_ = m;
        n_ = n;

        xMin_ = xMin;
        xMax_ = xMax;

        yMin_ = yMin;
        yMax_ = yMax;
    }

    Scalar xMin() const
    { return xMin_; }

    Scalar xMax() const
    { return xMax_; }

    Scalar yMin() const
    { return yMin_; }

    Scalar yMax() const
    { return yMax_; }

    unsigned numX() const
    { return m_; }

    unsigned numY() const
    { return n_; }

    Scalar iToX(unsigned i) const
    {
        assert(0 <= i && i < numX());

        return xMin() + i*(xMax() - xMin())/(numX() - 1);
    }

    Scalar jToY(unsigned j) const
    {
        assert(0 <= j && j < numY());

        return yMin() + j*(yMax() - yMin())/(numY() - 1);
    }

    template <class Evaluation>
    Evaluation xToI(const Evaluation& x) const
    { return (x - xMin())/(xMax() - xMin())*(numX() - 1); }

    template <class Evaluation>
    Evaluation yToJ(const Evaluation& y) const
    { return (y - yMin())/(yMax() - yMin())*(numY() - 1); }

    template <class Evaluation>
    bool applies(const Evaluation& x, const Evaluation& y) const
    {
        return
            xMin() <= x && x <= xMax() &&
            yMin() <= y && y <= yMax();
    }

    template <class Evaluation>
    Evaluation eval(const Evaluation& x, const Evaluation& y) const
    {
        typedef MathToolbox<Evaluation> Toolbox;

#ifndef NDEBUG
        if (!applies(x,y))
        {
            OPM_THROW(NumericalProblem,
                       "Attempt to get tabulated value for ("
                       << x << ", " << y
                       << ") on a table of extend "
                       << xMin() << " to " << xMax() << " times "
                       << yMin() << " to " << yMax());
        };
#endif

        Evaluation alpha = xToI(x);
        Evaluation beta = yToJ(y);

        unsigned i = std::max(0U, std::min(static_cast<unsigned>(numX()) - 2,
                                           static_cast<unsigned>(Toolbox::value(alpha))));
        unsigned j = std::max(0U, std::min(static_cast<unsigned>(numY()) - 2,
                                           static_cast<unsigned>(Toolbox::value(beta))));

        alpha -= i;
        beta -= j;

        // bi-linear interpolation
        const Evaluation& s1 = getSamplePoint(i, j)*(1.0 - alpha) + getSamplePoint(i + 1, j)*alpha;
        const Evaluation& s2 = getSamplePoint(i, j + 1)*(1.0 - alpha) + getSamplePoint(i + 1, j + 1)*alpha;
        return s1*(1.0 - beta) + s2*beta;
    }

    Scalar getSamplePoint(unsigned i, unsigned j) const
    {
        assert(0 <= i && i < m_);
        assert(0 <= j && j < n_);

        return samples_[j*m_ + i];
    }

    void setSamplePoint(unsigned i, unsigned j, Scalar value)
    {
        assert(0 <= i && i < m_);
        assert(0 <= j && j < n_);

        samples_[j*m_ + i] = value;
    }

private:
    // the vector which contains the values of the sample points
    // f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
    // instead!
    std::vector<Scalar> samples_;

    // the number of sample points in x direction
    unsigned m_;

    // the number of sample points in y direction
    unsigned n_;

    // the range of the tabulation on the x axis
    Scalar xMin_;
    Scalar xMax_;

    // the range of the tabulation on the y axis
    Scalar yMin_;
    Scalar yMax_;
};
} // namespace Opm

#endif