FieldMatVec.inc.hpp

#pragma once

#include <algorithm>
#include <limits>

#include <dune/functions/functionspacebases/flatvectorview.hh>

#include <amdis/common/Math.hpp>
#include <amdis/operations/Arithmetic.hpp>
#include <amdis/operations/MaxMin.hpp>

#ifndef DOXYGEN

namespace Dune {

namespace MatVec {

  template <class A>
  auto negate(A const& a)
  {
    return multiplies(a, -1);
  }

  // returns a+b
  template <class A, class B>
  auto plus(A const& a, B const& b)
  {
    if constexpr(IsNumber<A>::value && IsNumber<B>::value)
      return a + b;
    else {
      using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
      typename MakeMatVec<A,T>::type c(a);

      auto b_ = Dune::Functions::flatVectorView(b);
      auto c_ = Dune::Functions::flatVectorView(c);
      assert(int(b_.size()) == int(c_.size()));
      for(int i = 0; i < int(b_.size()); ++i)
        c_[i] += b_[i];

      return c;
    }
  }

  // returns a-b
  template <class A, class B>
  auto minus(A const& a, B const& b)
  {
    if constexpr(IsNumber<A>::value && IsNumber<B>::value)
      return a - b;
    else {
      using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
      typename MakeMatVec<A,T>::type c(a);

      auto b_ = Dune::Functions::flatVectorView(b);
      auto c_ = Dune::Functions::flatVectorView(c);
      assert(int(b_.size()) == int(c_.size()));
      for(int i = 0; i < int(b_.size()); ++i)
        c_[i] -= b_[i];

      return c;
    }
  }

  template <class A, class B,
    std::enable_if_t<IsNumber<A>::value && IsNumber<B>::value, int>>
  auto multiplies(A const& a, B const& b)
  {
    return a * b;
  }

  template <class A, class B,
    std::enable_if_t<IsNumber<A>::value != IsNumber<B>::value, int>>
  auto multiplies(A const& a, B const& b)
  {
    using T = std::common_type_t<typename FieldTraits<A>::field_type, typename FieldTraits<B>::field_type>;
    if constexpr(IsNumber<A>::value) {
      typename MakeMatVec<B,T>::type b_(b);
      return b_ *= a;
    } else {
      typename MakeMatVec<A,T>::type a_(a);
      return a_ *= b;
    }
  }

  template <class T, int N, class S>
  auto multiplies(FieldVector<T,N> const& a, FieldVector<S,N> const& b)
  {
    return a.dot(b);
  }


  template <class Mat, class Vec,
    std::enable_if_t<IsMatrix<Mat>::value && IsVector<Vec>::value, int>>
  auto multiplies(Mat const& mat, Vec const& vec)
  {
    static_assert(int(Mat::cols) == int(Vec::dimension), "");
    using T = std::common_type_t<typename FieldTraits<Vec>::field_type, typename FieldTraits<Mat>::field_type>;
    FieldVector<T,Mat::rows> y;
    mat.mv(vec, y);
    return y;
  }


  template <class Vec, class Mat,
    std::enable_if_t<IsVector<Vec>::value && IsMatrix<Mat>::value, int>>
  auto multiplies(Vec const& vec, Mat const& mat)
  {
    static_assert(int(Mat::rows) == int(Vec::dimension), "");
    using T = std::common_type_t<typename FieldTraits<Vec>::field_type, typename FieldTraits<Mat>::field_type>;
    FieldVector<T,Mat::cols> y;
    mat.mtv(vec, y);
    return y;
  }


  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)
  {
    FieldMatrix<std::common_type_t<T,S>,L,N> C;

    for (int i = 0; i < L; ++i) {
      for (int j = 0; j < N; ++j) {
        C[i][j] = 0;
        for (int k = 0; k < M; ++k)
          C[i][j] += a[i][k]*b[k][j];
      }
    }
    return C;
  }

  template <class T, int SIZE>
  class MatrixView<DiagonalMatrix<T,SIZE>>
  {
    using Matrix = DiagonalMatrix<T,SIZE>;
    using size_type = typename Matrix::size_type;

    struct RowView
    {
      T operator[](size_type c) const
      {
        assert(0 <= c && c < mat_->M());
        return c == r_ ? mat_->diagonal(r_) : T(0);
      }

      Matrix const* mat_;
      size_type r_;
    };

  public:
    MatrixView(Matrix const& mat)
      : mat_(mat)
    {}

    RowView operator[](size_type r) const
    {
      assert(0 <= r && r < mat_.N());
      return {&mat_,r};
    }

    size_type N() const
    {
      return mat_.N();
    }

    size_type M() const
    {
      return mat_.M();
    }

  private:
    Matrix const& mat_;
  };

} // end namespace MatVec

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

template <class T>
FieldVector<T, 2> cross(FieldVector<T, 2> const& a)
{
  return {{ a[1], -a[0] }};
}

template <class T>
FieldVector<T, 3> cross(FieldVector<T, 3> const& a, FieldVector<T, 3> const& b)
{
  return {{ a[1]*b[2] - a[2]*b[1],
            a[2]*b[0] - a[0]*b[2],
            a[0]*b[1] - a[1]*b[0] }};
}

