multitypeblockvector.hh

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- 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 <dune/common/std/type_traits.hh>
 
#include "istlexception.hh"
 
// forward declaration
namespace Dune {
  template < typename... Args >
  class MultiTypeBlockVector;
}
 
#include "gsetc.hh"
 
namespace Dune {
 
  template <typename... Args>
  struct FieldTraits< MultiTypeBlockVector<Args...> >
  {
    using field_type = typename MultiTypeBlockVector<Args...>::field_type;
    using real_type  = typename MultiTypeBlockVector<Args...>::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;
 
    using field_type = Std::detected_t<std::common_type_t, typename FieldTraits< std::decay_t<Args> >::field_type...>;
 
    using real_type = Std::detected_t<std::common_type_t, typename FieldTraits< std::decay_t<Args> >::real_type...>;
 
    // make sure that we have an std::common_type: using an assertion produces a more readable error message
    // than a compiler template instantiation error
    static_assert ( sizeof...(Args) == 0 or
                    not (std::is_same_v<field_type, Std::nonesuch> or std::is_same_v<real_type, Std::nonesuch>),
        "No std::common_type implemented for the given field_type/real_type of the Args. Please provide an implementation and check that a FieldTraits class is present for your type.");
 
 
    static constexpr size_type size()
    {
      return sizeof...(Args);
    }
 
    static constexpr size_type N()
    {
      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[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
    {
      return std::get<index>(*this);
    }
 
    template< size_type index >
    const typename std::tuple_element<index,TupleType>::type&
    operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) const
    {
      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, real_type(0), [&](auto&& a, auto&& entry) {
        return a + entry.one_norm();
      });
    }
 
    auto one_norm_real() const {
      using namespace Dune::Hybrid;
      return accumulate(*this, real_type(0), [&](auto&& a, auto&& entry) {
        return a + entry.one_norm_real();
      });
    }
 
    real_type two_norm2() const {
      using namespace Dune::Hybrid;
      return accumulate(*this, real_type(0), [&](auto&& a, auto&& entry) {
        return a + entry.two_norm2();
      });
    }
 
    real_type two_norm() const {return sqrt(this->two_norm2());}
 
    real_type infinity_norm() const
    {
      using namespace Dune::Hybrid;
      using std::max;
 
      real_type result = 0.0;
      // Compute max norm tracking appearing nan values
      // if the field type supports nan.
      if constexpr (HasNaN<field_type>()) {
        // 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);
      } else {
        using namespace Dune::Hybrid; // needed for icc, see issue #31
        forEach(*this, [&](auto&& entry) {
          result = max(entry.infinity_norm(), result);
        });
      }
      return result;
    }
 
    real_type infinity_norm_real() const
    {
      using namespace Dune::Hybrid;
      using std::max;
 
      real_type result = 0.0;
      // Compute max norm tracking appearing nan values
      // if the field type supports nan.
      if constexpr (HasNaN<field_type>()) {
        // 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);
      } else {
        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 Dune
 
namespace std
{
  template <size_t i, typename... Args>
  struct tuple_element<i,Dune::MultiTypeBlockVector<Args...> >
  {
    using type = typename std::tuple_element<i, std::tuple<Args...> >::type;
  };
 
  template <typename... Args>
  struct tuple_size<Dune::MultiTypeBlockVector<Args...> >
    : std::integral_constant<std::size_t, sizeof...(Args)>
  {};
}
 
#endif