#pragma once #include #include #include #include #include #include namespace std { template struct common_type, T> { using type = T; }; template struct common_type, T> { using type = T; }; } namespace Dune { namespace MatVec { /// Traits to detect fixed size containers like FieldVector and FieldMatrix /// @{ template struct IsMatrix : std::false_type {}; template struct IsMatrix> : std::true_type {}; template struct IsMatrix> : std::true_type {}; template struct IsVector : std::false_type {}; template struct IsVector> : std::true_type {}; template struct IsMatVec : std::integral_constant::value || IsVector::value> {}; /// @} /// Convert the field_Type of container to type S /// @{ template struct MakeMatVec { using type = A; }; template struct MakeMatVec,S> { using type = FieldMatrix; }; template struct MakeMatVec,S> { using type = DiagonalMatrix; }; template struct MakeMatVec,S> { using type = FieldVector; }; /// @} /// Convert pseudo-scalar to real scalar type /// @{ template decltype(auto) as_scalar(T&& t) { return FWD(t); } template T as_scalar(FieldVector const& t) { return t[0]; } template T as_scalar(FieldMatrix const& t) { return t[0][0]; } template T as_scalar(DiagonalMatrix const& t) { return t.diagonal(0); } /// @} /// Convert pseudo-vector to real vector type /// @{ template decltype(auto) as_vector(T&& t) { return FWD(t); } template FieldVector const& as_vector(FieldMatrix const& t) { return t[0]; } template FieldVector& as_vector(FieldMatrix& t) { return t[0]; } /// @} /// Convert pseudo-matrix to real matrix type with proper operator[][] /// @{ template decltype(auto) as_matrix(T&& t) { return FWD(t); } template class MatrixView; template MatrixView> as_matrix(DiagonalMatrix const& mat) { return {mat}; } /// @} // returns -a template auto negate(A const& a); // returns a+b template auto plus(A const& a, B const& b); // returns a-b template auto minus(A const& a, B const& b); // returns a*b template ::value && IsNumber::value, int> = 0> auto multiplies(A const& a, B const& b); template ::value != IsNumber::value, int> = 0> auto multiplies(A const& a, B const& b); template auto multiplies(FieldVector const& a, FieldVector const& b); template ::value && IsVector::value, int> = 0> auto multiplies(Mat const& mat, Vec const& vec); template ::value && IsMatrix::value, int> = 0> auto multiplies(Vec const& vec, Mat const& mat); template auto multiplies(FieldMatrix const& a, FieldMatrix const& b); // return a/b template auto divides(A a, B const& b) { return a /= b; } } // end namespace MatVec // some arithmetic operations with FieldVector and FieldMatrix template ::value, int> = 0> auto operator-(A const& a) { return MatVec::negate(MatVec::as_scalar(a)); } template ::value || MatVec::IsMatVec::value, int> = 0> auto operator+(A const& a, B const& b) { return MatVec::plus(MatVec::as_scalar(a), MatVec::as_scalar(b)); } template ::value || MatVec::IsMatVec::value, int> = 0> auto operator-(A const& a, B const& b) { return MatVec::minus(MatVec::as_scalar(a), MatVec::as_scalar(b)); } template ::value || MatVec::IsMatVec::value, int> = 0> auto operator*(A const& a, B const& b) { return MatVec::multiplies(MatVec::as_scalar(a), MatVec::as_scalar(b)); } template ::value || MatVec::IsMatVec::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 FieldVector cross(FieldVector const& a); /// Cross-product of two 3d-vectors template FieldVector cross(FieldVector const& a, FieldVector const& b); /// Dot product (vec1^T * vec2) template auto dot(FieldVector const& vec1, FieldVector const& vec2); template auto dot(FieldMatrix const& vec1, FieldMatrix const& vec2); // ---------------------------------------------------------------------------- /// Sum of vector entires. template T sum(FieldVector const& x); template T sum(FieldMatrix const& x); /// Dot-product with the vector itself template ::value, int> = 0> auto unary_dot(T const& x); template auto unary_dot(FieldVector const& x); template auto unary_dot(FieldMatrix const& x); /// Maximum over all vector entries template auto max(FieldVector const& x); template auto max(FieldMatrix const& x); /// Minimum over all vector entries template auto min(FieldVector const& x); template auto min(FieldMatrix const& x); /// Maximum of the absolute values of vector entries template auto abs_max(FieldVector const& x); template auto abs_max(FieldMatrix const& x); /// Minimum of the absolute values of vector entries template auto abs_min(FieldVector const& x); template auto abs_min(FieldMatrix const& x); // ---------------------------------------------------------------------------- /** \ingroup vector_norms * \brief The 1-norm of a vector = sum_i |x_i| **/ template auto one_norm(FieldVector const& x); template auto one_norm(FieldMatrix const& x); /** \ingroup vector_norms * \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 ) **/ template ::value, int> = 0> auto two_norm(T const& x); template auto two_norm(FieldVector const& x); template auto two_norm(FieldMatrix const& x); /** \ingroup vector_norms * \brief The p-norm of a vector = ( sum_i |x_i|^p )^(1/p) **/ template auto p_norm(FieldVector const& x); template auto p_norm(FieldMatrix const& x); /** \ingroup vector_norms * \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max **/ template auto infty_norm(FieldVector const& x); template auto infty_norm(FieldMatrix const& x); // ---------------------------------------------------------------------------- /// The euklidean distance between two vectors = |lhs-rhs|_2 template ::value, int> = 0> T distance(T const& lhs, T const& rhs); template T distance(FieldVector const& lhs, FieldVector const& rhs); // ---------------------------------------------------------------------------- /// Outer product (vec1 * vec2^T) template auto outer(FieldMatrix const& vec1, FieldMatrix const& vec2); // ---------------------------------------------------------------------------- template T det(FieldMatrix const& /*mat*/); /// Determinant of a 1x1 matrix template T det(FieldMatrix const& mat); /// Determinant of a 2x2 matrix template T det(FieldMatrix const& mat); /// Determinant of a 3x3 matrix template T det(FieldMatrix const& mat); /// Determinant of a NxN matrix template T det(FieldMatrix const& mat); /// Return the inverse of the matrix `mat` template auto inv(FieldMatrix mat); /// Solve the linear system A*x = b template void solve(FieldMatrix const& A, FieldVector& x, FieldVector const& b); /// Gramian determinant = sqrt( det( DT^T * DF ) ) template T gramian(FieldMatrix const& DF); /// Gramian determinant, specialization for 1 column matrices template T gramian(FieldMatrix const& DF); // ---------------------------------------------------------------------------- // some arithmetic operations with FieldMatrix template FieldMatrix trans(FieldMatrix const& A); template DiagonalMatrix const& trans(DiagonalMatrix const& A) { return A; } // ----------------------------------------------------------------------------- template FieldMatrix,M,N> multiplies_AtB(FieldMatrix const& A, FieldMatrix const& B); template FieldMatrix,M,N> multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B); template FieldMatrix& multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B, FieldMatrix& C); template FieldMatrix& multiplies_ABt(FieldMatrix const& A, DiagonalMatrix const& B, FieldMatrix& C); template FieldVector& multiplies_ABt(FieldMatrix const& A, DiagonalMatrix const& B, FieldVector& C); // ----------------------------------------------------------------------------- template T const& at(FieldMatrix const& vec, std::size_t i); template T const& at(FieldMatrix const& vec, std::size_t i); // necessary specialization to resolve ambiguous calls template T const& at(FieldMatrix const& vec, std::size_t i); template T const& at(FieldVector const& vec, std::size_t i); } // end namespace Dune namespace AMDiS { using Dune::FieldMatrix; using Dune::FieldVector; } #include "FieldMatVec.inc.hpp"