diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp
index 5ed9efc2c..7108e0390 100644
--- a/opm/simulators/wells/StandardWellEquations.cpp
+++ b/opm/simulators/wells/StandardWellEquations.cpp
@@ -35,6 +35,72 @@ StandardWellEquations(const ParallelWellInfo& parallel_well_info)
     invDuneD_.setBuildMode(DiagMatWell::row_wise);
 }
 
+template<class Scalar, int numEq>
+void StandardWellEquations<Scalar,numEq>::
+init(const int num_cells,
+     const int numWellEq,
+     const int numPerfs,
+     const std::vector<int>& cells)
+{
+    // setup sparsity pattern for the matrices
+    //[A C^T    [x    =  [ res
+    // B D] x_well]      res_well]
+    // set the size of the matrices
+    duneD_.setSize(1, 1, 1);
+    duneB_.setSize(1, num_cells, numPerfs);
+    duneC_.setSize(1, num_cells, numPerfs);
+
+    for (auto row = duneD_.createbegin(),
+              end = duneD_.createend(); row != end; ++row) {
+        // Add nonzeros for diagonal
+        row.insert(row.index());
+    }
+      // the block size is run-time determined now
+    duneD_[0][0].resize(numWellEq, numWellEq);
+
+    for (auto row = duneB_.createbegin(),
+              end = duneB_.createend(); row != end; ++row) {
+        for (int perf = 0 ; perf < numPerfs; ++perf) {
+            const int cell_idx = cells[perf];
+            row.insert(cell_idx);
+        }
+    }
+
+    for (int perf = 0 ; perf < numPerfs; ++perf) {
+        const int cell_idx = cells[perf];
+        // the block size is run-time determined now
+        duneB_[0][cell_idx].resize(numWellEq, numEq);
+    }
+
+         // make the C^T matrix
+    for (auto row = duneC_.createbegin(),
+              end = duneC_.createend(); row != end; ++row) {
+        for (int perf = 0; perf < numPerfs; ++perf) {
+            const int cell_idx = cells[perf];
+            row.insert(cell_idx);
+        }
+    }
+
+    for (int perf = 0; perf < numPerfs; ++perf) {
+        const int cell_idx = cells[perf];
+        duneC_[0][cell_idx].resize(numWellEq, numEq);
+    }
+
+    resWell_.resize(1);
+    // the block size of resWell_ is also run-time determined now
+    resWell_[0].resize(numWellEq);
+
+    // resize temporary class variables
+    Bx_.resize(duneB_.N());
+    for (unsigned i = 0; i < duneB_.N(); ++i) {
+        Bx_[i].resize(numWellEq);
+    }
+
+    invDrw_.resize(duneD_.N());
+    for (unsigned i = 0; i < duneD_.N(); ++i) {
+        invDrw_[i].resize(numWellEq);
+    }
+}
 
 template<class Scalar, int numEq>
 void StandardWellEquations<Scalar,numEq>::clear()
diff --git a/opm/simulators/wells/StandardWellEquations.hpp b/opm/simulators/wells/StandardWellEquations.hpp
index a7eeaa2f4..066aea6cf 100644
--- a/opm/simulators/wells/StandardWellEquations.hpp
+++ b/opm/simulators/wells/StandardWellEquations.hpp
@@ -60,6 +60,16 @@ public:
 
     StandardWellEquations(const ParallelWellInfo& parallel_well_info);
 
+    //! \brief Setup sparsity pattern for the matrices.
+    //! \param num_cells Total number of cells
+    //! \param numWellEq Number of well equations
+    //! \param numPerfs Number of perforations
+    //! \param cells Cell indices for perforations
+    void init(const int num_cells,
+              const int numWellEq,
+              const int numPerfs,
+              const std::vector<int>& cells);
+
     //! \brief Set all coefficients to 0.
     void clear();
 
diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp
index a5c5da6a3..374d08088 100644
--- a/opm/simulators/wells/StandardWellEval.cpp
+++ b/opm/simulators/wells/StandardWellEval.cpp
@@ -1040,63 +1040,8 @@ init(std::vector<double>& perf_depth,
     primary_variables_evaluation_.resize(numWellEq_, EvalWell{numWellEq_ + Indices::numEq, 0.0});
 
     // setup sparsity pattern for the matrices
-    //[A C^T    [x    =  [ res
-    // B D] x_well]      res_well]
-    // set the size of the matrices
-    this->linSys_.duneD_.setSize(1, 1, 1);
-    this->linSys_.duneB_.setSize(1, num_cells, baseif_.numPerfs());
-    this->linSys_.duneC_.setSize(1, num_cells, baseif_.numPerfs());
-
-    for (auto row = this->linSys_.duneD_.createbegin(),
-              end = this->linSys_.duneD_.createend(); row != end; ++row) {
-        // Add nonzeros for diagonal
-        row.insert(row.index());
-    }
-    // the block size is run-time determined now
-    this->linSys_.duneD_[0][0].resize(numWellEq_, numWellEq_);
-
-    for (auto row = this->linSys_.duneB_.createbegin(),
-              end = this->linSys_.duneB_.createend(); row != end; ++row) {
-        for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
-            const int cell_idx = baseif_.cells()[perf];
-            row.insert(cell_idx);
-        }
-    }
-
-    for (int perf = 0 ; perf < baseif_.numPerfs(); ++perf) {
-        const int cell_idx = baseif_.cells()[perf];
-         // the block size is run-time determined now
-         this->linSys_.duneB_[0][cell_idx].resize(numWellEq_, Indices::numEq);
-    }
-
-    // make the C^T matrix
-    for (auto row = this->linSys_.duneC_.createbegin(),
-              end = this->linSys_.duneC_.createend(); row != end; ++row) {
-        for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
-            const int cell_idx = baseif_.cells()[perf];
-            row.insert(cell_idx);
-        }
-    }
-
-    for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
-        const int cell_idx = baseif_.cells()[perf];
-        this->linSys_.duneC_[0][cell_idx].resize(numWellEq_, Indices::numEq);
-    }
-
-    this->linSys_.resWell_.resize(1);
-    // the block size of resWell_ is also run-time determined now
-    this->linSys_.resWell_[0].resize(numWellEq_);
-
-    // resize temporary class variables
-    this->linSys_.Bx_.resize(this->linSys_.duneB_.N());
-    for (unsigned i = 0; i < this->linSys_.duneB_.N(); ++i) {
-        this->linSys_.Bx_[i].resize(numWellEq_);
-    }
-
-    this->linSys_.invDrw_.resize(this->linSys_.duneD_.N());
-    for (unsigned i = 0; i < this->linSys_.duneD_.N(); ++i) {
-        this->linSys_.invDrw_[i].resize(numWellEq_);
-    }
+    this->linSys_.init(num_cells, this->numWellEq_,
+                       baseif_.numPerfs(), baseif_.cells());
 }
 
 template<class FluidSystem, class Indices, class Scalar>