SystemVector.hpp

#pragma once

#include <vector>
#include <string>
#include <memory>
#include <tuple>

#include <dune/amdis/Basic.hpp>
#include <dune/amdis/Log.hpp>
#include <dune/amdis/Loops.hpp>

#include <dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp>
#include <dune/amdis/linear_algebra/mtl/DOFVector.hpp>

// for explicit instantiation
#include <dune/amdis/ProblemStatBase.hpp>

namespace AMDiS
{
  /// \cond HIDDEN_SYMBOLS
  namespace Impl
  {
    // DOFVectors = std::tuple< DOFVector<FeSpace, ValueType>... >
    template <class DOFVectors>
    struct BuildDOFVectors
    {
      template <size_t I>
      using DOFVectorIdx = std::tuple_element_t<I, DOFVectors>;
      
      template <size_t... I, class FeSpaces, class MultiVector>
      static DOFVectors make(Seq<I...>, 
                              FeSpaces&& feSpaces, 
                              std::vector<std::string> const& names,
                              MultiVector&& multiVector)
      {
        return DOFVectors{DOFVectorIdx<I>(
          std::get<I>(std::forward<FeSpaces>(feSpaces)),
          names[I],
          multiVector[I] 
          )...};
      }
    };
        
    // construction method to construct a tuple of DOFVectors
    template <class DOFVectors, class FeSpaces, class MultiVector>
    DOFVectors build_dofvectors(FeSpaces&& feSpaces, 
                                    std::vector<std::string> const& names,
                                    MultiVector&& multiVector)
    {
      return BuildDOFVectors<DOFVectors>::make(
          MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
          std::forward<FeSpaces>(feSpaces),
          names,
          std::forward<MultiVector>(multiVector));
    }
    
  } // end namespace Impl
  /// \endcond
  
  
  /// \brief Container that repesents multiple data-Vectors of different value types.
  /** 
    *  Represents a std::array of \ref mtl::dense_vector + a tuple of corresponding 
    *  feSpaces. This I'th matrix combined with the I'th FeSpace 
    *  builds a \ref DOFVector  and can be return by the \ref operator().
    **/
  template <class FeSpaces, class ValueType = double>
  class SystemVector
  {
    using Self = SystemVector;
    
    template <class FeSpace>
    using DOFVectorGenerator = DOFVector<FeSpace, ValueType>;
    
  public:
    /// The number of blocks
    static constexpr size_t nComponents = std::tuple_size<FeSpaces>::value;
    AMDIS_STATIC_ASSERT( nComponents > 0 );
    
    /// The type of the block-vector
    using BaseVector = typename DOFVectorGenerator<std::tuple_element_t<0, FeSpaces>>::BaseVector;
    using MultiVector = BlockMTLVector<BaseVector, nComponents>;
    
    /// The type of the vector of DOFVectors
    using DOFVectors  = MakeTuple_t<DOFVectorGenerator, FeSpaces>;
    
    /// Constructor that constructs new base-vectors
    SystemVector(FeSpaces const& feSpaces, std::vector<std::string> const& names)
      : feSpaces(feSpaces)
      , names(names)
      , vector{}
      , dofvectors(Impl::build_dofvectors<DOFVectors>(feSpaces, names, vector))
    {
      compress();
    }
    

    /// Return the number of blocks
    static constexpr size_t size()
    {
      return std::tuple_size<FeSpaces>::value;
    }

    /// Return the number of blocks
    int count() const
    {
      return int(size());
    }
    
    /// Return a shared pointer to the I'th underlaying base vector
    template <size_t I>
    auto& getVector(const index_<I> = {})
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return vector[I];
    }
    
    /// Return a shared pointer to the I'th underlaying base vector
    template <size_t I>
    auto const& getVector(const index_<I> = {}) const
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return vector[I];
    }
    
    /// Return the underlying multi vector
    MultiVector&        getVector()       { return vector; }
    MultiVector const&  getVector() const { return vector; }
    
    /// Return the I'th finite element space
    template <size_t I = 0>
    auto const& getFeSpace(const index_<I> = {}) const
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return std::get<I>(feSpaces);
    }
    
    /// Return the name of the i'th DOFVector
    std::string getName(size_t i) const
    {
        return names[i];
    }
    
    /// Return the I'th DOFVector
    template <size_t I>
    auto& operator[](const index_<I>)
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return std::get<I>(dofvectors);
    }
    
    /// Return the I'th DOFVector
    template <size_t I>
    auto const& operator[](const index_<I>) const
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return std::get<I>(dofvectors);
    }
    
    /// Return the I'th DOFVector
    template <size_t I>
    auto& getDOFVector(const index_<I> = {})
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return std::get<I>(dofvectors);
    }
    
    /// Return the I'th DOFVector
    template <size_t I>
    auto const& getDOFVector(const index_<I> = {}) const
    {
      static_assert(I < size(), "Index out of range [0,SIZE)" );
      return std::get<I>(dofvectors);
    }
    
    /// Resize the I'th \ref vector to the sizes of the I'th \ref feSpaces.
    template <size_t I = 0>
    void compress(const index_<I> = {})
    {
      std::get<I>(dofvectors).compress();
    }
    
    /// Resize the \ref vectors to the sizes of the \ref feSpaces.
    void compress()
    {
      For<0, size()>::loop([this](const auto I) {
        std::get<I>(dofvectors).compress();
      });
    }
    
    /// Copy the Systemvector \p that element-wise
    void copy(Self const& that)
    {
      For<0, size()>::loop([this, &that](const auto I) {
        std::get<I>(dofvectors).copy(that[I]);
      });
    }
    
    /// Scale each DOFVector by the factor \p s.
    template <class Scalar>
    std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
    operator*=(Scalar s)
    {
      For<0, size()>::loop([this, s](const auto I) {
        vector[I] *= s;
      });
      return *this;
    }
    
    /// Sets each DOFVector to the scalar \p s.
    template <class Scalar>
    std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
    operator=(Scalar s)
    {
      For<0, size()>::loop([this, s](const auto I) {
        vector[I] = s;
      });
      return *this;
    }
      
  private:
    /// a tuple of feSpaces
    FeSpaces const& feSpaces;	
    
    /// The names of the DOFVectors
    std::vector<std::string> names;
    
    /// a tuple of base vectors, i.e. mtl::dense_vectors
    MultiVector vector;
    
    /// a tuple of dofvectors referencing \ref vector
    DOFVectors dofvectors;
  };
  
  
  /// Construct a systemvector from given feSpace-tuple and names-vector
  template <class FeSpacesType, class ValueType = double>
  SystemVector<FeSpacesType, ValueType>
  make_systemvector(FeSpacesType const& feSpaces, std::vector<std::string> names)
  {
    return {feSpaces, names};
  }
  
  /// Construct a systemvector from given feSpace-tuple. Create vector of names
  /// from base_name, i.e. {base_name_0, base_name_1, ...}
  template <class FeSpaces, class ValueType = double>
  SystemVector<FeSpaces, ValueType>
  make_systemvector(FeSpaces const& feSpaces, std::string base_name)
  {
    std::vector<std::string> names;
    for (size_t i = 0; i < std::tuple_size<FeSpaces>::value; ++i)
      names.push_back(base_name + "_" + std::to_string(i));
    return {feSpaces, names};
  }
    
  
#ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
    // explicit instantiation in SystemVector.cpp
    extern template class SystemVector<typename ProblemStatTraits<2,2,1,2>::FeSpaces>;
#endif

} // end namespace AMDiS