#pragma once #include #include #include #include #include #include #include #include namespace AMDiS { using Dune::FieldVector; using Dune::FieldMatrix; // some arithmetic operations with FieldVector template ) > FieldVector operator*(FieldVector v, S factor) { return v *= factor; } template ) > FieldVector operator*(S factor, FieldVector v) { return v *= factor; } template ) > FieldVector operator/(FieldVector v, S factor) { return v /= factor; } template FieldVector operator*(FieldVector const& v, FieldVector const& w) { return {v[0] * w[0]}; } template FieldVector operator*(FieldVector const& factor, FieldVector v) { return v *= factor[0]; } template FieldVector operator*(FieldVector v, FieldVector const& factor) { return v *= factor[0]; } // ---------------------------------------------------------------------------- /// Cross-product a 2d-vector = orthogonal vector template FieldVector cross(FieldVector const& a) { return {{ a[1], -a[0] }}; } /// Cross-product of two vectors (in 3d only) template FieldVector cross(FieldVector const& a, FieldVector 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] }}; } /// Dot product (vec1^T * vec2) template auto dot(FieldVector const& vec1, FieldVector const& vec2) { return vec1.dot(vec2); } template auto operator*(FieldVector v, FieldVector const& w) { static_assert(M == N, "Requires vectors of the same type!"); return v.dot(w); } // ---------------------------------------------------------------------------- namespace Impl { template T accumulate(FieldVector const& x, Operation op) { T result = 0; for (int i = 0; i < N; ++i) result = op(result, x[i]); return result; } template T accumulate(FieldMatrix const& x, Operation op) { T result = 0; for (int i = 0; i < N; ++i) result = op(result, x[0][i]); return result; } } // end namespace Impl /// Sum of vector entires. template T sum(FieldVector const& x) { return Impl::accumulate(x, Operation::Plus{}); } template T sum(FieldMatrix const& x) { return Impl::accumulate(x, Operation::Plus{}); } /// Dot-product with the vector itself template auto unary_dot(FieldVector const& x) { auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); }; return Impl::accumulate(x, op); } template auto unary_dot(FieldMatrix const& x) { auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); }; return Impl::accumulate(x, op); } /// Maximum over all vector entries template auto max(FieldVector const& x) { return Impl::accumulate(x, Operation::Max{}); } template auto max(FieldMatrix const& x) { return Impl::accumulate(x, Operation::Max{}); } /// Minimum over all vector entries template auto min(FieldVector const& x) { return Impl::accumulate(x, Operation::Min{}); } template auto min(FieldMatrix const& x) { return Impl::accumulate(x, Operation::Min{}); } /// Maximum of the absolute values of vector entries template auto abs_max(FieldVector const& x) { return Impl::accumulate(x, Operation::AbsMax{}); } template auto abs_max(FieldMatrix const& x) { return Impl::accumulate(x, Operation::AbsMax{}); } /// Minimum of the absolute values of vector entries template auto abs_min(FieldVector const& x) { return Impl::accumulate(x, Operation::AbsMin{}); } template auto abs_min(FieldMatrix const& x) { return Impl::accumulate(x, Operation::AbsMin{}); } // ---------------------------------------------------------------------------- /** \ingroup vector_norms * \brief The 1-norm of a vector = sum_i |x_i| **/ template auto one_norm(FieldVector const& x) { auto op = [](auto const& a, auto const& b) { return a + std::abs(b); }; return Impl::accumulate(x, op); } template auto one_norm(FieldMatrix const& x) { auto op = [](auto const& a, auto const& b) { return a + std::abs(b); }; return Impl::accumulate(x, op); } /** \ingroup vector_norms * \brief The euklidean 2-norm of a vector = sqrt( sum_i |x_i|^2 ) **/ template auto two_norm(FieldVector const& x) { return std::sqrt(unary_dot(x)); } template auto two_norm(FieldMatrix const& x) { return std::sqrt(unary_dot(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) { auto op = [](auto const& a, auto const& b) { return a + Math::pow

(std::abs(b)); }; return std::pow( Impl::accumulate(x, op), 1.0/p ); } template auto p_norm(FieldMatrix const& x) { auto op = [](auto const& a, auto const& b) { return a + Math::pow

(std::abs(b)); }; return std::pow( Impl::accumulate(x, op), 1.0/p ); } /** \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) { return abs_max(x); } template auto infty_norm(FieldMatrix const& x) { return abs_max(x); } // ---------------------------------------------------------------------------- /// The euklidean distance between two vectors = |lhs-rhs|_2 template T distance(FieldVector const& lhs, FieldVector const& rhs) { T result = 0; for (int i = 0; i < N; ++i) result += Math::sqr(lhs[i] - rhs[i]); return std::sqrt(result); } // ---------------------------------------------------------------------------- /// Outer product (vec1 * vec2^T) template auto outer(FieldMatrix const& vec1, FieldMatrix const& vec2) { using result_type = FieldMatrix() * std::declval() ), 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 T det(FieldMatrix const& /*mat*/) { return 0; } /// Determinant of a 1x1 matrix template T det(FieldMatrix const& mat) { return mat[0][0]; } /// Determinant of a 2x2 matrix template T det(FieldMatrix const& mat) { return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0]; } /// Determinant of a 3x3 matrix template T det(FieldMatrix 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]); } /// Determinant of a NxN matrix template T det(FieldMatrix const& mat) { return mat.determinant(); } /// Return the inverse of the matrix mat template auto inv(FieldMatrix mat) { mat.invert(); return mat; } /// Solve the linear system A*x = b template void solve(FieldMatrix const& A, FieldVector& x, FieldVector const& b) { A.solve(x, b); } /// Gramian determinant = sqrt( det( DT^T * DF ) ) template T gramian(FieldMatrix const& DF) { using std::sqrt; return sqrt( det(outer(DF, DF)) ); } /// Gramian determinant, specialization for 1 column matrices template T gramian(FieldMatrix const& DF) { using std::sqrt; return sqrt(dot(DF[0], DF[0])); } // ---------------------------------------------------------------------------- // some arithmetic operations with FieldMatrix template FieldMatrix trans(FieldMatrix const& A) { FieldMatrix At; for (int i = 0; i < M; ++i) for (int j = 0; j < N; ++j) At[j][i] = A[i][j]; return At; } template FieldMatrix operator*(T scalar, FieldMatrix A) { return A *= scalar; } template FieldMatrix operator*(FieldMatrix A, T scalar) { return A *= scalar; } template FieldMatrix operator/(FieldMatrix A, T scalar) { return A /= scalar; } template FieldMatrix operator+(FieldMatrix A, FieldMatrix const& B) { return A += B; } template FieldMatrix operator-(FieldMatrix A, FieldMatrix const& B) { return A -= B; } template FieldVector operator*(FieldMatrix const& mat, FieldVector const& vec) { return Dune::FMatrixHelp::mult(mat, vec); } template FieldMatrix multiplies(FieldMatrix const& A, FieldMatrix const& B) { return A.rightmultiplyany(B); } template FieldMatrix multiplies_AtB(FieldMatrix const& A, FieldMatrix const& B) { FieldMatrix 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 FieldMatrix multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B) { FieldMatrix 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 FieldMatrix& multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B, FieldMatrix& 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 FieldMatrix& multiplies_ABt(FieldMatrix const& A, Dune::DiagonalMatrix const& B, FieldMatrix& 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; } } // end namespace AMDiS