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 <dune/common/typetraits.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:
 
    using size_type = std::size_t;
 
    using TupleType::TupleType;
 
    typedef MultiTypeBlockVector<Args...> type;
 
    typedef double field_type;
 
    static constexpr size_type size()
    {
      return sizeof...(Args);
    }
 
    static constexpr size_type N()
    {
      return sizeof...(Args);
    }
 
    int count() const DUNE_DEPRECATED_MSG("Use method 'N' instead")
    {
      return sizeof...(Args);
    }
 
    size_type dim() const
    {
      size_type result = 0;
      Hybrid::forEach(std::make_index_sequence<N()>{},
                      [&](auto i){result += std::get<i>(*this).dim();});
 
      return result;
    }
 
    template< size_type index >
    typename std::tuple_element<index,TupleType>::type&
    operator[] ( const std::integral_constant< size_type, index > indexVariable )
    {
      DUNE_UNUSED_PARAMETER(indexVariable);
      return std::get<index>(*this);
    }
 
    template< size_type index >
    const typename std::tuple_element<index,TupleType>::type&
    operator[] ( const std::integral_constant< size_type, 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];
      });
    }
 
    template<class T,
             std::enable_if_t< IsNumber<T>::value, int> = 0>
    void operator*= (const T& w) {
      Hybrid::forEach(*this, [&](auto&& entry) {
        entry *= w;
      });
    }
 
    template<class T,
             std::enable_if_t< IsNumber<T>::value, int> = 0>
    void operator/= (const T& w) {
      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]);
      });
    }
 
    auto one_norm() const {
      using namespace Dune::Hybrid;
      return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
        return a + entry.one_norm();
      });
    }
 
    auto one_norm_real() const {
      using namespace Dune::Hybrid;
      return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
        return a + entry.one_norm_real();
      });
    }
 
    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(HasNaN<field_type>(), [&](auto&& id) {
        // This variable will preserve any nan value
        real_type nanTracker = 1.0;
        using namespace Dune::Hybrid; // needed for icc, see issue #31
        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) {
        using namespace Dune::Hybrid; // needed for icc, see issue #31
        forEach(*this, [&](auto&& entry) {
          result = max(entry.infinity_norm(), result);
        });
      });
      return result;
    }
 
    typename FieldTraits<field_type>::real_type infinity_norm_real() 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(HasNaN<field_type>(), [&](auto&& id) {
        // This variable will preserve any nan value
        real_type nanTracker = 1.0;
        using namespace Dune::Hybrid; // needed for icc, see issue #31
        forEach(*this, [&](auto&& entry) {
          real_type entryNorm = entry.infinity_norm_real();
          result = max(entryNorm, result);
          nanTracker += entryNorm;
        });
        // Incorporate possible nan value into result
        result *= (nanTracker / nanTracker);
      }, [&](auto&& id) {
        using namespace Dune::Hybrid; // needed for icc, see issue #31
        forEach(*this, [&](auto&& entry) {
          result = max(entry.infinity_norm_real(), 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