tabulated2dfunction.hh

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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_TABULATED_2D_FUNCTION_HH
#define DUMUX_TABULATED_2D_FUNCTION_HH
 
#include <vector>
#include <cassert>
#include <algorithm>
 
namespace Dumux {
 
template<class Scalar>
class Tabulated2DFunction
{
public:
    Tabulated2DFunction()
    { };
 
    Tabulated2DFunction(Scalar xMin, Scalar xMax, int m,
                        Scalar yMin, Scalar yMax, int n)
    {
        resize(xMin, xMax, m, yMin, yMax, n);
    };
 
    void resize(Scalar xMin, Scalar xMax, int m,
                Scalar yMin, Scalar yMax, int n)
    {
        samples_.resize(m*n);
 
        m_ = m;
        n_ = n;
 
        xMin_ = xMin;
        xMax_ = xMax;
 
        yMin_ = yMin;
        yMax_ = yMax;
    }
 
    Scalar iToX(int i) const
    {
        assert(0 <= i && i < m_);
 
        return xMin_ + i*(xMax_ - xMin_)/(m_ - 1);
    }
 
    Scalar jToY(int j) const
    {
        assert(0 <= j && j < n_);
 
        return yMin_ + j*(yMax_ - yMin_)/(n_ - 1);
    }
 
    Scalar xToI(Scalar x) const
    {
        return (x - xMin_)/(xMax_ - xMin_)*(m_ - 1);
    }
 
 
    Scalar yToJ(Scalar y) const
    {
        return (y - yMin_)/(yMax_ - yMin_)*(n_ - 1);
    }
 
 
    Scalar getSamplePoint(int i, int j) const
    {
        assert(0 <= i && i < m_);
        assert(0 <= j && j < n_);
 
        return samples_[j*m_ + i];
    }
 
    void setSamplePoint(int i, int j, Scalar value)
    {
        assert(0 <= i && i < m_);
        assert(0 <= j && j < n_);
 
        samples_[j*m_ + i] = value;
    }
 
    Scalar get(Scalar x, Scalar y) const
    {
        Scalar alpha = xToI(x);
        Scalar beta = yToJ(y);
 
        using std::clamp;
        int i = clamp(static_cast<int>(alpha), 0, m_ - 2);
        int j = clamp(static_cast<int>(beta),  0, n_ - 2);
 
        alpha -= i;
        beta -= j;
 
        // bi-linear interpolation
        Scalar s1 = getSamplePoint(i, j)*(1.0 - alpha) + getSamplePoint(i + 1, j)*alpha;
        Scalar s2 = getSamplePoint(i, j + 1)*(1.0 - alpha) + getSamplePoint(i + 1, j + 1)*alpha;
        return s1*(1.0 - beta) + s2*beta;
    }
 
    Scalar operator()(Scalar x, Scalar y) const
    { return get(x, y); }
 
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
    int m_;
 
    // the number of sample points in y direction
    int 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_;
 
};
 
} // end namespace Dumux
 
#endif