FieldMatVec.hpp

#pragma once

#include <type_traits>

#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>

#include <amdis/common/TypeTraits.hpp>

namespace std
{
  template <class T, int N>
  struct common_type<Dune::FieldVector<T,N>, T>
  {
    using type = T;
  };

  template <class T, int N, int M>
  struct common_type<Dune::FieldMatrix<T,N,M>, T>
  {
    using type = T;
  };
}

namespace Dune
{
  namespace MatVec
  {
    /// Traits to detect fixed size containers like FieldVector and FieldMatrix
    /// @{
    template <class T>
    struct IsMatrix : std::false_type {};

    template <class T, int M, int N>
    struct IsMatrix<FieldMatrix<T,M,N>> : std::true_type {};

    template <class T, int N>
    struct IsMatrix<DiagonalMatrix<T,N>> : std::true_type {};


    template <class T>
    struct IsVector : std::false_type {};

    template <class T, int N>
    struct IsVector<FieldVector<T,N>> : std::true_type {};


    template <class T>
    struct IsMatVec
      : std::integral_constant<bool, IsMatrix<T>::value || IsVector<T>::value> {};
    /// @}

    /// Convert the field_Type of container to type S
    /// @{
    template <class A, class S>
    struct MakeMatVec
    {
      using type = A;
    };

    template <class T, int M, int N, class S>
    struct MakeMatVec<FieldMatrix<T,M,N>,S>
    {
      using type = FieldMatrix<S,M,N>;
    };

    template <class T, int N, class S>
    struct MakeMatVec<DiagonalMatrix<T,N>,S>
    {
      using type = DiagonalMatrix<S,N>;
    };

    template <class T, int N, class S>
    struct MakeMatVec<FieldVector<T,N>,S>
    {
      using type = FieldVector<S,N>;
    };
    /// @}

    /// Convert pseudo-scalar to real scalar type
    /// @{
    template <class T>
    decltype(auto) as_scalar(T&& t)
    {
      return FWD(t);
    }

    template <class T>
    T as_scalar(FieldVector<T,1> const& t)
    {
      return t[0];
    }

    template <class T>
    T as_scalar(FieldMatrix<T,1,1> const& t)
    {
      return t[0][0];
    }

    template <class T>
    T as_scalar(DiagonalMatrix<T,1> const& t)
    {
      return t.diagonal(0);
    }
    /// @}

    /// Convert pseudo-vector to real vector type
    /// @{
    template <class T>
    decltype(auto) as_vector(T&& t)
    {
      return FWD(t);
    }

    template <class T, int N>
    FieldVector<T,N> const& as_vector(FieldMatrix<T,1,N> const& t)
    {
      return t[0];
    }

    template <class T, int N>
    FieldVector<T,N>& as_vector(FieldMatrix<T,1,N>& t)
    {
      return t[0];
    }
    /// @}


    /// Convert pseudo-matrix to real matrix type with proper operator[][]
    /// @{
    template <class T>
    decltype(auto) as_matrix(T&& t)
    {
      return FWD(t);
    }

    template <class Mat>
    class MatrixView;

    template <class T, int N>
    MatrixView<DiagonalMatrix<T,N>> as_matrix(DiagonalMatrix<T,N> const& mat)
    {
      return {mat};
    }
    /// @}

    // returns -a
    template <class A>
    auto negate(A const& a);

    // returns a+b
    template <class A, class B>
    auto plus(A const& a, B const& b);

    // returns a-b
    template <class A, class B>
    auto minus(A const& a, B const& b);

    // returns a*b
    template <class A, class B,
      std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int> = 0>
    auto multiplies(A const& a, B const& b);

    template <class A, class B,
      std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int> = 0>
    auto multiplies(A const& a, B const& b);

    template <class T, int N, class S>
    auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b);

    template <class Mat, class Vec,
      std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int> = 0>
    auto multiplies(Mat const& mat, Vec const& vec);

    template <class Vec, class Mat,
      std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int> = 0>
    auto multiplies(Vec const& vec, Mat const& mat);

    template <class T, int L, int M, int N, class S>
    auto multiplies(FieldMatrix<T,L,M> const& a, FieldMatrix<S,M,N> const& b);

