diff --git a/libsrc/csg/zrefine.cpp b/libsrc/csg/zrefine.cpp
index 7af6d71f..649bb50c 100644
--- a/libsrc/csg/zrefine.cpp
+++ b/libsrc/csg/zrefine.cpp
@@ -45,9 +45,10 @@ namespace netgen
   void MakePrismsSingEdge (Mesh & mesh, INDEX_2_HASHTABLE<int> & singedges)
   {
     // volume elements
-    for (int i = 1; i <= mesh.GetNE(); i++)
+    // for (int i = 1; i <= mesh.GetNE(); i++)
+    for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
       {
-	Element & el = mesh.VolumeElement(i);
+	Element & el = mesh.VolumeElement(ei);
 	if (el.GetType() != TET) continue;
 
 	for (int j = 1; j <= 3; j++)
@@ -76,9 +77,9 @@ namespace netgen
       }
 
     // surface elements
-    for (int i = 1; i <= mesh.GetNSE(); i++)
+    for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
       {
-	Element2d & el = mesh.SurfaceElement(i);
+	Element2d & el = mesh.SurfaceElement(sei);
 	if (el.GetType() != TRIG) continue;
 
 	for (int j = 1; j <= 3; j++)
@@ -110,14 +111,14 @@ namespace netgen
   */
   void MakePrismsClosePoints (Mesh & mesh)
   {
-    int i, j, k;
-    for (i = 1; i <= mesh.GetNE(); i++)
+    // int i, j, k;
+    for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
       {
-	Element & el = mesh.VolumeElement(i);
+	Element & el = mesh.VolumeElement(ei);
 	if (el.GetType() == TET)
 	  {
-	    for (j = 1; j <= 3; j++)
-	      for (k = j+1; k <= 4; k++)
+	    for (int j = 1; j <= 3; j++)
+	      for (int k = j+1; k <= 4; k++)
 		{
 		  INDEX_2 edge(el.PNum(j), el.PNum(k));
 		  edge.Sort();
@@ -145,7 +146,7 @@ namespace netgen
 	  {
 	    // pyramid, base face = 1,2,3,4
 	  
-	    for (j = 0; j <= 1; j++)
+	    for (int j = 0; j <= 1; j++)
 	      {
 		PointIndex pi1 = el.PNum( (j+0) % 4 + 1);
 		PointIndex pi2 = el.PNum( (j+1) % 4 + 1);
@@ -175,14 +176,14 @@ namespace netgen
 	  }
       }
   
-    for (i = 1; i <= mesh.GetNSE(); i++)
+    for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
       {
-	Element2d & el = mesh.SurfaceElement(i);
+	Element2d & el = mesh.SurfaceElement(sei);
 	if (el.GetType() != TRIG) continue;
 
-	for (j = 1; j <= 3; j++)
+	for (int j = 1; j <= 3; j++)
 	  {
-	    k = (j % 3) + 1;
+	    int k = (j % 3) + 1;
 	    INDEX_2 edge(el.PNum(j), el.PNum(k));
 	    edge.Sort();
 
@@ -244,7 +245,7 @@ namespace netgen
   void RefinePrisms (Mesh & mesh, const CSGeometry * geom, 
 		     ZRefinementOptions & opt)
   {
-    int i, j;
+    // int i, j;
     bool found, change;
     int cnt = 0;
 
@@ -264,8 +265,8 @@ namespace netgen
         auto & identpts =
           mesh.GetIdentifications().GetIdentifiedPoints ();
         
-	for (i = 1; i <= identpts.GetNBags(); i++)
-	  for (j = 1; j <= identpts.GetBagSize(i); j++)
+	for (int i = 1; i <= identpts.GetNBags(); i++)
+	  for (int j = 1; j <= identpts.GetBagSize(i); j++)
 	    {
 	      INDEX_3 pair;
 	      int dummy;
@@ -314,7 +315,7 @@ namespace netgen
 	found = 0;
 	// mark prisms due to close surface flags:
 	int oldsize = ref_uniform.Size();
-	for (i = 1; i <= oldsize; i++)
+	for (int i = 1; i <= oldsize; i++)
 	  {
 	    int pi1 = ref_uniform.Get(i).I1();
 	    int pi2 = ref_uniform.Get(i).I2();
@@ -340,7 +341,7 @@ namespace netgen
 		ref_uniform.Append (INDEX_3(pi2, npi, levels-1));
 	      }
 	  }
-	for (i = 1; i <= ref_singular.Size(); i++)
+	for (int i = 1; i <= ref_singular.Size(); i++)
 	  {
 	    int pi1 = ref_singular.Get(i).I1();
 	    int pi2 = ref_singular.Get(i).I2();
@@ -368,7 +369,7 @@ namespace netgen
 	      }
 	  }
 
-	for (i = 1; i <= ref_slices.Size(); i++)
+	for (int i = 1; i <= ref_slices.Size(); i++)
 	  {
 	    int pi1 = ref_slices.Get(i)[0];
 	    int pi2 = ref_slices.Get(i)[1];
@@ -414,13 +415,13 @@ namespace netgen
 
 
 
-	for (i = 1; i <= mesh.GetNE(); i++)
+	for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
 	  {
-	    Element & el = mesh.VolumeElement (i);
+	    Element & el = mesh.VolumeElement (ei);
 	    if (el.GetType() != PRISM)
 	      continue;
 
-	    for (j = 1; j <= 3; j++)
+	    for (int j = 1; j <= 3; j++)
 	      {
 		int pi1 = el.PNum(j);
 		int pi2 = el.PNum(j+3);
@@ -466,14 +467,14 @@ namespace netgen
 	  {
 	    PrintMessage (5, "start loop");
 	    change = 0;
-	    for (i = 1; i <= mesh.GetNE(); i++)
+	    for (ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
 	      {
-		Element & el = mesh.VolumeElement (i);
+		Element & el = mesh.VolumeElement (ei);
 		if (el.GetType() != PRISM)
 		  continue;
 	      
 		bool hasref = 0, hasnonref = 0;
-		for (j = 1; j <= 3; j++)
+		for (int j = 1; j <= 3; j++)
 		  {
 		    int pi1 = el.PNum(j);
 		    int pi2 = el.PNum(j+3);
@@ -492,7 +493,7 @@ namespace netgen
 		  {
 		    //		  cout << "el " << i << " in closure" << endl;
 		    change = 1;
-		    for (j = 1; j <= 3; j++)
+		    for (int j = 1; j <= 3; j++)
 		      {
 			int pi1 = el.PNum(j);
 			int pi2 = el.PNum(j+3);
@@ -519,7 +520,7 @@ namespace netgen
 
 	int oldns = mesh.GetNSeg();
 
-	for (i = 1; i <= oldns; i++)
+	for (int i = 1; i <= oldns; i++)
 	  {
 	    const Segment & el = mesh.LineSegment(i);
 
@@ -573,14 +574,14 @@ namespace netgen
 
 	// do refinement
 	int oldne = mesh.GetNE();
-	for (i = 1; i <= oldne; i++)
+	for (ElementIndex ei = 0; ei < oldne; ei++)
 	  {
-	    Element & el = mesh.VolumeElement (i);
+	    Element & el = mesh.VolumeElement (ei);
 	    if (el.GetNP() != 6)
 	      continue;
 
 	    int npi[3];
-	    for (j = 1; j <= 3; j++)
+	    for (int j = 1; j <= 3; j++)
 	      {
 		int pi1 = el.PNum(j);
 		int pi2 = el.PNum(j+3);
@@ -608,7 +609,7 @@ namespace netgen
 	    if (npi[0])
 	      {
 		Element nel1(6), nel2(6);
-		for (j = 1; j <= 3; j++)
+		for (int j = 1; j <= 3; j++)
 		  {
 		    nel1.PNum(j) = el.PNum(j);
 		    nel1.PNum(j+3) = npi[j-1];
@@ -617,7 +618,7 @@ namespace netgen
 		  }
 		nel1.SetIndex (el.GetIndex());
 		nel2.SetIndex (el.GetIndex());
-		mesh.VolumeElement (i) = nel1;
+		mesh.VolumeElement (ei) = nel1;
 		mesh.AddVolumeElement (nel2);
 	      }
 	  }
@@ -629,15 +630,15 @@ namespace netgen
 	// do surface elements
 	int oldnse = mesh.GetNSE();
 	//      cout << "oldnse = " << oldnse << endl;
-	for (i = 1; i <= oldnse; i++)
+	for (SurfaceElementIndex sei = 0; sei < oldnse; sei++)
 	  {
-	    Element2d & el = mesh.SurfaceElement (i);
+	    Element2d & el = mesh.SurfaceElement (sei);
 	    if (el.GetType() != QUAD)
 	      continue;
 
 	    int index = el.GetIndex();
 	    int npi[2];
-	    for (j = 1; j <= 2; j++)
+	    for (int j = 1; j <= 2; j++)
 	      {
 		int pi1, pi2;
 
@@ -670,7 +671,7 @@ namespace netgen
 	    if (npi[0])
 	      {
 		Element2d nel1(QUAD), nel2(QUAD);
-		for (j = 1; j <= 4; j++)
+		for (int j = 1; j <= 4; j++)
 		  {
 		    nel1.PNum(j) = el.PNum(j);
 		    nel2.PNum(j) = el.PNum(j);
@@ -691,7 +692,7 @@ namespace netgen
 		nel1.SetIndex (el.GetIndex());
 		nel2.SetIndex (el.GetIndex());
 
-		mesh.SurfaceElement (i) = nel1;
+		mesh.SurfaceElement (sei) = nel1;
 		mesh.AddSurfaceElement (nel2);
 
 		int si = mesh.GetFaceDescriptor (index).SurfNr();
@@ -717,9 +718,9 @@ namespace netgen
 
   void CombineSingularPrisms(Mesh& mesh)
   {
-    for(int i = 1; i<=mesh.GetNE(); i++)
+    for(ElementIndex ei = 0; ei < mesh.GetNE(); ei++)
       {
-        Element& el = mesh.VolumeElement(i);
+        Element& el = mesh.VolumeElement(ei);
         if(el.GetType() != PRISM)
           continue;
         if(el.PNum(3) == el.PNum(6))
diff --git a/libsrc/general/template.hpp b/libsrc/general/template.hpp
index 2083c93b..89103ea4 100644
--- a/libsrc/general/template.hpp
+++ b/libsrc/general/template.hpp
@@ -138,7 +138,7 @@ public:
   ///
 
 
-  INDEX_2 Sort ()
+  constexpr INDEX_2 Sort ()
   {
     if (i[0] > i[1]) 
       {
diff --git a/libsrc/interface/writediffpack.cpp b/libsrc/interface/writediffpack.cpp
index 9797dbed..54fe840d 100644
--- a/libsrc/interface/writediffpack.cpp
+++ b/libsrc/interface/writediffpack.cpp
@@ -40,7 +40,7 @@ void WriteDiffPackFormat (const Mesh & mesh,
       int nse = mesh.GetNSE();
       NgArray <int> BIname;
       NgArray <int> BCsinpoint;
-      int i, j, k, l;
+      // int i, j, k, l;
 
 
       outfile.precision(6);
@@ -58,18 +58,18 @@ void WriteDiffPackFormat (const Mesh & mesh,
 	"  Only one subdomain               : dpFALSE\n"
 	"  Lattice data                     ? 0\n\n\n\n";
       
-      for (i = 1; i <= nse; i++) 
+      for (int i = 1; i <= nse; i++) 
 	{
 	  int BI=mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex()).BCProperty();
 	  int nbi=BIname.Size();
 	  int found=0;
-	  for (j = 1; j <= nbi; j++)
+	  for (int j = 1; j <= nbi; j++)
 	    if(BI == BIname.Get(j)) found = 1;
 	  if( ! found ) BIname.Append(BI);	    	     
 	}
       
       outfile << "  " << BIname.Size() <<  " Boundary indicators:  ";
-      for (i =1 ; i <= BIname.Size(); i++)
+      for (int i =1 ; i <= BIname.Size(); i++)
 	outfile << BIname.Get(i) << " ";
       outfile << "\n\n\n";
       
@@ -92,7 +92,7 @@ void WriteDiffPackFormat (const Mesh & mesh,
 	}
 
 
-      for (i = 1; i <= np; i++)
+      for (int i = 1; i <= np; i++)
         {
           const Point3d & p = mesh.Point(i);
 
@@ -114,14 +114,14 @@ void WriteDiffPackFormat (const Mesh & mesh,
 	      NgFlatArray<SurfaceElementIndex> sels = point2sel[i];
 	      for (int jj = 0; jj < sels.Size(); jj++)
 		{
-		  for (k = 1; k <= mesh[sels[jj]].GetNP(); k++) 
+		  for (int k = 1; k <= mesh[sels[jj]].GetNP(); k++) 
 		    {
 		      if(mesh[sels[jj]].PNum(k)==i) 
 			{
 			  int BC=mesh.GetFaceDescriptor(mesh[sels[jj]].GetIndex()).BCProperty();
 			  int nbcsp=BCsinpoint.Size();
 			  int found = 0;
-			  for (l = 1; l <= nbcsp; l++)
+			  for (int l = 1; l <= nbcsp; l++)
 			    if(BC == BCsinpoint.Get(l)) found = 1;
 			  if( ! found ) BCsinpoint.Append(BC); 	    	     
 			}
@@ -129,7 +129,7 @@ void WriteDiffPackFormat (const Mesh & mesh,
 		}
 	      int nbcsp = BCsinpoint.Size();
 	      outfile << "[" << nbcsp << "] ";
-	      for (j = 1; j <= nbcsp; j++)
+	      for (int j = 1; j <= nbcsp; j++)
 		outfile << BCsinpoint.Get(j) << " ";
 	      outfile << "\n";
             }
@@ -146,7 +146,7 @@ void WriteDiffPackFormat (const Mesh & mesh,
 	"   - the global node numbers of the nodes in the element.\n"
 	"#\n";
 
-      for (i = 1; i <= ne; i++)
+      for (int i = 1; i <= ne; i++)
         {
           const Element & el = mesh.VolumeElement(i);
           outfile.width(5);
diff --git a/libsrc/interface/writeuser.cpp b/libsrc/interface/writeuser.cpp
index bfab7fbb..27370eb0 100644
--- a/libsrc/interface/writeuser.cpp
+++ b/libsrc/interface/writeuser.cpp
@@ -62,7 +62,7 @@ void WriteNeutralFormat (const Mesh & mesh,
   int ne = mesh.GetNE();
   int nse = mesh.GetNSE();
   int nseg = mesh.GetNSeg();
-  int i, j;
+  // int i, j;
 
   int inverttets = mparam.inverttets;
   int invertsurf = mparam.inverttrigs;
@@ -75,7 +75,7 @@ void WriteNeutralFormat (const Mesh & mesh,
 
   outfile << np << "\n";
 
-  for (i = 1; i <= np; i++)
+  for (int i = 1; i <= np; i++)
     {
       const Point3d & p = mesh.Point(i);
 
@@ -94,14 +94,14 @@ void WriteNeutralFormat (const Mesh & mesh,
   if (mesh.GetDimension() == 3)
     {
       outfile << ne << "\n";
-      for (i = 1; i <= ne; i++)
+      for (int i = 1; i <= ne; i++)
 	{
 	  Element el = mesh.VolumeElement(i);
 	  if (inverttets)
 	    el.Invert();
 	  outfile.width(4);
 	  outfile << el.GetIndex() << "  ";
-	  for (j = 1; j <= el.GetNP(); j++)
+	  for (int j = 1; j <= el.GetNP(); j++)
 	    {
 	      outfile << " ";
 	      outfile.width(8);
@@ -112,14 +112,14 @@ void WriteNeutralFormat (const Mesh & mesh,
     }
 
   outfile << nse << "\n";
-  for (i = 1; i <= nse; i++)
+  for (int i = 1; i <= nse; i++)
     {
       Element2d el = mesh.SurfaceElement(i);
       if (invertsurf)
 	el.Invert();
       outfile.width(4);
       outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << "    ";
-      for (j = 1; j <= el.GetNP(); j++)
+      for (int j = 1; j <= el.GetNP(); j++)
 	{
 	  outfile << " ";
 	  outfile.width(8);
diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index e55dd657..d7a4062d 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -1404,7 +1404,7 @@ namespace netgen
 
 
   void BTBisectIdentification (const MarkedIdentification & oldid,
-			       NgArray<PointIndex> & newp,
+			       Array<PointIndex> & newp,
 			       MarkedIdentification & newid1,
 			       MarkedIdentification & newid2)
   {
@@ -2157,9 +2157,10 @@ namespace netgen
 	      }
 	  }
 	
-	for (int i = 1; i <= nse; i++)
+	// for (int i = 1; i <= nse; i++)
+        for (SurfaceElementIndex sei = 0; sei < nse; sei++)
 	  {
-	    const Element2d & el = mesh.SurfaceElement(i);
+	    const Element2d & el = mesh[sei];
 	    if (el.GetType() == TRIG ||
 		el.GetType() == TRIG6)
 	      {
@@ -2907,7 +2908,7 @@ namespace netgen
     // INDEX_2_HASHTABLE<int> cutedges(10 + 5 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
     // INDEX_2_CLOSED_HASHTABLE<PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
     // ClosedHashTable<INDEX_2, PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
-    ClosedHashTable<PointIndices<2>, PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
+    ClosedHashTable<SortedPointIndices<2>, PointIndex> cutedges(10 + 9 * (mtets.Size()+mprisms.Size()+mtris.Size()+mquads.Size()));
 
     bool noprojection = false;
     NgProfiler::StopTimer (timer1a);
@@ -3125,14 +3126,15 @@ namespace netgen
 
 	    int cnttrig = 0;
 	    int cntquad = 0;
-	    for (int i = 1; i <= mesh.GetNSE(); i++)
+	    // for (int i = 1; i <= mesh.GetNSE(); i++)
+            for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++)
 	      {
-		if (mesh.SurfaceElement(i).GetType() == TRIG ||
-		    mesh.SurfaceElement(i).GetType() == TRIG6)
+		if (mesh[sei].GetType() == TRIG ||
+		    mesh[sei].GetType() == TRIG6)
 		  {
 		    cnttrig++;
 		    mtris.Elem(cnttrig).marked =
-		      mesh.SurfaceElement(i).TestRefinementFlag() ? (opt.onlyonce ? 1 : 2) : 0;
+		      mesh[sei].TestRefinementFlag() ? (opt.onlyonce ? 1 : 2) : 0;
 		    // mtris.Elem(cnttrig).marked = 0;
 		    if (mtris.Elem(cnttrig).marked)
 		      cntm++;
@@ -3142,7 +3144,7 @@ namespace netgen
 		    cntquad++;
                     // 2d: marked=2, 3d prisms: marked=1
 		    mquads.Elem(cntquad).marked =
-                        mesh.SurfaceElement(i).TestRefinementFlag() ? 4-mesh.GetDimension() : 0 ;
+                        mesh[sei].TestRefinementFlag() ? 4-mesh.GetDimension() : 0 ;
 		    // mquads.Elem(cntquad).marked = 0;
 		    if (mquads.Elem(cntquad).marked)
 		      cntm++;
@@ -3338,9 +3340,8 @@ namespace netgen
 		{
 		  MarkedTet oldtet = mtets[i];
                   
-		  PointIndices<2> edge(oldtet.pnums[oldtet.tetedge1],
-                                       oldtet.pnums[oldtet.tetedge2]);
-		  edge.Sort();
+		  SortedPointIndices<2> edge(oldtet.pnums[oldtet.tetedge1],
+                                             oldtet.pnums[oldtet.tetedge2]);
 
 		  PointIndex newp;
                   if (auto optnewp = cutedges.GetIfUsed(edge))
@@ -3377,12 +3378,10 @@ namespace netgen
 		    pi1++;
 		  int pi2 = 3-pi1-oldprism.markededge;
                   
-		  PointIndices<2> edge1(oldprism.pnums[pi1],
-                                        oldprism.pnums[pi2]);
-		  PointIndices<2> edge2(oldprism.pnums[pi1+3],
-                                        oldprism.pnums[pi2+3]);
-		  edge1.Sort();
-		  edge2.Sort();
+		  SortedPointIndices<2> edge1(oldprism.pnums[pi1],
+                                              oldprism.pnums[pi2]);
+		  SortedPointIndices<2> edge2(oldprism.pnums[pi1+3],
+                                              oldprism.pnums[pi2+3]);
 
 		  if (cutedges.Used (edge1))
 		    newp1 = cutedges.Get(edge1);
@@ -3414,7 +3413,7 @@ namespace netgen
 	      if (mids.Elem(i).marked)
 		{
 		  MarkedIdentification oldid,newid1,newid2;
-		  NgArray<PointIndex> newp;
+		  Array<PointIndex> newp;
 
 		  oldid = mids.Get(i);
 		  
diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp
index d1baf385..3c07a6b3 100644
--- a/libsrc/meshing/clusters.cpp
+++ b/libsrc/meshing/clusters.cpp
@@ -231,7 +231,9 @@ namespace netgen
       
 	for (int i = 1; i <= ne; i++)
 	  {
-	    const Element & el = mesh.VolumeElement(i);
+            ElementIndex ei(i-1);
+            
+	    const Element & el = mesh[ei];
 	    ELEMENT_TYPE typ = el.GetType();
 	  	  
 	    const int * clustertab = NULL;
@@ -279,8 +281,8 @@ namespace netgen
               {
                 // top.GetElementEdges (i, ednums);
                 // top.GetElementFaces (i, fanums);
-                auto ednums = top.GetEdges (ElementIndex(i-1));
-                auto fanums = top.GetFaces (ElementIndex(i-1));                
+                auto ednums = top.GetEdges (ei);
+                auto fanums = top.GetFaces (ei);
                 
                 int elnv = top.GetNVertices (typ);
                 int elned = ednums.Size();
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 9a102f16..dd783130 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -3303,47 +3303,53 @@ namespace netgen
   void Mesh :: FreeOpenElementsEnvironment (int layers)
   {
     static Timer timer("FreeOpenElementsEnvironment"); RegionTimer rt(timer);
-    int i, j, k;
-    PointIndex pi;
     const int large = 9999;
-    NgArray<int,PointIndex::BASE> dist(GetNP());
+    Array<int,PointIndex> dist(GetNP());
 
     dist = large;
 
     for (int i = 1; i <= GetNOpenElements(); i++)
       {
         const Element2d & face = OpenElement(i);
-        for (j = 0; j < face.GetNP(); j++)
+        for (int j = 0; j < face.GetNP(); j++)
           dist[face[j]] = 1;
       }
 
-    for (k = 1; k <= layers; k++)
+    for (int k = 1; k <= layers; k++)
+      /*
       for (i = 1; i <= GetNE(); i++)
         {
           const Element & el = VolumeElement(i);
+      */
+      for (auto & el : VolumeElements())
+        {
           if (el[0] == -1 || el.IsDeleted()) continue;
 
           int elmin = large;
-          for (j = 0; j < el.GetNP(); j++)
+          for (int j = 0; j < el.GetNP(); j++)
             if (dist[el[j]] < elmin)
               elmin = dist[el[j]];
-
+          
           if (elmin < large)
             {
-              for (j = 0; j < el.GetNP(); j++)
+              for (int j = 0; j < el.GetNP(); j++)
                 if (dist[el[j]] > elmin+1)
                   dist[el[j]] = elmin+1;
             }
         }
 
     int cntfree = 0;
-    for (i = 1; i <= GetNE(); i++)
+    /*
+    for (int i = 1; i <= GetNE(); i++)
       {
         Element & el = VolumeElement(i);
+    */
+    for (auto & el : VolumeElements())
+      {
         if (el[0] == -1 || el.IsDeleted()) continue;
 
         int elmin = large;
-        for (j = 0; j < el.GetNP(); j++)
+        for (int j = 0; j < el.GetNP(); j++)
           if (dist[el[j]] < elmin)
             elmin = dist[el[j]];
 
@@ -3357,7 +3363,7 @@ namespace netgen
     PrintMessage (5, "free: ", cntfree, ", fixed: ", GetNE()-cntfree);
     (*testout) << "free: " << cntfree << ", fixed: " << GetNE()-cntfree << endl;
 
-    for (pi = PointIndex::BASE; 
+    for (PointIndex pi = PointIndex::BASE; 
          pi < GetNP()+PointIndex::BASE; pi++)
       {
         if (dist[pi] > layers+1)
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 4700d447..46889487 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -344,6 +344,22 @@ namespace netgen
     PointIndex get() const { return PointIndex(INDEX_4::operator[](J)); }    
   };
 
+
+  template <int N>
+  class SortedPointIndices : public PointIndices<N>
+  {
+    using PointIndices<N>::Sort;
+  public:
+    constexpr SortedPointIndices (PointIndices<N> pnts)
+      : PointIndices<N>(pnts.Sort()) { } 
+    
+      template <typename ...Pnts>
+    constexpr SortedPointIndices (Pnts ...pnts)
+      : PointIndices<N>(pnts...)
+    { Sort(); }
+  };
+  
+  
 }
 
 namespace ngcore
@@ -357,6 +373,11 @@ namespace ngcore
   template <>
   constexpr inline netgen::PointIndices<3> InvalidHash<netgen::PointIndices<3>> ()
   { return netgen::PointIndices<3>{netgen::PointIndex::INVALID, netgen::PointIndex::INVALID, netgen::PointIndex::INVALID}; }
+
+  template <>
+  constexpr inline netgen::SortedPointIndices<2> InvalidHash<netgen::SortedPointIndices<2>> ()
+  { return InvalidHash<netgen::PointIndices<2>>(); }
+
 }
 
 
@@ -376,7 +397,9 @@ namespace netgen
     int i;
   public:
     ElementIndex () = default;
+    // private:
     constexpr ElementIndex (int ai) : i(ai) { ; }
+  public:
     ElementIndex & operator= (const ElementIndex & ai) { i = ai.i; return *this; }
     ElementIndex & operator= (int ai) { i = ai; return *this; }
     constexpr operator int () const { return i; }
diff --git a/libsrc/meshing/refine.cpp b/libsrc/meshing/refine.cpp
index 6c45efe7..b6b27369 100644
--- a/libsrc/meshing/refine.cpp
+++ b/libsrc/meshing/refine.cpp
@@ -790,14 +790,15 @@ namespace netgen
 
     int cnttrials = 10;
     int wrongels = 0;
-    for (int i = 1; i <= mesh.GetNE(); i++)
-      if (mesh.VolumeElement(i).Volume(mesh.Points()) < 0)
+
+    for (auto & el : mesh.VolumeElements())
+      if (el.Volume(mesh.Points()) < 0)
 	{
 	  wrongels++;
-	  mesh.VolumeElement(i).Flags().badel = 1;
+	  el.Flags().badel = 1;
 	}
       else
-	mesh.VolumeElement(i).Flags().badel = 0;
+	el.Flags().badel = 0;
 
     if (wrongels)
       {
diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp
index 7734f237..47d3211c 100644
--- a/libsrc/meshing/secondorder.cpp
+++ b/libsrc/meshing/secondorder.cpp
@@ -507,12 +507,15 @@ namespace netgen
 
 	TBitArray<PointIndex> boundp(np);
 	boundp.Clear();
+        /*
 	for (int i = 1; i <= mesh.GetNSE(); i++)
 	  {
 	    const Element2d & sel = mesh.SurfaceElement(i);
-	    for (int j = 1; j <= sel.GetNP(); j++)
-	      boundp.SetBit(sel.PNum(j));
-	  }
+        */
+        for (auto & sel : mesh.SurfaceElements())
+          for (int j = 1; j <= sel.GetNP(); j++)
+            boundp.SetBit(sel.PNum(j));
+        // }
 
 
 	(*testout) << "bpoints:" << endl;
diff --git a/libsrc/meshing/validate.cpp b/libsrc/meshing/validate.cpp
index f3c254a1..7e5dbe7d 100644
--- a/libsrc/meshing/validate.cpp
+++ b/libsrc/meshing/validate.cpp
@@ -198,19 +198,21 @@ namespace netgen
 
     NgArray<int> surfaceindex(np);
     surfaceindex = -1;
-    
+
+    /*
     for (int i = 1; i <= mesh.GetNSE(); i++)
       {
 	const Element2d & sel = mesh.SurfaceElement(i);
-	for (int j = 1; j <= sel.GetNP(); j++)
-	  if(!isedgepoint.Test(sel.PNum(j)))
-	    {
-	      isboundarypoint.SetBit(sel.PNum(j));
-	      surfaceindex[sel.PNum(j) - PointIndex::BASE] = 
-		mesh.GetFaceDescriptor(sel.GetIndex()).SurfNr();
-	    }
-      }
-
+    */
+    for (auto & sel : mesh.SurfaceElements())
+      for (int j = 1; j <= sel.GetNP(); j++)
+        if(!isedgepoint.Test(sel.PNum(j)))
+          {
+            isboundarypoint.SetBit(sel.PNum(j));
+            surfaceindex[sel.PNum(j) - PointIndex::BASE] = 
+              mesh.GetFaceDescriptor(sel.GetIndex()).SurfNr();
+          }
+    
 
 
     Validate(mesh,bad_elements,pure_badness,
@@ -306,9 +308,13 @@ namespace netgen
 	  {
 	    for(int i=0; i<nv.Size(); i++)
 	      *nv[i] = Vec<3>(0,0,0);
+            /*
 	    for (int i = 1; i <= mesh.GetNSE(); i++)
 	      {
 		const Element2d & sel = mesh.SurfaceElement(i);
+            */
+            for (auto & sel : mesh.SurfaceElements())
+              {
 		Vec<3> auxvec = Cross(mesh.Point(sel.PNum(2))-mesh.Point(sel.PNum(1)),
                                       mesh.Point(sel.PNum(3))-mesh.Point(sel.PNum(1)));
 		auxvec.Normalize();