diff --git a/CMakeLists.txt b/CMakeLists.txt
index 55c1d328a6331d742d743f94ddf941c4dfd72c9e..8fd26a40be807bf479348512a2b07c037c101a4a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -21,6 +21,7 @@ dune_enable_all_packages(MODULE_LIBRARIES amdis)
 
 include(AmdisMacros)
 add_subdirectory("src")
+add_subdirectory("test")
 add_subdirectory("dune")
 add_subdirectory("doc")
 add_subdirectory("cmake/modules")
diff --git a/dune/amdis/FileWriter.hpp b/dune/amdis/FileWriter.hpp
index 2dbeebf774a4a64a474046281a47ea34ebbf1d19..e0a16cfa74c6074f5b83de2b7c8e3e08a612b918 100644
--- a/dune/amdis/FileWriter.hpp
+++ b/dune/amdis/FileWriter.hpp
@@ -91,6 +91,8 @@ namespace AMDiS
     /// Implements \ref FileWriterInterface::writeFiles
     virtual void writeFiles(AdaptInfo& adaptInfo, bool force) override
     {
+      if (!filesystem::exists(dir_))
+        error_exit("Output directory '",dir_,"' does not exist!");
       if (vtkSeqWriter_)
         vtkSeqWriter_->write(adaptInfo.getTime(), mode_);
     }
diff --git a/dune/amdis/Mesh.hpp b/dune/amdis/Mesh.hpp
index ea09a268bc58867a1b676369b28aff67f818e4d3..9384f98a91a02385f3b3b455c3c8f37f2e984f7f 100644
--- a/dune/amdis/Mesh.hpp
+++ b/dune/amdis/Mesh.hpp
@@ -146,7 +146,7 @@ namespace AMDiS
 
       msg("L = ", L);
 
-      auto s = Dune::filledArray<std::size_t(dim)>(2); // number of cells on coarse mesh in each direction
+      auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
       Parameters::get(meshName + "->num cells", s);
 
       // TODO: add more parameters for yasp-grid (see constructor)
@@ -168,7 +168,7 @@ namespace AMDiS
       Parameters::get(meshName + "->min corner", lowerleft);
       Parameters::get(meshName + "->max corner", upperright);
 
-      auto s = Dune::filledArray<std::size_t(dim)>(2); // number of cells on coarse mesh in each direction
+      auto s = Dune::filledArray<std::size_t(dim)>(1); // number of cells on coarse mesh in each direction
       Parameters::get(meshName + "->num cells", s);
 
       // TODO: add more parameters for yasp-grid (see constructor)
diff --git a/dune/amdis/assembler/FirstOrderTestPartialTrial.hpp b/dune/amdis/assembler/FirstOrderTestPartialTrial.hpp
index 36b0649f7ee1549c5cbc40e5f608cdfd6d57d816..c7ebd529d07def559c23e22b97683b6668bd8488 100644
--- a/dune/amdis/assembler/FirstOrderTestPartialTrial.hpp
+++ b/dune/amdis/assembler/FirstOrderTestPartialTrial.hpp
@@ -59,6 +59,7 @@ namespace AMDiS
 
         // The transposed inverse Jacobian of the map from the reference element to the element
         const auto jacobian = context.geometry.jacobianInverseTransposed(local);
+        assert(jacobian.M() == Context::dim);
 
         // The multiplicative factor in the integral transformation formula
         double factor = context.integrationElement(quad[iq].position()) * quad[iq].weight();
diff --git a/dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp b/dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp
index 5a3dee5729a46e98e2b4243a4d54d2ecb924ae51..a56500a478990cff1e99fe7eb60c9fe064c98910 100644
--- a/dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp
+++ b/dune/amdis/assembler/ZeroOrderTestvecTrialvec.hpp
@@ -41,8 +41,8 @@ namespace AMDiS
                                 ElementMatrix& elementMatrix,
                                 RowNode const& rowNode,
                                 ColNode const& colNode,
-                                std::integral_constant<bool, sameFE>,
-                                std::integral_constant<bool, sameNode>)
+                                bool_t<sameFE>,
+                                bool_t<sameNode>)
     {
       static_assert(RowNode::isPower && ColNode::isPower,
         "Operator can be applied to Power-Nodes only.");
@@ -50,12 +50,10 @@ namespace AMDiS
       static_assert( (RowNode::CHILDREN == ColNode::CHILDREN), "" );
       // theoretically the restriction of the same nr of childs would not be necessary
 
-      using Category = ValueCategory_t<typename GridFct::Range>;
-      static_assert( std::is_same<Category, tag::scalar>::value || std::is_same<Category, tag::matrix>::value, "" );
-
+      using Category = ValueCategory_t<expr_value_type>;
       calcElementMatrix(context, quad, elementMatrix, rowNode, colNode,
-        std::integral_constant<bool, sameFE>{},
-        std::integral_constant<bool, sameNode>{},
+        bool_<sameFE>,
+        bool_<sameNode>,
         Category{});
     }
 
@@ -68,7 +66,7 @@ namespace AMDiS
                            ElementMatrix& elementMatrix,
                            RowNode const& rowNode,
                            ColNode const& colNode,
-                           std::integral_constant<bool, sameFE>,
+                           bool_t<sameFE>,
                            std::false_type /*sameNode*/,
                            tag::scalar)
     {
@@ -165,12 +163,12 @@ namespace AMDiS
                            ElementMatrix& elementMatrix,
                            RowNode const& rowNode,
                            ColNode const& colNode,
-                           std::integral_constant<bool, sameFE>,
-                           std::integral_constant<bool, sameNode>,
+                           bool_t<sameFE>,
+                           bool_t<sameNode>,
                            tag::matrix)
     {
       static const std::size_t CHILDREN = RowNode::CHILDREN;
-      static_assert( std::is_same<typename GridFct::Range, Dune::FieldMatrix<double, CHILDREN, CHILDREN>>::value, "" );
+      static_assert( std::is_same<expr_value_type, Dune::FieldMatrix<double, CHILDREN, CHILDREN>>::value, "" );
 
       auto const& rowLocalFE = rowNode.child(0).finiteElement();
       auto const& colLocalFE = colNode.child(0).finiteElement();
@@ -191,10 +189,10 @@ namespace AMDiS
           colLocalFE.localBasis().evaluateFunction(local, colShapeValues_);
 
         for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
-          double value = factor * (*rowShapeValues)[i];
+          double value0 = factor * (*rowShapeValues)[i];
 
           for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
-            value *= (*colShapeValues)[j];
+            double value = value0 * (*colShapeValues)[j];
             auto mat = exprValue * value;
 
             for (std::size_t k0 = 0; k0 < CHILDREN; ++k0) {
diff --git a/dune/amdis/common/CMakeLists.txt b/dune/amdis/common/CMakeLists.txt
index fdcc6f1d9c378e16f7bc4823b79c43c16f59d3e9..6c97462b8b72508ef8df20362167ce1e290c3b6a 100644
--- a/dune/amdis/common/CMakeLists.txt
+++ b/dune/amdis/common/CMakeLists.txt
@@ -11,7 +11,7 @@ install(FILES
     Loops.hpp
     Math.hpp
     Mpl.hpp
-    MultiTypeFieldVector.hpp
+    MultiTypeVector.hpp
     ScalarTypes.hpp
     Size.hpp
     TupleUtility.hpp
diff --git a/dune/amdis/common/Concepts.hpp b/dune/amdis/common/Concepts.hpp
index 9e31806000b2ad08e60ea494df8c74008d38e11f..846d2805fa45657213b4a6ed3eec7ea0fa4ae766 100644
--- a/dune/amdis/common/Concepts.hpp
+++ b/dune/amdis/common/Concepts.hpp
@@ -108,6 +108,19 @@ namespace AMDiS
           >;
       };
 
+      // idx[0]
+      struct MultiIndex
+      {
+        template <class MI>
+        auto requires_(MI&& idx) -> decltype(
+          Concepts::valid_expr(
+            idx[0],
+            idx.size(),
+            idx.max_size()
+            /* ,idx.resize() */
+          ));
+      };
+
     } // end namespace Definition
 #endif // DOXYGEN
 
@@ -166,6 +179,9 @@ namespace AMDiS
     template <class F, class... Args>
     constexpr bool Predicate = Functor<F, bool(Args...)>;
 
+    template <class MI>
+    constexpr bool MultiIndex = models<Definition::MultiIndex(MI)>;
+
     /** @} **/
 
   } // end namespace Concepts
