diff --git a/ProcessLib/LocalAssemblerTraits.h b/ProcessLib/LocalAssemblerTraits.h
new file mode 100644
index 0000000000000000000000000000000000000000..c2d8dfed11248ed024d3c34ee0b9e62e9e5049f3
--- /dev/null
+++ b/ProcessLib/LocalAssemblerTraits.h
@@ -0,0 +1,234 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_LOCALASSEMBLERTRAITS_H
+#define PROCESSLIB_LOCALASSEMBLERTRAITS_H
+
+namespace ProcessLib
+{
+
+namespace detail
+{
+
+/*! This class makes it possible to handle Eigen's block expressions transparently,
+ * both for fixed-size and dynamically allocated Eigen matrices.
+ *
+ * \tparam ShpPol the shape matrix policy used
+ * \tparam NIntPts the number of integration points for the FEM element
+ * \tparam NodalDOF the number of d.o.f. per node.
+ * If zero, this struct's methods will determine the size of
+ * the generated Eigen::block's at runtime. If non-zero, the
+ * size will be determined at compile-time.
+ * \tparam Dim global spatial dimension
+ */
+template<typename ShpPol, unsigned NIntPts, unsigned NodalDOF, unsigned Dim>
+struct LocalAssemblerTraitsFixed
+{
+private:
+ template<int N, int M>
+ using Matrix = typename ShpPol::template MatrixType<N, M>;
+ template<int N>
+ using Vector = typename ShpPol::template VectorType<N>;
+
+public:
+ using ShapeMatrices = typename ShpPol::ShapeMatrices;
+
+ //! Square matrix of the given space dimensionality
+ using MatrixDimDim = Matrix<Dim, Dim>;
+
+ // TODO That only works if all process variables are single component
+ // and use the same shape fucntions.
+ //! Local matrix for the given number of d.o.f.\ per node and number of
+ //! integration points
+ using LocalMatrix = Matrix<NIntPts*NodalDOF, NIntPts*NodalDOF>;
+ //! Local vector for the given number of d.o.f.\ per node and number of
+ //! integration points
+ using LocalVector = Vector<NIntPts*NodalDOF>;
+
+ //! Local vector for one component of one process variable.
+ //! The size is the number of nodes in the element.
+ using Vector1Comp = typename ShapeMatrices::ShapeType;
+
+ //! Laplace matrix for the given space dimensionality and number of d.o.f
+ //! per node.
+ using LaplaceMatrix = Matrix<Dim*NodalDOF, Dim*NodalDOF>;
+
+
+ //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF != 0,
+ decltype(std::declval<const Mat>().template block<Dim, Dim>(0u, 0u))
+ >::type
+ blockDimDim(Mat const& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows==Dim && ncols==Dim);
+ (void) nrows; (void) ncols;
+ return mat.template block<Dim, Dim>(top, left);
+ }
+ //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF == 0,
+ decltype(std::declval<const Mat>().block(0u, 0u, 0u, 0u))
+ >::type
+ blockDimDim(Mat const& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows == ncols);
+ return mat.block(top, left, nrows, ncols);
+ }
+ //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF != 0,
+ decltype(std::declval<Mat>().template block<Dim, Dim>(0u, 0u))
+ >::type
+ blockDimDim(Mat& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows==Dim && ncols==Dim);
+ (void) nrows; (void) ncols;
+ return mat.template block<Dim, Dim>(top, left);
+ }
+ //! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF == 0,
+ decltype(std::declval<Mat>().block(0u, 0u, 0u, 0u))
+ >::type
+ blockDimDim(Mat& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows == ncols);
+ return mat.block(top, left, nrows, ncols);
+ }
+
+ //! Get a block \c NIntPts x \c NIntPts whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF != 0,
+ decltype(std::declval<const Mat>().template block<NIntPts, NIntPts>(0u, 0u))
+ >::type
+ blockShpShp(Mat const& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows==NIntPts && ncols==NIntPts);
+ (void) nrows; (void) ncols;
+ return mat.template block<NIntPts, NIntPts>(top, left);
+ }
+ //! Get a block \c NIntPts x \c NIntPts whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF == 0,
+ decltype(std::declval<const Mat>().block(0u, 0u, 0u, 0u))
+ >::type
+ blockShpShp(Mat const& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows == ncols);
+ return mat.block(top, left, nrows, ncols);
+ }
+ //! Get a block \c NIntPts x \c NIntPts whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF != 0,
+ decltype(std::declval<Mat>().template block<NIntPts, NIntPts>(0u, 0u))
+ >::type
+ blockShpShp(Mat& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows==NIntPts && ncols==NIntPts);
+ (void) nrows; (void) ncols;
+ return mat.template block<NIntPts, NIntPts>(top, left);
+ }
+ //! Get a block \c NIntPts x \c NIntPts whose upper left corner is at \c top and \c left.
+ template<typename Mat>
+ static
+ typename std::enable_if<NodalDOF == 0,
+ decltype(std::declval<Mat>().block(0u, 0u, 0u, 0u))
+ >::type
+ blockShpShp(Mat& mat,
+ unsigned top, unsigned left, unsigned nrows, unsigned ncols)
+ {
+ assert(nrows == ncols);
+ return mat.block(top, left, nrows, ncols);
+ }
+
+ //! Get a block \c NIntPts x 1 starting at the \c top'th row.
+ template<typename Vec>
+ static
+ typename std::enable_if<NodalDOF != 0,
+ decltype(std::declval<const Vec>().template block<NIntPts, 1>(0u, 0u))
+ >::type
+ blockShp(Vec const& vec, unsigned top, unsigned nrows)
+ {
+ assert(nrows==NIntPts);
+ (void) nrows;
+ return vec.template block<NIntPts, 1>(top, 0);
+ }
+ //! Get a block \c NIntPts x 1 starting at the \c top'th row.
+ template<typename Vec>
+ static
+ typename std::enable_if<NodalDOF == 0,
+ decltype(std::declval<const Vec>().block(0u, 0u, 0u, 0u))
+ >::type
+ blockShp(Vec const& vec, unsigned top, unsigned nrows)
+ {
+ return vec.block(top, 0, nrows, 1);
+ }
+ //! Get a block \c NIntPts x 1 starting at the \c top'th row.
+ template<typename Vec>
+ static
+ typename std::enable_if<NodalDOF != 0,
+ decltype(std::declval<Vec>().template block<NIntPts, 1>(0u, 0u))
+ >::type
+ blockShp(Vec& vec, unsigned top, unsigned nrows)
+ {
+ assert(nrows==NIntPts);
+ (void) nrows;
+ return vec.template block<NIntPts, 1>(top, 0);
+ }
+ //! Get a block \c NIntPts x 1 starting at the \c top'th row.
+ template<typename Vec>
+ static
+ typename std::enable_if<NodalDOF == 0,
+ decltype(std::declval<Vec>().block(0u, 0u, 0u, 0u))
+ >::type
+ blockShp(Vec& vec, unsigned top, unsigned nrows)
+ {
+ return vec.block(top, 0, nrows, 1);
+ }
+};
+
+} // namespace detail
+
+
+#ifndef EIGEN_DYNAMIC_SHAPE_MATRICES
+
+template<typename ShpPol, unsigned NIntPts, unsigned NodalDOF, unsigned Dim>
+using LocalAssemblerTraits = detail::LocalAssemblerTraitsFixed<ShpPol, NIntPts, NodalDOF, Dim>;
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0,
+ "Inconsistent use of the macro OGS_EIGEN_DYNAMIC_SHAPE_MATRICES."
+ " Maybe you forgot to include some header file.");
+#else
+
+template<typename ShpPol, unsigned NIntPts, unsigned NodalDOF, unsigned Dim>
+using LocalAssemblerTraits = detail::LocalAssemblerTraitsFixed<ShapeMatrixPolicyType<void, 0>, 0, 0, 0>;
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1,
+ "Inconsistent use of the macro OGS_EIGEN_DYNAMIC_SHAPE_MATRICES."
+ " Maybe you forgot to include some header file.");
+#endif
+
+}
+
+#endif // PROCESSLIB_LOCALASSEMBLERTRAITS_H