spline.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_SPLINE_HH
#define DUMUX_SPLINE_HH
 
#include "fixedlengthspline_.hh"
#include "variablelengthspline_.hh"
#include "splinecommon_.hh"
 
namespace Dumux {
 
template<class Scalar, int numSamples = 2>
class Spline : public FixedLengthSpline_<Scalar, numSamples>
{
public:
    Spline()
    { };
 
    template <class ScalarArray>
    Spline(const ScalarArray &x,
           const ScalarArray &y)
    { this->setXYArrays(numSamples, x, y); }
 
    template <class PointArray>
    Spline(const PointArray &points)
    { this->setArrayOfPoints(numSamples, points); }
 
    template <class ScalarArray>
    Spline(const ScalarArray &x,
           const ScalarArray &y,
           Scalar m0,
           Scalar m1)
    { this->setXYArrays(numSamples, x, y, m0, m1); }
 
    template <class PointArray>
    Spline(const PointArray &points,
           Scalar m0,
           Scalar m1)
    { this->setArrayOfPoints(numSamples, points, m0, m1); }
};
 
template<class Scalar>
class Spline<Scalar, /*numSamples=*/-1> : public VariableLengthSpline_<Scalar>
{
public:
    Spline()
    { }
 
    template <class ScalarArrayX, class ScalarArrayY>
    Spline(int nSamples,
           const ScalarArrayX &x,
           const ScalarArrayY &y)
    { this->setXYArrays(nSamples, x, y); }
 
    template <class PointArray>
    Spline(int nSamples,
           const PointArray &points)
    { this->setArrayOfPoints(nSamples, points); }
 
    template <class ScalarContainer>
    Spline(const ScalarContainer &x,
           const ScalarContainer &y)
    { this->setXYContainers(x, y); }
 
    template <class PointContainer>
    Spline(const PointContainer &points)
    { this->setContainerOfPoints(points); }
 
    template <class ScalarArray>
    Spline(int nSamples,
           const ScalarArray &x,
           const ScalarArray &y,
           Scalar m0,
           Scalar m1)
    { this->setXYArrays(nSamples, x, y, m0, m1); }
 
    template <class PointArray>
    Spline(int nSamples,
           const PointArray &points,
           Scalar m0,
           Scalar m1)
    { this->setArrayOfPoints(nSamples, points, m0, m1); }
 
    template <class ScalarContainerX, class ScalarContainerY>
    Spline(const ScalarContainerX &x,
           const ScalarContainerY &y,
           Scalar m0,
           Scalar m1)
    { this->setXYContainers(x, y, m0, m1); }
 
    template <class PointContainer>
    Spline(const PointContainer &points,
           Scalar m0,
           Scalar m1)
    { this->setContainerOfPoints(points, m0, m1); }
};
 
template<class Scalar>
class Spline<Scalar, /*numSamples=*/0>
// Splines with zero sampling points do not make sense!
{ private: Spline() { }; };
 
template<class Scalar>
class Spline<Scalar, /*numSamples=*/1>
// Splines with one sampling point do not make sense!
{ private: Spline() { }; };
 
template<class Scalar>
class Spline<Scalar, 2> : public SplineCommon_<Scalar, Spline<Scalar, 2> >
{
    friend class  SplineCommon_<Scalar, Spline<Scalar, 2> >;
    using Vector = Dune::FieldVector<Scalar, 2>;
    using Matrix = Dune::FieldMatrix<Scalar, 2, 2>;
 
public:
    Spline()
    {};
 
    template <class ScalarArrayX, class ScalarArrayY>
    Spline(const ScalarArrayX &x,
           const ScalarArrayY &y,
           Scalar m0, Scalar m1)
    { setXYArrays(2, x, y, m0, m1); }
 
    template <class PointArray>
    Spline(const PointArray &points,
           Scalar m0,
           Scalar m1)
    { this->setArrayOfPoints(2, points, m0, m1); }
 
    Spline(Scalar x0, Scalar x1,
           Scalar y0, Scalar y1,
           Scalar m0, Scalar m1)
    {
        set(x0, x1,
            y0, y1,
            m0, m1);
    };
 
    int numSamples() const
    { return 2; }
 
    void set(Scalar x0, Scalar x1,
             Scalar y0, Scalar y1,
             Scalar m0, Scalar m1)
    {
        Matrix M(numSamples());
        Vector d;
        assignXY_(x0, x1, y0, y1);
        this->makeFullSystem_(M, d, m0, m1);
 
        // solve for the moments
        M.solve(m_, d);
    }
 
    template <class ScalarContainer>
    void setXYArrays(int nSamples,
                     const ScalarContainer &x,
                     const ScalarContainer &y,
                     Scalar m0, Scalar m1)
    {
        assert(nSamples == 2);
        set(x[0], x[1], y[0], y[1], m0, m1);
    }
 
    template <class ScalarContainerX, class ScalarContainerY>
    void setXYContainers(const ScalarContainerX &x,
                         const ScalarContainerY &y,
                         Scalar m0, Scalar m1)
    {
        assert(x.size() == y.size());
        assert(x.size() == 2);
 
        Matrix M(numSamples());
        Vector d;
 
        typename ScalarContainerX::const_iterator xIt0 = x.begin();
        typename ScalarContainerX::const_iterator xIt1 = xIt0;
        ++xIt1;
        typename ScalarContainerY::const_iterator yIt0 = y.begin();
        typename ScalarContainerY::const_iterator yIt1 = yIt0;
        ++yIt1;
        set(*xIt0, *xIt1, *yIt0, *yIt1);
    }
 
    template <class PointArray>
    void setArrayOfPoints(int nSamples,
                          const PointArray &points,
                          Scalar m0,
                          Scalar m1)
    {
        assert(nSamples == 2);
 
        set(points[0][0],
            points[1][0],
            points[0][1],
            points[1][1],
            m0, m1);
    }
 
    template <class PointContainer>
    void setContainerOfPoints(const PointContainer &points,
                              Scalar m0,
                              Scalar m1)
    {
        assert(points.size() == 2);
 
        Matrix M;
        Vector d;
        typename PointContainer::const_iterator it0 = points.begin();
        typename PointContainer::const_iterator it1 = it0;
        ++it1;
 
        set((*it0)[0],
            (*it0)[1],
            (*it1)[0],
            (*it1)[1],
            m0, m1);
    }
 
    template <class TupleContainer>
    void setContainerOfTuples(const TupleContainer &tuples,
                              Scalar m0,
                              Scalar m1)
    {
        assert(tuples.size() == 2);
 
        typename TupleContainer::const_iterator it0 = tuples.begin();
        typename TupleContainer::const_iterator it1 = it0;
        ++it1;
 
        set(std::get<0>(*it0),
            std::get<1>(*it0),
            std::get<0>(*it1),
            std::get<1>(*it1),
            m0, m1);
    }
 
protected:
    void assignXY_(Scalar x0, Scalar x1,
                   Scalar y0, Scalar y1)
    {
        if (x0 > x1) {
            xPos_[0] = x1;
            xPos_[1] = x0;
            yPos_[0] = y1;
            yPos_[1] = y0;
        }
        else {
            xPos_[0] = x0;
            xPos_[1] = x1;
            yPos_[0] = y0;
            yPos_[1] = y1;
        }
    };
 
    Scalar x_(int i) const
    { return xPos_[i]; }
 
    Scalar y_(int i) const
    { return yPos_[i]; }
 
    Scalar moment_(int i) const
    { return m_[i]; }
 
    Vector xPos_;
    Vector yPos_;
    Vector m_;
};
 
} // end namespace Dumux
 
#endif