diff --git a/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h b/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h
index 7771993294a74c0b32eff9836791af65cb14734e..eeb5ee23787a6e716fe830ace6048e1bd2186bf1 100644
--- a/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h
+++ b/NumLib/Fem/CoordinatesMapping/ShapeMatrices-impl.h
@@ -103,24 +103,6 @@ inline void ShapeMatrices<T_N, T_DNDR, T_J, T_DNDX>::setZero()
     detail::setZero(*this, detail::ShapeDataFieldType<T_SHAPE_MATRIX_TYPE>());
 }
 
-template <class T_N, class T_DNDR, class T_J, class T_DNDX>
-void ShapeMatrices<T_N, T_DNDR, T_J, T_DNDX>::resize(std::size_t const dim, std::size_t n_nodes)
-{
-    resize(dim, dim, n_nodes);
-}
-
-template <class T_N, class T_DNDR, class T_J, class T_DNDX>
-void ShapeMatrices<T_N, T_DNDR, T_J, T_DNDX>::resize(std::size_t const local_dim, std::size_t const global_dim, std::size_t n_nodes)
-{
-    N.resize(n_nodes);
-    dNdr.resize(local_dim, n_nodes);
-    J.resize(local_dim, local_dim);
-    invJ.resize(local_dim, local_dim);
-    dNdx.resize(global_dim, n_nodes);
-
-    setZero();
-}
-
 template <class T_N, class T_DNDR, class T_J, class T_DNDX>
 void ShapeMatrices<T_N, T_DNDR, T_J, T_DNDX>::write(std::ostream& out) const
 {
diff --git a/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h b/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
index 4e8363e04cbc25bc39e22287dc6fcb8d12d0bc84..a467ccec811207b154f2408e3df3613a702bc49b 100644
--- a/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
+++ b/NumLib/Fem/CoordinatesMapping/ShapeMatrices.h
@@ -58,11 +58,14 @@ struct ShapeMatrices
     JacobianType invJ;  ///< Inverse matrix of the Jacobian
     DxShapeType dNdx;   ///< Matrix of gradient of shape functions in physical coordinates, dN(r)/dx
 
-    /** The default constructor is used by fixed-size (at compile-time)
-     * matrix/vector types where no resizing of the matrices/vectors is required
-     * in this constructor.
+    /** Not default constructible, dimensions always must be given.
+     *
+     * The default constructor has been deleted explicitly, because with
+     * dynamically allocated matrices it is rather easy to forget the
+     * required <tt>resize()</tt> call. Note: the <tt>resize()</tt> member
+     * is also deleted now.
      */
-    ShapeMatrices() = default;
+    ShapeMatrices() = delete;
 
     /**
      * Initialize matrices and vectors
@@ -82,21 +85,9 @@ struct ShapeMatrices
           invJ(local_dim, local_dim),
           dNdx(global_dim, n_nodes)
     {
-        this->setZero();
+        setZero();
     }
 
-    ~ShapeMatrices() {}
-
-    /// Reinitialize the matrices and vectors. When using dynamic size matrices
-    /// and the default constructor of this class memory allocation must be
-    /// performed calling this function.
-    /// For fixed size matrices no memory reallocation happens and the
-    /// matrix/vector sizes must be the same as at construction (given by the
-    /// template parameters).
-    void resize(std::size_t const dim, std::size_t n_nodes);
-
-    void resize(std::size_t const local_dim, std::size_t const global_dim, std::size_t n_nodes);
-
     /// reset all data with zero
     void setZero();
 
diff --git a/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp b/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bf65581b0c9e033f2c14bbdfca536f9f492b5a47
--- /dev/null
+++ b/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp
@@ -0,0 +1,92 @@
+/**
+ * \file ShapeStaticConsts.cpp
+ *
+ * \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
+ *
+ * The purpose of this file is to prevent linker errors.
+ *
+ * The classes defined in the included headers each have <tt>static const</tt> members
+ * \c DIM and \c NPOINTS. The compiler is allowed to optimize those symbols out.
+ * Therefore the linker later on might not find those symbols anymore and raises errors.
+ * The <tt>static const</tt> members are explicitly instantiated here in order to rule
+ * such errors out.
+ *
+ * Helper to generate that code:
+ *
+ * \code{.unparsed}
+ * find . \( -name '*.h' -and -not -name '*-impl.h' \) -printf '%f\n' \
+ * | sort \
+ * | while read f; do
+ *     c="${f%.h}";
+ *     echo "const std::size_t $c::DIM;";
+ *     echo "const std::size_t $c::NPOINTS;";
+ * done >ShapeStaticConsts.cpp
+ * \endcode
+ *
+ */
+
+#include "ShapeHex20.h"
+#include "ShapeHex8.h"
+#include "ShapeLine2.h"
+#include "ShapeLine3.h"
+#include "ShapePoint1.h"
+#include "ShapePrism15.h"
+#include "ShapePrism6.h"
+#include "ShapePyra13.h"
+#include "ShapePyra5.h"
+#include "ShapeQuad4.h"
+#include "ShapeQuad8.h"
+#include "ShapeQuad9.h"
+#include "ShapeTet10.h"
+#include "ShapeTet4.h"
+#include "ShapeTri3.h"
+#include "ShapeTri6.h"
+
+namespace NumLib
+{
+
+const std::size_t ShapeHex20::DIM;
+const std::size_t ShapeHex20::NPOINTS;
+const std::size_t ShapeHex8::DIM;
+const std::size_t ShapeHex8::NPOINTS;
+
+const std::size_t ShapeLine2::DIM;
+const std::size_t ShapeLine2::NPOINTS;
+const std::size_t ShapeLine3::DIM;
+const std::size_t ShapeLine3::NPOINTS;
+
+const std::size_t ShapePoint1::DIM;
+const std::size_t ShapePoint1::NPOINTS;
+
+const std::size_t ShapePrism15::DIM;
+const std::size_t ShapePrism15::NPOINTS;
+const std::size_t ShapePrism6::DIM;
+const std::size_t ShapePrism6::NPOINTS;
+
+const std::size_t ShapePyra13::DIM;
+const std::size_t ShapePyra13::NPOINTS;
+const std::size_t ShapePyra5::DIM;
+const std::size_t ShapePyra5::NPOINTS;
+
+const std::size_t ShapeQuad4::DIM;
+const std::size_t ShapeQuad4::NPOINTS;
+const std::size_t ShapeQuad8::DIM;
+const std::size_t ShapeQuad8::NPOINTS;
+const std::size_t ShapeQuad9::DIM;
+const std::size_t ShapeQuad9::NPOINTS;
+
+const std::size_t ShapeTet10::DIM;
+const std::size_t ShapeTet10::NPOINTS;
+const std::size_t ShapeTet4::DIM;
+const std::size_t ShapeTet4::NPOINTS;
+
+const std::size_t ShapeTri3::DIM;
+const std::size_t ShapeTri3::NPOINTS;
+const std::size_t ShapeTri6::DIM;
+const std::size_t ShapeTri6::NPOINTS;
+
+}
diff --git a/ProcessLib/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlowFEM.h
index b2cca22b1a84a39212c615df0f07239d3ef7f88b..748f13f728655a1ea9a445e0e7b34d1c0fe2207b 100644
--- a/ProcessLib/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlowFEM.h
@@ -73,9 +73,10 @@ public:
         IntegrationMethod_ integration_method(_integration_order);
         std::size_t const n_integration_points = integration_method.getNPoints();
 
-        _shape_matrices.resize(n_integration_points);
+        _shape_matrices.reserve(n_integration_points);
         for (std::size_t ip(0); ip < n_integration_points; ip++) {
-            _shape_matrices[ip].resize(ShapeFunction::DIM, ShapeFunction::NPOINTS);
+            _shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
+                                         ShapeFunction::NPOINTS);
             fe.computeShapeFunctions(
                     integration_method.getWeightedPoint(ip).getCoords(),
                     _shape_matrices[ip]);
diff --git a/ProcessLib/NeumannBcAssembler.h b/ProcessLib/NeumannBcAssembler.h
index e4fd38aa35410dc821e3c97d0d000700785cc208..064288059a4f8390b7134f85f8805d8a305b8f06 100644
--- a/ProcessLib/NeumannBcAssembler.h
+++ b/ProcessLib/NeumannBcAssembler.h
@@ -70,9 +70,10 @@ public:
         IntegrationMethod_ integration_method(_integration_order);
         std::size_t const n_integration_points = integration_method.getNPoints();
 
-        _shape_matrices.resize(n_integration_points);
+        _shape_matrices.reserve(n_integration_points);
         for (std::size_t ip(0); ip < n_integration_points; ip++) {
-            _shape_matrices[ip].resize(ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+            _shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
+                                         ShapeFunction::NPOINTS);
             fe.computeShapeFunctions(
                     integration_method.getWeightedPoint(ip).getCoords(),
                     _shape_matrices[ip]);
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
index 1d2335e32a4802d766ec3d4eaca6599130a05a02..dd43128f41867a126e6a72b594343461cddad83c 100644
--- a/Tests/NumLib/TestFunctionInterpolation.cpp
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -67,7 +67,8 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
     const unsigned integration_order = 2;
     IntegrationMethod integration_method(integration_order);
 
-    ShapeMatricesType::ShapeMatrices shape_matrix;
+    ShapeMatricesType::ShapeMatrices shape_matrix(
+        ShapeFunction::DIM, ShapeFunction::DIM, ShapeFunction::NPOINTS);
 
     finite_element.computeShapeFunctions(
             integration_method.getWeightedPoint(0).getCoords(),