diff --git a/dune/amdis/common/FieldMatVec.hpp b/dune/amdis/common/FieldMatVec.hpp
index 1f0602bd291c4ca2753b9a15e197db22092d64da..fff29c748182dfa4832cfb7fbd2ea11cefb690a1 100644
--- a/dune/amdis/common/FieldMatVec.hpp
+++ b/dune/amdis/common/FieldMatVec.hpp
@@ -2,6 +2,7 @@
 
 #include <algorithm>
 
+#include <dune/common/diagonalmatrix.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 
@@ -371,20 +372,23 @@ namespace AMDiS
   }
 
 
-  template <class T, int M, int N>
-  FieldMatrix<T,M,N> operator*(T scalar, FieldMatrix<T, M, N> 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;
   }
 
-  template <class T, int M, int N>
-  FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, T 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, T 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;
   }
diff --git a/dune/amdis/common/FieldTraits.hpp b/dune/amdis/common/FieldTraits.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bd16c80e92568b0ceeeeb4bf49405096a06eb759
--- /dev/null
+++ b/dune/amdis/common/FieldTraits.hpp
@@ -0,0 +1,16 @@
+#pragma once
+
+#include <dune/common/ftraits.hh>
+
+#include <dune/amdis/common/Utility.hpp>
+
+namespace AMDiS
+{
+  template <class... T>
+  struct CommonFieldTraits
+  {
+    using field_type = CommonType_t<typename Dune::FieldTraits<T>::field_type...>;
+    using real_type = CommonType_t<typename Dune::FieldTraits<T>::real_type...>;
+  };
+
+} // end namespace AMDiS
diff --git a/dune/amdis/common/Mpl.hpp b/dune/amdis/common/Mpl.hpp
index b3282eac8b3b64c37614f53266aa853e7bd40c45..4c4f65b28d7b6eacaed24f6dc48244bad34ecada 100644
--- a/dune/amdis/common/Mpl.hpp
+++ b/dune/amdis/common/Mpl.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 // std c++ headers
+#include <tuple>
 #include <type_traits>
 #include <utility>
 
@@ -152,8 +153,8 @@ namespace AMDiS
   constexpr bool_t<!B> operator!(bool_t<B>) { return {}; }
 
 
-  template <class T, T A, T B>
-  using is_equal = std::conditional_t<A == B, std::true_type, std::false_type>;
+  template <class T, T value0, T... values>
+  using is_equal = all_of_t<(value0 == values)...>;
 
   template <class T, class... Ts>
   using is_one_of = or_t<std::is_same<T, Ts>::value...>;
