diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc
index 39faa08522c5e967a822168bfe38e04ac3b85357..b12b46c03f6c6f3876414303a4a4f4a9baada820 100644
--- a/examples/navier_stokes.cc
+++ b/examples/navier_stokes.cc
@@ -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>
diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index 8f9b0c904f7678e76ddefbf2f2eee8d095a2e16f..03492346bd8b0e507998ca66b8be429a9b1f9a94 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/src/amdis/common/CMakeLists.txt
@@ -6,6 +6,7 @@ install(FILES
     Concepts.hpp
     ConceptsBase.hpp
     FieldMatVec.hpp
+    FieldMatVec.inc.hpp
     FieldTraits.hpp
     IndexSeq.hpp
     Literals.hpp
diff --git a/src/amdis/common/ConceptsBase.hpp b/src/amdis/common/ConceptsBase.hpp
index acb6861afe597fd3112b4706763c2fa603fa0517..7df39159cfb8163ffa54034b26a62945ce94a7e9 100644
--- a/src/amdis/common/ConceptsBase.hpp
+++ b/src/amdis/common/ConceptsBase.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
 
diff --git a/src/amdis/common/FieldMatVec.hpp b/src/amdis/common/FieldMatVec.hpp
index 3a62f5242d5a008b55d5096fddf8afe44a7132de..b8fac3891e22d708fb6fba10616bd42cf935897d 100644
--- a/src/amdis/common/FieldMatVec.hpp
+++ b/src/amdis/common/FieldMatVec.hpp
@@ -1,17 +1,11 @@
 #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"
diff --git a/src/amdis/common/FieldMatVec.inc.hpp b/src/amdis/common/FieldMatVec.inc.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8353b618d2725aef034be1640605f2c59a985dff
--- /dev/null
+++ b/src/amdis/common/FieldMatVec.inc.hpp
@@ -0,0 +1,473 @@
+#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)
+{
+  return v *= 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;
+}
+
+template <class T, int N, class S,
+  REQUIRES_(Concepts::Arithmetic<S>) >
+FieldVector<T,N> operator/(FieldVector<T,N> v, S factor)
+{
+  return v /= factor;
+}
+
+template <class T>
+FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w)
+{
+  return {v[0] * w[0]};
+}
+
+template <class T, int N>
+FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v)
+{
+  return v *= factor[0];
+}
+
+template <class T, int N>
+FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor)
+{
+  return v *= factor[0];
+}
+
+// ----------------------------------------------------------------------------
+
+/// Cross-product a 2d-vector = orthogonal vector
+template <class T>
+FieldVector<T, 2> cross(FieldVector<T, 2> const& a)
+{
+  return {{ a[1], -a[0] }};
+}
+
+/// Cross-product of two vectors (in 3d only)
+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] }};
+}
+
+/// 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);
+}
+
+template <class T, int N, int M,
+  REQUIRES_( N!=1 && M!=1 )>
+auto operator*(FieldVector<T,N> const& v, FieldVector<T,M> const& w)
+{
+  static_assert(M == N, "Requires vectors of the same type!");
+  return v.dot(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{});
+}
+
+template <class T, int N>
+T sum(FieldMatrix<T, 1, N> const& x)
+{
+  return Impl::accumulate(x, T(0), Operation::Plus{});
+}
+
+
+/// 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);
+}
+
+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);
+}
+
+/// 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{});
+}
+
+template <class T, int N>
+auto max(FieldMatrix<T, 1, N> const& x)
+{
+  return Impl::accumulate(x, std::numeric_limits<T>::lowest(), Operation::Max{});
+}
+
+/// 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{});
+}
+
+template <class T, int N>
+auto min(FieldMatrix<T, 1, N> const& x)
+{
+  return Impl::accumulate(x, std::numeric_limits<T>::max(), Operation::Min{});
+}
+
+/// 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{});
+}
+
+template <class T, int N>
+auto abs_max(FieldMatrix<T, 1, N> const& x)
+{
+  return Impl::accumulate(x, T(0), Operation::AbsMax{});
+}
+
+/// 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{});
+}
+
+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{});
+}
+
+// ----------------------------------------------------------------------------
+
+/** \ingroup vector_norms
+  *  \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);
+}
+
+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);
+}
+
+/** \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));
+}
+
+template <class T, int N>
+auto two_norm(FieldMatrix<T, 1, N> 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 <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 );
+}
+
+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 );
+}
+
+/** \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);
+}
+
+template <class T, int N>
+auto infty_norm(FieldMatrix<T, 1, N> const& x)
+{
+  return abs_max(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);
+}
+
+// ----------------------------------------------------------------------------
+
+/// 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;
+}
+
+// ----------------------------------------------------------------------------
+
+template <class T>
+T det(FieldMatrix<T, 0, 0> const& /*mat*/)
+{
+  return 0;
+}
+
+/// Determinant of a 1x1 matrix
+template <class T>
+T det(FieldMatrix<T, 1, 1> const& mat)
+{
+  return mat[0][0];
+}
+
+/// 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];
+}
+
+/// 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]);
+}
+
+/// Determinant of a NxN matrix
+template <class T,  int N>
+T det(FieldMatrix<T, N, N> const& mat)
+{
+  return mat.determinant();
+}
+
+/// Return the inverse of the matrix `mat`
+template <class T, int N>
+auto inv(FieldMatrix<T, N, N> mat)
+{
+  mat.invert();
+  return 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);
+}
+
+
+/// 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)) );
+}
+
+/// 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]));
+}
+
+// ----------------------------------------------------------------------------
+// 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 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;
+}
+
+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;
+}
+
+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;
+}
+
+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;
+}
+
+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;
+}
+
+
+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);
+}
+
+
+
+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);
+}
+
+
+
+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;
+}
+
+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;
+}
+
+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;
+}
+
+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;
+}
+
+} // end namespace AMDiS
+
+#endif
diff --git a/src/amdis/gridfunctions/ConstantGridFunction.hpp b/src/amdis/gridfunctions/ConstantGridFunction.hpp
index 6e7539f44d7acad1215e9cad05ed256e118b5bc3..9856df9619c1b4c682d00d6a41bf42af0db9a065 100644
--- a/src/amdis/gridfunctions/ConstantGridFunction.hpp
+++ b/src/amdis/gridfunctions/ConstantGridFunction.hpp
@@ -47,9 +47,6 @@ namespace AMDiS
       return 0;
     }
 
