diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp
index eaced993..68cac826 100644
--- a/libsrc/csg/edgeflw.cpp
+++ b/libsrc/csg/edgeflw.cpp
@@ -656,9 +656,10 @@ namespace netgen
       }
     
 
-    for (int i = 1; i <= nseg; i++)
+    // for (int i = 1; i <= nseg; i++)
+    for (Segment & seg : mesh.LineSegments())
       {
-	Segment & seg = mesh.LineSegment (i);
+	// Segment & seg = mesh.LineSegment (i);
 	if (seg.edgenr >= 1 && seg.edgenr <= cntedge)
 	  {
 	    if (osedges.Get(seg.edgenr) != -1)
diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp
index 7e656b17..c19862b0 100644
--- a/libsrc/meshing/clusters.cpp
+++ b/libsrc/meshing/clusters.cpp
@@ -96,7 +96,7 @@ namespace netgen
              
              nnums.SetSize(elnv+elned+elnfa+1);
              for (int j = 1; j <= elnv; j++)
-               nnums.Elem(j) = el.PNum(j);
+               nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
              for (int j = 1; j <= elned; j++)
                nnums.Elem(elnv+j) = nv+ednums.Elem(j);
              for (int j = 1; j <= elnfa; j++)
@@ -108,7 +108,6 @@ namespace netgen
            }
        });
     
-
     NgProfiler::StopTimer(timer1);
     NgProfiler::StartTimer(timer2);      
 
@@ -125,7 +124,7 @@ namespace netgen
 	  
 	nnums.SetSize(elnv+elned+1);
 	for (int j = 1; j <= elnv; j++)
-	  nnums.Elem(j) = el.PNum(j);
+	  nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
 	for (int j = 1; j <= elned; j++)
 	  nnums.Elem(elnv+j) = nv+ednums.Elem(j);
 	nnums.Elem(elnv+elned+1) = fanum;
@@ -197,7 +196,7 @@ namespace netgen
 	  
 	    nnums.SetSize(elnv+elned+elnfa+1);
 	    for (int j = 1; j <= elnv; j++)
-	      nnums.Elem(j) = el.PNum(j);
+	      nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
 	    for (int j = 1; j <= elned; j++)
 	      nnums.Elem(elnv+j) = nv+ednums.Elem(j);
 	    for (int j = 1; j <= elnfa; j++)
@@ -220,23 +219,23 @@ namespace netgen
 		break;
 	      case TET:
 	      case TET10:
-		if (cluster_reps.Get(el.PNum(1)) == 
-		    cluster_reps.Get(el.PNum(2)))
+		if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == 
+		    cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE))
 		  clustertab = tet_cluster12;
-		else if (cluster_reps.Get(el.PNum(1)) == 
-			 cluster_reps.Get(el.PNum(3)))
+		else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == 
+			 cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
 		  clustertab = tet_cluster13;
-		else if (cluster_reps.Get(el.PNum(1)) == 
-			 cluster_reps.Get(el.PNum(4)))
+		else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == 
+			 cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
 		  clustertab = tet_cluster14;
-		else if (cluster_reps.Get(el.PNum(2)) == 
-			 cluster_reps.Get(el.PNum(3)))
+		else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == 
+			 cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
 		  clustertab = tet_cluster23;
-		else if (cluster_reps.Get(el.PNum(2)) == 
-			 cluster_reps.Get(el.PNum(4)))
+		else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == 
+			 cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
 		  clustertab = tet_cluster24;
-		else if (cluster_reps.Get(el.PNum(3)) == 
-			 cluster_reps.Get(el.PNum(4)))
+		else if (cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE) == 
+			 cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
 		  clustertab = tet_cluster34;
 
 		else
diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp
index ee4e7c41..7aa5b43d 100644
--- a/libsrc/meshing/improve3.cpp
+++ b/libsrc/meshing/improve3.cpp
@@ -274,16 +274,8 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh,
 void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 				     OPTIMIZEGOAL goal)
 {
-  int j, k, l;
-  Point3d p1, p2, pnew;
-
-  ElementIndex ei;
-  SurfaceElementIndex sei;
-  PointIndex pi1, pi2;
-
   double bad1, bad2, badmax, badlimit;
 
-
   int cnt = 0;
   int np = mesh.GetNP();
   int ne = mesh.GetNE();
@@ -291,7 +283,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
   TABLE<ElementIndex,PointIndex::BASE> elementsonnode(np); 
   Array<ElementIndex> hasbothpoints;
 
-  BitArray origpoint(np), boundp(np);
+  BitArray origpoint(np+1), boundp(np+1);  // big enough for 0 and 1-based
   origpoint.Set();
 
   Array<double> elerrs(ne);
@@ -301,7 +293,6 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
   const char * savetask = multithread.task;
   multithread.task = "Split Improve";
 
-
   PrintMessage (3, "SplitImprove");
   (*testout)  << "start SplitImprove" << "\n";
 
@@ -311,10 +302,11 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 
   bad1 = 0;
   badmax = 0;
-  for (ei = 0; ei < ne; ei++)
+  for (ElementIndex ei = 0; ei < ne; ei++)
     {
-      if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
+      if(mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
       	continue;
+      
       elerrs[ei] = CalcBad (mesh.Points(), mesh[ei], 0);
       bad1 += elerrs[ei];
       if (elerrs[ei] > badmax) badmax = elerrs[ei];
@@ -325,41 +317,40 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 
 
   boundp.Clear();
-  for (sei = 0; sei < mesh.GetNSE(); sei++)
-    for (j = 0; j < 3; j++)
-      boundp.Set (mesh[sei][j]);
-
+  for (auto & el : mesh.SurfaceElements())
+    for (PointIndex pi : el.PNums())
+      boundp.Set (pi);
+  
   if (goal == OPT_QUALITY)
     {
       bad1 = CalcTotalBad (mesh.Points(), mesh.VolumeElements());
       (*testout) << "Total badness = " << bad1 << endl;
     }
 
-  for (ei = 0; ei < ne; ei++)
-    for (j = 0; j < mesh[ei].GetNP(); j++)
-      elementsonnode.Add (mesh[ei][j], ei);
-
+  for (ElementIndex ei : mesh.VolumeElements().Range())
+    for (PointIndex pi : mesh[ei].PNums())
+      elementsonnode.Add (pi, ei);
 
   mesh.MarkIllegalElements();
   if (goal == OPT_QUALITY || goal == OPT_LEGAL)
     {
       int cntill = 0;
-      for (ei = 0; ei < ne; ei++)
+      for (ElementIndex ei = 0; ei < ne; ei++)
 	{
-	    //	  if (!LegalTet (volelements.Get(i)))
-	    if (mesh[ei].flags.illegal)
-	      {
-		cntill++;
-		illegaltet.Set (ei+1);
-	      }
+          //	  if (!LegalTet (volelements.Get(i)))
+          if (mesh[ei].flags.illegal)
+            {
+              cntill++;
+              illegaltet.Set (ei);
+            }
 	}
-      //      (*mycout) << cntill << " illegal tets" << endl;
     }
 
-
-  for (ei = 0; ei < ne; ei++)
+  for (ElementIndex ei : mesh.VolumeElements().Range())    
     {
-      if(mesh.GetDimension()==3 && mp.only3D_domain_nr && mp.only3D_domain_nr != mesh.VolumeElement(ei).GetIndex())
+      Element & elem = mesh[ei];
+      
+      if(mp.only3D_domain_nr && mp.only3D_domain_nr != elem.GetIndex())
       	continue;
       if (multithread.terminate)
 	break;
@@ -368,82 +359,74 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 
       bool ltestmode = 0;
 
-
-      if (elerrs[ei] < badlimit && !illegaltet.Test(ei+1)) continue;
+      if (elerrs[ei] < badlimit && !illegaltet.Test(ei)) continue;
 
       if ((goal == OPT_LEGAL) &&
-	  !illegaltet.Test(ei+1) &&
-	  CalcBad (mesh.Points(), mesh[ei], 0) < 1e3) 
+	  !illegaltet.Test(ei) &&
+	  CalcBad (mesh.Points(), elem, 0) < 1e3) 
 	continue;
 
-      
-      Element & elem = mesh[ei];
-
       if (ltestmode)
 	{
 	  (*testout) << "test el " << ei << endl;
-	  for (j = 0; j < 4; j++)
+	  for (int j = 0; j < 4; j++)
 	    (*testout) << elem[j] << " ";
 	  (*testout) << endl;
 	}
 
-
-      for (j = 0; j < 6; j++)
+      for (int j = 0; j < 6; j++)
 	{
-
 	  static const int tetedges[6][2] =
 	    { { 0, 1 }, { 0, 2 }, { 0, 3 },
 	      { 1, 2 }, { 1, 3 }, { 2, 3 } };
 
-	  pi1 = elem[tetedges[j][0]];
-	  pi2 = elem[tetedges[j][1]];
+	  PointIndex pi1 = elem[tetedges[j][0]];
+	  PointIndex pi2 = elem[tetedges[j][1]];
 
 	  if (pi2 < pi1) Swap (pi1, pi2);
-	  if (pi2 > elementsonnode.Size()) continue;
+	  if (pi2 >= elementsonnode.Size()+PointIndex::BASE) continue; // old number of points
 
 	  if (!origpoint.Test(pi1) || !origpoint.Test(pi2))
 	    continue;
-
   
 	  INDEX_2 i2(pi1, pi2);
 	  i2.Sort();
 
 	  if (mesh.BoundaryEdge (pi1, pi2)) continue;
 
-	  if (edgetested.Used (i2) && !illegaltet.Test(ei+1)) continue;
+	  if (edgetested.Used (i2) && !illegaltet.Test(ei)) continue;
 	  edgetested.Set (i2, 1);
 
 	  hasbothpoints.SetSize (0);
-	  for (k = 1; k <= elementsonnode.EntrySize(pi1); k++)
+          /*
+	  for (int k = 1; k <= elementsonnode.EntrySize(pi1); k++)
 	    {
-	      bool has1 = 0, has2 = 0;
-
 	      ElementIndex elnr = elementsonnode.Get(pi1, k);
-	      Element & el = mesh[elnr];
+          */
+          for (ElementIndex ei : elementsonnode[pi1])
+            {
+	      Element & el = mesh[ei];
+              bool has1 = el.PNums().Contains(pi1);
+              bool has2 = el.PNums().Contains(pi2);
 
-	      for (l = 0; l < el.GetNP(); l++)
-		{
-		  if (el[l] == pi1) has1 = 1;
-		  if (el[l] == pi2) has2 = 1;
-		}
 	      if (has1 && has2) 
-		if (!hasbothpoints.Contains (elnr))
-		  hasbothpoints.Append (elnr);
+		if (!hasbothpoints.Contains (ei))
+		  hasbothpoints.Append (ei);
 	    }
 	  
 	  bad1 = 0;
-	  for (k = 0; k < hasbothpoints.Size(); k++)
-	    bad1 += CalcBad (mesh.Points(), mesh[hasbothpoints[k]], 0);
-	  
 
+          for (ElementIndex ei : hasbothpoints)
+	    bad1 += CalcBad (mesh.Points(), mesh[ei], 0);            
+          
 	  bool puretet = 1;
-	  for (k = 0; k < hasbothpoints.Size(); k++)
-	    if (mesh[hasbothpoints[k]].GetType() != TET)
+	  for (ElementIndex ei : hasbothpoints)
+	    if (mesh[ei].GetType() != TET)
 	      puretet = 0;
 	  if (!puretet) continue;
 
-	  p1 = mesh[pi1];
-	  p2 = mesh[pi2];
+	  Point3d p1 = mesh[pi1];
+	  Point3d p2 = mesh[pi2];
 
 	  /*
 	    pnew = Center (p1, p2);
@@ -463,13 +446,12 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 	    points.Elem(pi2) = p2;
 	  */
 
-
 	  locfaces.SetSize (0);
-	  for (k = 0; k < hasbothpoints.Size(); k++)
+	  for (ElementIndex ei : hasbothpoints)
 	    {
-	      const Element & el = mesh[hasbothpoints[k]];
-
-	      for (l = 0; l < 4; l++)
+	      const Element & el = mesh[ei];
+              
+	      for (int l = 0; l < 4; l++)
 		if (el[l] == pi1 || el[l] == pi2)
 		  {
 		    INDEX_3 i3;
@@ -486,7 +468,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 	  par.maxit_linsearch = 50;
 	  par.maxit_bfgs = 20;
 
-	  pnew = Center (p1, p2);
+	  Point3d pnew = Center (p1, p2);
 	  Vector px(3);
 	  px(0) = pnew.X();
 	  px(1) = pnew.Y();
@@ -501,11 +483,10 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 	  pnew.Y() = px(1);
 	  pnew.Z() = px(2);
 
-
-	  int hpinew = mesh.AddPoint (pnew);
+	  PointIndex hpinew = mesh.AddPoint (pnew);
 	  //	  ptyps.Append (INNERPOINT);
 
-	  for (k = 0; k < hasbothpoints.Size(); k++)
+	  for (int k = 0; k < hasbothpoints.Size(); k++)
 	    {
 	      Element & oldel = mesh[hasbothpoints[k]];
 	      Element newel1 = oldel;
@@ -515,7 +496,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 	      newel1.flags.illegal_valid = 0;
 	      newel2.flags.illegal_valid = 0;
 	      
-	      for (l = 0; l < 4; l++)
+	      for (int l = 0; l < 4; l++)
 		{
 		  if (newel1[l] == pi2) newel1[l] = hpinew;
 		  if (newel2[l] == pi1) newel2[l] = hpinew;
@@ -536,15 +517,15 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 
 	      PointIndex pinew = mesh.AddPoint (pnew);
 	      
-	      for (k = 0; k < hasbothpoints.Size(); k++)
+	      for (ElementIndex ei : hasbothpoints)
 		{
-		  Element & oldel = mesh[hasbothpoints[k]];
+		  Element & oldel = mesh[ei];
 		  Element newel = oldel;
-
+                  
 		  newel.flags.illegal_valid = 0;
 		  oldel.flags.illegal_valid = 0;
 
-		  for (l = 0; l < 4; l++)
+		  for (int l = 0; l < 4; l++)
 		    {
 		      origpoint.Clear (oldel[l]);
 
@@ -554,7 +535,7 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 		  mesh.AddVolumeElement (newel);
 		}
 	      
-	      j = 10;
+	      j = 10;  // end j-loop
 	    }
 	}
     }
@@ -572,11 +553,9 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh,
 
       int cntill = 0;
       ne = mesh.GetNE();
-      for (ei = 0; ei < ne; ei++)
-	{
-	  if (!mesh.LegalTet (mesh[ei]))
-	    cntill++;
-	}
+      for (ElementIndex ei = 0; ei < ne; ei++)
+        if (!mesh.LegalTet (mesh[ei]))
+          cntill++;
       //      cout << cntill << " illegal tets" << endl;
     }
 
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 7475e512..d76248a4 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -5367,30 +5367,9 @@ namespace netgen
   int Mesh :: MarkIllegalElements ()
   {
     int cnt = 0;
-    int i;
-
-    for (i = 1; i <= GetNE(); i++)
-      {
-        LegalTet (VolumeElement(i));
-
-        /*
-          Element & el = VolumeElement(i);
-          int leg1 = LegalTet (el);
-          el.flags.illegal_valid = 0;
-          int leg2 = LegalTet (el);
-
-          if (leg1 != leg2) 
-          {
-          cerr << "legal differs!!" << endl;
-          (*testout) << "legal differs" << endl;
-          (*testout) << "elnr = " << i << ", el = " << el
-          << " leg1 = " << leg1 << ", leg2 = " << leg2 << endl;
-          }
-
-          //      el.flags.illegal = !LegalTet (el);
-          */
-        cnt += VolumeElement(i).Illegal();
-      }
+    for (auto & el : VolumeElements())
+      if (!LegalTet (el))
+        cnt++;
     return cnt;
   }
 
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 9b747ac7..83565105 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -24,8 +24,8 @@ namespace netgen
   {
   public:
     typedef ::netgen::T_POINTS T_POINTS;
-    typedef Array<Element> T_VOLELEMENTS;
-    typedef Array<Element2d> T_SURFELEMENTS;
+    typedef Array<Element, 0, ElementIndex> T_VOLELEMENTS;
+    typedef Array<Element2d, 0, SurfaceElementIndex> T_SURFELEMENTS;
 
   private:
     /// point coordinates
@@ -192,24 +192,24 @@ namespace netgen
     void ClearSurfaceElements();
 
     ///
-	DLL_HEADER void ClearVolumeElements()
+    DLL_HEADER void ClearVolumeElements()
     {
       volelements.SetSize(0); 
       timestamp = NextTimeStamp();
     }
 
     ///
-	DLL_HEADER void ClearSegments()
+    DLL_HEADER void ClearSegments()
     { 
       segments.SetSize(0); 
       timestamp = NextTimeStamp();
     }
-  
+    
     ///
     bool TestOk () const;
 
     void SetAllocSize(int nnodes, int nsegs, int nsel, int nel);
-
+    
 
     DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer = 1);
     DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type);
@@ -236,13 +236,17 @@ namespace netgen
       segments.Elem(segnr)[0].Invalidate();
       segments.Elem(segnr)[1].Invalidate();
     }
+    /*
     void FullDeleteSegment (int segnr)  // von wem ist das ???
     {
       segments.Delete(segnr-PointIndex::BASE);
     }
+    */
 
     int GetNSeg () const { return segments.Size(); }
+    [[deprecated("Use LineSegment(SegmentIndex) instead of int !")]]                
     Segment & LineSegment(int i) { return segments.Elem(i); }
+    [[deprecated("Use LineSegment(SegmentIndex) instead of int !")]]                    
     const Segment & LineSegment(int i) const { return segments.Get(i); }
 
     Segment & LineSegment(SegmentIndex si) { return segments[si]; }
@@ -256,6 +260,8 @@ namespace netgen
     Array<Element0d> pointelements;  // only via python interface
 
     DLL_HEADER SurfaceElementIndex AddSurfaceElement (const Element2d & el);
+    
+    [[deprecated("Use DeleteSurfaceElement(SurfaceElementIndex) instead of int !")]]
     void DeleteSurfaceElement (int eli)
     { 
       surfelements.Elem(eli).Delete();
@@ -267,11 +273,16 @@ namespace netgen
 
     void DeleteSurfaceElement (SurfaceElementIndex eli)
     {
-      DeleteSurfaceElement (int(eli)+1);
+      for (auto & p : surfelements[eli].PNums()) p.Invalidate();
+      surfelements[eli].Delete();
+      timestamp = NextTimeStamp();
     }
 
     int GetNSE () const { return surfelements.Size(); }
+
+    [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]    
     Element2d & SurfaceElement(int i) { return surfelements.Elem(i); }
+    [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]]        
     const Element2d & SurfaceElement(int i) const { return surfelements.Get(i); }
     Element2d & SurfaceElement(SurfaceElementIndex i) { return surfelements[i]; }
     const Element2d & SurfaceElement(SurfaceElementIndex i) const { return surfelements[i]; }
@@ -309,9 +320,8 @@ namespace netgen
     ELEMENTTYPE ElementType (ElementIndex i) const 
     { return (volelements[i].flags.fixed) ? FIXEDELEMENT : FREEELEMENT; }
 
-    const T_VOLELEMENTS & VolumeElements() const { return volelements; }
-    T_VOLELEMENTS & VolumeElements() { return volelements; }
-
+    const auto & VolumeElements() const { return volelements; }
+    auto & VolumeElements() { return volelements; }
 
     ///
     DLL_HEADER double ElementError (int eli, const MeshingParameters & mp) const;
@@ -321,17 +331,13 @@ namespace netgen
     ///
     void ClearLockedPoints ();
 
-    const Array<PointIndex> & LockedPoints() const
-    { return lockedpoints; }
+    const auto & LockedPoints() const { return lockedpoints; }
 
     /// Returns number of domains
     int GetNDomains() const;
-
     ///
-    int GetDimension() const 
-    { return dimension; }
-    void SetDimension(int dim)
-    { dimension = dim; }
+    int GetDimension() const { return dimension; }
+    void SetDimension (int dim) { dimension = dim; }
 
     /// sets internal tables
     void CalcSurfacesOfNode ();
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 062c3e99..9af43005 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -708,7 +708,7 @@ namespace netgen
 	  break;
 
 	PrintMessage (5, nillegal, " illegal tets");
-	optmesh.SplitImprove (mesh3d, OPT_LEGAL);
+        optmesh.SplitImprove (mesh3d, OPT_LEGAL);
 
 	mesh3d.MarkIllegalElements();  // test
 	optmesh.SwapImprove (mesh3d, OPT_LEGAL);
diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp
index c71f3251..5f12c32d 100644
--- a/libsrc/meshing/meshing3.cpp
+++ b/libsrc/meshing/meshing3.cpp
@@ -658,8 +658,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind);
 	    }
 
-          cout << "bbb121212" << endl;
-          
 	  for (int i = 1; i <= locelements.Size(); i++)
 	    {
 	      Point3d * hp1, * hp2, * hp3, * hp4;
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index ec467c99..a38bcad0 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -156,8 +156,10 @@ namespace netgen
     ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; }
     ElementIndex & operator= (int ai) { i = ai; return *this; }
     operator int () const { return i; }
-    ElementIndex & operator++ (int) { i++; return *this; }
-    ElementIndex & operator-- (int) { i--; return *this; }
+    ElementIndex operator++ (int) { return ElementIndex(i++); }    
+    ElementIndex operator-- (int) { return ElementIndex(i--); }
+    ElementIndex & operator++ () { ++i; return *this; }
+    ElementIndex & operator-- () { --i; return *this; }
   };
 
   inline istream & operator>> (istream & ist, ElementIndex & pi)
