diff --git a/mpp b/mpp
index a343ae441ef32bdd9faba1b52e6e60ca1674d3fd..ee33eb4836a746c7a74b83248660aad95a1f4e16 160000
--- a/mpp
+++ b/mpp
@@ -1 +1 @@
-Subproject commit a343ae441ef32bdd9faba1b52e6e60ca1674d3fd
+Subproject commit ee33eb4836a746c7a74b83248660aad95a1f4e16
diff --git a/src/CardMechIProblem.hpp b/src/CardMechIProblem.hpp
index 1a53c04f484b9478c0291b98f86703201e63297f..a1d474e728555841c15766f36d3f5d60db057def 100644
--- a/src/CardMechIProblem.hpp
+++ b/src/CardMechIProblem.hpp
@@ -56,7 +56,7 @@ public:
 
   const Mesh &GetMesh(int level) const { return (*meshes)[level]; }
 
-  const Mesh &GetMesh(LevelPair level) const { return (*meshes)[level]; }
+  const Mesh &GetMesh(MeshIndex level) const { return (*meshes)[level]; }
 
   virtual std::string Name() const = 0;
 };
diff --git a/src/CardiacData.hpp b/src/CardiacData.hpp
index 148f2ed333fbd2e7de7ede8d31f30c1e07f292be..1bbaca31d046f7ba9e098128a18244c75e9016d4 100644
--- a/src/CardiacData.hpp
+++ b/src/CardiacData.hpp
@@ -114,19 +114,19 @@ void PlotFibreOrientation(const Vector &U);
 
 static double Amplitude(const DataContainer &data) { return data[2]; }
 
-static double Amplitude(const Vector &vec, row &r, int j = 0) { return vec(r, 2 * r.n() + j); }
+static double Amplitude(const Vector &vec, row &r, int j = 0) { return vec(r, 2 * r.NumberOfDofs() + j); }
 
 static double VectorAmplitude(const Vector &vec, row &r, int j = 0) {
   return vec(r,
-             2 * (r.n() / 3) +
+             2 * (r.NumberOfDofs() / 3) +
              j);
 }
 
 static double Duration(const DataContainer &data) { return data[1]; }
 
-static double Duration(const Vector &vec, row &r, int j = 0) { return vec(r, r.n() + j); }
+static double Duration(const Vector &vec, row &r, int j = 0) { return vec(r, r.NumberOfDofs() + j); }
 
-static double VectorDuration(const Vector &vec, row &r, int j = 0) { return vec(r, r.n() / 3 + j); }
+static double VectorDuration(const Vector &vec, row &r, int j = 0) { return vec(r, r.NumberOfDofs() / 3 + j); }
 
 static double Excitation(const DataContainer &data) { return data[0]; }
 
