Commit 45abb132 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/clang_errors' into 'develop'

several issues found by clang corrected, especiall related to FieldMatVec…

See merge request !11
parents f099e85c 949c8fad
Pipeline #1035 failed with stage
in 9 seconds
......@@ -3,6 +3,7 @@
#include <cmath>
#include <amdis/AMDiS.hpp>
#include <amdis/common/FieldMatVec.hpp>
#include <amdis/AdaptInstationary.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/ProblemInstat.hpp>
......
......@@ -6,6 +6,7 @@ install(FILES
Concepts.hpp
ConceptsBase.hpp
FieldMatVec.hpp
FieldMatVec.inc.hpp
FieldTraits.hpp
IndexSeq.hpp
Literals.hpp
......
......@@ -4,9 +4,11 @@
#ifdef DOXYGEN
#define REQUIRES(...)
#define REQUIRES_(...)
#define CONCEPT constexpr
#else
#define REQUIRES(...) std::enable_if_t<__VA_ARGS__ , int>* = nullptr
#define REQUIRES_(...) std::enable_if_t<__VA_ARGS__ , int>*
#define CONCEPT constexpr
#endif
......
#pragma once
#include <algorithm>
#include <limits>
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <amdis/common/Math.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/ScalarTypes.hpp>
#include <amdis/operations/Arithmetic.hpp>
#include <amdis/operations/MaxMin.hpp>
namespace AMDiS
{
......@@ -22,178 +16,87 @@ namespace AMDiS
template <class T, int N, class S,
REQUIRES(Concepts::Arithmetic<S>) >
FieldVector<T,N> operator*(FieldVector<T,N> v, S factor)
{
return v *= factor;
}
FieldVector<T,N> operator*(FieldVector<T,N> v, S factor);
template <class S, class T, int N,
REQUIRES(Concepts::Arithmetic<S>) >
FieldVector<T,N> operator*(S factor, FieldVector<T,N> v)
{
return v *= factor;
}
FieldVector<T,N> operator*(S factor, FieldVector<T,N> v);
template <class T, int N, class S,
REQUIRES(Concepts::Arithmetic<S>) >
FieldVector<T,N> operator/(FieldVector<T,N> v, S factor)
{
return v /= factor;
}
FieldVector<T,N> operator/(FieldVector<T,N> v, S factor);
template <class T>
FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w)
{
return {v[0] * w[0]};
}
FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w);
template <class T, int N>
FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v)
{
return v *= factor[0];
}
FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v);
template <class T, int N>
FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor)
{
return v *= factor[0];
}
FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor);
// ----------------------------------------------------------------------------
/// Cross-product a 2d-vector = orthogonal vector
template <class T>
FieldVector<T, 2> cross(FieldVector<T, 2> const& a)
{
return {{ a[1], -a[0] }};
}
FieldVector<T, 2> cross(FieldVector<T, 2> const& a);
/// Cross-product of two vectors (in 3d only)
/// Cross-product of two 3d-vectors
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] }};
}
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)
{
return vec1.dot(vec2);
}
auto dot(FieldVector<T,N> const& vec1, FieldVector<S,N> const& vec2);
template <class T, int N, int M,
REQUIRES( N!=1 && M!=1 )>
auto operator*(FieldVector<T,N> v, FieldVector<T,M> const& w)
{
static_assert(M == N, "Requires vectors of the same type!");
return v.dot(w);
}
auto operator*(FieldVector<T,N> const& v, FieldVector<T,M> const& w);
// ----------------------------------------------------------------------------
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
/// Sum of vector entires.
template <class T, int N>
T sum(FieldVector<T, N> const& x)
{
return Impl::accumulate(x, T(0), Operation::Plus{});
}
T sum(FieldVector<T, N> const& x);
template <class T, int N>
T sum(FieldMatrix<T, 1, N> const& x)
{
return Impl::accumulate(x, T(0), Operation::Plus{});
}
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)
{
auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
return Impl::accumulate(x, T(0), op);
}
auto unary_dot(FieldVector<T, N> const& x);
template <class T, int N>
auto unary_dot(FieldMatrix<T, 1, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + Math::sqr(std::abs(b)); };
return Impl::accumulate(x, T(0), op);
}
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)
{
return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
}
auto max(FieldVector<T, N> const& x);
template <class T, int N>
auto max(FieldMatrix<T, 1, N> const& x)
{
return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
}
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)
{
return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
}
auto min(FieldVector<T, N> const& x);
template <class T, int N>
auto min(FieldMatrix<T, 1, N> const& x)
{
return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
}
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)
{
return Impl::accumulate(x, T(0), Operation::AbsMax{});
}
auto abs_max(FieldVector<T, N> const& x);
template <class T, int N>
auto abs_max(FieldMatrix<T, 1, N> const& x)
{
return Impl::accumulate(x, T(0), Operation::AbsMax{});
}
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)
{
return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::AbsMin{});
}
auto abs_min(FieldVector<T, N> const& x);
template <class T, int N>
auto abs_min(FieldMatrix<T, 1, N> const& x)
{
return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::AbsMin{});
}
auto abs_min(FieldMatrix<T, 1, N> const& x);
// ----------------------------------------------------------------------------
......@@ -201,278 +104,138 @@ namespace AMDiS
* \brief The 1-norm of a vector = sum_i |x_i|
**/
template <class T, int N>
auto one_norm(FieldVector<T, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
return Impl::accumulate(x, T(0), op);
}
auto one_norm(FieldVector<T, N> const& x);
template <class T, int N>
auto one_norm(FieldMatrix<T, 1, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + std::abs(b); };
return Impl::accumulate(x, T(0), op);
}
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)
{
return std::sqrt(unary_dot(x));
}
auto two_norm(FieldVector<T, N> const& x);
template <class T, int N>
auto two_norm(FieldMatrix<T, 1, N> const& x)
{
return std::sqrt(unary_dot(x));
}
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)
{
auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
}
auto p_norm(FieldVector<T, N> const& x);
template <int p, class T, int N>
auto p_norm(FieldMatrix<T, 1, N> const& x)
{
auto op = [](auto const& a, auto const& b) { return a + Math::pow<p>(std::abs(b)); };
return std::pow( Impl::accumulate(x, T(0), op), 1.0/p );
}
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)
{
return abs_max(x);
}
auto infty_norm(FieldVector<T, N> const& x);
template <class T, int N>
auto infty_norm(FieldMatrix<T, 1, N> const& x)
{
return abs_max(x);
}
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)
{
T result = 0;
for (int i = 0; i < N; ++i)
result += Math::sqr(lhs[i] - rhs[i]);
return std::sqrt(result);
}
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)
{
using result_type = FieldMatrix<decltype( 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;
}
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*/)
{
return 0;
}
T det(FieldMatrix<T, 0, 0> const& /*mat*/);
/// Determinant of a 1x1 matrix
template <class T>
T det(FieldMatrix<T, 1, 1> const& mat)
{
return mat[0][0];
}
T det(FieldMatrix<T, 1, 1> const& mat);
/// Determinant of a 2x2 matrix
template <class T>
T det(FieldMatrix<T, 2, 2> const& mat)
{
return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
}
T det(FieldMatrix<T, 2, 2> const& mat);
/// Determinant of a 3x3 matrix
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]);
}
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 mat.determinant();
}
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)
{
mat.invert();
return mat;
}
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)
{
A.solve(x, b);
}
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)
{
using std::sqrt;
return sqrt( det(outer(DF, DF)) );
}
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)
{
using std::sqrt;
return sqrt(dot(DF[0], DF[0]));
}
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)
{
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;
}
FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A);
template <class T, int M, int N, class S,
REQUIRES(Concepts::Arithmetic<S>) >
FieldMatrix<T,M,N> operator*(S scalar, FieldMatrix<T, M, N> A)
{
return A *= scalar;
}
FieldMatrix<T,M,N> operator*(S scalar, FieldMatrix<T, M, N> A);
template <class T, int M, int N, class S,
REQUIRES(Concepts::Arithmetic<S>) >
FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, S scalar)
{
return A *= scalar;
}
FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, S scalar);
template <class T, int M, int N, class S,
REQUIRES(Concepts::Arithmetic<S>) >
FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar)
{
return A /= scalar;
}
FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, S scalar);
template <class T, int M, int N>
FieldMatrix<T,M,N> operator+(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
{
return A += B;
}
FieldMatrix<T,M,N> operator+(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B);
template <class T, int M, int N>
FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
{
return A -= B;
}
FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B);
template <class T, int N, int M>
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec)
{
return Dune::FMatrixHelp::mult(mat, vec);
}
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec);
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A, FieldMatrix<T, L, N> const& B)
{
return A.rightmultiplyany(B);
}
FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A, FieldMatrix<T, L, N> const& B);
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B)
{
FieldMatrix<T,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;
}
FieldMatrix<T,M,N> multiplies_AtB(FieldMatrix<T, L, M> const& A, FieldMatrix<T, N, L> const& B);
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B)
{
FieldMatrix<T,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;
}
FieldMatrix<T,M,N> multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B);
template <class T, int M, int N, int L>
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B, FieldMatrix<T,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;
}
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, L> const& A, FieldMatrix<T, N, L> const& B, FieldMatrix<T,M,N>& C);
template <class T, int M, int N>
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, Dune::DiagonalMatrix<T, N> const& B, FieldMatrix<T,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;
}
FieldMatrix<T,M,N>& multiplies_ABt(FieldMatrix<T, M, N> const& A, Dune::DiagonalMatrix<T, N> const& B, FieldMatrix<T,M,N>& C);
} // end namespace AMDiS
#include "FieldMatVec.inc.hpp"
#pragma once
#include <algorithm>
#include <limits>
#include <amdis/common/Math.hpp>
#include <amdis/operations/Arithmetic.hpp>
#include <amdis/operations/MaxMin.hpp>
#ifndef DOXYGEN
namespace AMDiS {