@@ -181,9 +183,11 @@ namespace netgen
     { i = ai.i; return *this; }
     SurfaceElementIndex & operator= (int ai) { i = ai; return *this; }
     operator int () const { return i; }
-    SurfaceElementIndex & operator++ (int) { i++; return *this; }
+    SurfaceElementIndex operator++ (int) { SurfaceElementIndex hi(*this); i++; return hi; }
+    SurfaceElementIndex operator-- (int) { SurfaceElementIndex hi(*this); i--; return hi; }
+    SurfaceElementIndex & operator++ () { ++i; return *this; }
+    SurfaceElementIndex & operator-- () { --i; return *this; }
     SurfaceElementIndex & operator+= (int inc) { i+=inc; return *this; }
-    SurfaceElementIndex & operator-- (int) { i--; return *this; }
   };
 
   inline istream & operator>> (istream & ist, SurfaceElementIndex & pi)
@@ -389,6 +393,8 @@ namespace netgen
 
     FlatArray<const PointIndex> PNums () const 
     { return FlatArray<const PointIndex> (np, &pnum[0]); }
+    FlatArray<PointIndex> PNums ()
+    { return FlatArray<PointIndex> (np, &pnum[0]); }
     
     ///
     PointIndex & PNum (int i) { return pnum[i-1]; }