template <class T, class S, int N>
auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2)
{
  return vec1.dot(vec2);
}

template <class T, class S, int N>
auto dot(FieldMatrix<T,1,N> const& vec1, FieldMatrix<S,1,N> const& vec2)
{
  return vec1[0].dot(vec2[0]);
}

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

namespace Impl
{
  template <class T, int N, class Operation>
  T accumulate(FieldVector<T, N> const& x, T init, Operation op)
  {
    for (int i = 0; i < N; ++i)
      init = op(init, x[i]);
    return init;
  }

  template <class T, int N, class Operation>
  T accumulate(FieldMatrix<T, 1, N> const& x, T init, Operation op)
  {
    for (int i = 0; i < N; ++i)
      init = op(init, x[0][i]);
    return init;
  }

} // end namespace Impl

template <class T, int N>
T sum(FieldVector<T, N> const& x)
{
  return Impl::accumulate(x, T(0), AMDiS::Operation::Plus{});
}

template <class T, int N>
T sum(FieldMatrix<T, 1, N> const& x)
{
  return Impl::accumulate(x, T(0), AMDiS::Operation::Plus{});
}


template <class T,
  std::enable_if_t<Dune::IsNumber<T>::value, int> >
auto unary_dot(T const& x)
{
  using std::abs;
  return AMDiS::Math::sqr(abs(x));
}

template <class T, int N>
auto unary_dot(FieldVector<T, N> const& x)
{
  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
  return Impl::accumulate(x, T(0), op);
}

template <class T, int N>
auto unary_dot(FieldMatrix<T, 1, N> const& x)
{
  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::sqr(abs(b)); };
  return Impl::accumulate(x, T(0), op);
}

template <class T, int N>
auto max(FieldVector<T, N> const& x)
{
  return Impl::accumulate(x, std::numeric_limits<T>::lowest(), AMDiS::Operation::Max{});
}

template <class T, int N>
auto max(FieldMatrix<T, 1, N> const& x)
{
  return Impl::accumulate(x, std::numeric_limits<T>::lowest(), AMDiS::Operation::Max{});
}

template <class T, int N>
auto min(FieldVector<T, N> const& x)
{
  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::Min{});
}

template <class T, int N>
auto min(FieldMatrix<T, 1, N> const& x)
{
  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::Min{});
}

template <class T, int N>
auto abs_max(FieldVector<T, N> const& x)
{
  return Impl::accumulate(x, T(0), AMDiS::Operation::AbsMax{});
}

template <class T, int N>
auto abs_max(FieldMatrix<T, 1, N> const& x)
{
  return Impl::accumulate(x, T(0), AMDiS::Operation::AbsMax{});
}

template <class T, int N>
auto abs_min(FieldVector<T, N> const& x)
{
  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::AbsMin{});
}

template <class T, int N>
auto abs_min(FieldMatrix<T, 1, N> const& x)
{
  return Impl::accumulate(x, std::numeric_limits<T>::max(), AMDiS::Operation::AbsMin{});
}

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

template <class T, int N>
auto one_norm(FieldVector<T, N> const& x)
{
  auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
  return Impl::accumulate(x, T(0), op);
}

template <class T, int N>
auto one_norm(FieldMatrix<T, 1, N> const& x)
{
  auto op = [](auto const& a, auto const& b) { using std::abs; return a + abs(b); };
  return Impl::accumulate(x, T(0), op);
}

template <class T,
  std::enable_if_t<Dune::IsNumber<T>::value, int> >
auto two_norm(T const& x)
{
  using std::abs;
  return abs(x);
}

template <class T, int N>
auto two_norm(FieldVector<T, N> const& x)
{
  using std::sqrt;
  return sqrt(unary_dot(x));
}

template <class T, int N>
auto two_norm(FieldMatrix<T, 1, N> const& x)
{
  using std::sqrt;
  return sqrt(unary_dot(x));
}

template <int p, class T, int N>
auto p_norm(FieldVector<T, N> const& x)
{
  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
  return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
}

template <int p, class T, int N>
auto p_norm(FieldMatrix<T, 1, N> const& x)
{
  auto op = [](auto const& a, auto const& b) { using std::abs; return a + AMDiS::Math::pow<p>(abs(b)); };
  return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
}

template <class T, int N>
auto infty_norm(FieldVector<T, N> const& x)
{
  return abs_max(x);
}

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

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

template <class T,
  std::enable_if_t<Dune::IsNumber<T>::value, int> >
T distance(T const& lhs, T const& rhs)
{
  using std::abs;
  return abs(lhs - rhs);
}

template <class T, int N>
T distance(FieldVector<T, N> const& lhs, FieldVector<T, N> const& rhs)
{
  using std::sqrt;
  T result = 0;
  for (int i = 0; i < N; ++i)
    result += AMDiS::Math::sqr(lhs[i] - rhs[i]);
  return sqrt(result);
}

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

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)
{
  using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
  result_type mat;
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < M; ++j)
      mat[i][j] = vec1[i].dot(vec2[j]);
  return mat;
}

