ModelArray.hpp

 
#ifndef MODELARRAY_HPP
#define MODELARRAY_HPP
 
#include <Eigen/Core>
#include <cstddef>
#include <map>
#include <string>
#include <utility>
#include <vector>
 
namespace Nextsim {
 
/*
 * Set the storage order to row major. This matches with DGVector when there is
 * more than one DG component. If there is only one DG component (the finite
 * element component), then the order of the data in the buffer is the same,
 * and data can be safely transferred between the buffers underlying the two
 * types. The physics code will never directly touch the spatially varying
 * components, so the choice of storage order should not matter.
 */
const static Eigen::StorageOptions majority = Eigen::RowMajor;
 
class ModelArray {
public:
    // Forward defines make Eclipse less red and squiggly
    enum class Type;
    enum class Dimension;
 
#include "include/ModelArrayDetails.hpp"
 
    static const int N_DEFINED_DIMENSIONS = static_cast<int>(Dimension::COUNT);
 
    struct DimensionSpec {
        std::string name;
        size_t length;
    };
    typedef std::map<Type, std::vector<Dimension>> TypeDimensions;
 
    static TypeDimensions typeDimensions;
    static std::map<Dimension, DimensionSpec> definedDimensions;
    static const std::map<Type, std::string> typeNames;
    // The dimension that defines the components of each ModelArray type, if any
    static const std::map<Type, Dimension> componentMap;
 
    typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, majority> DataType;
 
    typedef DataType::RowXpr Component;
    typedef DataType::ConstRowXpr ConstComponent;
 
    ModelArray();
    ModelArray(const Type type);
    ModelArray(const ModelArray&);
    virtual ~ModelArray() {};
 
    ModelArray& operator=(const ModelArray&);
    ModelArray& operator=(const double& val);
 
    // ModelArray arithmetic
    ModelArray& operator+=(const ModelArray& b)
    {
        m_data += b.m_data;
        return *this;
    }
    ModelArray& operator-=(const ModelArray& b)
    {
        m_data -= b.m_data;
        return *this;
    }
    ModelArray& operator*=(const ModelArray& b)
    {
        m_data *= b.m_data;
        return *this;
    }
    ModelArray& operator/=(const ModelArray& b)
    {
        m_data /= b.m_data;
        return *this;
    }
 
    ModelArray& operator+=(double b)
    {
        m_data += b;
        return *this;
    }
    ModelArray& operator-=(double b)
    {
        m_data -= b;
        return *this;
    }
    ModelArray& operator*=(double b)
    {
        m_data *= b;
        return *this;
    }
    ModelArray& operator/=(double b)
    {
        m_data /= b;
        return *this;
    }
 
    ModelArray operator+(const ModelArray&) const;
    ModelArray operator-(const ModelArray&) const;
    ModelArray operator*(const ModelArray&) const;
    ModelArray operator/(const ModelArray&) const;
    // Returns a ModelArray containing the element-wise negation of this.
    ModelArray operator-() const;
 
    ModelArray operator+(const double&) const;
    ModelArray operator-(const double&) const;
    ModelArray operator*(const double&) const;
    ModelArray operator/(const double&) const;
 
    ModelArray max(double max) const;
    ModelArray min(double min) const;
    ModelArray max(const ModelArray& maxArr) const;
    ModelArray min(const ModelArray& minArr) const;
 
    ModelArray& clampAbove(double max);
    ModelArray& clampBelow(double min);
    ModelArray& clampAbove(const ModelArray& maxArr);
    ModelArray& clampBelow(const ModelArray& minArr);
 
    typedef std::vector<size_t> MultiDim;
 
    size_t nDimensions() const { return nDimensions(type); }
    static size_t nDimensions(Type type) { return m_dims.at(type).size(); }
    const MultiDim& dimensions() const { return dimensions(type); }
    static const MultiDim& dimensions(Type type) { return m_dims.at(type); }
    size_t size() const { return size(type); }
    static size_t size(Type type) { return m_sz.at(type); }
    size_t trueSize() const { return m_data.rows(); }
 
    const double* getData() const { return m_data.data(); }
 
    const DataType& data() const { return m_data; }
    Type getType() const { return type; }
 
    static void setDimensions(Type, const MultiDim&);
    void setDimensions(const MultiDim& dims)
    {
        setDimensions(type, dims);
        resize();
    }
 
    static void setDimension(Dimension dim, size_t length);
 
    void resize()
    {
        if (size() != trueSize()) {
            if (hasDoF(type)) {
                m_data.resize(m_sz.at(type), definedDimensions.at(componentMap.at(type)).length);
            } else {
                m_data.resize(m_sz.at(type), Eigen::NoChange);
            }
        }
    }
 
    void setData(double value);
    void setData(const double* pData);
    void setData(const DataType& data);
    void setData(const ModelArray& source);
 
private:
    // Fast special case for 1-d indexing
    template <typename T, typename I> static inline T indexr(const T* dims, I first)
    {
        return static_cast<T>(first);
    }
 
