Commit 949c8fad authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

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

several issues found by clang corrected, especiall related to FieldMatVec definitions and localFunction/derivative of constantGridFunction
parent f099e85c
......@@ -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 {
// some arithmetic operations with FieldVector
template <class T, int N, class S,
REQUIRES_(Concepts::Arithmetic<S>) >
FieldVector<T,N> operator*(FieldVector<T,N> v, S factor)
{