VectorFacade.hpp

#pragma once
 
#include <string>
#include <type_traits>
#include <utility>
 
#include <dune/common/classname.hh>
#include <dune/common/deprecated.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/functions/functionspacebases/concepts.hh>
 
#include <amdis/Output.hpp>
#include <amdis/algorithm/ForEach.hpp>
#include <amdis/algorithm/InnerProduct.hpp>
#include <amdis/algorithm/Transform.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/FakeContainer.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/operations/Assigner.hpp>
#include <amdis/typetree/MultiIndex.hpp>
 
namespace AMDiS
{
  // States the vector can be in. Is changed on various insertion or gathering methods.
  enum class VectorState
  {
    unknown = 0,
    synchronized = 1,
    insert_values = 2,
    add_values = 3
  };
 
 
  template <class T, template <class> class VectorImpl>
  class VectorFacade
  {
    using Self = VectorFacade;
    using Impl = VectorImpl<T>;
 
  private:
    template <class A>
    using VectorStateOf_t = std::conditional_t<std::is_same_v<A,Assigner::plus_assign>,
      std::integral_constant<VectorState, VectorState::add_values>,
      std::integral_constant<VectorState, VectorState::insert_values>>;
 
  public:
    using size_type = typename Impl::size_type;
    using value_type = typename Impl::value_type;
 
    template <class GlobalBasis,
      class = std::void_t<decltype(std::declval<GlobalBasis>().indexDistribution())> >
    VectorFacade(GlobalBasis const& basis)
      : impl_(basis.indexDistribution())
    {
      resizeZero(basis);
    }
 
    Impl const& impl() const { return impl_; }
    Impl&       impl()       { return impl_; }
 
 
    template <class V>
    using HasLocalSize = decltype(std::declval<V>().localSize());
 
    template <class V>
    using HasGlobalSize = decltype(std::declval<V>().globalSize());
 
    std::size_t localSize() const
    {
      if constexpr (Dune::Std::is_detected<HasLocalSize,Impl>::value)
        return impl_.localSize();
      else
        return impl_.size();
    }
 
    std::size_t globalSize() const
    {
      if constexpr (Dune::Std::is_detected<HasGlobalSize,Impl>::value)
        return impl_.globalSize();
      else
        return impl_.size();
    }
 
    template <class Basis>
    void resize(Basis const& basis)
    {
      init(basis, false);
    }
 
    template <class Basis>
    void resizeZero(Basis const& basis)
    {
      init(basis, true);
    }
 
    template <class Basis>
    void init(Basis const& basis, bool clear)
    {
      impl_.init(basis, clear);
      state_ = clear ? VectorState::synchronized : VectorState::unknown;
    }
 
    void finish()
    {
      impl_.finish();
      state_ = VectorState::unknown;
    }
 
 
    template <class Index, class Assign = Assigner::plus_assign,
      REQUIRES(Concepts::MultiIndex<Index>)>
    void insert(Index const& idx, typename Impl::value_type const& value, Assign assign = {})
    {
      test_exit_dbg(state_ == VectorStateOf_t<Assign>::value ||
                    state_ == VectorState::unknown ||
                    state_ == VectorState::synchronized,
        "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
 
      impl_.insert(idx, value, assign);
 
      // set the state to insert_values or add_values
      state_ = VectorStateOf_t<Assign>::value;
    }
 
    template <class Index,
      REQUIRES(Concepts::MultiIndex<Index>)>
    void set(Index const& idx, typename Impl::value_type const& value)
    {
      insert(idx, value, Assigner::assign{});
    }
 
    template <class Index,
      REQUIRES(Concepts::MultiIndex<Index>)>
    void add(Index const& idx, typename Impl::value_type const& value)
    {
      insert(idx, value, Assigner::plus_assign{});
    }
 
    template <class Index,
      REQUIRES(Concepts::MultiIndex<Index>)>
    typename Impl::value_type get(Index const& idx) const
    {
      const_cast<Self*>(this)->synchronize();
      return impl_.at(idx);
    }
 
 
 
    template <class LocalView, class Node, class Buffer,
      REQUIRES(Concepts::LocalView<LocalView>),
      REQUIRES(Concepts::BasisNode<Node>)>
    void gather(LocalView const& localView, Node const& node, Buffer& buffer) const
    {
      test_exit(state_ == VectorState::unknown ||
                state_ == VectorState::synchronized,
        "Vector in invalid state {} for gather operations.", to_string(state_));
 
      const_cast<Self*>(this)->synchronize();
 
      buffer.resize(node.size());
      impl_.gather(nodeIndices(localView, node), buffer.begin());
    }
 
    // [[expects: localView is bound to an element]]
    template <class LocalView, class Buffer,
      REQUIRES(Concepts::LocalView<LocalView>)>
    void gather(LocalView const& localView, Buffer& buffer) const
    {
      gather(localView, localView.tree(), buffer);
    }
 
    template <class Index, class Buffer,
      REQUIRES(Concepts::MultiIndex<Index>)>
    void gather(std::vector<Index> const& indices, Buffer& buffer)
    {
      test_exit(state_ == VectorState::unknown ||
                state_ == VectorState::synchronized,
        "Vector in invalid state {} for gather operations.", to_string(state_));
 
      const_cast<Self*>(this)->synchronize();
 
      buffer.resize(indices.size());
      impl_.gather(indices, buffer.begin());
    }
 
 
    template <class LocalView, class Node, class NodeVector, class MaskRange, class Assign,
      REQUIRES(Concepts::LocalView<LocalView>),
      REQUIRES(Concepts::BasisNode<Node>)>
    void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector,
                 MaskRange const& mask, Assign assign)
    {
      test_exit(state_ == VectorStateOf_t<Assign>::value ||
                state_ == VectorState::unknown ||
                state_ == VectorState::synchronized,
        "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());
 
      assert(localVector.size() == node.size());
 
      // create a range of DOF indices on the node
      impl_.scatter(nodeIndices(localView, node), localVector, mask, assign);
 
      // set the state to insert_values or add_values
      state_ = VectorStateOf_t<Assign>::value;
    }
 
    // [[expects: localView is bound to an element]]
    // [[expects: node is in localView.tree()]]
    template <class LocalView, class Node, class NodeVector, class Assign,
      REQUIRES(Concepts::LocalView<LocalView>),
      REQUIRES(Concepts::BasisNode<Node>)>
    void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector, Assign assign)
    {
      scatter(localView, node, localVector, FakeContainer<bool,true>{}, assign);
    }
 