@@ -477,8 +483,7 @@ namespace netgen
     void Delete ()
     {
       deleted = 1;
-      for (PointIndex & p : pnum)
-        p.Invalidate(); 
+      for (PointIndex & p : pnum) p.Invalidate(); 
     }
     
     bool IsDeleted () const 
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 8ac969c1..10a19466 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -304,8 +304,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
 
   
 
-  ExportArray<Element>(m);
-  ExportArray<Element2d>(m);
+  ExportArray<Element,0,ElementIndex>(m);
+  ExportArray<Element2d,0,SurfaceElementIndex>(m);
   ExportArray<Segment>(m);
   ExportArray<Element0d>(m);
   ExportArray<MeshPoint,PointIndex::BASE,PointIndex>(m);
@@ -384,11 +384,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension)
 
     .def("Elements3D", 
-         static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
+         static_cast<Array<Element,0,ElementIndex>&(Mesh::*)()> (&Mesh::VolumeElements),
          py::return_value_policy::reference)
 
     .def("Elements2D", 
-         static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements),
+         static_cast<Array<Element2d,0,SurfaceElementIndex>&(Mesh::*)()> (&Mesh::SurfaceElements),
          py::return_value_policy::reference)
 
     .def("Elements1D", 
diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp
index 686b7190..4e8a5ae7 100644
--- a/libsrc/meshing/refine.cpp
+++ b/libsrc/meshing/refine.cpp
@@ -383,11 +383,11 @@ namespace netgen
 	     swap (pnums.Elem(3), pnums.Elem(4));
 
 	     for (int j = 0; j < 6; j++)
-	     {
-	       INDEX_2 i2;
-	       i2.I1() = pnums.Get(betw[j][0]);
-	       i2.I2() = pnums.Get(betw[j][1]);
-	       i2.Sort();
+               {
+                 PointIndex pi1 = pnums.Get(betw[j][0]);
+                 PointIndex pi2 = pnums.Get(betw[j][1]);
+                 INDEX_2 i2 (pi1, pi2);
+                 i2.Sort();
 
 	       /*
 	       if (between.Used(i2))
@@ -405,8 +405,8 @@ namespace netgen
 	       if (!pointset[pinew])
 		 {
 		   pointset[pinew] = true;
-		   mesh.Point(pinew) = Center(mesh.Point(i2.I1()),
-					      mesh.Point(i2.I2()));
+		   mesh.Point(pinew) = Center(mesh.Point(pi1),
+					      mesh.Point(pi2));
 		 }
 	    }
 
diff --git a/libsrc/visualization/vispar.hpp b/libsrc/visualization/vispar.hpp
index 1661ac3b..94026440 100644
--- a/libsrc/visualization/vispar.hpp
+++ b/libsrc/visualization/vispar.hpp
@@ -76,7 +76,7 @@ public:
   int drawcurveprojedge;
   
 
-  int centerpoint;
+  PointIndex centerpoint;
   int drawelement;
 
   // stl:
diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp
index cd366aa5..08c2f6de 100644
--- a/libsrc/visualization/vsmesh.cpp
+++ b/libsrc/visualization/vsmesh.cpp
@@ -440,12 +440,12 @@ namespace netgen
 	char buf[30];
 
 	if (vispar.drawpointnumbers)
-	  for (int i = 1; i <= mesh->GetNP(); i++)
+	  for (PointIndex pi : mesh->Points().Range())
             {
-	      const Point3d & p = mesh->Point(i);
+	      const Point3d & p = mesh->Point(pi);
 	      glRasterPos3d (p.X(), p.Y(), p.Z());
 
-	      sprintf (buf, "%d", i);
+	      sprintf (buf, "%d", int(pi));
 
 	      // glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
 	      MyOpenGLText (buf);
@@ -663,12 +663,12 @@ namespace netgen
 
 
 
-	for (int i = 1; i <= mesh->GetNE(); i++)
+	for (ElementIndex ei : mesh->VolumeElements().Range())
 	  {
-            if (mesh->VolumeElement(i).flags.badel)
+            if (mesh->VolumeElement(ei).flags.badel)
 	      {
 		// copy to be thread-safe
-		Element el = mesh->VolumeElement (i);
+		Element el = mesh->VolumeElement (ei);
 		if ( (el.GetNP() == 4) || (el.GetNP() == 10))
 		  {
 		    glBegin (GL_LINES);
@@ -747,17 +747,16 @@ namespace netgen
 	      }
 	  }
 
-
-	for (int i = 1; i <= mesh->GetNSE(); i++)
+        for (SurfaceElementIndex sei : mesh->SurfaceElements().Range())
 	  {
-            Element2d el = mesh->SurfaceElement(i);
+            Element2d el = mesh->SurfaceElement(sei); // copy to be thread-safe
             if (!el.BadElement())
 	      continue;
 
-            int drawel = 1;
+            bool drawel = true;
             for (int j = 1; j <= el.GetNP(); j++)
-	      if (!el.PNum(j))
-		drawel = 0;
+	      if (!el.PNum(j).IsValid())
+		drawel = false;
 
             if (!drawel)
 	      continue;