diff --git a/dune/amdis/common/MultiTypeFieldVector.hpp b/dune/amdis/common/MultiTypeFieldVector.hpp
deleted file mode 100644
index 86fc5442fa74322b580ed67f8aba39278201b6de..0000000000000000000000000000000000000000
--- a/dune/amdis/common/MultiTypeFieldVector.hpp
+++ /dev/null
@@ -1,143 +0,0 @@
-#pragma once
-
-#include <tuple>
-
-#include <dune/common/ftraits.hh>
-
-#include <dune/amdis/common/Concepts.hpp>
-#include <dune/amdis/common/Mpl.hpp>
-#include <dune/amdis/common/Size.hpp>
-
-namespace AMDiS
-{
-  // forward declaration
-  template <class FV0, class... FV>
-  class MultiTypeFieldVector;
-}
-
-namespace Dune
-{
-  template <class FV0, class... FV>
-  struct FieldTraits<AMDiS::MultiTypeFieldVector<FV0,FV...>>
-  {
-    using field_type = typename Dune::FieldTraits<FV0>::field_type;
-    using real_type = typename Dune::FieldTraits<FV0>::real_type;
-  };
-}
-
-
-namespace AMDiS
-{
-  template <class FV0, class... FV>
-  class MultiTypeFieldVector
-      : public std::tuple<FV0,FV...>
-  {
-    using Self = MultiTypeFieldVector;
-    using Super = std::tuple<FV0,FV...>;
-
-  public:
-    using field_type = typename Dune::FieldTraits<FV0>::field_type;
-    using real_type = typename Dune::FieldTraits<FV0>::real_type;
-    using size_type = std::size_t;
-
-    enum {
-      dimension = std::tuple_size<Super>::value
-    };
-
-    template <class FV0_, class... FV_,
-      REQUIRES( Concepts::Similar<Types<FV0,FV...>, Types<FV0_,FV_...>> )>
-    MultiTypeFieldVector(FV0_&& fv0, FV_&&... fv)
-      : Super(std::forward<FV0_>(fv0), std::forward<FV_>(fv)...)
-    {}
-
-    /// Default construction of tuple of FieldVectors
-    MultiTypeFieldVector() = default;
-
-    /// Construct tuple by initializing all tuple elements with a constant value
-    MultiTypeFieldVector(real_type value)
-    {
-      *this = value;
-    }
-
-    /// Assignment of real number to all tuple elements
-    MultiTypeFieldVector& operator=(real_type value)
-    {
-      forEach(*this, [value](auto& fv) { fv = value; });
-      return *this;
-    }
-
-    // Compound assignment operator +=
-    MultiTypeFieldVector& operator+=(MultiTypeFieldVector const& that)
-    {
-      forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
-      return *this;
-    }
-
-    // Compound assignment operator -=
-    MultiTypeFieldVector& operator-=(MultiTypeFieldVector const& that)
-    {
-      forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
-      return *this;
-    }
-
-    // Scaling of all tuple elements by a constant value
-    MultiTypeFieldVector& operator*=(real_type value)
-    {
-      forEach(*this, [value](auto& fv) { fv *= value; });
-      return *this;
-    }
-
-    // Scaling of all tuple elements by the inverse of a constant value
-    MultiTypeFieldVector& operator/=(real_type value)
-    {
-      forEach(*this, [value](auto& fv) { fv /= value; });
-      return *this;
-    }
-
-    /// Const access to the tuple elements
-    template <std::size_t I>
-    decltype(auto) operator[](const index_t<I>&) const
-    {
-      return std::get<I>(*this);
-    }
-
-    /// Mutable access to the tuple elements
-    template <std::size_t I>
-    decltype(auto) operator[](const index_t<I>&)
-    {
-      return std::get<I>(*this);
-    }
-
-    /// Return number of elements of the tuple
-    static constexpr std::size_t size()
-    {
-      return std::tuple_size<Super>::value;
-    }
-
-    // /// Return number of elements in a flat vector
-    // static constexpr std::size_t flat_size()
-    // {
-    //   return Size<Self>;
-    // }
-
-    // /// Create flat FieldVector
-    // operator Dune::FieldVector<field_type, flat_size()>() const
-    // {
-    //   Dune::FieldVector<field_type, flat_size()> result;
-
-    //   int shift = 0;
-    //   forEach(*this, [&result,&shift](auto& fv)
-    //   {
-    //     static const int sub_size = Size<decltype(fv)>;
-    //     Dune::FieldVector<field_type, sub_size> sub = fv;
-    //     for (int i = 0; i < sub_size; ++i)
-    //       result[shift + i] = sub[i];
-    //     shift += sub_size;
-    //   });
-
-    //   return result;
-    // }
-
-  };
-
-} // end namesüace AMDiS
diff --git a/dune/amdis/common/MultiTypeMatrix.hpp b/dune/amdis/common/MultiTypeMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..66e0b8e06d7b6b2631812a4b7dadf509a43618c9
--- /dev/null
+++ b/dune/amdis/common/MultiTypeMatrix.hpp
@@ -0,0 +1,131 @@
+#pragma once
+
+#include <tuple>
+
+#include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/common/FieldTraits.hpp>
+#include <dune/amdis/common/Loops.hpp>
+#include <dune/amdis/common/Mpl.hpp>
+#include <dune/amdis/common/MultiTypeVector.hpp>
+#include <dune/amdis/common/Size.hpp>
+#include <dune/amdis/utility/MultiIndex.hpp>
+
+namespace AMDiS
+{
+  // forward declaration
+  template <class... Rows>
+  class MultiTypeMatrix;
+}
+
+namespace Dune
+{
+  template <class... Rows>
+  struct FieldTraits<AMDiS::MultiTypeMatrix<Rows...>>
+  {
+    using field_type = typename AMDiS::CommonFieldTraits<Rows...>::field_type;
+    using real_type = typename AMDiS::CommonFieldTraits<Rows...>::real_type;
+  };
+}
+
+
+namespace AMDiS
+{
+  // Rows should be of type MultiTypeVector
+  template <class... Rows>
+  class MultiTypeMatrix
+      : public std::tuple<Rows...>
+  {
+    using Self = MultiTypeMatrix;
+    using Super = std::tuple<Rows...>;
+
+    static_assert(is_equal<int, Rows::dimension...>::value,
+      "All columns must have the same length.");
+
+  public:
+    using field_type = typename Dune::FieldTraits<Self>::field_type;
+    using real_type = typename Dune::FieldTraits<Self>::real_type;
+    using size_type = std::size_t;
+
+    enum {
+      rows = std::tuple_size<Super>::value,
+      cols = Math::max(Rows::dimension...)
+    };
+
+    template <class... Rows_,
+      REQUIRES( Concepts::Similar<Types<Rows...>, Types<Rows_...>> )>
+    MultiTypeMatrix(Rows_&&... rows)
+      : Super(std::forward<Rows_>(rows)...)
+    {}
+
+    /// Default construction of tuple of FieldVectors
+    MultiTypeMatrix() = default;
+
+    /// Construct tuple by initializing all tuple elements with a constant value
+    MultiTypeMatrix(real_type value)
+    {
+      *this = value;
+    }
+
+    /// Assignment of real number to all tuple elements
+    MultiTypeMatrix& operator=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv = value; });
+      return *this;
+    }
+
+    // Compound assignment operator +=
+    MultiTypeMatrix& operator+=(MultiTypeMatrix const& that)
+    {
+      forEach(range_<0,rows>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
+      return *this;
+    }
+
+    // Compound assignment operator -=
+    MultiTypeMatrix& operator-=(MultiTypeMatrix const& that)
+    {
+      forEach(range_<0,rows>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
+      return *this;
+    }
+
+    // Scaling of all tuple elements by a constant value
+    MultiTypeMatrix& operator*=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv *= value; });
+      return *this;
+    }
+
+    // Scaling of all tuple elements by the inverse of a constant value
+    MultiTypeMatrix& operator/=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv /= value; });
+      return *this;
+    }
+
+    /// Const access to the tuple elements
+    template <std::size_t I, std::size_t J>
+    decltype(auto) operator()(const index_t<I>&, const index_t<J>&) const
+    {
+      return std::get<J>(std::get<I>(*this));
+    }
+
+    /// Mutable access to the tuple elements
+    template <std::size_t I, std::size_t J>
+    decltype(auto) operator()(const index_t<I>&, const index_t<J>&)
+    {
+      return std::get<J>(std::get<I>(*this));
+    }
+
+    /// Return number of elements of the tuple
+    static constexpr std::size_t num_rows()
+    {
+      return rows;
+    }
+
+    /// Return number of elements of the tuple
+    static constexpr std::size_t num_cols()
+    {
+      return cols;
+    }
+  };
+
+} // end namespace AMDiS
diff --git a/dune/amdis/common/MultiTypeVector.hpp b/dune/amdis/common/MultiTypeVector.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f31823ee6c6e4852e2b9302e9fb49b9a64ce733
--- /dev/null
+++ b/dune/amdis/common/MultiTypeVector.hpp
@@ -0,0 +1,136 @@
+#pragma once
+
+#include <tuple>
+
+#include <dune/functions/common/indexaccess.hh>
+
+#include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/common/FieldTraits.hpp>
+#include <dune/amdis/common/Loops.hpp>
+#include <dune/amdis/common/Mpl.hpp>
+#include <dune/amdis/common/Size.hpp>
+
+namespace AMDiS
+{
+  // forward declaration
+  template <class... FV>
+  class MultiTypeVector;
+}
+
+namespace Dune
+{
+  template <class... FV>
+  struct FieldTraits<AMDiS::MultiTypeVector<FV...>>
+  {
+    using field_type = typename AMDiS::CommonFieldTraits<FV...>::field_type;
+    using real_type = typename AMDiS::CommonFieldTraits<FV...>::real_type;
+  };
+}
+
+
+namespace AMDiS
+{
+  template <class... FV>
+  class MultiTypeVector
+      : public std::tuple<FV...>
+  {
+    using Self = MultiTypeVector;
+    using Super = std::tuple<FV...>;
+
+  public:
+    using field_type = typename Dune::FieldTraits<Self>::field_type;
+    using real_type = typename Dune::FieldTraits<Self>::real_type;
+    using size_type = std::size_t;
+
+    enum {
+      dimension = std::tuple_size<Super>::value
+    };
+
+    template <class... FV_,
+      REQUIRES( Concepts::Similar<Types<FV...>, Types<FV_...>> )>
+    MultiTypeVector(FV_&&... fv)
+      : Super(std::forward<FV_>(fv)...)
+    {}
+
+    /// Default construction of tuple of FieldVectors
+    MultiTypeVector() = default;
+
+    /// Construct tuple by initializing all tuple elements with a constant value
+    explicit MultiTypeVector(real_type value)
+    {
+      *this = value;
+    }
+
+    /// Assignment of real number to all tuple elements
+    MultiTypeVector& operator=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv = value; });
+      return *this;
+    }
+
+    // Compound assignment operator +=
+    MultiTypeVector& operator+=(MultiTypeVector const& that)
+    {
+      forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] += that[_i]; });
+      return *this;
+    }
+
+    // Compound assignment operator -=
+    MultiTypeVector& operator-=(MultiTypeVector const& that)
+    {
+      forEach(range_<0,dimension>, [&that,this](auto const _i) { (*this)[_i] -= that[_i]; });
+      return *this;
+    }
+
+    // Scaling of all tuple elements by a constant value
+    MultiTypeVector& operator*=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv *= value; });
+      return *this;
+    }
+
+    // Scaling of all tuple elements by the inverse of a constant value
+    MultiTypeVector& operator/=(real_type value)
+    {
+      forEach(*this, [value](auto& fv) { fv /= value; });
+      return *this;
+    }
+
+    /// Const access to the tuple elements
+    template <std::size_t I>
+    decltype(auto) operator[](const index_t<I>&) const
+    {
+      return std::get<I>(*this);
+    }
+
+    /// Mutable access to the tuple elements
+    template <std::size_t I>
+    decltype(auto) operator[](const index_t<I>&)
+    {
+      return std::get<I>(*this);
+    }
+
+    /// Const access to the vector using multi-indices
+    template <class Index,
+      REQUIRES( Concepts::MultiIndex<Index> )>
+    decltype(auto) operator[](Index const& index) const
+    {
+      return Dune::Functions::hybridMultiIndexAccess<field_type const&>(*this, index);
+    }
+
+    /// Mutable access to the vector using multi-indices
+    template <class Index,
+      REQUIRES( Concepts::MultiIndex<Index> )>
+    decltype(auto) operator[](Index const& index)
+    {
+      return Dune::Functions::hybridMultiIndexAccess<field_type&>(*this, index);
+    }
+
+    /// Return number of elements of the tuple
+    static constexpr std::size_t size()
+    {
+      return dimension;
+    }
+  };
+
+} // end namespace AMDiS
diff --git a/dune/amdis/common/TypeDefs.hpp b/dune/amdis/common/TypeDefs.hpp
index 3f6aeffe784ddc86eb741b14143e4a3d586b0c5e..1bafb5a11a3061249762a9ab054c1c1afee02763 100644
--- a/dune/amdis/common/TypeDefs.hpp
+++ b/dune/amdis/common/TypeDefs.hpp
@@ -2,9 +2,7 @@
 
 #include <type_traits>
 
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <dune/amdis/linear_algebra/Mtl.hpp>
 
 namespace AMDiS
 {
diff --git a/dune/amdis/common/Utility.hpp b/dune/amdis/common/Utility.hpp
index 6908b2f986bc51aa65dff82754e9af4426e2a378..dd897b88aeb70d0833899abd705c9c5909490c7a 100644
--- a/dune/amdis/common/Utility.hpp
+++ b/dune/amdis/common/Utility.hpp
@@ -25,6 +25,25 @@ namespace AMDiS
   using Underlying_t = typename Impl::UnderlyingType<T>::type;
 
 
+  namespace Impl
+  {
+    template <class T0, class... Ts>
+    struct CommonType
+    {
+      using type = std::common_type_t<T0, typename CommonType<Ts...>::type>;
+    };
+
+    template <class T0>
+    struct CommonType<T0>
+    {
+      using type = T0;
+    };
+  }
+
+  template <class... T>
+  using CommonType_t = typename Impl::CommonType<T...>::type;
+
+
   // ---------------------------------------------------------------------------
 
 
diff --git a/dune/amdis/io/FileWriterInterface.hpp b/dune/amdis/io/FileWriterInterface.hpp
index a21f5fc2aebcce39aec649d1a56b22c2a33bdca7..3601272c9b291215a1150c702e173642e9923af7 100644
--- a/dune/amdis/io/FileWriterInterface.hpp
+++ b/dune/amdis/io/FileWriterInterface.hpp
@@ -22,9 +22,6 @@ namespace AMDiS
       Parameters::get(base + "->output directory", dir_);
       name_ = "solution";
       Parameters::get(base + "->name", name_);
-
-      if (!filesystem::exists(dir_))
-        error_exit("Output directory '",dir_,"' does not exist!");
     }
 
 
