multitypeblockvector.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
#define DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
 
#if HAVE_DUNE_BOOST
#ifdef HAVE_BOOST_FUSION
 
#include <cmath>
#include <iostream>
 
#include <dune/common/dotproduct.hh>
#include <dune/common/ftraits.hh>
 
#include "istlexception.hh"
 
#include <boost/fusion/sequence.hpp>
#include <boost/fusion/container.hpp>
#include <boost/fusion/iterator.hpp>
#include <boost/typeof/typeof.hpp>
#include <boost/fusion/algorithm.hpp>
 
namespace mpl=boost::mpl;
namespace fusion=boost::fusion;
 
// forward decl
namespace Dune {
  template<typename T1, typename T2=fusion::void_, typename T3=fusion::void_, typename T4=fusion::void_,
      typename T5=fusion::void_, typename T6=fusion::void_, typename T7=fusion::void_,
      typename T8=fusion::void_, typename T9=fusion::void_>
  class MultiTypeBlockVector;
}
 
#include "gsetc.hh"
 
namespace Dune {
 
  template<int current_element, int remaining_elements, typename TVec>
  class MultiTypeBlockVector_Print {
  public:
    static void print(const TVec& v) {
      std::cout << "\t(" << current_element << "):\n" << fusion::at_c<current_element>(v) << "\n";
      MultiTypeBlockVector_Print<current_element+1,remaining_elements-1,TVec>::print(v);   //next element
    }
  };
  template<int current_element, typename TVec>            //recursion end (remaining_elements=0)
  class MultiTypeBlockVector_Print<current_element,0,TVec> {
  public:
    static void print(const TVec& v) {std::cout << "\n";}
  };
 
 
 
  template<int count, typename T1, typename T2>
  class MultiTypeBlockVector_Ident {
  public:
 
    static void equalize(T1& a, const T2& b) {
      fusion::at_c<count-1>(a) = b;           //equalize current elements
      MultiTypeBlockVector_Ident<count-1,T1,T2>::equalize(a,b);    //next elements
    }
  };
  template<typename T1, typename T2>                      //recursion end (count=0)
  class MultiTypeBlockVector_Ident<0,T1,T2> {public: static void equalize (T1& a, const T2& b) {} };
 
 
 
 
 
 
  template<int count, typename T>
  class MultiTypeBlockVector_Add {
  public:
 
    static void add (T& a, const T& b) {    //add vector elements
      fusion::at_c<(count-1)>(a) += fusion::at_c<(count-1)>(b);
      MultiTypeBlockVector_Add<count-1,T>::add(a,b);
    }
 
    static void sub (T& a, const T& b) {    //sub vector elements
      fusion::at_c<(count-1)>(a) -= fusion::at_c<(count-1)>(b);
      MultiTypeBlockVector_Add<count-1,T>::sub(a,b);
    }
  };
  template<typename T>                                    //recursion end; specialization for count=0
  class MultiTypeBlockVector_Add<0,T> {public: static void add (T& a, const T& b) {} static void sub (T& a, const T& b) {} };
 
 
 
  template<int count, typename TVec, typename Ta>
  class MultiTypeBlockVector_AXPY {
  public:
 
    static void axpy(TVec& x, const Ta& a, const TVec& y) {
      fusion::at_c<(count-1)>(x).axpy(a,fusion::at_c<(count-1)>(y));
      MultiTypeBlockVector_AXPY<count-1,TVec,Ta>::axpy(x,a,y);
    }
  };
  template<typename TVec, typename Ta>                    //specialization for count=0
  class MultiTypeBlockVector_AXPY<0,TVec,Ta> {public: static void axpy (TVec& x, const Ta& a, const TVec& y) {} };
 
 
  template<int count, typename TVec, typename Ta>
  class MultiTypeBlockVector_Mulscal {
  public:
 
    static void mul(TVec& x, const Ta& a) {
      fusion::at_c<(count-1)>(x) *= a;
      MultiTypeBlockVector_Mulscal<count-1,TVec,Ta>::mul(x,a);
    }
  };
  template<typename TVec, typename Ta>                    //specialization for count=0
  class MultiTypeBlockVector_Mulscal<0,TVec,Ta> {public: static void mul (TVec& x, const Ta& a) {} };
 
 
 