template <class T, class S, int N, int M>
auto outer(FieldVector<T,N> const& vec1, FieldVector<S,M> const& vec2)
{
  using result_type = FieldMatrix<TYPEOF( std::declval<T>() * std::declval<S>() ), N, M>;
  result_type mat;
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < M; ++j)
      mat[i][j] = vec1[i] * vec2[j];
  return mat;
}

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

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

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

template <class T>
T det(FieldMatrix<T, 2, 2> const& mat)
{
  return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
}

template <class T>
T det(FieldMatrix<T, 3, 3> const& mat)
{
  return mat[0][0]*mat[1][1]*mat[2][2] + mat[0][1]*mat[1][2]*mat[2][0] + mat[0][2]*mat[1][0]*mat[2][1]
      - (mat[0][2]*mat[1][1]*mat[2][0] + mat[0][1]*mat[1][0]*mat[2][2] + mat[0][0]*mat[1][2]*mat[2][1]);
}

template <class T,  int N>
T det(FieldMatrix<T, N, N> const& mat)
{
  return mat.determinant();
}

template <class T, int N>
auto inv(FieldMatrix<T, N, N> mat)
{
  mat.invert();
  return mat;
}

template <class T, int N>
void solve(FieldMatrix<T, N, N> const& A,  FieldVector<T, N>& x,  FieldVector<T, N> const& b)
{
  A.solve(x, b);
}


template <class T, int N, int M>
T gramian(FieldMatrix<T,N,M> const& DF)
{
  using std::sqrt;
  return sqrt( det(outer(DF, DF)) );
}

template <class T, int M>
T gramian(FieldMatrix<T, 1, M> const& DF)
{
  using std::sqrt;
  return sqrt(dot(DF[0], DF[0]));
}

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

template <class T, int M, int N>
FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
{
  FieldMatrix<T,N,M> At;
  for (int i = 0; i < M; ++i)
    for (int j = 0; j < N; ++j)
      At[j][i] = A[i][j];

  return At;
}


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)
{
  FieldMatrix<std::common_type_t<T1,T2>,M,N> C;

  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      C[m][n] = 0;
      for (int l = 0; l < L; ++l)
        C[m][n] += A[l][m] * B[n][l];
    }
  }
  return C;
}

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)
{
  FieldMatrix<std::common_type_t<T1,T2>,M,N> C;
  return multiplies_ABt(A,B,C);
}

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)
{
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      C[m][n] = 0;
      for (int l = 0; l < L; ++l)
        C[m][n] += A[m][l] * B[n][l];
    }
  }
  return C;
}

template <class T1, class T2, class T3, int N, int L>
FieldVector<T3,N>& multiplies_ABt(FieldMatrix<T1, 1, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C)
{
  for (int n = 0; n < N; ++n) {
    C[n] = 0;
    for (int l = 0; l < L; ++l)
      C[n] += A[0][l] * B[n][l];
  }
  return C;
}

template <class T1, class T2, class T3, int N, int L>
FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldVector<T3,N>& C)
{
  for (int n = 0; n < N; ++n) {
    C[n] = 0;
    for (int l = 0; l < L; ++l)
      C[n] += A[l] * B[n][l];
  }
  return C;
}

template <class T1, class T2, class T3, int N, int L>
FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, L> const& A,  FieldMatrix<T2, N, L> const& B, FieldMatrix<T3,1,N>& C)
{
  for (int n = 0; n < N; ++n) {
    C[0][n] = 0;
    for (int l = 0; l < L; ++l)
      C[0][n] += A[l] * B[n][l];
  }
  return 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)
{
  for (int m = 0; m < M; ++m) {
    for (int n = 0; n < N; ++n) {
      C[m][n] = A[m][n] * B.diagonal(n);
    }
  }
  return 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)
{
  for (int n = 0; n < N; ++n) {
    C[n] = A[0][n] * B.diagonal(n);
  }
  return C;
}

template <class T1, class T2, class T3, int N>
FieldVector<T3,N>& multiplies_ABt(FieldVector<T1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldVector<T3,N>& C)
{
  for (int n = 0; n < N; ++n) {
    C[n] = A[n] * B.diagonal(n);
  }
  return C;
}

template <class T1, class T2, class T3, int N>
FieldMatrix<T3,1,N>& multiplies_ABt(FieldVector<T1, N> const& A,  DiagonalMatrix<T2, N> const& B, FieldMatrix<T3,1,N>& C)
{
  for (int n = 0; n < N; ++n) {
    C[0][n] = A[n] * B.diagonal(n);
  }
  return C;
}


template <class T, int N>
T const& at(FieldMatrix<T,N,1> const& vec, std::size_t i)
{
  return vec[i][0];
}

template <class T, int M>
T const& at(FieldMatrix<T,1,M> const& vec, std::size_t i)
{
  return vec[0][i];
}

template <class T>
T const& at(FieldMatrix<T,1,1> const& vec, std::size_t i)
{
  return vec[0][i];
}

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

} // end namespace AMDiS

#endif