diff --git a/src/CardiacInterpolation.cpp b/src/CardiacInterpolation.cpp
index c1f89bed82a00173c0f470ae6537b55912d529a4..7256d6645bac1ed8e7ae3123472f23969ebddb0a 100644
--- a/src/CardiacInterpolation.cpp
+++ b/src/CardiacInterpolation.cpp
@@ -15,7 +15,7 @@ bool compare_excitation(const DataContainer &a, const DataContainer &b) {
 
 
 void fillLagrangeData(Vector &vec, row r, DataContainer d) {
-  int singleSize = r.n() / d.size();
+  int singleSize = r.NumberOfDofs() / d.size();
   if (singleSize==1){
     for (int i = 0; i < d.size(); ++i) {
       vec(r, i) = d[i];
@@ -29,12 +29,12 @@ void fillLagrangeData(Vector &vec, row r, DataContainer d) {
     vec(r,5) =0.0;
 
   }else{
-    for(int i = 0;i<r.n();++i){
+    for(int i = 0;i<r.NumberOfDofs();++i){
       vec(r,i) = d[i%d.size()];
     }
 
   }
-  /*int singleSize = r.n() / d.size();
+  /*int singleSize = r.NumberOfDofs() / d.size();
   for (int i = 0; i < d.size(); ++i) {
     for (int j = 0; j < singleSize; ++j)
       if (j == 1 && i == 0) {
diff --git a/src/coupled/problem/CoupledProblem.cpp b/src/coupled/problem/CoupledProblem.cpp
index 9657a51485abc3eb68c52d34161f37fba724e435..bfd47c218d11908cccc9f2b301a7781ae061e223 100644
--- a/src/coupled/problem/CoupledProblem.cpp
+++ b/src/coupled/problem/CoupledProblem.cpp
@@ -84,7 +84,7 @@ const Meshes &CoupledProblem::GetElphyMeshes() const { return *elphyMeshes; }
 
 const Mesh &CoupledProblem::GetElphyMesh(int level) const { return (*elphyMeshes)[level]; }
 
-const Mesh &CoupledProblem::GetElphyMesh(LevelPair level) const { return (*elphyMeshes)[level]; }
+const Mesh &CoupledProblem::GetElphyMesh(MeshIndex level) const { return (*elphyMeshes)[level]; }
 
 const CoarseGeometry &CoupledProblem::MechDomain() const { return *mechGeo; }
 
@@ -92,7 +92,7 @@ const Meshes &CoupledProblem::GetMechMeshes() const { return *mechMeshes; }
 
 const Mesh &CoupledProblem::GetMechMesh(int level) const { return (*mechMeshes)[level]; }
 
-const Mesh &CoupledProblem::GetMechMesh(LevelPair level) const { return (*mechMeshes)[level]; }
+const Mesh &CoupledProblem::GetMechMesh(MeshIndex level) const { return (*mechMeshes)[level]; }
 
 void CoupledProblem::InitializeEvaluationPoints() {
 
diff --git a/src/coupled/problem/CoupledProblem.hpp b/src/coupled/problem/CoupledProblem.hpp
index 58278945271fd900096e3b6d414115af28445035..e12614e1d67fac1c70c5c6362241dca73424d868 100644
--- a/src/coupled/problem/CoupledProblem.hpp
+++ b/src/coupled/problem/CoupledProblem.hpp
@@ -81,7 +81,7 @@ public:
 
   const Mesh &GetElphyMesh(int level) const;
 
-  const Mesh &GetElphyMesh(LevelPair level) const;
+  const Mesh &GetElphyMesh(MeshIndex level) const;
 
 
   const CoarseGeometry &MechDomain() const;;
@@ -90,7 +90,7 @@ public:
 
   const Mesh &GetMechMesh(int level) const;
 
-  const Mesh &GetMechMesh(LevelPair level) const;
+  const Mesh &GetMechMesh(MeshIndex level) const;
 
   virtual std::string Name() const override = 0;
   virtual void InitializeEvaluationPoints() override;
diff --git a/src/electrophysiology/problem/DiffusionProblem.cpp b/src/electrophysiology/problem/DiffusionProblem.cpp
index 82d53f826f0ac808c38135a06a0dab030c5a0fb2..6a0bb43b58f2a6e35113d4f6ffc885fba54b9d55 100644
--- a/src/electrophysiology/problem/DiffusionProblem.cpp
+++ b/src/electrophysiology/problem/DiffusionProblem.cpp
@@ -21,7 +21,7 @@ vector<double> DiffusionProblem::EvaluationResults(const Vector &solution) const
 
   Vector error(solution);
   for (row r = solution.rows(); r != solution.rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
       error(r, j) -= Potential(0.5, r());
     }
   }
diff --git a/src/electrophysiology/solvers/LinearImplicitSolver.cpp b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
index da62aad48f23645769864528160fb94a35d9a78c..178c7c52d6929e7652181c3a59237f0936a4f9b5 100644
--- a/src/electrophysiology/solvers/LinearImplicitSolver.cpp
+++ b/src/electrophysiology/solvers/LinearImplicitSolver.cpp
@@ -93,7 +93,7 @@ void LinearImplicitSolver::SolveGating(double t, double dt) {
   std::vector<double> vcw(M.Size());
   int position = M.GatingIndex();
   for (row r = (*gating)[0].rows(); r != (*gating)[0].rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
       for (int i = 0; i < M.Size(); ++i) {
         vcw[i] = (*gating)[i](r, j);
       }
@@ -114,7 +114,7 @@ void LinearImplicitSolver::SolveGatingOnCells(double t,double dt){
   std::vector<double> vcwPerCell(M.Size());
   int position = M.GatingIndex();
   for (row r = (*vcw_c)[0].rows(); r != (*vcw_c)[0].rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
       for (int i = 0; i < M.Size(); ++i) {
         vcwPerCell[i] = (*vcw_c)[i](r(), j);
       }
@@ -172,7 +172,7 @@ void LinearImplicitSolver::updateValues(Vectors &values) {
 
 void LinearImplicitSolver::updateMaxMin() {
   for (row r = (*maxgating).rows(); r != (*maxgating).rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
       for (int i = 0; i < M.Size(); i++) {
         if ((*gating)[i](r, j) < (*mingating)[i](r, j)) {
           (*mingating)[i](r, j) = (*gating)[i](r, j);
@@ -199,7 +199,7 @@ void LinearImplicitSolver::PrintMaxMinGating() {
     minGatingOverVertices[i] = 10000000;
   }
   for (row r = (*maxgating).rows(); r != (*maxgating).rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
       for (int i = 0; i < M.Size(); i++) {
         if ((*mingating)[i](r, j) < minGatingOverVertices[i]) {
           minGatingOverVertices[i] = (*mingating)[i](r, j);
@@ -233,7 +233,7 @@ void LinearImplicitSolver::SolveConcentration(IElphyAssemble &A, double dt, cons
     G[i] = std::vector<double>(M.Size());
   }
   for (row r = (*gating)[0].rows(); r != (*gating)[0].rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
       for (int i = 0; i < M.Size(); ++i) {
         vcw[i] = (*gating)[i](r, j);
       }
@@ -271,7 +271,7 @@ void LinearImplicitSolver::SolveConcentrationOnCells(IElphyAssemble &A, const Ve
       rate = 1 / A.GetElphyProblem().Jacobian_rate(*c);
     }
     //row r = (*vcw_c)[0].find_row(c());
-    //for (int j = 0; j < r.n(); ++j) {
+    //for (int j = 0; j < r.NumberOfDofs(); ++j) {
       for (int i = 0; i < M.Size(); ++i) {
         vcwPerCell[i] = (*vcw_c)[i](c(), 0);
       }
@@ -410,7 +410,7 @@ void SemiImplicitSolverOnCells:: Initialize(IElphyAssemble &A, Vector &V){
   Vector excitation( *potential);
   Vector duration( *potential);
   for (row r = (*potential).rows(); r != (*potential).rows_end(); ++r) {
-    for (int j = 0; j < r.n(); ++j) {
+    for (int j = 0; j < r.NumberOfDofs(); ++j) {
        amplitude(r,j)=Amplitude(A.ExcitationData(), r, j);
        excitation(r,j)=Excitation(A.ExcitationData(), r, j);
        duration(r,j)=Duration(A.ExcitationData(), r, j);
@@ -456,7 +456,7 @@ void SemiImplicitSolverOnCells::updateValues(Vectors &values, Vector &gamma_f_c_
   values[0] = *potential;
   if (M.Type() == TENSION){
       for (row r = (*vcw_c)[0].rows(); r != (*vcw_c)[0].rows_end(); ++r) {
-        for (int j = 0; j < r.n(); ++j) {
+        for (int j = 0; j < r.NumberOfDofs(); ++j) {
           gamma_f_c_pot(r(),j) = (*vcw_c)[M.ElphySize()](r(),j);
         }
       }