diff --git a/libsrc/meshing/adfront3.cpp b/libsrc/meshing/adfront3.cpp
index f1622f68..9f2925ef 100644
--- a/libsrc/meshing/adfront3.cpp
+++ b/libsrc/meshing/adfront3.cpp
@@ -182,9 +182,13 @@ void AdFront3 :: DeleteFace (INDEX fi)
 {
   nff--;
 
+  /*
   for (int i = 1; i <= faces.Get(fi).Face().GetNP(); i++)
     {
       PointIndex pi = faces.Get(fi).Face().PNum(i);
+  */
+  for (PointIndex pi : faces.Get(fi).Face().PNums())
+    {
       points[pi].RemoveFace();
       if (!points[pi].Valid())
 	delpointl.Append (pi);
@@ -397,16 +401,10 @@ void AdFront3 :: RebuildInternalTables ()
 
 
 
-  int negvol = 0;
-  /*
-  for (int i = PointIndex::BASE; 
-       i < clvol.Size()+PointIndex::BASE; i++)
-  */
+  bool negvol = false;
   for (auto i : clvol.Range())
-    {
-      if (clvol[i] < 0)
-	negvol = 1;
-    }
+    if (clvol[i] < 0)
+      negvol = true;
   
   if (negvol)
     {
@@ -427,8 +425,6 @@ void AdFront3 :: RebuildInternalTables ()
 
 int AdFront3 :: SelectBaseElement ()
 {
-  int i, hi, fstind;
-
   /*
   static int minval = -1;
   static int lasti = 0;
@@ -453,12 +449,12 @@ int AdFront3 :: SelectBaseElement ()
     }
     */
   
-  fstind = 0;
+  int fstind = 0;
 
-  for (i = lasti+1; i <= faces.Size() && !fstind; i++)
+  for (int i = lasti+1; i <= faces.Size() && !fstind; i++)
     if (faces.Elem(i).Valid())
       {
-	hi = faces.Get(i).QualClass() +
+	int hi = faces.Get(i).QualClass() +
 	  points[faces.Get(i).Face().PNum(1)].FrontNr() +
 	  points[faces.Get(i).Face().PNum(2)].FrontNr() +
 	  points[faces.Get(i).Face().PNum(3)].FrontNr();
@@ -474,10 +470,10 @@ int AdFront3 :: SelectBaseElement ()
   if (!fstind)
     {
       minval = INT_MAX;
-      for (i = 1; i <= faces.Size(); i++)
+      for (int i = 1; i <= faces.Size(); i++)
 	if (faces.Elem(i).Valid())
 	  {
-	    hi = faces.Get(i).QualClass() +
+	    int hi = faces.Get(i).QualClass() +
 	      points[faces.Get(i).Face().PNum(1)].FrontNr() +
 	      points[faces.Get(i).Face().PNum(2)].FrontNr() +
 	      points[faces.Get(i).Face().PNum(3)].FrontNr();
@@ -499,9 +495,9 @@ int AdFront3 :: SelectBaseElement ()
 
 int AdFront3 :: GetLocals (int fstind,
 			   Array<Point3d, PointIndex> & locpoints,
-			   NgArray<MiniElement2d> & locfaces,   // local index
+			   Array<MiniElement2d> & locfaces,   // local index
 			   Array<PointIndex, PointIndex> & pindex,
-			   NgArray<INDEX> & findex,
+			   Array<INDEX> & findex,
 			   INDEX_2_HASHTABLE<int> & getconnectedpairs,
 			   float xh,
 			   float relh,
@@ -518,7 +514,7 @@ int AdFront3 :: GetLocals (int fstind,
       hashcreated=1;
     }
 
-  INDEX i, j;
+  INDEX i;
   PointIndex pstind;
   Point3d midp, p0;
 
@@ -599,28 +595,37 @@ int AdFront3 :: GetLocals (int fstind,
 
 
   invpindex.SetSize (points.Size());
+  /*
   for (i = 1; i <= locfaces.Size(); i++)
     for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
       {
 	PointIndex pi = locfaces.Get(i).PNum(j);
         invpindex[pi] = PointIndex::INVALID;
       }
+  */
+  for (auto & f : locfaces)
+    for (int j = 1; j <= f.GetNP(); j++)
+      {
+	PointIndex pi = f.PNum(j);
+        invpindex[pi] = PointIndex::INVALID;
+      }
 
-  for (i = 1; i <= locfaces.Size(); i++)
+  // for (i = 1; i <= locfaces.Size(); i++)
+  for (auto & f : locfaces)
     {
-      for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
+      // for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
+      for (int j = 1; j <= f.GetNP(); j++)
 	{
-	  PointIndex pi = locfaces.Get(i).PNum(j);
+	  // PointIndex pi = locfaces.Get(i).PNum(j);
+          PointIndex pi = f.PNum(j);
 	  if (!invpindex[pi].IsValid())
 	    {
 	      pindex.Append (pi);
               locpoints.Append (points[pi].P());
 	      invpindex[pi] = pindex.Size()-1+IndexBASE<PointIndex>();
             }
-          // locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P());
-          // }
-	  // else
-          locfaces.Elem(i).PNum(j) = invpindex[pi];
+          // locfaces.Elem(i).PNum(j) = invpindex[pi];
+          f.PNum(j) = invpindex[pi];          
 	}
     }
 
@@ -675,10 +680,10 @@ int AdFront3 :: GetLocals (int fstind,
 
 // returns all points connected with fi
 void AdFront3 :: GetGroup (int fi,
-			   NgArray<MeshPoint, PointIndex::BASE> & grouppoints,
-			   NgArray<MiniElement2d> & groupelements,
+			   Array<MeshPoint, PointIndex> & grouppoints,
+			   Array<MiniElement2d> & groupelements,
 			   Array<PointIndex, PointIndex> & pindex,
-			   NgArray<INDEX> & findex) 
+			   Array<INDEX> & findex) 
 {
   // static NgArray<char> pingroup;
   int changed;
@@ -796,8 +801,8 @@ void AdFront3 :: SetStartFront (int /* baseelnp */)
     */
 }
 
-  bool AdFront3 :: PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
-                                  const NgArray<MiniElement2d> &groupfaces) const
+bool AdFront3 :: PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
+                                  const Array<MiniElement2d> &groupfaces) const
 {
   for(auto pi : Range(points))
     {
diff --git a/libsrc/meshing/adfront3.hpp b/libsrc/meshing/adfront3.hpp
index fb61a628..b0acff20 100644
--- a/libsrc/meshing/adfront3.hpp
+++ b/libsrc/meshing/adfront3.hpp
@@ -95,7 +95,8 @@ public:
   const PointIndex PNum (int i) const { return pnum[i-1]; }
   PointIndex & PNum (int i) { return pnum[i-1]; }
   const PointIndex PNumMod (int i) const { return pnum[(i-1)%np]; }
-  auto PNums() const { return NgFlatArray<const PointIndex> (np, &pnum[0]); }
+  auto PNums() { return FlatArray<PointIndex> (np, &pnum[0]); }
+  auto PNums() const { return FlatArray<const PointIndex> (np, &pnum[0]); }
   void Delete () { deleted = true; for (PointIndex & p : pnum) p.Invalidate(); }
   bool IsDeleted () const { return deleted; }
 };
@@ -238,9 +239,9 @@ public:
   ///
   int GetNF() const
   { return nff; }
-  ///
+  /// 1-based
   const MiniElement2d & GetFace (int i) const
-  { return faces.Get(i).Face(); }
+  { return faces[i-1].Face(); }
   const auto & Faces() const { return faces; }
   ///
   void Print () const;
@@ -265,7 +266,7 @@ public:
 			     NgArray<int> & ifaces) const;
 
   bool PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
-                        const NgArray<MiniElement2d>& groupfaces) const;
+                        const Array<MiniElement2d>& groupfaces) const;
 
   ///
   void GetFaceBoundingBox (int i, Box3d & box) const;
@@ -273,9 +274,9 @@ public:
   ///
   int GetLocals (int baseelement,
 		 Array<Point3d, PointIndex> & locpoints,
-                 NgArray<MiniElement2d> & locfaces,   // local index
+                 Array<MiniElement2d> & locfaces,   // local index
                  Array<PointIndex, PointIndex> & pindex,
-                 NgArray<INDEX> & findex,
+                 Array<INDEX> & findex,
 		 INDEX_2_HASHTABLE<int> & connectedpairs,
                  float xh,
 		 float relh,
@@ -283,10 +284,10 @@ public:
   
   ///
   void GetGroup (int fi,
-                 NgArray<MeshPoint, PointIndex::BASE> & grouppoints,
-                 NgArray<MiniElement2d> & groupelements,
+                 Array<MeshPoint, PointIndex> & grouppoints,
+                 Array<MiniElement2d> & groupelements,
                  Array<PointIndex, PointIndex> & pindex,
-                 NgArray<INDEX> & findex);
+                 Array<INDEX> & findex);
 
   ///
   void DeleteFace (INDEX fi);
@@ -298,11 +299,11 @@ public:
   INDEX AddConnectedPair (PointIndices<2> pair);
   ///
   void IncrementClass (INDEX fi)
-  { faces.Elem(fi).IncrementQualClass(); }
+  { faces[fi-1].IncrementQualClass(); }
 
   ///
   void ResetClass (INDEX fi)
-  { faces.Elem(fi).ResetQualClass(); }
+  { faces[fi-1].ResetQualClass(); }
 
   ///
   void SetStartFront (int baseelnp = 0);
diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index c5c8a99c..36293b6c 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -527,12 +527,12 @@ namespace netgen
 		  {
 		    PointIndices<2> e1(el2d[i], el2d[(i+1) % el2d.GetNP()]);
 		    e1.Sort();
-		    INDEX_2 e2;
+		    PointIndices<2> e2;
 		    
 		    for(k = 0; k < idmaps.Size(); k++)
 		      {
-			e2.I1() = (*idmaps[k])[e1.I1()];
-			e2.I2() = (*idmaps[k])[e1.I2()];
+			e2[0] = (*idmaps[k])[e1[0]];
+			e2[1] = (*idmaps[k])[e1[1]];
 			
 			if(e2.I1() == 0 || e2.I2() == 0 ||
 			   e1.I1() == e2.I1() || e1.I2() == e2.I2())
@@ -985,7 +985,7 @@ namespace netgen
 	if(j == 0 || mi.pnums[j+mi.np] < min2)
 	  min2 = mi.pnums[j+mi.np];
 
-	identified = (mi.pnums[j+mi.np] != 0 && mi.pnums[j+mi.np] != mi.pnums[j]);
+	identified = (mi.pnums[j+mi.np].IsValid() && mi.pnums[j+mi.np] != mi.pnums[j]);
       }
 
     identified = identified && (min1 < min2);
@@ -1929,13 +1929,14 @@ namespace netgen
 
     ist >> size;
     mtets.SetSize(size);
+    constexpr auto PI0 = IndexBASE<PointIndex>();
     for(int i=0; i<size; i++)
       {
         ist >> mtets[i];
-        if(mtets[i].pnums[0] > mesh.GetNV() || 
-           mtets[i].pnums[1] > mesh.GetNV() || 
-           mtets[i].pnums[2] > mesh.GetNV() || 
-           mtets[i].pnums[3] > mesh.GetNV())
+        if(mtets[i].pnums[0] >= PI0+mesh.GetNV() || 
+           mtets[i].pnums[1] >= PI0+mesh.GetNV() || 
+           mtets[i].pnums[2] >= PI0+mesh.GetNV() || 
+           mtets[i].pnums[3] >= PI0+mesh.GetNV())
           return false;
       }
 
diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp
index b8f442be..3eaaa508 100644
--- a/libsrc/meshing/boundarylayer.cpp
+++ b/libsrc/meshing/boundarylayer.cpp
@@ -146,7 +146,7 @@ namespace netgen
           auto & seg = mesh[segi];
           if(seg.edgenr != edgenr+1)
               continue;
-          PointIndex other = seg[0]+seg[1]-pi;
+          PointIndex other = seg[0]-pi+seg[1];
           if(!pts.Contains(other))
             pts.Append(other);
       }
@@ -1353,7 +1353,8 @@ namespace netgen
                             nel[1] = p2;
                             nel[2] = p3;
                             nel[3] = p4;
-                            nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1];
+                            // nel[4] = el[0] + el[1] + el[2] + el[3] - fixed[0] - moved[0] - moved[1];
+                            nel[4] = el[0]-fixed[0] + el[1]-moved[0] + el[2]-moved[1] + el[3];
                             if(Cross(mesh[p2]-mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0)
                               Swap(nel[1], nel[3]);
                             nel.SetIndex(el.GetIndex());
@@ -1399,7 +1400,8 @@ namespace netgen
                         nel[1] = p2;
                         nel[2] = p3;
                         nel[3] = p4;
-                        nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - fixed[1] - moved[0] - moved[1];
+                        // nel[4] = el[0] + el[1] + el[2] + el[3] + el[4] - fixed[0] - fixed[1] - moved[0] - moved[1];
+                        nel[4] = el[0]-fixed[0] + el[1]-fixed[1] + el[2]-moved[0] + el[3]-moved[1] + el[4];
                         if(Cross(mesh[p2] - mesh[p1], mesh[p4]-mesh[p1]) * (mesh[nel[4]]-mesh[nel[1]]) > 0)
                           Swap(nel[1], nel[3]);
                         nel.SetIndex(el.GetIndex());
diff --git a/libsrc/meshing/classifyhpel.hpp b/libsrc/meshing/classifyhpel.hpp
index c7fbb3ba..def1c67c 100644
--- a/libsrc/meshing/classifyhpel.hpp
+++ b/libsrc/meshing/classifyhpel.hpp
@@ -10,6 +10,7 @@ HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
   HPREF_ELEMENT_TYPE type = HP_NONE; 
   
   int debug = 0;
+  /*
   for (int j = 0;j < 4; j++)
     {
       if (el.pnums[j] == 444) debug++;
@@ -18,7 +19,7 @@ HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges
       if (el.pnums[j] == 281) debug++;
     }
   if (debug < 4) debug = 0;
-  
+  */
 
   // *testout << "new el" << endl;
   
diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp
index b64fbc3f..845b43d3 100644
--- a/libsrc/meshing/clusters.cpp
+++ b/libsrc/meshing/clusters.cpp
@@ -23,6 +23,7 @@ namespace netgen
     // static int timer2 = NgProfiler::CreateTimer ("clusters2");
     // static int timer3 = NgProfiler::CreateTimer ("clusters3");
     RegionTimer reg (timer);
+    constexpr auto PI0 = IndexBASE<PointIndex>();
 
     const MeshTopology & top = mesh.GetTopology();
 
@@ -103,7 +104,7 @@ namespace netgen
              
              nnums.SetSize(elnv+elned+elnfa+1);
              for (int j = 0; j < elnv; j++)
-               nnums[j] = el[j]+1-PointIndex::BASE;
+               nnums[j] = el[j]+1-PI0;
              for (int j = 0; j < elned; j++)
                nnums[elnv+j] = nv+ednums[j]+1;
              for (int j = 0; j < elnfa; j++)
@@ -131,7 +132,7 @@ namespace netgen
 	  
 	nnums.SetSize(elnv+elned+1);
 	for (int j = 1; j <= elnv; j++)
-	  nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
+	  nnums.Elem(j) = el.PNum(j)+1-PI0;
 	for (int j = 1; j <= elned; j++)
 	  nnums.Elem(elnv+j) = nv+ednums.Elem(j);
 	nnums.Elem(elnv+elned+1) = fanum;
@@ -162,7 +163,7 @@ namespace netgen
              
              nnums.SetSize(elnv+elned+1);
              for (int j = 0; j < elnv; j++)
-               nnums[j] = el[j]+1-PointIndex::BASE;
+               nnums[j] = el[j]+1-PI0;
              for (int j = 0; j < elned; j++)
                nnums[elnv+j] = nv+ednums[j]+1;
              nnums[elnv+elned] = fanum;             
@@ -219,7 +220,6 @@ namespace netgen
       { 2, 3, 1, 1, 4, 5, 1, 6, 4, 5, 5, 4, 7, 7, 7 };
 
     // int cnt = 0;
-
     do
       {
         static Timer t("update cluster, identify");
@@ -247,23 +247,23 @@ namespace netgen
 		break;
 	      case TET:
 	      case TET10:
-		if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == 
-		    cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE))
+		if (cluster_reps.Get(el.PNum(1)+1-PI0) == 
+		    cluster_reps.Get(el.PNum(2)+1-PI0))
 		  clustertab = tet_cluster12;
-		else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == 
-			 cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
+		else if (cluster_reps.Get(el.PNum(1)+1-PI0) == 
+			 cluster_reps.Get(el.PNum(3)+1-PI0))
 		  clustertab = tet_cluster13;
-		else if (cluster_reps.Get(el.PNum(1)+1-PointIndex::BASE) == 
-			 cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
+		else if (cluster_reps.Get(el.PNum(1)+1-PI0) == 
+			 cluster_reps.Get(el.PNum(4)+1-PI0))
 		  clustertab = tet_cluster14;
-		else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == 
-			 cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE))
+		else if (cluster_reps.Get(el.PNum(2)+1-PI0) == 
+			 cluster_reps.Get(el.PNum(3)+1-PI0))
 		  clustertab = tet_cluster23;
-		else if (cluster_reps.Get(el.PNum(2)+1-PointIndex::BASE) == 
-			 cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
+		else if (cluster_reps.Get(el.PNum(2)+1-PI0) == 
+			 cluster_reps.Get(el.PNum(4)+1-PI0))
 		  clustertab = tet_cluster24;
-		else if (cluster_reps.Get(el.PNum(3)+1-PointIndex::BASE) == 
-			 cluster_reps.Get(el.PNum(4)+1-PointIndex::BASE))
+		else if (cluster_reps.Get(el.PNum(3)+1-PI0) == 
+			 cluster_reps.Get(el.PNum(4)+1-PI0))
 		  clustertab = tet_cluster34;
 
 		else
@@ -286,7 +286,7 @@ namespace netgen
                 
                 nnums.SetSize(elnv+elned+elnfa+1);
                 for (int j = 0; j < elnv; j++)
-                  nnums[j] = el[j]+1-PointIndex::BASE;
+                  nnums[j] = el[j]+1-IndexBASE<PointIndex>();
                 for (int j = 0; j < elned; j++)
                   nnums[elnv+j] = nv+ednums[j]+1;
                 for (int j = 0; j < elnfa; j++)
diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp
index c5b5fa18..a3be91a8 100644
--- a/libsrc/meshing/curvedelems.cpp
+++ b/libsrc/meshing/curvedelems.cpp
@@ -1194,9 +1194,10 @@ namespace netgen
 	    // if (el.GetType() == TRIG && order >= 3)
 	    if (top.GetFaceType(facenr+1) == TRIG && order >= 3)
 	      {
-		NgArrayMem<int, 3> verts(3);
-		top.GetFaceVertices (facenr+1, verts);
-
+                // NgArrayMem<PointIndex, 3> verts(3);
+                // top.GetFaceVertices (facenr+1, verts);
+                auto verts = top.GetFaceVertices(facenr);
+                                                 
 		int fnums[] = { 0, 1, 2 };
 		/*
 		if (el[fnums[0]] > el[fnums[1]]) swap (fnums[0], fnums[1]);
diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp
index a6540cf2..afdf8fa1 100644
--- a/libsrc/meshing/delaunay.cpp
+++ b/libsrc/meshing/delaunay.cpp
@@ -34,12 +34,10 @@ namespace netgen
     int NB(int i) const { return nb[i]; }
 
 
-    int FaceNr (INDEX_3 & face) const  // which face nr is it ?
+    int FaceNr (const PointIndices<3> & face) const  // which face nr is it ?
     {
       for (int i = 0; i < 3; i++)
-	if (pnums[i] != face.I1() && 
-	    pnums[i] != face.I2() && 
-	    pnums[i] != face.I3())
+	if (pnums[i] != face[0] && pnums[i] != face[1] && pnums[i] != face[2])
 	  return i;
       return 3;
     }
@@ -673,7 +671,7 @@ namespace netgen
     IndexSet closesphere(mesh.GetNP());
 
     // "random" reordering of points  (speeds a factor 3 - 5 !!!)
-    NgArray<PointIndex, PointIndex::BASE, PointIndex> mixed(np);
+    Array<PointIndex, PointIndex> mixed(np);
     // int prims[] = { 11, 13, 17, 19, 23, 29, 31, 37 };
     // int prims[] = { 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
     int prims[] = { 211, 223, 227, 229, 233, 239, 241, 251, 257, 263 }; 
@@ -919,7 +917,7 @@ namespace netgen
                 {
                     auto & tri_other = mesh.OpenElement(i_other);
                     PointIndex pi2 = tri[(edge+2)%3];
-                    PointIndex pi3 = tri_other[0]+tri_other[1]+tri_other[2] - pi0 - pi1;
+                    PointIndex pi3 = tri_other[0]-pi0+tri_other[1]-pi1+tri_other[2];
                     if(pi2>pi3)
                         Swap(pi2, pi3);
 
diff --git a/libsrc/meshing/delaunay2d.cpp b/libsrc/meshing/delaunay2d.cpp
index c6674719..787fd7cc 100644
--- a/libsrc/meshing/delaunay2d.cpp
+++ b/libsrc/meshing/delaunay2d.cpp
@@ -819,7 +819,7 @@ namespace netgen
       for(auto sei : els)
       {
         auto & el = tempmesh[sei];
-        PointIndex pi2 = el[0]+el[1]+el[2] - seg[0] - seg[1];
+        PointIndex pi2 = el[0]-seg[0]+el[1]-seg[1]+el[2];
         bool is_left = ::netgen::Area(P2(tempmesh[seg[0]]), P2(tempmesh[seg[1]]), P2(tempmesh[pi2]))>0.0;
         POSITION pos;
 
diff --git a/libsrc/meshing/geomsearch.cpp b/libsrc/meshing/geomsearch.cpp
index 2b3fabba..65717753 100644
--- a/libsrc/meshing/geomsearch.cpp
+++ b/libsrc/meshing/geomsearch.cpp
@@ -190,7 +190,7 @@ namespace netgen
     MinCoords(maxextreal,maxp);
 
 
-    int cluster = faces->Get(fstind).Cluster();
+    PointIndex cluster = faces->Get(fstind).Cluster();
   
     int sx = int((minp.X()-minext.X())/elemsize.X()+1.);
     int ex = int((maxp.X()-minext.X())/elemsize.X()+1.);
diff --git a/libsrc/meshing/improve2gen.cpp b/libsrc/meshing/improve2gen.cpp
index 02fa5fd6..6a5f1371 100644
--- a/libsrc/meshing/improve2gen.cpp
+++ b/libsrc/meshing/improve2gen.cpp
@@ -226,7 +226,7 @@ namespace netgen
         */
         for (const Element2d & el : rule.oldels)
           for (PointIndex pi : el.PNums())
-            rule.incelsonnode[pi-PointIndex::BASE]--;
+            rule.incelsonnode[pi-IndexBASE<PointIndex>()]--;
 
 	for (int j = 1; j <= rule.newels.Size(); j++)
 	  {
diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp
index f0f566ef..10175c91 100644
--- a/libsrc/meshing/improve3.cpp
+++ b/libsrc/meshing/improve3.cpp
@@ -1585,11 +1585,11 @@ void MeshOptimize3d :: SwapImproveSurface (
 	  bool found = false;
 	  for(int k=0; !found && k<used_idmaps->Size(); k++)
 	    {
-	      if(pi2 < (*used_idmaps)[k]->Size() + PointIndex::BASE)
+	      if(pi2 < (*used_idmaps)[k]->Size() + IndexBASE<PointIndex>())
 		{
 		  pi1other = (*(*used_idmaps)[k])[pi1];
 		  pi2other = (*(*used_idmaps)[k])[pi2];
-		  found = (pi1other != 0 && pi2other != 0 && pi1other != pi1 && pi2other != pi2);
+		  found = (pi1other.IsValid() && pi2other.IsValid() && pi1other != pi1 && pi2other != pi2);
 		  if(found)
 		    idnum = k;
 		}
diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp
index 2fe2440a..9d03f5b9 100644
--- a/libsrc/meshing/meshing3.cpp
+++ b/libsrc/meshing/meshing3.cpp
@@ -182,10 +182,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 
 
   Array<Point3d, PointIndex> locpoints;      // local points
-  NgArray<MiniElement2d> locfaces;                   // local faces
+  Array<MiniElement2d> locfaces;                   // local faces
   Array<PointIndex, PointIndex> pindex;      // mapping from local to front point numbering
   Array<int, PointIndex> allowpoint;         // point is allowed ?
-  NgArray<INDEX> findex;                             // mapping from local to front face numbering
+  Array<INDEX> findex;                             // mapping from local to front face numbering
   //INDEX_2_HASHTABLE<int> connectedpairs(100);    // connecgted pairs for prism meshing
 
   Array<Point3d, PointIndex> plainpoints;    // points in reference coordinates
@@ -195,7 +195,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
   int j, oldnp, oldnf;
   int found;
   referencetransform trans;
-  int rotind;
   Point3d inp;
   float err;
 
@@ -211,10 +210,10 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 
   
   // for star-shaped domain meshing
-  NgArray<MeshPoint, PointIndex::BASE> grouppoints;      
-  NgArray<MiniElement2d> groupfaces;
+  Array<MeshPoint, PointIndex> grouppoints;      
+  Array<MiniElement2d> groupfaces;
   Array<PointIndex, PointIndex> grouppindex;
-  NgArray<INDEX> groupfindex;
+  Array<INDEX> groupfindex;
   
   
   float minerr;
@@ -313,7 +312,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 
       if (loktestmode)
 	{
-	  (*testout) << "baseel = " << baseelem << ", ind = " << findex.Get(1) << endl;
+	  (*testout) << "baseel = " << baseelem << ", ind = " << findex[0] << endl;
           int pi1 = pindex[locfaces[0].PNum(1)];
           int pi2 = pindex[locfaces[0].PNum(2)];
           int pi3 = pindex[locfaces[0].PNum(3)];
@@ -368,10 +367,11 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 			       grouppindex, groupfindex);
 
 	  bool onlytri = 1;
-	  for (auto i : groupfaces.Range())
-	    if (groupfaces[i].GetNP() != 3) 
-	      onlytri = 0;
-	  
+          for (auto & f : groupfaces)
+            if (f.GetNP() != 3) 
+              onlytri = 0;
+            
+              
 	  if (onlytri && groupfaces.Size() <= 20 + 2*stat.qualclass &&
 	      FindInnerPoint (grouppoints, groupfaces, inp) &&
               !adfront->PointInsideGroup(grouppindex, groupfaces))
@@ -379,11 +379,11 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      (*testout) << "inner point found" << endl;
 
 	      for(int i = 1; i <= groupfaces.Size(); i++)
-		adfront -> DeleteFace (groupfindex.Get(i));
+		adfront -> DeleteFace (groupfindex[i-1]);
 	      
 	      for(int i = 1; i <= groupfaces.Size(); i++)
 		for (j = 1; j <= locfaces.Size(); j++)
-		  if (findex.Get(j) == groupfindex.Get(i))
+		  if (findex[j-1] == groupfindex[i-1])
 		    delfaces.Append (j);
 	      
 	      
@@ -396,13 +396,13 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      newel.SetNP(4);
 	      newel.PNum(4) = npi;
 	      
-	      for(int i = 1; i <= groupfaces.Size(); i++)
+	      for(int i = 0; i < groupfaces.Size(); i++)
 		{
 		  for (j = 1; j <= 3; j++)
 		    {
 		      newel.PNum(j) = 
 			adfront->GetGlobalIndex 
-			(grouppindex[groupfaces.Get(i).PNum(j)]);
+			(grouppindex[groupfaces[i].PNum(j)]);
 		    }
 		  mesh.AddVolumeElement (newel);
 		}
@@ -437,7 +437,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 
       bool impossible = 1;
 
-      for (rotind = 1; rotind <= locfaces[0].GetNP(); rotind++)
+      for (int rotind = 1; rotind <= locfaces[0].GetNP(); rotind++)
 	{
 	  // set transformatino to reference coordinates
 
@@ -687,26 +687,26 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      stat.cntelem++;
 	    }
 
-	  for(int i = oldnf+1; i <= locfaces.Size(); i++)
+	  for(int i = oldnf; i < locfaces.Size(); i++)
 	    {
-	      for (j = 1; j <= locfaces.Get(i).GetNP(); j++)
-		locfaces.Elem(i).PNum(j) = 
-		  pindex[locfaces.Get(i).PNum(j)];
+	      for (j = 1; j <= locfaces[i].GetNP(); j++)
+		locfaces[i].PNum(j) = 
+		  pindex[locfaces[i].PNum(j)];
 	      // (*testout) << "add face " << locfaces.Get(i) << endl;
-	      adfront->AddFace (locfaces.Get(i));
+	      adfront->AddFace (locfaces[i]);
 	    }
 	  
 	  for(int i = 1; i <= delfaces.Size(); i++)
-	    adfront->DeleteFace (findex.Get(delfaces.Get(i)));
+	    adfront->DeleteFace (findex[delfaces.Get(i)-1]);
 	}
       else
 	{
-	  adfront->IncrementClass (findex.Get(1));
+	  adfront->IncrementClass (findex[0]);
 	  if (impossible && mp.check_impossible)
 	    {
 	      (*testout) << "skip face since it is impossible" << endl;
 	      for (j = 0; j < 100; j++)
-		adfront->IncrementClass (findex.Get(1));
+		adfront->IncrementClass (findex[0]);
 	    }
 	}
 
@@ -938,7 +938,7 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh)
 
   for(int i = 1; i <= n; i++)
     {
-      pointnr.Elem(i) = 0;
+      pointnr.Elem(i) = PointIndex::INVALID;
       frontpointnr.Elem(i) = 0;
     }
   
@@ -954,7 +954,7 @@ void Meshing3 :: BlockFill (Mesh & mesh, double gh)
 		  for (j3 = i3; j3 <= i3+1; j3++)
 		    {
 		      j = j3 + (j2-1) * n3 + (j1-1) * n2 * n3;
-		      if (pointnr.Get(j) == 0)
+		      if (!pointnr.Get(j).IsValid())
 			{
 			  Point3d hp(xmin + (j1-1) * gh, 
 				     ymin + (j2-1) * gh, 
diff --git a/libsrc/meshing/meshing3.hpp b/libsrc/meshing/meshing3.hpp
index 85aa4a3f..2289826b 100644
--- a/libsrc/meshing/meshing3.hpp
+++ b/libsrc/meshing/meshing3.hpp
@@ -47,7 +47,7 @@ public:
   ///
   int ApplyRules (Array<Point3d, PointIndex> & lpoints,
                   Array<int, PointIndex> & allowpoint,
-		  NgArray<MiniElement2d> & lfaces, INDEX lfacesplit,
+		  Array<MiniElement2d> & lfaces, INDEX lfacesplit,
 		  INDEX_2_HASHTABLE<int> & connectedpairs,
 		  NgArray<Element> & elements,
 		  NgArray<INDEX> & delfaces, int tolerance, 
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index cb750166..37f6736a 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -221,7 +221,6 @@ namespace netgen
     friend istream & operator>> (istream &, PointIndex &);    
     friend ostream & operator<< (ostream &, const PointIndex &);
     template <int N> friend class PointIndices;
-
     /*
     friend PointIndex operator+ (PointIndex, int);
     friend PointIndex operator+ (PointIndex, size_t);    
@@ -235,7 +234,6 @@ namespace netgen
     friend bool operator== (PointIndex a, PointIndex b);
     friend bool operator!= (PointIndex a, PointIndex b);
     */
-    
   public:
     constexpr PointIndex (t_invalid inv) : i(PointIndex::BASE-1) { ; }
     // PointIndex & operator= (const PointIndex &ai) { i = ai.i; return *this; }
@@ -250,6 +248,7 @@ namespace netgen
     PointIndex operator+= (int add) { i += add; return *this; }
     void Invalidate() { i = PointIndex::BASE-1; }
     bool IsValid() const { return i != PointIndex::BASE-1; }
+    // operator bool() const { return IsValid(); }
 #ifdef BASE0
     static constexpr size_t BASE = 0;
 #else
@@ -258,7 +257,7 @@ namespace netgen
 
     void DoArchive (Archive & ar) { ar & i; }
   };
-  
+
   /*
   inline PointIndex operator+ (PointIndex pi, int i) { return PointIndex(pi.i+i); }
   inline PointIndex operator+ (PointIndex pi, size_t i) { return PointIndex(pi.i+i); }  
@@ -308,6 +307,19 @@ namespace netgen
     template <size_t J>
     PointIndex get() const { return PointIndex(INDEX_2::operator[](J)); }    
   };
+  template <> class PointIndices<3> : public INDEX_3
+  {
+  public:
+    PointIndices () = default;
+    PointIndices (INDEX_3 i3) : INDEX_3(i3) { ; }
+    PointIndices (PointIndex i1, PointIndex i2, PointIndex i3) : INDEX_3(i1,i2,i3) { ; } 
+    PointIndex operator[] (int i) const { return PointIndex(INDEX_3::operator[](i)); }
+    PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_3::operator[](i)); }
+    using INDEX_3::Sort;
+    static PointIndices Sort(PointIndex i1, PointIndex i2, PointIndex i3) { return INDEX_3::Sort(i1, i2, i3); }
+    template <size_t J>
+    PointIndex get() const { return PointIndex(INDEX_3::operator[](J)); }    
+  };
 }
 
 namespace std
diff --git a/libsrc/meshing/ruler3.cpp b/libsrc/meshing/ruler3.cpp
index 188e6a25..20a58f17 100644
--- a/libsrc/meshing/ruler3.cpp
+++ b/libsrc/meshing/ruler3.cpp
@@ -51,7 +51,7 @@ int Meshing3 :: ApplyRules
 (
  Array<Point3d, PointIndex> & lpoints,     // in: local points, out: old+new local points
  Array<int, PointIndex> & allowpoint,      // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means
- NgArray<MiniElement2d> & lfaces,    // in: local faces, out: old+new local faces
+ Array<MiniElement2d> & lfaces,    // in: local faces, out: old+new local faces
  INDEX lfacesplit,	       // for local faces in outer radius
  INDEX_2_HASHTABLE<int> & connectedpairs,  // connected pairs for prism-meshing
  NgArray<Element> & elements,    // out: new elements
@@ -314,7 +314,7 @@ int Meshing3 :: ApplyRules
 		  
 		  if (fnearness.Get(locfi) > rule->GetFNearness (nfok) ||
 		      fused.Get(locfi) ||
-		      actfnp != lfaces.Get(locfi).GetNP() )
+		      actfnp != lfaces[locfi-1].GetNP() )
 		    {
 		      // face not feasible in any rotation
 
@@ -325,7 +325,7 @@ int Meshing3 :: ApplyRules
 		      
 		      ok = 1;
 		      
-		      locface = &lfaces.Get(locfi);
+		      locface = &lfaces[locfi-1];
 
 		      
 		      // reference point already mapped differently ?
@@ -528,7 +528,7 @@ int Meshing3 :: ApplyRules
 			    (*testout) << pi << " ";
 			  (*testout) << endl;
 			  snprintf (problems.Elem(ri), 255, "mapping found");
-			  (*testout) << rule->GetNP(1) << " = " << lfaces.Get(1).GetNP() << endl;
+			  (*testout) << rule->GetNP(1) << " = " << lfaces[0].GetNP() << endl;
 			}
 		      
 		      ok = 1;
@@ -691,7 +691,7 @@ int Meshing3 :: ApplyRules
 			  if (!fused.Get(i))
 			    { 
 			      int triin;
-			      const MiniElement2d & lfacei = lfaces.Get(i);
+			      const MiniElement2d & lfacei = lfaces[i-1];
 
 			      if (!triboxes.Elem(i).Intersect (rule->fzbox))
 				triin = 0;
@@ -744,19 +744,19 @@ int Meshing3 :: ApplyRules
 
 				  if (loktestmode)
 				    {
-				      (*testout) << "El with " << lfaces.Get(i).GetNP() << " points in freezone: "
-						 << lfaces.Get(i).PNum(1) << " - " 
-						 << lfaces.Get(i).PNum(2) << " - "
-						 << lfaces.Get(i).PNum(3) << " - "
-						 << lfaces.Get(i).PNum(4) << endl;
-				      for (int lj = 1; lj <= lfaces.Get(i).GetNP(); lj++)
-					(*testout) << lpoints[lfaces.Get(i).PNum(lj)] << " ";
+				      (*testout) << "El with " << lfaces[i-1].GetNP() << " points in freezone: "
+						 << lfaces[i-1].PNum(1) << " - " 
+						 << lfaces[i-1].PNum(2) << " - "
+						 << lfaces[i-1].PNum(3) << " - "
+						 << lfaces[i-1].PNum(4) << endl;
+				      for (int lj = 1; lj <= lfaces[i-1].GetNP(); lj++)
+					(*testout) << lpoints[lfaces[i-1].PNum(lj)] << " ";
 
 				      (*testout) << endl;
 
 				      sprintf (problems.Elem(ri), "triangle (%d, %d, %d) in Freezone",
-					       lfaces.Get(i).PNum(1), lfaces.Get(i).PNum(2),
-					       lfaces.Get(i).PNum(3));
+					       lfaces[i-1].PNum(1), lfaces[i-1].PNum(2),
+					       lfaces[i-1].PNum(3));
 				    }
 #else
 				  if (loktestmode)
@@ -789,9 +789,9 @@ int Meshing3 :: ApplyRules
 						     << endl;
 
 				      snprintf (problems.Elem(ri), 255, "triangle (%d, %d, %d) in Freezone",
-					       int(lfaces.Get(i).PNum(1)), 
-					       int(lfaces.Get(i).PNum(2)),
-					       int(lfaces.Get(i).PNum(3)));
+					       int(lfaces[i-1].PNum(1)), 
+					       int(lfaces[i-1].PNum(2)),
+					       int(lfaces[i-1].PNum(3)));
 				    }	
 
 				  hc = 0;
@@ -802,9 +802,9 @@ int Meshing3 :: ApplyRules
 					  rule->GetPointNr(k, 3) <= rule->GetNOldP())
 					{
 					  for (int j = 1; j <= 3; j++)
-					    if (lfaces.Get(i).PNumMod(j  ) == pmap.Get(rule->GetPointNr(k, 1)) &&
-						lfaces.Get(i).PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) &&
-						lfaces.Get(i).PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2)))
+					    if (lfaces[i-1].PNumMod(j  ) == pmap.Get(rule->GetPointNr(k, 1)) &&
+						lfaces[i-1].PNumMod(j+1) == pmap.Get(rule->GetPointNr(k, 3)) &&
+						lfaces[i-1].PNumMod(j+2) == pmap.Get(rule->GetPointNr(k, 2)))
 					      {
 						fmapi.Elem(k) = i;
 						hc = 1;
@@ -827,14 +827,14 @@ int Meshing3 :: ApplyRules
 				      if (loktestmode)
 					{
 					  (*testout) << "Triangle in freezone: "
-						     << lfaces.Get(i).PNum(1) << " - " 
-						     << lfaces.Get(i).PNum(2) << " - "
-						     << lfaces.Get(i).PNum(3) << endl;
+						     << lfaces[i-1].PNum(1) << " - " 
+						     << lfaces[i-1].PNum(2) << " - "
+						     << lfaces[i-1].PNum(3) << endl;
 
 					  snprintf (problems.Elem(ri), 255, "triangle (%d, %d, %d) in Freezone",
-						   int (lfaces.Get(i).PNum(1)), 
-						   int (lfaces.Get(i).PNum(2)),
-						   int (lfaces.Get(i).PNum(3)));
+						   int (lfaces[i-1].PNum(1)), 
+						   int (lfaces[i-1].PNum(2)),
+						   int (lfaces[i-1].PNum(3)));
 					}
 				      ok = 0;
 				    }
diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp
index ccf47ce4..761158cf 100644
--- a/libsrc/meshing/topology.hpp
+++ b/libsrc/meshing/topology.hpp
@@ -137,6 +137,8 @@ public:
   
   DLL_HEADER void GetFaceVertices (int fnr, NgArray<int> & vertices) const;
   DLL_HEADER void GetFaceVertices (int fnr, int * vertices) const;
+  auto GetFaceVertices (int fnr) const
+  { return FlatArray (face2vert[fnr][3].IsValid() ? 4 : 3, &face2vert[fnr][0]); }
   [[deprecated("use GetEdgeVertices -> tupe(v0,v1) instead")]]                            
   DLL_HEADER void GetEdgeVertices (int enr, int & v1, int & v2) const;
   [[deprecated("use GetEdgeVertices -> tupe(v0,v1) instead")]]