-    template <class R_, class D_, class L_, class T_>
-    friend auto derivative(ConstantLocalFunction<R_(D_), L_, T_> const& lf);
-
   private:
     T value_;
   };
@@ -104,22 +101,25 @@ namespace AMDiS
       return value_;
     }
 
-    /// Return the LocalFunction of the \ref ConstantGridFunction. \relates ConstantGridFunction
-    friend LocalFunction localFunction(ConstantGridFunction const& gf)
-    {
-      return LocalFunction{gf.value_};
-    }
-
     EntitySet const& entitySet() const
     {
       return entitySet_;
     }
 
+    T const& value() const { return value_; }
+
   private:
     T value_;
     EntitySet entitySet_;
   };
 
+  /// Return the LocalFunction of the \ref ConstantGridFunction. \relates ConstantGridFunction
+  template <class T, class GV>
+  typename ConstantGridFunction<T,GV>::LocalFunction localFunction(ConstantGridFunction<T,GV> const& gf)
+  {
+    return {gf.value()};
+  }
+
   /** @} **/
 
 
diff --git a/src/amdis/gridfunctions/DerivativeGridFunction.hpp b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
index bdfc164d8fbc000ee9de62087a7237cc258e87a9..f5ae0601133aa38abfcfafe94977afec481a2495 100644
--- a/src/amdis/gridfunctions/DerivativeGridFunction.hpp
+++ b/src/amdis/gridfunctions/DerivativeGridFunction.hpp
@@ -76,10 +76,11 @@ namespace AMDiS
 
 #ifndef DOXYGEN
   template <class GridFct,
-    REQUIRES(not Concepts::HasDerivative<GridFct> &&
-                Concepts::HasLocalFunctionDerivative<GridFct>)>
+            class LocalFct = decltype(localFunction(std::declval<GridFct>())),
+    REQUIRES(not Concepts::HasDerivative<GridFct>)>
   auto derivative(GridFct const& gridFct)
   {
+    static_assert(Concepts::HasDerivative<LocalFct>, "derivative(LocalFunction) not defined!");
     return DerivativeGridFunction<GridFct>{gridFct};
   }
 
diff --git a/test/ConceptsTest.cpp b/test/ConceptsTest.cpp
index 118f25bc0076e00b2922f430000ffa812bc30516..40d14bc335ac9eeda076ead3bf960efa7d5ac0e3 100644
--- a/test/ConceptsTest.cpp
+++ b/test/ConceptsTest.cpp
@@ -4,8 +4,8 @@
 #include <dune/common/reservedvector.hh>
 #include <dune/functions/functionspacebases/flatmultiindex.hh>
 
-#include <amdis/common/Concepts.hpp>
 #include <amdis/common/FieldMatVec.hpp>
+#include <amdis/common/Concepts.hpp>
 
 using namespace AMDiS;
 
diff --git a/test/FilesystemTest.cpp b/test/FilesystemTest.cpp
index 8d54d868cc7e7b3b12941e9abeecc2bcbd73b150..d4a7fa8bd3707958335dc9e889a8c962e3d1c97a 100644
--- a/test/FilesystemTest.cpp
+++ b/test/FilesystemTest.cpp
@@ -14,9 +14,9 @@ void test0()
   AMDIS_TEST_EQ( path("a/b/c.txt").filename(), path("c.txt") );
   AMDIS_TEST_EQ( path("a/b/c.txt").stem(), path("c") );
 
-  AMDIS_TEST_EQ( path("a/b/c.txt").extension(), path(".txt") );
-  AMDIS_TEST_EQ( path("a/b/c.").extension(), path(".") );
-  AMDIS_TEST_EQ( path(".txt").extension(), path(".txt") );
+  AMDIS_TEST_EQ( path("a/b/c.txt").extension().string(), ".txt" );
+  AMDIS_TEST_EQ( path("a/b/c.").extension().string(), "." );
+  AMDIS_TEST_EQ( path(".txt").extension().string(), ".txt" );
 
   AMDIS_TEST( !path("a/b/c").is_absolute() );
   AMDIS_TEST( !path("a\\b\\c").is_absolute() );
@@ -50,7 +50,7 @@ void test1()
   AMDIS_TEST( path("/tmp").is_directory());
   AMDIS_TEST( !path("/tmp").is_file());
   AMDIS_TEST( exists( path("/tmp") ) );
-  
+
   AMDIS_TEST( path("/etc/fstab").is_file() );
   AMDIS_TEST( !path("/etc/fstab").is_directory() );
   AMDIS_TEST( exists("/etc/fstab") );