diff --git a/dune/amdis/linear_algebra/DOFVectorBase.hpp b/dune/amdis/linear_algebra/DOFVectorBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2b3902eed018517804ee296c1ba3f36288269f73
--- /dev/null
+++ b/dune/amdis/linear_algebra/DOFVectorBase.hpp
@@ -0,0 +1,30 @@
+#pragma once
+
+namespace AMDiS
+{
+  template <class GlobalBasisType>
+  class DOFVectorBase
+  {
+  public:
+    using GlobalBasis = GlobalBasisType;
+    using SizeType = std::size_t;
+
+    DOFVectorBase(GlobalBasisType const& basis)
+      : basis_(basis)
+    {}
+
+    SizeType size() const
+    {
+      return basis_.dimension();
+    }
+
+    GlobalBasis const& basis() const
+    {
+      return GlobalBasis;
+    }
+
+  private:
+    GlobalBasis const& basis_;
+  };
+
+} // end namespace AMDiS
diff --git a/dune/amdis/linear_algebra/Mtl.hpp b/dune/amdis/linear_algebra/Mtl.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d89ea427fd21f3cd57d2926100c19586478b532
--- /dev/null
+++ b/dune/amdis/linear_algebra/Mtl.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <dune/common/ftraits.hh>
+
+#include <boost/numeric/mtl/matrix/compressed2D.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+
+namespace Dune
+{
+  template <class Elt, class Parameters>
+  struct FieldTraits<mtl::compressed2D<Elt,Parameters>>
+  {
+    using field_type = typename FieldTraits<Elt>::field_type;
+    using real_type = typename FieldTraits<Elt>::real_type;
+  };
+
+  template <class Value, class Parameters>
+  struct FieldTraits<mtl::dense2D<Value,Parameters>>
+  {
+    using field_type = typename FieldTraits<Value>::field_type;
+    using real_type = typename FieldTraits<Value>::real_type;
+  };
+
+  template <class Value, class Parameters>
+  struct FieldTraits<mtl::dense_vector<Value,Parameters>>
+  {
+    using field_type = typename FieldTraits<Value>::field_type;
+    using real_type = typename FieldTraits<Value>::real_type;
+  };
+
+} // end namespace Dune
diff --git a/dune/amdis/operations/Basic.hpp b/dune/amdis/operations/Basic.hpp
index e37aed4bc31f5b525a0ede071f69ad3042d22fac..bcc5d0f8f00ff0867abd3b8368312afde9f86b0c 100644
--- a/dune/amdis/operations/Basic.hpp
+++ b/dune/amdis/operations/Basic.hpp
@@ -2,7 +2,9 @@
 
 #include <algorithm>
 