    // Fast special case for 2-d indexing
    template <typename T, typename I> static inline T indexr(const T* dims, I first, I second)
    {
        return first + second * dims[0];
    }
 
    // Indices as separate function parameters
    template <typename T, typename I, typename... Args>
    static inline T indexr(const T* dims, I first, Args... args)
    {
        std::initializer_list<I> loc { first, args... };
        return indexrHelper(dims, loc);
    }
 
    // Indices as a Dimensions object
    template <typename T> static T indexr(const T* dims, const ModelArray::MultiDim& loc)
    {
        return indexrHelper(dims, loc);
    }
 
    // Generic index generator that will work on any container
    template <typename T, typename C> static T indexrHelper(const T* dims, const C& loc)
    {
        size_t ndims = loc.size();
        T stride = 1;
        T ii = 0;
        auto iloc = begin(loc);
        for (size_t dim = 0; dim < ndims; ++dim) {
            ii += stride * (*iloc++);
            stride *= dims[dim];
        }
        return ii;
    }
 
public:
    const double& operator[](size_t i) const { return m_data(i, 0); }
    const double& operator[](const MultiDim& dims) const;
 
    template <typename... Args> const double& operator()(Args... args) const
    {
        return (*this)[indexr(dimensions().data(), args...)];
    }
 
    double& operator[](size_t i) { return const_cast<double&>(std::as_const(*this)(i)); }
    double& operator[](const MultiDim&);
    template <typename... Args> double& operator()(Args... args)
    {
        return const_cast<double&>(std::as_const(*this)(args...));
    }
 
    static void setNComponents(std::map<Type, size_t> cMap);
 
    static void setNComponents(Type type, size_t nComp)
    {
        if (hasDoF(type)) {
            definedDimensions.at(componentMap.at(type)).length = nComp;
        }
    }
 
    inline void setNComponents(size_t nComp) { setNComponents(type, nComp); }
 
    Component components(size_t i) { return m_data.row(i); }
 
    const ConstComponent components(size_t i) const { return m_data.row(i); }
 
    Component components(const MultiDim& loc);
    const ConstComponent components(const MultiDim& loc) const;
 
    double& zIndexAndLayer(size_t hIndex, size_t layer)
    {
        return this->operator[](zLayerIndex(hIndex, layer));
    }
    const double& zIndexAndLayer(size_t hIndex, size_t layer) const
    {
        return this->operator[](zLayerIndex(hIndex, layer));
    }
 
    size_t indexFromLocation(const MultiDim& loc) const { return indexFromLocation(type, loc); }
 
    static size_t indexFromLocation(Type type, const MultiDim& loc);
 
    MultiDim locationFromIndex(size_t index) const { return locationFromIndex(type, index); }
 
    static MultiDim locationFromIndex(Type type, size_t index);
 
    inline size_t nComponents() const { return nComponents(type); }
    inline static size_t nComponents(const Type type)
    {
        return (hasDoF(type)) ? definedDimensions.at(componentMap.at(type)).length : 1;
    }
    inline bool hasDoF() const { return hasDoF(type); }
    static bool hasDoF(const Type type);
 
protected:
    Type type;
 
    size_t zLayerIndex(size_t hIndex, size_t layer) const
    {
        return hIndex + layer * dimensions()[0] * dimensions()[1];
    }
 
private:
    static bool areMapsInvalid;
    static void validateMaps();
    class SizeMap {
    public:
        SizeMap();
        size_t& at(const Type& type) { return m_sizes.at(type); }
        const size_t& at(const Type& type) const { return m_sizes.at(type); }
 
        size_t& operator[](const Type& type) { return m_sizes[type]; }
        size_t& operator[](Type&& type) { return m_sizes[type]; }
 
        size_t size() const noexcept { return m_sizes.size(); }
 
        void validate();
 
    protected:
        std::map<Type, size_t> m_sizes;
    };
    static SizeMap m_sz;
 
    class DimensionMap {
    public:
        DimensionMap();
        MultiDim& at(const Type& type) { return m_dimensions.at(type); }
        const MultiDim& at(const Type& type) const { return m_dimensions.at(type); }
 
        MultiDim& operator[](const Type& type) { return m_dimensions[type]; }
        MultiDim& operator[](Type&& type) { return m_dimensions[type]; }
 
        size_t size() const noexcept { return m_dimensions.size(); }
 
        void validate();
 
    private:
        std::map<Type, MultiDim> m_dimensions;
    };
    static DimensionMap m_dims;
    DataType m_data;
};
 
#include "include/ModelArrayTypedefs.hpp"
 
// ModelArray arithmetic with doubles
ModelArray operator+(const double&, const ModelArray&);
ModelArray operator-(const double&, const ModelArray&);
ModelArray operator*(const double&, const ModelArray&);
ModelArray operator/(const double&, const ModelArray&);
} /* namespace Nextsim */
 
#endif /* MODELARRAY_HPP */