PiecewiseConstantInterpolation.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <algorithm>
#include <vector>
 
#include "BaseLib/Logging.h"
 
namespace MathLib
{
template <typename T = double>
class PiecewiseConstantInterpolation
{
public:
    PiecewiseConstantInterpolation(
        std::vector<T> const& supporting_points,
        std::vector<double> const& values_at_supp_pnts)
        : supp_pnts_(supporting_points),
          values_at_supp_pnts_(values_at_supp_pnts)
    {
        if (supp_pnts_.size() != values_at_supp_pnts_.size())
        {
            OGS_FATAL(
                "Inconsistent data given to PiecewiseConstantInterpolation, "
                "number of given supporting points is {}, number of given "
                "values is {}.",
                supp_pnts_.size(), values_at_supp_pnts_.size());
        }
        if (supp_pnts_.empty())
        {
            ERR("PiecewiseConstantInterpolation: passed empty vector.");
        }
    }

    double value(double const pnt_to_interpolate) const
    {
        if (pnt_to_interpolate <= supp_pnts_.front())
        {
            return values_at_supp_pnts_.front();
        }
 
        if (supp_pnts_.back() <= pnt_to_interpolate)
        {
            return values_at_supp_pnts_.back();
        }
 
        auto const& it(std::upper_bound(supp_pnts_.begin(), supp_pnts_.end(),
                                        pnt_to_interpolate));
        // Here access the iterator it without checking is okay since the
        // corner cases are checked above.
        auto const interval_idx = std::distance(supp_pnts_.begin(), it) - 1;
 
        return values_at_supp_pnts_[interval_idx];
    }
 
private:
    std::vector<T> const& supp_pnts_;
    std::vector<double> const& values_at_supp_pnts_;
};
}  // end namespace MathLib