diff --git a/MathLib/LinAlg/MatrixVectorTraits.cpp b/MathLib/LinAlg/MatrixVectorTraits.cpp
index cba84eff4fbedb4c5d540581a82f0f3d1729c876..d1036ed4347d240bd323dd387d49741e5f10e2c3 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.cpp
+++ b/MathLib/LinAlg/MatrixVectorTraits.cpp
@@ -7,6 +7,7 @@
  *
  */
 
+#include "AssemblerLib/LocalToGlobalIndexMap.h"
 #include "MatrixVectorTraits.h"
 
 #ifdef OGS_USE_EIGEN
@@ -32,7 +33,10 @@ std::unique_ptr<Eigen::MatrixXd>
 MatrixVectorTraits<Eigen::MatrixXd>::
 newInstance(MatrixSpecifications const& spec)
 {
-    return std::unique_ptr<Eigen::MatrixXd>(new Eigen::MatrixXd(spec.nrows, spec.ncols));
+    auto const nrows = spec.dof_table ? spec.dof_table->dofSize() : spec.nrows;
+    auto const ncols = spec.dof_table ? nrows : spec.ncols;
+
+    return std::unique_ptr<Eigen::MatrixXd>(new Eigen::MatrixXd(nrows, ncols));
 }
 
 std::unique_ptr<Eigen::VectorXd>
@@ -53,7 +57,9 @@ std::unique_ptr<Eigen::VectorXd>
 MatrixVectorTraits<Eigen::VectorXd>::
 newInstance(MatrixSpecifications const& spec)
 {
-    return std::unique_ptr<Eigen::VectorXd>(new Eigen::VectorXd(spec.nrows));
+    auto const nrows = spec.dof_table ? spec.dof_table->dofSize() : spec.nrows;
+
+    return std::unique_ptr<Eigen::VectorXd>(new Eigen::VectorXd(nrows));
 }
 
 } // namespace MathLib
@@ -63,6 +69,8 @@ newInstance(MatrixSpecifications const& spec)
 
 #ifdef USE_PETSC
 
+#include "MeshLib/NodePartitionedMesh.h"
+
 namespace MathLib
 {
 
@@ -70,7 +78,7 @@ std::unique_ptr<PETScMatrix>
 MatrixVectorTraits<PETScMatrix>::
 newInstance()
 {
-    return std::unique_ptr<PETScMatrix>(new PETScMatrix(0, 0)); // TODO default constructor
+    return std::unique_ptr<PETScMatrix>(new PETScMatrix);
 }
 
 std::unique_ptr<PETScMatrix>
@@ -84,7 +92,32 @@ std::unique_ptr<PETScMatrix>
 MatrixVectorTraits<PETScMatrix>::
 newInstance(MatrixSpecifications const& spec)
 {
-    return std::unique_ptr<PETScMatrix>(new PETScMatrix(spec.nrows)); // TODO sparsity pattern
+    auto const nrows = spec.dof_table ? spec.dof_table->dofSizeLocal() : spec.nrows;
+    auto const ncols = spec.dof_table ? nrows : spec.ncols;
+
+    // TODO I guess it is not hard to make AssemblerLib::computeSparsityPattern()
+    //      work also for NodePartitionedMesh'es.
+    assert(((!spec.sparsity_pattern) || spec.sparsity_pattern->size() == 0u) &&
+           "The sparsity pattern is not used with PETSc so I rather crash than"
+           " silently accept a precomputed sparsity pattern.");
+
+    if (spec.mesh)
+    {
+        assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(spec.mesh));
+        auto const& mesh = *static_cast<MeshLib::NodePartitionedMesh const*>(spec.mesh);
+
+        auto const max_nonzeroes = mesh.getMaximumNConnectedNodesToNode();
+
+        PETScMatrixOption mat_opt;
+        mat_opt.d_nz = max_nonzeroes;
+        mat_opt.o_nz = max_nonzeroes;
+        mat_opt.is_global_size = false;
+        return std::unique_ptr<PETScMatrix>(
+            new PETScMatrix(nrows, ncols, mat_opt));
+    }
+    else
+        return std::unique_ptr<PETScMatrix>(
+            new PETScMatrix(nrows, ncols));
 }
 
 std::unique_ptr<PETScVector>
@@ -105,7 +138,16 @@ std::unique_ptr<PETScVector>
 MatrixVectorTraits<PETScVector>::
 newInstance(MatrixSpecifications const& spec)
 {
-    return std::unique_ptr<PETScVector>(new PETScVector(spec.nrows));
+    auto const is_global_size = false;
+
+    if (spec.dof_table) {
+        auto const& dt = *spec.dof_table;
+        return std::unique_ptr<PETScVector>(
+            new PETScVector(dt.dofSizeLocal(), dt.getGhostIndices(), is_global_size));
+    } else {
+        return std::unique_ptr<PETScVector>(
+            new PETScVector(spec.nrows, is_global_size));
+    }
 }
 
 } // namespace MathLib
@@ -134,7 +176,14 @@ std::unique_ptr<EigenMatrix>
 MatrixVectorTraits<EigenMatrix>::
 newInstance(MatrixSpecifications const& spec)
 {
-    return std::unique_ptr<EigenMatrix>(new EigenMatrix(spec.nrows)); // TODO sparsity pattern
+    auto const nrows = spec.dof_table ? spec.dof_table->dofSize() : spec.nrows;
+
+    auto A = std::unique_ptr<EigenMatrix>(new EigenMatrix(nrows));
+
+    if (spec.sparsity_pattern)
+        setMatrixSparsity(*A, *spec.sparsity_pattern);
+
+    return A;
 }
 
 std::unique_ptr<EigenVector>
@@ -155,7 +204,9 @@ std::unique_ptr<EigenVector>
 MatrixVectorTraits<EigenVector>::
 newInstance(MatrixSpecifications const& spec)
 {
-    return std::unique_ptr<EigenVector>(new EigenVector(spec.nrows));
+    auto const nrows = spec.dof_table ? spec.dof_table->dofSize() : spec.nrows;
+
+    return std::unique_ptr<EigenVector>(new EigenVector(nrows));
 }
 
 } // namespace MathLib