[NL] Specialize storage of Nx1 and 0xM matrices for eigen.

Currently Eigen-3.2.5 can store Nx1 matrices only
in column major form.

0xM fixed Eigen matrices lead to 0-size arrays which
cause compilation errors (e.g. on MSVC++)

NumLib/Fem/ShapeMatrixPolicy.h

@@ -13,17 +13,54 @@
 #include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
 
 #ifdef OGS_USE_EIGEN
-#include <Eigen/Eigen>
+#include <Eigen/Dense>
+
+namespace detail
+{
+    /// Forwards the Eigen::Matrix type for general N and M.
+    /// There is a partial specialization for M = 1 to store the matrix in
+    /// column major storage order.
+    template <int N, int M>
+    struct EigenMatrixType
+    {
+        using type = Eigen::Matrix<double, N, M, Eigen::RowMajor>;
+    };
+
+    /// Specialization for Nx1 matrices which can be stored only in column major
+    /// form in Eigen-3.2.5.
+    template <int N>
+    struct EigenMatrixType<N, 1>
+    {
+        using type = Eigen::Matrix<double, N, 1, Eigen::ColMajor>;
+    };
+
+    /// Specialization for 0xM matrices. Using fixed size Eigen matrices here
+    /// would lead to zero sized arrays, which cause compilation errors on
+    /// some compilers.
+    template <int M>
+    struct EigenMatrixType<0, M>
+    {
+        using type = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+    };
+
+    template <>
+    struct EigenMatrixType<0, 1>
+    {
+        using type = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
+    };
+
+}   // detail
 
 /// An implementation of ShapeMatrixPolicy using fixed size (compile-time) eigen
 /// matrices and vectors.
 template <typename ShapeFunction, unsigned GlobalDim>
 struct EigenFixedShapeMatrixPolicy
 {
-    template <int N, int M>
-    using _MatrixType = Eigen::Matrix<double, N, M, Eigen::RowMajor>;
     template <int N>
-    using _VectorType = Eigen::Matrix<double, N, 1>;
+    using _VectorType = typename ::detail::EigenMatrixType<N, 1>::type;
+
+    template <int N, int M>
+    using _MatrixType = typename ::detail::EigenMatrixType<N, M>::type;
 
     using NodalMatrixType = _MatrixType<ShapeFunction::NPOINTS, ShapeFunction::NPOINTS>;
     using NodalVectorType = _VectorType<ShapeFunction::NPOINTS>;