  template<int count, typename TVec>
  class MultiTypeBlockVector_Mul {
  public:
    static typename TVec::field_type mul(const TVec& x, const TVec& y) { return (fusion::at_c<count-1>(x) * fusion::at_c<count-1>(y)) + MultiTypeBlockVector_Mul<count-1,TVec>::mul(x,y); }
    static typename TVec::field_type dot(const TVec& x, const TVec& y) { return (Dune::dot(fusion::at_c<count-1>(x),fusion::at_c<count-1>(y))) + MultiTypeBlockVector_Mul<count-1,TVec>::dot(x,y); }
  };
  template<typename TVec>
  class MultiTypeBlockVector_Mul<0,TVec> {
  public:
    static typename TVec::field_type mul(const TVec& x, const TVec& y) {return 0;}
    static typename TVec::field_type dot(const TVec& x, const TVec& y) {return 0;}
  };
 
 
 
 
 
  template<int count, typename T>
  class MultiTypeBlockVector_Norm {
  public:
    typedef typename T::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
 
    static real_type result (const T& a) {             //result = sum of all elements' 2-norms
      return fusion::at_c<count-1>(a).two_norm2() + MultiTypeBlockVector_Norm<count-1,T>::result(a);
    }
  };
 
  template<typename T>                                    //recursion end: no more vector elements to add...
  class MultiTypeBlockVector_Norm<0,T> {
  public:
    typedef typename T::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
    static real_type result (const T& a) {return 0.0;}
  };
 
  template<typename T1, typename T2, typename T3, typename T4,
      typename T5, typename T6, typename T7, typename T8, typename T9>
  class MultiTypeBlockVector : public fusion::vector<T1, T2, T3, T4, T5, T6, T7, T8, T9> {
 
  public:
 
    typedef MultiTypeBlockVector<T1, T2, T3, T4, T5, T6, T7, T8, T9> type;
 
    typedef typename T1::field_type field_type;
 
    int count() {return mpl::size<type>::value;}
 
    template< std::size_t index >
    auto
    operator[] ( const std::integral_constant< std::size_t, index > indexVariable )
    -> decltype(fusion::at_c<index>(*this))
    {
      DUNE_UNUSED_PARAMETER(indexVariable);
      return fusion::at_c<index>(*this);
    }
 
    template< std::size_t index >
    auto
    operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const
    -> decltype(fusion::at_c<index>(*this))
    {
      DUNE_UNUSED_PARAMETER(indexVariable);
      return fusion::at_c<index>(*this);
    }
 
    template<typename T>
    void operator= (const T& newval) {MultiTypeBlockVector_Ident<mpl::size<type>::value,type,T>::equalize(*this, newval); }
 
    void operator+= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::add(*this,newv);}
 
    void operator-= (const type& newv) {MultiTypeBlockVector_Add<mpl::size<type>::value,type>::sub(*this,newv);}
 
    void operator*= (const int& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const int>::mul(*this,w);}
    void operator*= (const float& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const float>::mul(*this,w);}
    void operator*= (const double& w) {MultiTypeBlockVector_Mulscal<mpl::size<type>::value,type,const double>::mul(*this,w);}
 
    field_type operator* (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::mul(*this,newv);}
    field_type dot (const type& newv) const {return MultiTypeBlockVector_Mul<mpl::size<type>::value,type>::dot(*this,newv);}
 
    typename FieldTraits<field_type>::real_type two_norm2() const {return MultiTypeBlockVector_Norm<mpl::size<type>::value,type>::result(*this);}
 
    typename FieldTraits<field_type>::real_type two_norm() const {return sqrt(this->two_norm2());}
 
    template<typename Ta>
    void axpy (const Ta& a, const type& y) {
      MultiTypeBlockVector_AXPY<mpl::size<type>::value,type,Ta>::axpy(*this,a,y);
    }
 
  };
 
 
 
  template<typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9>
  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9>& v) {
    MultiTypeBlockVector_Print<0,mpl::size<MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::value,MultiTypeBlockVector<T1,T2,T3,T4,T5,T6,T7,T8,T9> >::print(v);
    return s;
  }
 
 
 
} // end namespace
 
#endif // end HAVE_BOOST_FUSION
#endif // end HAVE_DUNE_BOOST
 
#endif