diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index 818958f2..8043cb87 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -3087,13 +3087,14 @@ namespace netgen
 	      {
 		ElementIndex cnttet = IndexBASE<ElementIndex>();
 		int cntprism = 0;
-		for (int i = 1; i <= mesh.GetNE(); i++)
+		// for (int i = 1; i <= mesh.GetNE(); i++)
+                for (auto ei : mesh.VolumeElements().Range())
 		  {
-		    if (mesh.VolumeElement(i).GetType() == TET ||
-			mesh.VolumeElement(i).GetType() == TET10)
+		    if (mesh.VolumeElement(ei).GetType() == TET ||
+			mesh.VolumeElement(ei).GetType() == TET10)
 		      {
 			mtets[cnttet].marked =
-			  (opt.onlyonce ? 3 : 1) * mesh.VolumeElement(i).TestRefinementFlag();
+			  (opt.onlyonce ? 3 : 1) * mesh.VolumeElement(ei).TestRefinementFlag();
 			if (mtets[cnttet].marked)
 			  cntm++;
 			cnttet++;
@@ -3102,7 +3103,7 @@ namespace netgen
 		      {
 			cntprism++;
 			mprisms.Elem(cntprism).marked =
-			  2 * mesh.VolumeElement(i).TestRefinementFlag();
+			  2 * mesh.VolumeElement(ei).TestRefinementFlag();
 			if (mprisms.Elem(cntprism).marked)
 			  cntm++;
 		      }
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 0aeb3f87..90619969 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -994,7 +994,7 @@ namespace netgen
       NgArray<int> & segment_weights){ }
 #endif
 
-    NgArray<int> vol_partition;
+    Array<int, ElementIndex> vol_partition;
     NgArray<int> surf_partition;
     NgArray<int> seg_partition;
 
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 0341131f..44040fe3 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -362,7 +362,8 @@ namespace netgen
            
            meshing.GenerateMesh (mesh, mpquad);
            
-           for (int i = oldne + 1; i <= mesh.GetNE(); i++)
+           // for (int i = oldne + 1; i <= mesh.GetNE(); i++)
+           for (ElementIndex i : mesh.VolumeElements().Range().Modify(oldne, 0))
              mesh.VolumeElement(i).SetIndex (domain);
            
            (*testout) 
@@ -409,7 +410,8 @@ namespace netgen
 
       md.meshing->Delaunay (mesh, domain, mp);
 
-      for (int i = oldne + 1; i <= mesh.GetNE(); i++)
+      // for (int i = oldne + 1; i <= mesh.GetNE(); i++)
+      for (ElementIndex i : mesh.VolumeElements().Range().Modify(oldne, 0))
          mesh.VolumeElement(i).SetIndex (domain);
 
       PrintMessage (3, mesh.GetNP(), " points, ",
diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp
index 8770e07c..9c5f5556 100644
--- a/libsrc/meshing/refine.cpp
+++ b/libsrc/meshing/refine.cpp
@@ -854,9 +854,10 @@ namespace netgen
 
 		TBitArray<PointIndex> free (mesh.GetNP()), fhelp(mesh.GetNP());
 		free.Clear();
-		for (int i = 1; i <= mesh.GetNE(); i++)
+		// for (int i = 1; i <= mesh.GetNE(); i++)
+                for (ElementIndex ei : mesh.VolumeElements().Range())
 		  {
-		    const Element & el = mesh.VolumeElement(i);
+		    const Element & el = mesh.VolumeElement(ei);
 		    if (el.Volume(mesh.Points()) < 0)
 		      for (int j = 1; j <= el.GetNP(); j++)
 			free.SetBit (el.PNum(j));
diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp
index 47d3211c..5527b2a8 100644
--- a/libsrc/meshing/secondorder.cpp
+++ b/libsrc/meshing/secondorder.cpp
@@ -207,9 +207,10 @@ namespace netgen
 
 
     // refine volume elements
-    for (int i = 1; i <= mesh.GetNE(); i++)
+    // for (int i = 1; i <= mesh.GetNE(); i++)
+    for (ElementIndex ei : mesh.VolumeElements().Range())
       {
-	const Element & el = mesh.VolumeElement(i);
+	const Element & el = mesh.VolumeElement(ei);
 	int onp = 0;
 
 	Element newel(TET);
@@ -342,7 +343,7 @@ namespace netgen
 	      }
 	  }
 
-	mesh.VolumeElement (i) = newel;
+	mesh.VolumeElement (ei) = newel;
       }
 
 
@@ -430,9 +431,10 @@ namespace netgen
     for (int i = 1; i <= np; i++)
       parents.Elem(i) = INDEX_2(0,0);
 
-    for (int i = 1; i <= ne; i++)
+    // for (int i = 1; i <= ne; i++)
+    for (ElementIndex ei : mesh.VolumeElements().Range())
       {
-	const Element & el = mesh.VolumeElement(i);
+	const Element & el = mesh[ei];
 	if (el.GetType() == TET10)
 	  {
 	    static int betweentab[6][3] =
@@ -469,14 +471,15 @@ namespace netgen
 
     int cnttrials = 100;
     int wrongels = 0;
-    for (int i = 1; i <= ne; i++)
-      if (mesh.VolumeElement(i).CalcJacobianBadness (mesh.Points()) > 1e10)
+    // for (int i = 1; i <= ne; i++)
+    for (ElementIndex ei : mesh.VolumeElements().Range())
+      if (mesh.VolumeElement(ei).CalcJacobianBadness (mesh.Points()) > 1e10)
 	{
 	  wrongels++;
-	  mesh.VolumeElement(i).Flags().badel = 1;
+	  mesh.VolumeElement(ei).Flags().badel = 1;
 	}
       else
-	mesh.VolumeElement(i).Flags().badel = 0;
+	mesh.VolumeElement(ei).Flags().badel = 0;
 
     double facok = 0;
     double factry;