    // [[expects: localView is bound to an element]]
    template <class LocalView, class LocalVector, class Assign,
      REQUIRES(Concepts::LocalView<LocalView>)>
    void scatter(LocalView const& localView, LocalVector const& localVector, Assign assign)
    {
      scatter(localView, localView.tree(), localVector, assign);
    }
 
    template <class Index, class Buffer,
      REQUIRES(Concepts::MultiIndex<Index>)>
    void scatter(std::vector<Index> const& indices, Buffer const& values)
    {
      impl_.scatter(indices, values, FakeContainer<bool,true>{}, Assigner::assign{});
 
      // set the state to insert_values or add_values
      state_ = VectorStateOf_t<Assigner::assign>::value;
    }
 
 
 
    template <class S,
      std::enable_if_t<Dune::IsNumber<S>::value, int> = 0>
    Self& axpy (const S& alpha, const Self& x)
    {
      impl_.axpy(alpha, x.impl_);
      return *this;
    }
 
    template <class S,
      std::enable_if_t<Dune::IsNumber<S>::value, int> = 0>
    Self& aypx (const S& alpha, const Self& x)
    {
      impl_.aypx(alpha, x.impl_);
      return *this;
    }
 
    template <class S,
      std::enable_if_t<Dune::IsNumber<S>::value, int> = 0>
    Self& operator= (const S& alpha)
    {
      impl_.set(alpha);
      return *this;
    }
 
    Self& operator+= (const Self& x)
    {
      return this->axpy(1.0, x);
    }
 
    Self& operator-= (const Self& x)
    {
      return this->axpy(-1.0, x);
    }
 
    template <class S,
      std::enable_if_t<Dune::IsNumber<S>::value, int> = 0>
    Self& operator*= (const S& alpha)
    {
      impl_.scale(alpha);
      return *this;
    }
 
    template <class S,
      std::enable_if_t<Dune::IsNumber<S>::value, int> = 0>
    Self& operator/= (const S& alpha)
    {
      impl_.scale(S(1)/alpha);
      return *this;
    }
 
 
 
 
    template <class LocalInd, class Func>
    void forEach(LocalInd const& localInd, Func&& func)
    {
      synchronize();
      impl_.forEach(localInd, FWD(func));
    }
 
 
    template <class LocalInd, class Func>
    void forEach(LocalInd const& localInd, Func&& func) const
    {
      const_cast<Self*>(this)->synchronize();
      impl_.forEach(localInd, FWD(func));
    }
 
 
    VectorState state () const
    {
      return state_;
    }
 
  private:
    // synchronizes ghost values. Does not touch owned values
    void synchronize()
    {
      if (state_ != VectorState::synchronized)
        impl_.synchronize();
 
      state_ = VectorState::synchronized;
    }
 
    // print the vector state to string for debugging purposes
    static std::string to_string(VectorState state)
    {
      switch (state) {
        case VectorState::synchronized: return "synchronized";
        case VectorState::insert_values: return "insert_values";
        case VectorState::add_values: return "add_values";
        default: return "unknown";
      }
    }
 
  private:
    Impl impl_;
 
    VectorState state_ = VectorState::unknown;
  };
 
 
  namespace Recursive
  {
    template <class T, template <class> class Impl>
    struct ForEach<VectorFacade<T,Impl>>
    {
      template <class Vec, class UnaryFunction>
      static void impl (Vec&& vec, UnaryFunction f)
      {
        // NOTE: maybe needs synchronization
        Recursive::forEach(vec.impl(), f);
      }
    };
 
    template <class T, template <class> class Impl>
    struct Transform<VectorFacade<T,Impl>>
    {
      template <class Operation, class... Ts>
      static void impl (VectorFacade<T,Impl>& vecOut, Operation op, VectorFacade<Ts,Impl> const&... vecIn)
      {
        // NOTE: maybe needs synchronization
        Recursive::transform(vecOut.impl(), op, vecIn.impl()...);
      }
    };
 
    template <class S, template <class> class Impl>
    struct InnerProduct<VectorFacade<S,Impl>>
    {
      template <class In1, class In2, class T, class BinOp1, class BinOp2>
      static T impl (In1 const& in1, In2 const& in2, T init, BinOp1 op1, BinOp2 op2)
      {
        // NOTE: maybe needs synchronization
        return Recursive::innerProduct(in1.impl(), in2.impl(), std::move(init), op1, op2);
      }
    };
 
  } // end namespace Recursive
 
 
  namespace Impl
  {
    template <class T, class Facade>
    struct VectorTypeImpl;
 
    template <class T, class S, template <class> class Impl>
    struct VectorTypeImpl<T,VectorFacade<S,Impl>>
    {
      using type = VectorFacade<T,Impl>;
    };
 
  } // end namespace Impl
 
  template <class T, class Facade>
  using VectorType_t = typename Impl::VectorTypeImpl<T,Facade>::type;
 
} // end namespace AMDiS