fvector.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_COMMON_FVECTOR_HH
#define DUNE_COMMON_FVECTOR_HH

#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <cstdlib>
#include <cstring>
#include <type_traits>
#include <utility>
#include <initializer_list>

#include <dune/common/boundschecking.hh>
#include <dune/common/densevector.hh>
#include <dune/common/filledarray.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/math.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/typeutilities.hh>
#include <dune/common/concepts/number.hh>
#include <dune/common/std/algorithm.hh>
#include <dune/common/std/compare.hh>

namespace Dune {

  template< class K, int SIZE > class FieldVector;
  template< class K, int SIZE >
  struct DenseMatVecTraits< FieldVector<K,SIZE> >
  {
    typedef FieldVector<K,SIZE> derived_type;
    typedef std::array<K,SIZE> container_type;
    typedef K value_type;
    typedef typename container_type::size_type size_type;
  };

  template< class K, int SIZE >
  struct FieldTraits< FieldVector<K,SIZE> >
  {
    typedef typename FieldTraits<K>::field_type field_type;
    typedef typename FieldTraits<K>::real_type real_type;
  };

  template<typename C, int SIZE>
  struct IsFieldVectorSizeCorrect
  {
    constexpr static bool value = true;
  };

  template<typename T, int SIZE>
  struct IsFieldVectorSizeCorrect<FieldVector<T,SIZE>,SIZE>
  {
    constexpr static bool value = true;
  };

  template<typename T, int SIZE, int SIZE1>
  struct IsFieldVectorSizeCorrect<FieldVector<T,SIZE1>,SIZE>
  {
    constexpr static bool value = false;
  };


  template< class K, int SIZE >
  class FieldVector :
    public DenseVector< FieldVector<K,SIZE> >
  {
    using Base = DenseVector< FieldVector<K,SIZE> >;

    std::array<K,SIZE> _data;

  public:

    static constexpr int dimension = SIZE;

    using size_type = typename Base::size_type;

    using value_type = typename Base::value_type;

    using reference = value_type&;

    using const_reference = const value_type&;

  public:

    constexpr FieldVector ()
        noexcept(std::is_nothrow_default_constructible_v<K>)
      : _data{}
    {}

    explicit(SIZE != 1)
    constexpr FieldVector (const value_type& value) noexcept
      : _data{filledArray<SIZE>(value)}
    {}

    template<Concept::Number S>
      requires (std::constructible_from<K,S>)
    explicit(SIZE != 1)
    constexpr FieldVector (const S& scalar)
        noexcept(std::is_nothrow_constructible_v<K,S>)
      : _data{filledArray<SIZE,K>(K(scalar))}
    {}

    constexpr FieldVector (const std::initializer_list<K>& l)
      : _data{}
    {
      assert(l.size() == size());
      for (size_type i = 0; i < size(); ++i)
        _data[i] = std::data(l)[i];
    }

    template<class V>
      requires (IsFieldVectorSizeCorrect<V,SIZE>::value &&
        std::is_assignable_v<K&, decltype(std::declval<const V&>()[0])>)
    constexpr FieldVector (const DenseVector<V>& x)
    {
      assert(x.size() == size());
      for (size_type i = 0; i < size(); ++i)
        _data[i] = x[i];
    }

    template<class OtherK>
      requires (std::is_assignable_v<K&, const OtherK&>)
    explicit constexpr FieldVector (const FieldVector<OtherK, SIZE>& x)
        noexcept(std::is_nothrow_assignable_v<K&, const OtherK&>)
    {
      for (size_type i = 0; i < size(); ++i)
        _data[i] = x[i];
    }

    constexpr FieldVector (const FieldVector&) = default;


    template<class V>
      requires (IsFieldVectorSizeCorrect<V,SIZE>::value &&
        std::is_assignable_v<K&, decltype(std::declval<const V&>()[0])>)
    constexpr FieldVector& operator= (const DenseVector<V>& x)
    {
      assert(x.size() == size());
      for (size_type i = 0; i < size(); ++i)
        _data[i] = x[i];
      return *this;
    }

    template<Concept::Number S>
      requires std::constructible_from<K,S>
    constexpr FieldVector& operator= (const S& scalar)
        noexcept(std::is_nothrow_constructible_v<K,S>)
    {
      _data.fill(K(scalar));
      return *this;
    }

    template<class OtherK>
      requires (std::is_assignable_v<K&, const OtherK&>)
    constexpr FieldVector& operator= (const FieldVector<OtherK, SIZE>& x)
        noexcept(std::is_nothrow_assignable_v<K&, const OtherK&>)
    {
      for (size_type i = 0; i < size(); ++i)
        _data[i] = x[i];
      return *this;
    }

    constexpr FieldVector& operator= (const FieldVector&) = default;



    static constexpr size_type size () noexcept { return SIZE; }




    constexpr reference operator[] (size_type i)
    {
      DUNE_ASSERT_BOUNDS(i < size());
      return _data[i];
    }

    constexpr const_reference operator[] (size_type i) const
    {
      DUNE_ASSERT_BOUNDS(i < size());
      return _data[i];
    }

    constexpr K* data () noexcept
    {
      return _data.data();
    }

    constexpr const K* data () const noexcept
    {
      return _data.data();
    }

    constexpr operator const_reference () const noexcept
        requires(SIZE == 1)
    {
      return _data[0];
    }

    constexpr operator reference () noexcept
        requires(SIZE == 1)
    {
      return _data[0];
    }




    template<Concept::Number S>
    friend constexpr bool operator== (const FieldVector& a, const S& b) noexcept
        requires(SIZE == 1)
    {
      return a._data[0] == b;
    }

    template<Concept::Number S>
    friend constexpr bool operator== (const S& a, const FieldVector& b) noexcept
        requires(SIZE == 1)
    {
      return a == b._data[0];
    }

    template<class T>
      requires (Std::three_way_comparable_with<K,T>)
    friend constexpr auto operator<=> (const FieldVector& a, const FieldVector<T,SIZE>& b) noexcept
    {
#if __cpp_lib_three_way_comparison
      return a._data <=> b._data;
#else
      return Std::lexicographical_compare_three_way(a.begin(), a.end(), b.begin(), b.end());
#endif
    }

    template<Concept::Number S>
    friend constexpr auto operator<=> (const FieldVector& a, const S& b) noexcept
        requires(SIZE == 1)
    {
      return a._data[0] <=> b;
    }

    template<Concept::Number S>
    friend constexpr auto operator<=> (const S& a, const FieldVector& b) noexcept
        requires(SIZE == 1)
    {
      return a <=> b._data[0];
    }




    template<Concept::Number S>
    friend constexpr auto operator* (const FieldVector& a, const S& b) noexcept
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      FieldVector<ResultValueType,dimension> result;
      for (size_type i = 0; i < size(); ++i)
        result[i] = a[i] * b;
      return result;
    }