+#include <dune/amdis/common/IndexSeq.hpp>
 #include <dune/amdis/common/Math.hpp>
+#include <dune/amdis/common/Mpl.hpp>
 #include <dune/amdis/common/Concepts.hpp>
 #include <dune/amdis/common/ScalarTypes.hpp>
 
diff --git a/dune/amdis/utility/MultiIndex.hpp b/dune/amdis/utility/MultiIndex.hpp
index ebc86148806f2cfc5259f857300cf8d7de18b347..83a3ce04505d58b411edbb90861c810b0ff5ef4b 100644
--- a/dune/amdis/utility/MultiIndex.hpp
+++ b/dune/amdis/utility/MultiIndex.hpp
@@ -1,8 +1,11 @@
 #pragma once
 
 #include <dune/common/reservedvector.hh>
+#include <dune/functions/common/indexaccess.hh>
 #include <dune/functions/functionspacebases/flatmultiindex.hh>
 
+#include <dune/amdis/common/Mpl.hpp>
+
 namespace AMDiS
 {
   std::size_t flatMultiIndex(std::size_t idx) { return idx; }
@@ -11,4 +14,45 @@ namespace AMDiS
 
   std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx) { return idx[0]; }
 
+
+
+  template <int rank, class Index>
+  struct MultiIndexResolver
+  {
+    MultiIndexResolver(Index const& index)
+      : index_(index)
+    {}
+
+    template <class C>
+    decltype(auto) operator()(C&& c)
+    {
+      return resolve(std::forward<C>(c), int_<rank>);
+    }
+
+  private:
+    template <class C, int r>
+    decltype(auto) resolve(C&& c, int_t<r>)
+    {
+      using namespace Dune::Indices;
+      auto&& subIndex = Dune::Functions::shiftedMultiIndex(index_);
+      auto&& subIndexResolver = MultiIndexResolver<r-1, decltype(subIndex)>(subIndex);
+      return Dune::Functions::hybridIndexAccess(c, index_[_0], subIndexResolver);
+    }
+
+    template <class C>
+    decltype(auto) resolve(C&& c, int_t<0>)
+    {
+      return std::forward<C>(c);
+    }
+
+    Index const& index_;
+  };
+
+  template <int rank, class Container, class MultiIndex>
+  decltype(auto) multiIndexAccess(Container&& c, MultiIndex const& index)
+  {
+    MultiIndexResolver<rank, MultiIndex> multiIndexResolver(index);
+    return multiIndexResolver(c);
+  }
+
 } // end namespace AMDiS
diff --git a/dune/amdis/utility/RangeType.hpp b/dune/amdis/utility/RangeType.hpp
index 1ce5c7516a04513b4688e80cfa4b1d6b0ac68ed3..857af1e89de830fdf47de101eb256a255d3f1a36 100644
--- a/dune/amdis/utility/RangeType.hpp
+++ b/dune/amdis/utility/RangeType.hpp
@@ -3,7 +3,7 @@
 #include <type_traits>
 
 #include <dune/amdis/common/IndexSeq.hpp>