    // return a/b
    template <class A, class B>
    auto divides(A a, B const& b)
    {
      return a /= b;
    }

  } // end namespace MatVec

  // some arithmetic operations with FieldVector and FieldMatrix

  template <class A,
    std::enable_if_t<MatVec::IsMatVec<A>::value, int> = 0>
  auto operator-(A const& a)
  {
    return MatVec::negate(MatVec::as_scalar(a));
  }

  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator+(A const& a, B const& b)
  {
    return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b));
  }

  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator-(A const& a, B const& b)
  {
    return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b));
  }

  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator*(A const& a, B const& b)
  {
    return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b));
  }

  template <class A, class B,
    std::enable_if_t<MatVec::IsMatVec<A>::value || MatVec::IsMatVec<B>::value, int> = 0>
  auto operator/(A const& a, B const& b)
  {
    return MatVec::divides(MatVec::as_scalar(a), MatVec::as_scalar(b));
  }

  // ----------------------------------------------------------------------------

  /// Cross-product a 2d-vector = orthogonal vector
  template <class T>
  FieldVector<T, 2> cross(FieldVector<T, 2> const& a);

  /// Cross-product of two 3d-vectors
  template <class T>
  FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b);

  /// Dot product (vec1^T * vec2)
  template <class T, class S, int N>
  auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2);

  template <class T, class S, int N>
  auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2);

  // ----------------------------------------------------------------------------

  /// Sum of vector entires.
  template <class T, int N>
  T sum(FieldVector<T, N> const& x);

  template <class T, int N>
  T sum(FieldMatrix<T, 1, N> const& x);


  /// Dot-product with the vector itself
  template <class T, int N>
  auto unary_dot(FieldVector<T, N> const& x);

  template <class T, int N>
  auto unary_dot(FieldMatrix<T, 1, N> const& x);

  /// Maximum over all vector entries
  template <class T, int N>
  auto max(FieldVector<T, N> const& x);

  template <class T, int N>
  auto max(FieldMatrix<T, 1, N> const& x);

  /// Minimum over all vector entries
  template <class T, int N>
  auto min(FieldVector<T, N> const& x);

  template <class T, int N>
  auto min(FieldMatrix<T, 1, N> const& x);

  /// Maximum of the absolute values of vector entries
  template <class T, int N>
  auto abs_max(FieldVector<T, N> const& x);

  template <class T, int N>
  auto abs_max(FieldMatrix<T, 1, N> const& x);

  /// Minimum of the absolute values of vector entries
  template <class T, int N>
  auto abs_min(FieldVector<T, N> const& x);

  template <class T, int N>
  auto abs_min(FieldMatrix<T, 1, N> const& x);

  // ----------------------------------------------------------------------------

  /** \ingroup vector_norms
   *  \brief The 1-norm of a vector = sum_i |x_i|
   **/
  template <class T, int N>
  auto one_norm(FieldVector<T, N> const& x);

  template <class T, int N>
  auto one_norm(FieldMatrix<T, 1, N> const& x);

  /** \ingroup vector_norms
   *  \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 )
   **/
  template <class T, int N>
  auto two_norm(FieldVector<T, N> const& x);

  template <class T, int N>
  auto two_norm(FieldMatrix<T, 1, N> const& x);

  /** \ingroup vector_norms
   *  \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p)
   **/
  template <int p, class T, int N>
  auto p_norm(FieldVector<T, N> const& x);

  template <int p, class T, int N>
  auto p_norm(FieldMatrix<T, 1, N> const& x);

  /** \ingroup vector_norms
   *  \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max
   **/
  template <class T, int N>
  auto infty_norm(FieldVector<T, N> const& x);

  template <class T, int N>
  auto infty_norm(FieldMatrix<T, 1, N> const& x);

  // ----------------------------------------------------------------------------

  /// The euklidean distance between two vectors = |lhs-rhs|_2
  template <class T, int N>
  T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs);

  // ----------------------------------------------------------------------------

  /// Outer product (vec1 * vec2^T)
  template <class T, class S, int N, int M, int K>
  auto outer(FieldMatrix<T,N,K> const& vec1, FieldMatrix<S,M,K> const& vec2);

  // ----------------------------------------------------------------------------

  template <class T>
  T det(FieldMatrix<T, 0, 0> const& /*mat*/);

  /// Determinant of a 1x1 matrix
  template <class T>
  T det(FieldMatrix<T, 1, 1> const& mat);

  /// Determinant of a 2x2 matrix
  template <class T>
  T det(FieldMatrix<T, 2, 2> const& mat);

  /// Determinant of a 3x3 matrix
  template <class T>
  T det(FieldMatrix<T, 3, 3> const& mat);

  /// Determinant of a NxN matrix
  template <class T,  int N>
  T det(FieldMatrix<T, N, N> const& mat);


  /// Return the inverse of the matrix `mat`
  template <class T, int N>
  auto inv(FieldMatrix<T, N, N> mat);

  /// Solve the linear system A*x = b
  template <class T, int N>
  void solve(FieldMatrix<T, N, N> const& A,  FieldVector<T, N>& x,  FieldVector<T, N> const& b);


  /// Gramian determinant = sqrt( det( DT^T * DF ) )
  template <class T, int N, int M>
  T gramian(FieldMatrix<T,N,M> const& DF);

  /// Gramian determinant, specialization for 1 column matrices
  template <class T, int M>
  T gramian(FieldMatrix<T, 1, M> const& DF);

  // ----------------------------------------------------------------------------
  // some arithmetic operations with FieldMatrix

  template <class T, int M, int N>
  FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);

  template <class T, int N>
  DiagonalMatrix<T,N> const& trans(DiagonalMatrix<T,N> const& A)
  {
    return A;
  }

  // -----------------------------------------------------------------------------

  template <class T1, class T2, int M, int N, int L>
  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_AtB(FieldMatrix<T1, L, M> const& A,  FieldMatrix<T2, N, L> const& B);

  template <class T1, class T2, int M, int N, int L>
  FieldMatrix<std::common_type_t<T1,T2>,M,N> multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B);

  template <class T1, class T2, class T3, int M, int N, int L>
  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,M,N>& C);

  template <class T1, class T2, class T3, int M, int N>
  FieldMatrix<T3,M,N>& multiplies_ABt(FieldMatrix<T1, M, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,M,N>& C);

  template <class T1, class T2, class T3, int N>
  FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C);

  // -----------------------------------------------------------------------------

  template <class T, int N>
  T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i);

  template <class T, int M>
  T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i);

  // necessary specialization to resolve ambiguous calls
  template <class T>
  T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i);

  template <class T, int N>
  T const& at(FieldVector<T,N> const& vec, std::size_t i);

} // end namespace Dune

namespace AMDiS
{
  using Dune::FieldMatrix;
  using Dune::FieldVector;
}

#include "FieldMatVec.inc.hpp"