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
 
#include <cmath>
#include <iostream>
#include <tuple>
 
#include <dune/common/dotproduct.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/hybridutilities.hh>
 
#include "istlexception.hh"
 
// forward declaration
namespace Dune {
  template < typename... Args >
  class MultiTypeBlockVector;
}
 
#include "gsetc.hh"
 
namespace Dune {
 
  template <typename Arg0, typename... Args>
  struct FieldTraits< MultiTypeBlockVector<Arg0, Args...> >
  {
    typedef typename FieldTraits<Arg0>::field_type field_type;
    typedef typename FieldTraits<Arg0>::real_type real_type;
  };
  template < typename... Args >
  class MultiTypeBlockVector
  : public std::tuple<Args...>
  {
    typedef std::tuple<Args...> tupleType;
  public:
 
    typedef MultiTypeBlockVector<Args...> type;
 
    typedef double field_type;
 
    static constexpr std::size_t size()
    {
      return sizeof...(Args);
    }
 
    int count()
    {
      return sizeof...(Args);
    }
 
    template< std::size_t index >
    typename std::tuple_element<index,tupleType>::type&
    operator[] ( const std::integral_constant< std::size_t, index > indexVariable )
    {
      DUNE_UNUSED_PARAMETER(indexVariable);
      return std::get<index>(*this);
    }
 
    template< std::size_t index >
    const typename std::tuple_element<index,tupleType>::type&
    operator[] ( const std::integral_constant< std::size_t, index > indexVariable ) const
    {
      DUNE_UNUSED_PARAMETER(indexVariable);
      return std::get<index>(*this);
    }
 
    template<typename T>
    void operator= (const T& newval) {
      Dune::Hybrid::forEach(*this, [&](auto&& entry) {
        entry = newval;
      });
    }
 
    void operator+= (const type& newv) {
      using namespace Dune::Hybrid;
      forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
        (*this)[i] += newv[i];
      });
    }
 
    void operator-= (const type& newv) {
      using namespace Dune::Hybrid;
      forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
        (*this)[i] -= newv[i];
      });
    }
 
    // Once we have the IsNumber traits class the following
    // three implementations could be replaced by:
    //
    //    template<class T,
    //      std::enable_if_t< IsNumber<T>::value, int> = 0>
    //    void operator*= (const T& w) {
    //      Dune::Hybrid::forEach(*this, [&](auto&& entry) {
    //        entry *= w;
    //      });
    //    }
 
    void operator*= (const int& w) {
      Dune::Hybrid::forEach(*this, [&](auto&& entry) {
        entry *= w;
      });
    }
 
    void operator*= (const float& w) {
      Dune::Hybrid::forEach(*this, [&](auto&& entry) {
        entry *= w;
      });
    }
 
    void operator*= (const double& w) {
      Dune::Hybrid::forEach(*this, [&](auto&& entry) {
        entry *= w;
      });
    }
 
    field_type operator* (const type& newv) const {
      using namespace Dune::Hybrid;
      return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
        return a + (*this)[i]*newv[i];
      });
    }
 
    field_type dot (const type& newv) const {
      using namespace Dune::Hybrid;
      return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
        return a + (*this)[i].dot(newv[i]);
      });
    }
 
    typename FieldTraits<field_type>::real_type two_norm2() const {
      using namespace Dune::Hybrid;
      return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
        return a + entry.two_norm2();
      });
    }
 
    typename FieldTraits<field_type>::real_type two_norm() const {return sqrt(this->two_norm2());}
 
    typename FieldTraits<field_type>::real_type infinity_norm() const
    {
      using namespace Dune::Hybrid;
      using std::max;
      using real_type = typename FieldTraits<field_type>::real_type;
 
      real_type result = 0.0;
      // Compute max norm tracking appearing nan values
      // if the field type supports nan.
      ifElse(has_nan<field_type>(), [&](auto&& id) {
        // This variable will preserve any nan value
        real_type nanTracker = 1.0;
        forEach(*this, [&](auto&& entry) {
          real_type entryNorm = entry.infinity_norm();
          result = max(entryNorm, result);
          nanTracker += entryNorm;
        });
        // Incorporate possible nan value into result
        result *= (nanTracker / nanTracker);
      }, [&](auto&& id) {
        forEach(*this, [&](auto&& entry) {
          result = max(entry.infinity_norm(), result);
        });
      });
      return result;
    }
 
    template<typename Ta>
    void axpy (const Ta& a, const type& y) {
      using namespace Dune::Hybrid;
      forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
        (*this)[i].axpy(a, y[i]);
      });
    }
 
  };
 
 
 
  template <typename... Args>
  std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<Args...>& v) {
    using namespace Dune::Hybrid;
    forEach(integralRange(Dune::Hybrid::size(v)), [&](auto&& i) {
      s << "\t(" << i << "):\n" << v[i] << "\n";
    });
    return s;
  }
 
 
 
} // end namespace
 
#endif