-#include <dune/amdis/common/MultiTypeFieldVector.hpp>
+#include <dune/amdis/common/MultiTypeVector.hpp>
 
 namespace AMDiS
 {
@@ -62,7 +62,7 @@ namespace AMDiS
       {
         template <std::size_t J>
         using ChildNode = typename Node::template Child<J>::type;
-        using type = MultiTypeFieldVector<RangeType_t<ChildNode<I>>...>;
+        using type = MultiTypeVector<RangeType_t<ChildNode<I>>...>;
       };
 
       using type = typename RangeTypeGenerator<MakeSeq_t<Node::CHILDREN>>::type;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f60187bf477173f08c553e66990023d7bb11be58
--- /dev/null
+++ b/test/CMakeLists.txt
@@ -0,0 +1,19 @@
+
+dune_add_test(SOURCES ConceptsTest.cpp
+  LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES ExpressionsTest.cpp
+  LINK_LIBRARIES amdis
+  CMD_ARGS "${CMAKE_SOURCE_DIR}/init/ellipt.dat.2d")
+
+dune_add_test(SOURCES FieldMatVecTest.cpp
+  LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES FilesystemTest.cpp
+  LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES MultiTypeVectorTest.cpp
+  LINK_LIBRARIES amdis)
+
+dune_add_test(SOURCES StringTest.cpp
+  LINK_LIBRARIES amdis)
diff --git a/test/ConceptsTest.cpp b/test/ConceptsTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3ce8b5e6e978861285d045d84fab665003529720
--- /dev/null
+++ b/test/ConceptsTest.cpp
@@ -0,0 +1,39 @@
+#include <array>
+#include <vector>
+
+#include <dune/common/reservedvector.hh>
+#include <dune/functions/functionspacebases/flatmultiindex.hh>
+
+#include <dune/amdis/common/Concepts.hpp>
+#include <dune/amdis/common/FieldMatVec.hpp>
+
+using namespace AMDiS;
+
+int main()
+{
+  using V = Dune::FieldVector<double,2>;
+  using M = Dune::FieldMatrix<double,2,2>;
+
+  // arithmetic concepts
+  static_assert( Concepts::Addable<V,V>, "" );
+  static_assert( Concepts::Addable<M,M>, "" );
+
+  static_assert( Concepts::Multiplicable<V,double>, "" );
+  static_assert( Concepts::Multiplicable<M,double>, "" );
+  // static_assert( Concepts::Multiplicable<M,V>, "" );
+
+  // Concepts::Callable
+  auto f = [](V const&) { return 0.0; };
+  auto b = [](V const&) { return true; };
+
+  static_assert( Concepts::Callable<decltype(f),V>, "" );
+  static_assert( Concepts::Functor<decltype(f),double(V)>, "" );
+  static_assert( Concepts::Functor<decltype(b),bool(V)>, "" );
+  static_assert( Concepts::Predicate<decltype(b),V>, "" );
+
+  // Concepts::MultiIndex
+  static_assert( Concepts::MultiIndex<Dune::ReservedVector<std::size_t,2>>, "" );
+  static_assert( Concepts::MultiIndex<Dune::Functions::FlatMultiIndex<std::size_t>>, "" );
+  static_assert( Concepts::MultiIndex<std::array<std::size_t,2>>, "" );
+  static_assert( Concepts::MultiIndex<std::vector<std::size_t>>, "" );
+}
diff --git a/test/ExpressionsTest.cpp b/test/ExpressionsTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..80db279b37f3eae1cddd45ba1ef83812253303a3
--- /dev/null
+++ b/test/ExpressionsTest.cpp
@@ -0,0 +1,108 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <iostream>
+
+#include <dune/amdis/AMDiS.hpp>
+#include <dune/amdis/ProblemStat.hpp>
+#include <dune/amdis/Operators.hpp>
+#include <dune/amdis/common/Literals.hpp>
+#include <dune/amdis/common/FieldMatVec.hpp>
+#include <dune/amdis/gridfunctions/Integrate.hpp>
+
+
+using namespace AMDiS;
+
+using ElliptParam   = YaspGridBasis<2, 2>;
+using ElliptProblem = ProblemStat<ElliptParam>;
+
+int main(int argc, char** argv)
+{
+  AMDiS::init(argc, argv);
+
+  using namespace Dune::Indices;
+
+  ElliptProblem prob("ellipt");
+  prob.initialize(INIT_ALL);
+
+  auto u = prob.getSolution(_0);
+
+  // eval a functor at global coordinates (at quadrature points)
+  auto expr1 = evalAtQP([](Dune::FieldVector<double, 2> const& x) { return x[0] + x[1]; });
+  auto expr2 = [](Dune::FieldVector<double, 2> const& x) { return x[0] + x[1]; };
+  auto expr3 = [](auto const& x) { return x[0] + x[1]; };
+
+  // constant values at quadrature points
+  auto expr4 = 1.0;
+  auto expr5 = Dune::FieldVector<double, 2>{1.0, 2.0};
+  auto expr6 = std::ref(expr4);
+
+  // Coordinate vector and component
+  auto expr7 = X();
+  auto expr8 = X(0);
+
+  // Evaluation of the DOFVector (component) at quadrature points
+  auto expr9 = prob.getSolution(_0);
+
+  // ---------------------------------
+
+  // derivative of expressions
+
+  auto diff4 = gradientAtQP(expr4);
+  // auto diff5 = gradientAtQP(expr5);
+  auto diff6 = gradientAtQP(expr6);
+
+  // auto diff7 = gradientAtQP(expr7);
+  auto diff8 = gradientAtQP(expr8);
+
+  auto diff9 = gradientAtQP(expr9);
+
+  // ---------------------------------
+
+  u.interpolate(expr1);
+  u.interpolate(expr2);
+  u.interpolate(expr3);
+
+  u.interpolate(expr4);
+  u.interpolate(expr6);
+
+  u.interpolate(expr8);
+
+  u.interpolate(expr9);
+
+  // ---------------------------------
+
+  // operations with expressions
+
+  auto op1 = expr1 + expr2;
+  auto op2 = expr1 * expr4;
+  auto op3 = two_norm(expr5);
+  auto op4 = min(expr6, expr8);
+  auto op5 = one_norm(expr7);
+
+  auto op6 = invokeAtQP([](double v) { return 2*v; }, op5);
+
+
+  u.interpolate(two_norm(diff4));
+  u.interpolate(two_norm(diff6));
+
+  u.interpolate(two_norm(diff8));
+
+  u.interpolate(two_norm(diff9));
+
+  // ---------------------------------
+
+  // integration of expressions
+
+  auto gv = u.basis().gridView();
+
+  DUNE_UNUSED auto int1 = integrate(op1, gv, 5);
+  DUNE_UNUSED auto int2 = integrate(op2, gv, 5);
+  DUNE_UNUSED auto int3 = integrate(op3, gv);
+  DUNE_UNUSED auto int4 = integrate(op4, gv, 5);
+  DUNE_UNUSED auto int5 = integrate(op5, gv, 5);
+  DUNE_UNUSED auto int6 = integrate(op6, gv, 5);
+
+  AMDiS::finalize();
+  return 0;
+}
diff --git a/test/FieldMatVecTest.cpp b/test/FieldMatVecTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bf890851752e4716068c44dce85db8c7248b4332
--- /dev/null
+++ b/test/FieldMatVecTest.cpp
@@ -0,0 +1,107 @@
+#include <cmath>
+
+#include <dune/amdis/common/FieldMatVec.hpp>
+#include <dune/amdis/common/Math.hpp>
+
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+// -----------------------------------------------------------------------------
+void test0()
+{
+  using V = FieldVector<double, 3>;
+  V a{1.0, 2.0, 3.0};
+  V b{2.0, 3.0, 4.0};
+
+  AMDIS_TEST_EQ( a+b, (V{3.0, 5.0, 7.0}) );
+  AMDIS_TEST_EQ( a-b, (V{-1.0, -1.0, -1.0}) );
+
+  AMDIS_TEST_EQ( a*2,   (V{2.0, 4.0, 6.0}) );
+  AMDIS_TEST_EQ( 2.0*a, (V{2.0, 4.0, 6.0}) );
+
+  AMDIS_TEST_EQ( b/2.0, (V{1.0, 1.5, 2.0}) );
+
+  AMDIS_TEST_EQ( a + 2*b, (V{5.0, 8.0, 11.0}) );
+}
+
+// -----------------------------------------------------------------------------
+void test1()
+{
+  using V = FieldVector<double, 3>;
+  V a{1.0, 2.0, 3.0};
+  V b{2.0, 3.0, 4.0};
+
+  AMDIS_TEST_EQ( unary_dot(a), 14.0 );
+  AMDIS_TEST_EQ( dot(a, b), 20.0 );
+  AMDIS_TEST_EQ( dot(a, b), dot(b, a) );
+
+  AMDIS_TEST_APPROX( two_norm(a), std::sqrt(14.0) );
+  AMDIS_TEST_APPROX( distance(a, b), std::sqrt(3.0) );
+
+  AMDIS_TEST_EQ( cross(a, b), (V{-1.0, 2.0, -1.0}) );
+}
+
+// -----------------------------------------------------------------------------
+void test2()
+{
+  using V = FieldVector<double, 3>;
+  using M = FieldMatrix<double, 3, 3>;
+
+  M A{ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} };
+  M B{ {1.0, 1.0, 1.0}, {2.0, 1.0, 0.0}, {3.0, 1.0, -1.0} };
+  M C{ {0.5, 2.0, -1.0}, {1.0, -1.0, 1.0}, {-0.25, 0.5, 4.0} };
+  M D{ {-0.25, 1.0, -1.0}, {2.0, -2.0, 8.0}, {-0.125, -0.125, 2.0} };
+
+  AMDIS_TEST_EQ( A+B,   (M{ {2.0, 1.0, 1.0}, { 2.0, 2.0, 0.0}, { 3.0, 1.0, 0.0} }) );
+  AMDIS_TEST_EQ( A-B,   (M{ {0.0,-1.0,-1.0}, {-2.0, 0.0, 0.0}, {-3.0,-1.0, 2.0} }) );
+
+  AMDIS_TEST_EQ( A*2,   (M{ {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0} }) );
+  AMDIS_TEST_EQ( 2.0*A, (M{ {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0} }) );
+
+  AMDIS_TEST_EQ( A/2.0, (M{ {0.5, 0.0, 0.0}, {0.0, 0.5, 0.0}, {0.0, 0.0, 0.5} }) );
+
+
+  V a{1.0, 2.0, 3.0};
+  V b{2.0, 3.0, 4.0};
+
+  // AMDIS_TEST_EQ( outer(a, b), (M{ {2.0, 3.0, 4.0}, {4.0, 6.0, 8.0}, {6.0, 9.0, 12.0} }) );
+  // AMDIS_TEST_EQ( diagonal(a), (M{ {1.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 3.0} }) );
+
+  // AMDIS_TEST_EQ( diagonal(A), (V{1.0, 1.0, 1.0}) );
+
+  // some tests for the determinant calculation
+  // ------------------------------------------
+  AMDIS_TEST_EQ( det(A), 1.0 );
+  AMDIS_TEST_EQ( det(B), 0.0 );
+
+  // AMDIS_TEST_EQ( det(C*D), det(C)*det(D) );
+  AMDIS_TEST_EQ( det(trans(C)), det(C) );
+  AMDIS_TEST_EQ( det(2.0*D), Math::pow<3>(2.0)*det(D) );
+
+  // some tests for the trace
+  // ------------------------
+  // AMDIS_TEST_EQ( trace(D), trace(trans(D)) );
+  // AMDIS_TEST_EQ( trace(2.0*C + 0.5*D), 2.0*trace(C) + 0.5*trace(D) );
+  // AMDIS_TEST_EQ( trace(C*D), trace(D*C) );
+  // AMDIS_TEST_EQ( trace(C), sum(diagonal(C)) );
+
+  // some tests for the norm
+  // -----------------------
+  // AMDIS_TEST_APPROX( sqr(frobenius_norm(B)), trace(hermitian(B)*B) );
+  // AMDIS_TEST( frobenius_norm(C*D) <= frobenius_norm(B)*frobenius_norm(C) );
+
+  // test the inverse
+  // ----------------
+  // AMDIS_TEST_EQ( inv(C), adj(C)/det(C) );
+}
+
+
+int main()
+{
+  test0();
+  test1();
+  test2();
+
+  return report_errors();
+}
diff --git a/test/FilesystemTest.cpp b/test/FilesystemTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0686769f31070fad6b38ae296025c54b15e3b5de
--- /dev/null
+++ b/test/FilesystemTest.cpp
@@ -0,0 +1,56 @@
+#include <dune/amdis/utility/Filesystem.hpp>
+
+#include "Tests.hpp"
+
+
+using namespace AMDiS;
+
+// ---------------------------------------------------------------------------------------
+void test0()
+{
+  using namespace AMDiS::filesystem;
+  AMDIS_TEST_EQ( path("a/b/c").parent_path(), path("a/b") );
+  AMDIS_TEST_EQ( path("a/b/c.txt").parent_path(), path("a/b") );
+  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( !path("a/b/c").is_absolute() );
+  AMDIS_TEST( !path("a\\b\\c").is_absolute() );
+#ifdef _WIN32
+  AMDIS_TEST( path("a:\\b\\c").is_absolute() );
+  AMDIS_TEST( !path("a:\\b\\c").is_relative() );
+#else
+  AMDIS_TEST( path("/a/b/c").is_absolute() );
+  AMDIS_TEST( !path("/a/b/c").is_relative() );
+  AMDIS_TEST( path("/").is_absolute() );
+#endif
+
+  AMDIS_TEST( path("a/b/c").is_relative() );
+
+  AMDIS_TEST_EQ( path("a/b/c.txt").remove_filename(), path("a/b") );
+  AMDIS_TEST_EQ( path("a/b/c").string(), "a/b/c" );
+
+  AMDIS_TEST_EQ( path("a/b/../c").string(), "a/c" );
+  AMDIS_TEST_EQ( path("a/b/./c").string(), "a/b/c" );
+  AMDIS_TEST_EQ( path("../a/b/c").string(), "../a/b/c" );
+  AMDIS_TEST_EQ( path("./a/b/c").string(), "a/b/c" );
+}
+
+void test1()
+{
+  using namespace AMDiS::filesystem;
+  AMDIS_TEST( exists( path("/tmp") ) );
+  AMDIS_TEST( exists("/etc/fstab") );
+}
+
+int main()
+{
+  test0();
+  test1();
+
+  return report_errors();
+}
diff --git a/test/MultiTypeVectorTest.cpp b/test/MultiTypeVectorTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dbfe6ba85008bc54c3d39ef70dc1ee5f4165603c
--- /dev/null
+++ b/test/MultiTypeVectorTest.cpp
@@ -0,0 +1,52 @@
+#include <dune/amdis/common/FieldMatVec.hpp>
+#include <dune/amdis/common/MultiTypeVector.hpp>
+
+using namespace AMDiS;
+
+int main()
+{
+  using V = MultiTypeVector<FieldVector<double,2>, FieldVector<double,3>>;
+
+  // constructors
+  V v0{}; // default constructor
+
+  FieldVector<double,2> f2 = {1.0, 2.0};
+  FieldVector<double,3> f3 = {1.0, 2.0, 3.0};
+
+  V v1{f2,f3}; // constructor with FV arguments
+  V v2{v1}; // copy constructor
+  v0 = v1; // copy assignment
+
+  V v3{0.0}; // construction from constant
+  v3 = 1.0; // assignment of constant
+
+  {
+    V v0_tmp;
+    V v1_tmp{std::move(v0_tmp)}; // move constructor
+    v1 = std::move(v1_tmp); // move assignment
+  }
+
+  using namespace Dune::Indices;
+
+  // element access
+  v1[_0] = f2; // mutable block access
+  v1[_0][0] = 2.0; // mutable element access
+  DUNE_UNUSED auto x = v1[_1][1]; // const access
+
+  // compound assignment operators
+  v1[_1] += f3;
+  v1 += v2;
+
+  v1 *= 2.0;
+  v1 /= 2.0;
+
+  // multi-index access
+  std::array<std::size_t,2> index{1,1};
+  v1[index] = 2.0;
+
+  using FV = double; //FieldVector<double,1>;
+  using W0 = MultiTypeVector<FV,FV>;
+  using W = MultiTypeVector<W0,W0>;
+  W w;
+  w[index] = 3.0;
+}
\ No newline at end of file
diff --git a/test/StringTest.cpp b/test/StringTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dabdb755eb639f813ead124af8a02566af1bb672
--- /dev/null
+++ b/test/StringTest.cpp
@@ -0,0 +1,52 @@
+#include <string>
+#include <vector>
+
+#include <dune/amdis/utility/String.hpp>
+#include "Tests.hpp"
+
+using namespace AMDiS;
+
+// ---------------------------------------------------------------------------------------
+void test0()
+{
+  AMDIS_TEST_EQ( to_upper("Hallo Welt"), std::string("HALLO WELT") );
+  AMDIS_TEST_EQ( to_lower("Hallo Welt"), std::string("hallo welt") );
+
+  AMDIS_TEST_EQ( trim_copy("  Hallo Welt  "), std::string("Hallo Welt") );
+}
+
+// ---------------------------------------------------------------------------------------
+// test the split() function
+void test1()
+{
+  std::vector<std::string> subs;
+  std::string text = "Dies ist ein langer Text";
+  split(text.begin(), text.end(), ' ', [&subs](auto it0, auto it1) { subs.push_back(std::string(it0, it1)); });
+  AMDIS_TEST_EQ( subs.size(), 5 );
+  AMDIS_TEST_EQ( subs[0], std::string("Dies") );
+
+  int n = 0;
+  std::string sep = " x";
+  split(text.begin(), text.end(), sep.begin(), sep.end(), [&n](auto,auto) { ++n; });
+  AMDIS_TEST_EQ( n, 6 );
+
+  n = 0;
+  std::string text0 = "";
+  split(text0.begin(), text0.end(), ' ', [&n](auto,auto) { ++n; });
+  AMDIS_TEST_EQ( n, 0 );
+
+  std::string text1 = "Hallo";
+  split(text1.begin(), text1.end(), ' ', [](auto it0, auto it1)
+  {
+    AMDIS_TEST_EQ( std::string(it0, it1), std::string("Hallo") );
+  });
+}
+
+
+int main()
+{
+  test0();
+  test1();
+
+  return report_errors();
+}
diff --git a/test/Tests.hpp b/test/Tests.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c44226e5ea51d3e2f445aefa95d878df2c88e169
--- /dev/null
+++ b/test/Tests.hpp
@@ -0,0 +1,153 @@
+#pragma once
+
+#include <iostream>
+
+// inspired by boost/core/lightweight_test.hpp - lightweight test library
+
+#define AMDIS_TEST(expr) \
+  ::AMDiS::Impl::_test(expr, #expr, __FILE__, __LINE__)
+
+#define AMDIS_TEST_EQ(expr, value) \
+  ::AMDiS::Impl::_test_eq(expr, value, #expr, #value, __FILE__, __LINE__)
+
+#define AMDIS_TEST_NE(expr, value) \
+  ::AMDiS::Impl::_test_ne(expr, value, #expr, #value, __FILE__, __LINE__)
+
+#define AMDIS_TEST_TOL 1.e-10
+#define AMDIS_TEST_APPROX(expr, value) \
+  ::AMDiS::Impl::_test_approx(expr, value, #expr, #value, __FILE__, __LINE__)
+
+#ifdef DEC_NO_THROW
+  #define AMDIS_TEST_THROWS(expr)
+#else
+  #define AMDIS_TEST_THROWS(expr) \
+    try { \
+      expr; \
+      ::AMDiS::Impl::_test_throws(#expr, __FILE__, __LINE__); \
+    } catch(...) {}
+#endif
+
+namespace AMDiS
+{
+  namespace Impl
+  {
+    /// global counter for the number of detected errors
+    inline int& num_errors()
+    {
+      static int num = 0;
+      return num;
+    }
+
+  } // end namespace Impl
+
+
+  /// Returns 0 if no errors are detected, otherwise returns 1. Must be called at the end
+  /// of main().
+  inline int report_errors()
+  {
+    int errors = Impl::num_errors();
+
+    if (errors == 0) {
+      std::cout  << "No errors detected.\n";
+      return 0;
+    } else {
+      std::cerr << errors << " error(s) detected.\n";
+      return 1;
+    }
+  }
+
+
+  namespace Impl
+  {
+    /// Tests whether an expression evaluates to true
+    inline void _test(bool success,
+                      char const* expr, char const* file, size_t line)
+    {
+      if (!success) {
+        std::cerr << file << ":" << line
+                 << "  TEST( " << expr << " ) failed\n";
+        num_errors()++;
+      }
+    }
+
+    /// Tests whether the value of an expression is equal to an expected value
+    template <class T1, class T2>
+    inline void _test_eq(T1 const& expr_value, T2 const& value,
+                         char const* expr_str, char const* value_str, char const* file, size_t line)
+    {
+      using T = std::common_type_t<T1,T2>;
+      if (T(expr_value) != T(value)) {
+        std::cerr << file << ":" << line
+                 << "  TEST( " << expr_str << " == " << value_str << " ) failed:  " << expr_value << " != " << value << "\n";
+        num_errors()++;
+      }
+    }
+
+    /// Tests whether the value of an expression is not equal to an expected value
+    template <class T1, class T2>
+    inline void _test_ne(T1 const& expr_value, T2 const& value,
+                         char const* expr_str, char const* value_str, char const* file, size_t line)
+    {
+      using T = std::common_type_t<T1,T2>;
+      if (T(expr_value) == T(value)) {
+        std::cerr << file << ":" << line
+                 << "  TEST( " << expr_str << " != " << value_str << " ) failed:  " << expr_value << " == " << value << "\n";
+        num_errors()++;
+      }
+    }
+
+    template <class T>
+    constexpr T _abs(T const& x) { return x < T(0) ? -x : x; }
+
+    /// Tests whether the value of an expression is approximately equal to an expected value
+    // implementation for scalars
+    template <class T1, class T2,
+      std::enable_if_t<std::is_arithmetic<T1>::value && std::is_arithmetic<T2>::value, int>* = nullptr>
+    inline void _test_approx(T1 const& expr_value, T2 const& value,
+                             char const* expr_str, char const* value_str, char const* file, size_t line)
+    {
+      if (_abs(expr_value - value) > AMDIS_TEST_TOL) {
+        std::cerr << file << ":" << line
+                 << "  TEST( " << expr_str << " ~= " << value_str << " ) failed:  " << expr_value << " != " << value << "\n";
+        num_errors()++;
+      }
+    }
+
+    // implementation for pair
+    template <class T11, class T12, class T21, class T22>
+    inline void _test_approx(std::pair<T11,T12> const& expr_value, std::pair<T21,T22> const& value,
+                             char const* expr_str, char const* value_str, char const* file, size_t line)
+    {
+      if (_abs(expr_value.first - value.first) > AMDIS_TEST_TOL ||
+          _abs(expr_value.second - value.second) > AMDIS_TEST_TOL)
+      {
+        std::cerr << file << ":" << line
+                 << "  TEST( " << expr_str << " ~= " << value_str << " ) failed:  (" << expr_value.first << "," << expr_value.second << ") != (" << value.first << "," << value.second << ")\n";
+        num_errors()++;
+      }
+    }
+
+    // implementation for ranges
+    template <class T1, class T2,
+      class = decltype(((void)std::declval<T1>().cbegin(), (void)std::declval<T2>().cbegin())) >
+    inline void _test_approx(T1 const& expr_range, T2 const& range,
+                             char const* expr_str, char const* value_str, char const* file, size_t line)
+    {
+      auto it1 = expr_range.cbegin();
+      auto it2 = range.cbegin();
+
+      for(; it1 != expr_range.cend(); ++it1, ++it2)
+        _test_approx(*it1, *it2, expr_str, value_str, file, line);
+    }
+
+    /// Tests whether an expression throws an exception
+    template <class T1, class T2>
+    inline void _test_throws(char const* expr, char const* file, size_t line)
+    {
+      std::cerr << file << ":" << line
+               << "  EXPR( " << expr << " ) should throw\n";
+      num_errors()++;
+    }
+  } // end namspace Impl
+
+} // end namespace Dec