    template<Concept::Number S>
    friend constexpr auto operator* (const S& a, const FieldVector& b) noexcept
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      FieldVector<ResultValueType,dimension> result;
      for (size_type i = 0; i < size(); ++i)
        result[i] = a * b[i];
      return result;
    }

    template<Concept::Number S>
    friend constexpr auto operator/ (const FieldVector& a, const S& b) noexcept
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      FieldVector<ResultValueType,dimension> result;
      for (size_type i = 0; i < size(); ++i)
        result[i] = a[i] / b;
      return result;
    }

    template<Concept::Number S>
    friend constexpr FieldVector operator/ (const S& a, const FieldVector& b) noexcept
        requires(SIZE == 1)
    {
      return FieldVector{a / b[0]};
    }

    template<Concept::Number S>
    friend constexpr auto operator+ (const FieldVector& a, const S& b) noexcept
        requires(SIZE == 1)
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      return FieldVector<ResultValueType,dimension>{a[0] + b};
    }

    template<Concept::Number S>
    friend constexpr auto operator+ (const S& a, const FieldVector& b) noexcept
        requires(SIZE == 1)
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      return FieldVector<ResultValueType,dimension>{a + b[0]};
    }

    template<Concept::Number S>
    friend constexpr auto operator- (const FieldVector& a, const S& b) noexcept
        requires(SIZE == 1)
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      return FieldVector<ResultValueType,dimension>{a[0] - b};
    }

    template<Concept::Number S>
    friend constexpr auto operator- (const S& a, const FieldVector& b) noexcept
        requires(SIZE == 1)
    {
      using ResultValueType = typename PromotionTraits<K,S>::PromotedType;
      return FieldVector<ResultValueType,dimension>{a - b[0]};
    }

  };

  template<class K, int SIZE>
  std::istream& operator>> (std::istream& in, FieldVector<K, SIZE>& v)
  {
    FieldVector<K, SIZE> w;
    for (int i = 0; i < SIZE; ++i)
      in >> w[i];
    if (in)
      v = w;
    return in;
  }

  /* Overloads for common classification functions */
  namespace MathOverloads {

    template<class K, int SIZE>
    auto isFinite (const FieldVector<K,SIZE>& b, PriorityTag<2>, ADLTag)
    {
      bool out = true;
      for (int i = 0; i < SIZE; ++i) {
        out &= Dune::isFinite(b[i]);
      }
      return out;
    }

    template<class K, int SIZE>
    bool isInf (const FieldVector<K,SIZE>& b, PriorityTag<2>, ADLTag)
    {
      bool out = false;
      for (int i = 0; i < SIZE; ++i) {
        out |= Dune::isInf(b[i]);
      }
      return out;
    }

    template<class K, int SIZE,
      std::enable_if_t<HasNaN<K>::value, int> = 0>
    bool isNaN (const FieldVector<K,SIZE>& b, PriorityTag<2>, ADLTag)
    {
      bool out = false;
      for (int i = 0; i < SIZE; ++i) {
        out |= Dune::isNaN(b[i]);
      }
      return out;
    }

    template<class K,
      std::enable_if_t<HasNaN<K>::value, int> = 0>
    bool isUnordered (const FieldVector<K,1>& b, const FieldVector<K,1>& c,
                      PriorityTag<2>, ADLTag)
    {
      return Dune::isUnordered(b[0],c[0]);
    }

  } // end namespace MathOverloads

} // end namespace Dune

#endif // DUNE_COMMON_FVECTOR_HH