PAvgDynamicSourceData.hpp

/*
  Copyright 2023 Equinor 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 PAVE_DYNAMIC_SOURCE_DATA_HPP
#define PAVE_DYNAMIC_SOURCE_DATA_HPP
 
#include <algorithm>
#include <cstddef>
#include <optional>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <vector>
 
namespace Opm {

template <typename Scalar>
class PAvgDynamicSourceData
{
public:
    template <typename T>
    class SourceDataSpan
    {
    private:
        friend class PAvgDynamicSourceData<Scalar>;
 
    public:
        enum class Item
        {
            Pressure,           //< Dynamic pressure value
            MixtureDensity,     //< Dynamic mixture density
            PoreVol,            //< Dynamic pore volume
            Depth,              //< Constant depth location
 
            // ----------------------------------------
 
            Last_Do_Not_Use, //< Simplifies item count
        };
 
        using ElmT = std::remove_cv_t<T>;

        [[nodiscard]] constexpr ElmT operator[](const Item i) const
        {
            return this->begin_[this->index(i)];
        }

        template <typename Ret = SourceDataSpan&>
        constexpr std::enable_if_t<! std::is_const_v<T>, Ret>
        set(const Item i, const ElmT value)
        {
            this->begin_[this->index(i)] = value;
            return *this;
        }

        template <typename U, typename Ret = SourceDataSpan&>
        constexpr std::enable_if_t<! std::is_const_v<T>, Ret>
        operator=(const SourceDataSpan<U> src)
        {
            std::copy_n(src.begin_, NumItems, this->begin_);
            return *this;
        }
 
    private:
        static constexpr auto NumItems =
            static_cast<std::size_t>(Item::Last_Do_Not_Use);

        T* begin_{nullptr};

        explicit SourceDataSpan(T* begin)
            : begin_{begin}
        {}

        constexpr std::size_t index(const Item i) const
        {
            const auto ix = static_cast<std::size_t>(i);
            if (ix >= NumItems) {
                throw std::invalid_argument {
                    "Index out of bounds"
                };
            }
 
            return ix;
        }
    };

    explicit PAvgDynamicSourceData(const std::vector<std::size_t>& sourceLocations);

    virtual ~PAvgDynamicSourceData() {}

    [[nodiscard]] SourceDataSpan<Scalar>
    operator[](const std::size_t source);

    [[nodiscard]] SourceDataSpan<const Scalar>
    operator[](const std::size_t source) const;
 
protected:
    std::vector<Scalar> src_{};

    [[nodiscard]] SourceDataSpan<Scalar>
    sourceTerm(const std::size_t ix, std::vector<Scalar>& src);

    void reconstruct(const std::vector<std::size_t>& sourceLocations);

    static constexpr std::size_t numSpanItems() noexcept
    {
        return SourceDataSpan<const Scalar>::NumItems;
    }
 
private:
    std::unordered_map<std::size_t, typename std::vector<Scalar>::size_type> ix_{};

    void buildLocationMapping(const std::vector<std::size_t>& sourceLocations);

    [[nodiscard]] std::optional<typename std::vector<Scalar>::size_type>
    index(const std::size_t source) const;

    [[nodiscard]] virtual typename std::vector<Scalar>::size_type
    storageIndex(typename std::vector<Scalar>::size_type elemIndex) const
    {
        return elemIndex;
    }
};
 
} // namespace Opm
 
#endif // PAVE_DYNAMIC_SOURCE_DATA_HPP