diff --git a/libsrc/meshing/adfront2.cpp b/libsrc/meshing/adfront2.cpp
index 390e0879..5a9b6cab 100644
--- a/libsrc/meshing/adfront2.cpp
+++ b/libsrc/meshing/adfront2.cpp
@@ -152,15 +152,15 @@ namespace netgen
 
     if (allflines)
       {
-	if (allflines->Used (INDEX_2 (GetGlobalIndex (pi1), 
-				      GetGlobalIndex (pi2))))
+	if (allflines->Used (PointIndices<2>(GetGlobalIndex (pi1), 
+                                             GetGlobalIndex (pi2))))
 	  {
 	    cerr << "ERROR Adfront2::AddLine: line exists" << endl;
 	    (*testout) << "ERROR Adfront2::AddLine: line exists" << endl;
 	  }
 
-	allflines->Set (INDEX_2 (GetGlobalIndex (pi1), 
-				 GetGlobalIndex (pi2)), 1);
+	allflines->Set (PointIndices<2>(GetGlobalIndex (pi1), 
+                                        GetGlobalIndex (pi2)), 1);
       }
 
     return li;
@@ -194,8 +194,8 @@ namespace netgen
 
     if (allflines)
       {
-	allflines->Set (INDEX_2 (GetGlobalIndex (lines[li].L().I1()),
-				 GetGlobalIndex (lines[li].L().I2())), 2);
+	allflines->Set (PointIndices<2>(GetGlobalIndex (lines[li].L().I1()),
+                                        GetGlobalIndex (lines[li].L().I2())), 2);
       }
 
     lines[li].Invalidate();
diff --git a/libsrc/meshing/adfront3.cpp b/libsrc/meshing/adfront3.cpp
index 31370930..bcf9a67b 100644
--- a/libsrc/meshing/adfront3.cpp
+++ b/libsrc/meshing/adfront3.cpp
@@ -13,7 +13,7 @@ FrontPoint3 :: FrontPoint3 ()
   globalindex.Invalidate(); //  = -1;
   nfacetopoint = 0; 
   frontnr = 1000; 
-  cluster = 0;
+  cluster = PointIndex::INVALID;
 }
 
 
@@ -23,7 +23,7 @@ FrontPoint3 :: FrontPoint3 (const Point<3> & ap, PointIndex agi)
   globalindex = agi;
   nfacetopoint = 0; 
   frontnr = 1000; 
-  cluster = 0;
+  cluster = PointIndex::INVALID;
 }
 
 
@@ -35,7 +35,7 @@ FrontFace :: FrontFace ()
   qualclass = 1; 
   oldfront = 0; 
   hashvalue = 0;
-  cluster = 0;
+  cluster = PointIndex::INVALID;
 }
 
 FrontFace :: FrontFace (const MiniElement2d & af)
@@ -71,7 +71,7 @@ AdFront3 :: AdFront3 ()
     hashtable.Init(&points, &faces);
 
   facetree = NULL;
-  connectedpairs = NULL;
+  // connectedpairs = NULL;
 
   rebuildcounter = -1;
   lasti = 0;
@@ -82,7 +82,7 @@ AdFront3 :: AdFront3 ()
 AdFront3 :: ~AdFront3 ()
 {
   delete facetree;
-  delete connectedpairs;
+  // delete connectedpairs;
 }
 
 void AdFront3 :: GetPoints (NgArray<Point<3> > & apoints) const
@@ -153,10 +153,10 @@ INDEX AdFront3 :: AddFace (const MiniElement2d & aface)
     }
 
 
-  int cluster = 0;
+  PointIndex cluster = PointIndex::INVALID;
   for (i = 1; i <= aface.GetNP(); i++)
     {
-      if (points[aface.PNum(i)].cluster)
+      if (points[aface.PNum(i)].cluster.IsValid())
 	cluster = points[aface.PNum(i)].cluster;
     }
   for (i = 1; i <= aface.GetNP(); i++)
@@ -213,13 +213,13 @@ void AdFront3 :: DeleteFace (INDEX fi)
 }
 
 
-INDEX AdFront3 :: AddConnectedPair (const INDEX_2 & apair)
+INDEX AdFront3 :: AddConnectedPair (PointIndices<2> apair)
 {
   if (!connectedpairs)
-    connectedpairs = new TABLE<int, PointIndex::BASE> (GetNP());
+    connectedpairs = make_unique<DynamicTable<PointIndex, PointIndex>> (GetNP());
 
-  connectedpairs->Add (apair.I1(), apair.I2());
-  connectedpairs->Add (apair.I2(), apair.I1());
+  connectedpairs->Add (apair[0], apair[1]);
+  connectedpairs->Add (apair[1], apair[0]);
 
   return 0;
 }
@@ -231,11 +231,11 @@ void AdFront3 :: CreateTrees ()
   PointIndex pi;
   Point3d pmin, pmax;
 
-  for (pi = PointIndex::BASE; 
-       pi < GetNP()+PointIndex::BASE; pi++)
+  for (pi = IndexBASE<PointIndex>(); 
+       pi < GetNP()+IndexBASE<PointIndex>(); pi++)
     {
       const Point<3> & p = GetPoint(pi);
-      if (pi == PointIndex::BASE)
+      if (pi == IndexBASE<PointIndex>())
 	{
 	  pmin = p;
 	  pmax = p;
@@ -323,12 +323,12 @@ void AdFront3 :: RebuildInternalTables ()
 	{
 	  const MiniElement2d & el = faces.Get(i).Face();
 
-	  int mini = points[el.PNum(1)].cluster;
-	  int maxi = mini;
+	  PointIndex mini = points[el.PNum(1)].cluster;
+	  PointIndex maxi = mini;
 	  
 	  for (int j = 2; j <= 3; j++)
 	    {
-	      int ci = points[el.PNum(j)].cluster;
+	      PointIndex ci = points[el.PNum(j)].cluster;
 	      if (ci < mini) mini = ci;
 	      if (ci > maxi) maxi = ci;
 	    }
@@ -366,7 +366,7 @@ void AdFront3 :: RebuildInternalTables ()
       cntcl++;
   */
   
-  NgArray<double, PointIndex::BASE> clvol (np);
+  Array<double, PointIndex> clvol (np);
   clvol = 0.0;
 
   for (int i = 1; i <= faces.Size(); i++)
@@ -398,8 +398,11 @@ void AdFront3 :: RebuildInternalTables ()
 
 
   int negvol = 0;
+  /*
   for (int i = PointIndex::BASE; 
        i < clvol.Size()+PointIndex::BASE; i++)
+  */
+  for (auto i : clvol.Range())
     {
       if (clvol[i] < 0)
 	negvol = 1;
@@ -408,10 +411,10 @@ void AdFront3 :: RebuildInternalTables ()
   if (negvol)
     {
       for (int i = 1; i <= faces.Size(); i++)
-	faces.Elem(i).cluster = 1;
+	faces.Elem(i).cluster = IndexBASE<PointIndex>();
       // for (PointIndex pi = points.Begin(); pi < points.End(); pi++)
       for (PointIndex pi : points.Range())
-	points[pi].cluster = 1;
+	points[pi].cluster = IndexBASE<PointIndex>();
     }
 
   if (hashon) 
@@ -495,9 +498,9 @@ int AdFront3 :: SelectBaseElement ()
 
 
 int AdFront3 :: GetLocals (int fstind,
-			   NgArray<Point3d, PointIndex::BASE> & locpoints,
+			   Array<Point3d, PointIndex> & locpoints,
 			   NgArray<MiniElement2d> & locfaces,   // local index
-			   NgArray<PointIndex, PointIndex::BASE> & pindex,
+			   Array<PointIndex, PointIndex> & pindex,
 			   NgArray<INDEX> & findex,
 			   INDEX_2_HASHTABLE<int> & getconnectedpairs,
 			   float xh,
@@ -529,7 +532,7 @@ int AdFront3 :: GetLocals (int fstind,
   locfaces3.SetSize(0);
   findex2.SetSize(0);
 
-  int cluster = faces.Get(fstind).cluster;
+  PointIndex cluster = faces.Get(fstind).cluster;
 
   pstind = faces.Get(fstind).Face().PNum(1);
   p0 = points[pstind].P();
@@ -612,7 +615,7 @@ int AdFront3 :: GetLocals (int fstind,
 	    {
 	      pindex.Append (pi);
               locpoints.Append (points[pi].P());
-	      invpindex[pi] = pindex.Size()-1+PointIndex::BASE;
+	      invpindex[pi] = pindex.Size()-1+IndexBASE<PointIndex>();
             }
           // locfaces.Elem(i).PNum(j) = locpoints.Append (points[pi].P());
           // }
@@ -625,22 +628,25 @@ int AdFront3 :: GetLocals (int fstind,
 
   if (connectedpairs)
     {
-      for (i = 1; i <= locpoints.Size(); i++)
+      // for (i = 1; i <= locpoints.Size(); i++)
+      for (auto i : locpoints.Range())
 	{
-	  int pind = pindex.Get(i);
-	  if (pind >= 1 && pind <= connectedpairs->Size ())
+	  PointIndex pind = pindex[i]; // .Get(i);
+	  // if (pind.IsValid() && pind <= connectedpairs->Size ())
+          if (connectedpairs->Range().Contains(pind))
 	    {
 	      for (j = 1; j <= connectedpairs->EntrySize(pind); j++)
 		{
-		  int oi = connectedpairs->Get(pind, j);
-		  int other = invpindex.Get(oi);
-		  if (other >= 1 && other <= pindex.Size() && 
-		      pindex.Get(other) == oi)
+		  PointIndex oi = connectedpairs->Get(pind, j);
+		  PointIndex other = invpindex[oi];
+		  // if (other >= 1 && other <= pindex.Size() &&
+                  if (pindex.Range().Contains(other) &&
+		      pindex[other] == oi)
 		    {
 		      // INDEX_2 coned(i, other);
 		      // coned.Sort();
 		      // (*testout) << "connected: " << locpoints.Get(i) << "-" << locpoints.Get(other) << endl;
-		      getconnectedpairs.Set (INDEX_2::Sort (i, other), 1);
+		      getconnectedpairs.Set (PointIndices<2>::Sort (i, other), 1);
 		    }
 		}
 	    }
@@ -669,7 +675,7 @@ int AdFront3 :: GetLocals (int fstind,
 void AdFront3 :: GetGroup (int fi,
 			   NgArray<MeshPoint, PointIndex::BASE> & grouppoints,
 			   NgArray<MiniElement2d> & groupelements,
-			   NgArray<PointIndex, PointIndex::BASE> & pindex,
+			   Array<PointIndex, PointIndex> & pindex,
 			   NgArray<INDEX> & findex) 
 {
   // static NgArray<char> pingroup;
@@ -734,8 +740,9 @@ void AdFront3 :: GetGroup (int fi,
     if (points[pi].Valid())
       {
 	grouppoints.Append (points[pi].P());
-	pindex.Append (pi);
-	invpindex[pi] = pindex.Size();
+        pindex.Append (pi);
+        // invpindex[pi] = pindex.Size();
+	invpindex[pi] = pindex.Size()-1 + IndexBASE<PointIndex>();
       }
 
   for (int i = 1; i <= faces.Size(); i++)
@@ -787,7 +794,7 @@ void AdFront3 :: SetStartFront (int /* baseelnp */)
     */
 }
 
-bool AdFront3 :: PointInsideGroup(const NgArray<PointIndex, PointIndex::BASE> &grouppindex,
+  bool AdFront3 :: PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
                                   const NgArray<MiniElement2d> &groupfaces) const
 {
   for(auto pi : Range(points))
@@ -797,7 +804,7 @@ bool AdFront3 :: PointInsideGroup(const NgArray<PointIndex, PointIndex::BASE> &g
       for(const auto& f : groupfaces)
         {
           for(auto i : Range(3))
-            if(grouppindex.Get(f.PNum(i+1)) == pi)
+            if(grouppindex[f.PNum(i+1)] == pi)
               {
                 found = true;
                 break;
@@ -814,9 +821,9 @@ bool AdFront3 :: PointInsideGroup(const NgArray<PointIndex, PointIndex::BASE> &g
       int count = 0;
       for(const auto& f : groupfaces)
         {
-          const auto& p1 = points[grouppindex.Get(f.PNum(1))].P();
-          auto v1 = points[grouppindex.Get(f.PNum(2))].P() - p1;
-          auto v2 = points[grouppindex.Get(f.PNum(3))].P() - p1;
+          const auto& p1 = points[grouppindex[f.PNum(1)]].P();
+          auto v1 = points[grouppindex[f.PNum(2)]].P() - p1;
+          auto v2 = points[grouppindex[f.PNum(3)]].P() - p1;
           for(auto i : Range(3))
             {
               a(i,0) = v1[i];
diff --git a/libsrc/meshing/adfront3.hpp b/libsrc/meshing/adfront3.hpp
index 2e209d4a..fb61a628 100644
--- a/libsrc/meshing/adfront3.hpp
+++ b/libsrc/meshing/adfront3.hpp
@@ -31,7 +31,7 @@ class FrontPoint3
   /// distance to original boundary
   int frontnr;
   /// 
-  int cluster;
+  PointIndex cluster;
 public:
   ///
   FrontPoint3 ();
@@ -125,7 +125,7 @@ private:
   ///
   int hashvalue;
   ///
-  int cluster;
+  PointIndex cluster;
   
 public:
   ///
@@ -172,7 +172,7 @@ public:
   ///
   friend class AdFront3;
 
-  int Cluster () const { return cluster; }
+  PointIndex Cluster () const { return cluster; }
 };  
 
 
@@ -182,14 +182,16 @@ public:
 class AdFront3
 {
   ///
-  NgArray<FrontPoint3, PointIndex::BASE, PointIndex> points;
+  // NgArray<FrontPoint3, PointIndex::BASE, PointIndex> points;
+  Array<FrontPoint3, PointIndex> points;
   ///
   NgArray<FrontFace> faces;
   ///
-  NgArray<PointIndex> delpointl;
+  Array<PointIndex> delpointl;
   
   /// which points are connected to pi ?
-  TABLE<int, PointIndex::BASE> * connectedpairs;
+  // TABLE<PointIndex, PointIndex::BASE> * connectedpairs;
+  unique_ptr<DynamicTable<PointIndex, PointIndex>> connectedpairs;
   
   /// number of total front faces;
   int nff;
@@ -214,8 +216,8 @@ class AdFront3
   int lasti;
   /// minimal selection-value of baseelements
   int minval;
-  NgArray<PointIndex, PointIndex::BASE, PointIndex> invpindex;
-  NgArray<char, PointIndex::BASE> pingroup;
+  Array<PointIndex, PointIndex> invpindex;
+  Array<char, PointIndex> pingroup;
   
   ///
   class BoxTree<3> * facetree;
@@ -262,7 +264,7 @@ public:
   void GetIntersectingFaces (const Point<3> & pmin, const Point<3> & pmax, 
 			     NgArray<int> & ifaces) const;
 
-  bool PointInsideGroup(const NgArray<PointIndex, PointIndex::BASE> &grouppindex,
+  bool PointInsideGroup(const Array<PointIndex, PointIndex> &grouppindex,
                         const NgArray<MiniElement2d>& groupfaces) const;
 
   ///
@@ -270,9 +272,9 @@ public:
 
   ///
   int GetLocals (int baseelement,
-		 NgArray<Point3d, PointIndex::BASE> & locpoints,
+		 Array<Point3d, PointIndex> & locpoints,
                  NgArray<MiniElement2d> & locfaces,   // local index
-                 NgArray<PointIndex, PointIndex::BASE> & pindex,
+                 Array<PointIndex, PointIndex> & pindex,
                  NgArray<INDEX> & findex,
 		 INDEX_2_HASHTABLE<int> & connectedpairs,
                  float xh,
@@ -283,7 +285,7 @@ public:
   void GetGroup (int fi,
                  NgArray<MeshPoint, PointIndex::BASE> & grouppoints,
                  NgArray<MiniElement2d> & groupelements,
-                 NgArray<PointIndex, PointIndex::BASE> & pindex,
+                 Array<PointIndex, PointIndex> & pindex,
                  NgArray<INDEX> & findex);
 
   ///
@@ -293,7 +295,7 @@ public:
   ///
   INDEX AddFace (const MiniElement2d & e);
   ///
-  INDEX AddConnectedPair (const INDEX_2 & pair);
+  INDEX AddConnectedPair (PointIndices<2> pair);
   ///
   void IncrementClass (INDEX fi)
   { faces.Elem(fi).IncrementQualClass(); }
diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index 409750f3..3e659695 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -2174,7 +2174,7 @@ namespace netgen
                 auto seg = mesh[j];
                 for (auto map : idmaps)
                   {
-                    if (seg[0] > 0 && seg[1] > 0 && (*map)[seg[0]] && (*map)[seg[1]])
+                    if (seg[0].IsValid() && seg[1].IsValid() && (*map)[seg[0]] && (*map)[seg[1]])
                       {
                         MarkedIdentification mi;
                         mi.np = 2;
@@ -4013,7 +4013,7 @@ namespace netgen
 	do_repair = false;
 	for(int ii=1; ii<=mesh.GetNP(); ii++)
 	  {
-	    if(isnewpoint.Test(ii) && mesh.mlbetweennodes[ii][0] > 0)
+	    if(isnewpoint.Test(ii) && mesh.mlbetweennodes[ii][0].IsValid())
 	      {
 		mesh.Point(ii) = Center(mesh.Point(mesh.mlbetweennodes[ii][0]),
                                         mesh.Point(mesh.mlbetweennodes[ii][1]));
diff --git a/libsrc/meshing/geomsearch.cpp b/libsrc/meshing/geomsearch.cpp
index 88e9ebe9..2b3fabba 100644
--- a/libsrc/meshing/geomsearch.cpp
+++ b/libsrc/meshing/geomsearch.cpp
@@ -20,7 +20,7 @@ namespace netgen
       } 
   }
 
-  void GeomSearch3d :: Init (NgArray <FrontPoint3,PointIndex::BASE, PointIndex> *pointsi, NgArray <FrontFace> *facesi)
+  void GeomSearch3d :: Init (Array <FrontPoint3,PointIndex> *pointsi, NgArray <FrontFace> *facesi)
   {
     points = pointsi;
     faces = facesi;
diff --git a/libsrc/meshing/geomsearch.hpp b/libsrc/meshing/geomsearch.hpp
index fc817b5b..858107de 100644
--- a/libsrc/meshing/geomsearch.hpp
+++ b/libsrc/meshing/geomsearch.hpp
@@ -27,7 +27,7 @@ public:
   virtual ~GeomSearch3d();
 
   ///
-  void Init (NgArray <FrontPoint3,PointIndex::BASE, PointIndex> *pointsi, NgArray <FrontFace> *facesi);
+  void Init (Array <FrontPoint3,PointIndex> *pointsi, NgArray <FrontFace> *facesi);
 
   ///get elements max extension
   void ElemMaxExt(Point3d& minp, Point3d& maxp, const MiniElement2d& elem);
@@ -52,7 +52,7 @@ public:
 private:
   
   NgArray <FrontFace> *faces; // Pointers to Arrays in Adfront
-  NgArray <FrontPoint3,PointIndex::BASE, PointIndex> *points;
+  Array <FrontPoint3,PointIndex> *points;
 
   NgArray <NgArray <int>*> hashtable;
 
diff --git a/libsrc/meshing/improve2.hpp b/libsrc/meshing/improve2.hpp
index 539b0055..1cd74707 100644
--- a/libsrc/meshing/improve2.hpp
+++ b/libsrc/meshing/improve2.hpp
@@ -62,7 +62,7 @@ void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & element
         }
         QuickSort(local_edges);
 
-        auto edge_prev = std::make_tuple<PointIndex, PointIndex>(-1,-1);
+        auto edge_prev = std::make_tuple<PointIndex, PointIndex>(PointIndex::INVALID, PointIndex::INVALID);
 
         for(auto edge : local_edges)
             if(edge != edge_prev)
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 4be51c87..c3a241cd 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -564,7 +564,7 @@ namespace netgen
 
     bool IsSegment (PointIndex pi1, PointIndex pi2) const
     {
-      INDEX_2 i2 (pi1, pi2);
+      PointIndices<2> i2 (pi1, pi2);
       i2.Sort();
       return segmentht->Used (i2);
     }
diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp
index 35197cfd..2fe2440a 100644
--- a/libsrc/meshing/meshing3.cpp
+++ b/libsrc/meshing/meshing3.cpp
@@ -181,14 +181,14 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
   // NgProfiler::RegionTimer reg (meshing3_timer);
 
 
-  NgArray<Point3d, PointIndex::BASE> locpoints;      // local points
+  Array<Point3d, PointIndex> locpoints;      // local points
   NgArray<MiniElement2d> locfaces;                   // local faces
-  NgArray<PointIndex, PointIndex::BASE> pindex;      // mapping from local to front point numbering
-  NgArray<int, PointIndex::BASE> allowpoint;         // point is allowed ?
+  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
   //INDEX_2_HASHTABLE<int> connectedpairs(100);    // connecgted pairs for prism meshing
 
-  NgArray<Point3d, PointIndex::BASE> plainpoints;    // points in reference coordinates
+  Array<Point3d, PointIndex> plainpoints;    // points in reference coordinates
   NgArray<int> delpoints, delfaces;   // points and lines to be deleted
   NgArray<Element> locelements;       // new generated elements
 
@@ -213,7 +213,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
   // for star-shaped domain meshing
   NgArray<MeshPoint, PointIndex::BASE> grouppoints;      
   NgArray<MiniElement2d> groupfaces;
-  NgArray<PointIndex, PointIndex::BASE> grouppindex;
+  Array<PointIndex, PointIndex> grouppindex;
   NgArray<INDEX> groupfindex;
   
   
@@ -342,10 +342,11 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 
       allowpoint.SetSize(locpoints.Size());
       if (uselocalh && stat.qualclass <= 3)
-	for(int i = 1; i <= allowpoint.Size(); i++)
+	//for(int i = 1; i <= allowpoint.Size(); i++)
+        for(auto i : allowpoint.Range())
 	  {
-	    allowpoint.Elem(i) =
-	      (mesh.GetH (locpoints.Get(i)) > 0.4 * hshould / mp.sloppy) ? 2 : 1;
+	    allowpoint[i] =
+	      (mesh.GetH (locpoints[i]) > 0.4 * hshould / mp.sloppy) ? 2 : 1;
 	  }
       else
 	allowpoint = 2;
@@ -401,7 +402,7 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 		    {
 		      newel.PNum(j) = 
 			adfront->GetGlobalIndex 
-			(grouppindex.Get(groupfaces.Get(i).PNum(j)));
+			(grouppindex[groupfaces.Get(i).PNum(j)]);
 		    }
 		  mesh.AddVolumeElement (newel);
 		}
@@ -505,9 +506,13 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	  if (found) stat.cntsucc++;
 
 	  locpoints.SetSize (plainpoints.Size());
+          /*
 	  for (int i = oldnp+1; i <= plainpoints.Size(); i++)
 	    trans.FromPlain (plainpoints.Elem(i), locpoints.Elem(i));
-	  
+          */
+	  for (auto i : plainpoints.Range().Modify(oldnp,0))
+	    trans.FromPlain (plainpoints[i], locpoints[i]);
+          
 
 
 	  // avoid meshing from large to small mesh-size
@@ -565,14 +570,15 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      if (testmode)
 		{
 		  (*testout) << "found is active, 3" << endl;
-		  for(int i = 1; i <= plainpoints.Size(); i++)
+		  //for(int i = 1; i <= plainpoints.Size(); i++)
+                  for(auto i : plainpoints.Range())
 		    {
 		      (*testout) << "p";
-		      if (i <= pindex.Size())
-			(*testout) << pindex.Get(i) << ": ";
+		      if (i < pindex.Range().Next())
+			(*testout) << pindex[i] << ": ";
 		      else
 			(*testout) << "new: ";
-		      (*testout) << plainpoints.Get(i) << endl;
+		      (*testout) << plainpoints[i] << endl;
 		    }
 		}
 	      
@@ -582,20 +588,24 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	      minerr = err;
 	      
 	      tempnewpoints.SetSize (0);
-	      for(int i = oldnp+1; i <= locpoints.Size(); i++)
-		tempnewpoints.Append (locpoints.Get(i));
+	      // for(int i = oldnp+1; i <= locpoints.Size(); i++)
+              for (auto i : locpoints.Range().Modify(oldnp,0))
+		tempnewpoints.Append (locpoints[i]);
 	      
 	      tempnewfaces.SetSize (0);
-	      for(int i = oldnf+1; i <= locfaces.Size(); i++)
-		tempnewfaces.Append (locfaces.Get(i));
+	      // for(int i = oldnf+1; i <= locfaces.Size(); i++)
+              for (auto i : locfaces.Range().Modify(oldnf,0))              
+		tempnewfaces.Append (locfaces[i]);
 	      
 	      tempdelfaces.SetSize (0);
-	      for(int i = 1; i <= delfaces.Size(); i++)
-		tempdelfaces.Append (delfaces.Get(i));
+	      // for(int i = 1; i <= delfaces.Size(); i++)
+              for (auto i : delfaces.Range())
+		tempdelfaces.Append (delfaces[i]);
 	      
 	      templocelements.SetSize (0);
-	      for(int i = 1; i <= locelements.Size(); i++)
-		templocelements.Append (locelements.Get(i));
+	      // for(int i = 1; i <= locelements.Size(); i++)
+              for (auto i : locelements.Range())
+		templocelements.Append (locelements[i]);
 
 	      /*
 	      optother =
@@ -637,14 +647,14 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 	  if (loktestmode)
 	    {
 	      (*testout) << "apply rule" << endl;
-	      for(int i = 1; i <= locpoints.Size(); i++)
+	      for (auto i : locpoints.Range())
 		{
 		  (*testout) << "p";
-		  if (i <= pindex.Size())
-		    (*testout) << pindex.Get(i) << ": ";
+		  if (pindex.Range().Contains(i))
+		    (*testout) << pindex[i] << ": ";
 		  else
 		    (*testout) << "new: ";
-		  (*testout) << locpoints.Get(i) << endl;
+		  (*testout) << locpoints[i] << endl;
 		}
 	    }
 
@@ -652,10 +662,11 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp)
 
 	  pindex.SetSize(locpoints.Size());
 
-	  for (int i = oldnp+1; i <= locpoints.Size(); i++)
+	  // for (int i = oldnp+1; i <= locpoints.Size(); i++)
+          for (auto i : locpoints.Range().Modify(oldnp,0))
 	    {
-	      PointIndex globind = mesh.AddPoint (locpoints.Get(i));
-	      pindex.Elem(i) = adfront -> AddPoint (locpoints.Get(i), globind);
+	      PointIndex globind = mesh.AddPoint (locpoints[i]);
+	      pindex[i] = adfront -> AddPoint (locpoints[i], globind);
 	    }
 
 	  for (int i = 1; i <= locelements.Size(); i++)
diff --git a/libsrc/meshing/meshing3.hpp b/libsrc/meshing/meshing3.hpp
index cc28a276..85aa4a3f 100644
--- a/libsrc/meshing/meshing3.hpp
+++ b/libsrc/meshing/meshing3.hpp
@@ -45,8 +45,8 @@ public:
   MESHING3_RESULT GenerateMesh (Mesh & mesh, const MeshingParameters & mp);
   
   ///
-  int ApplyRules (NgArray<Point3d, PointIndex::BASE> & lpoints,
-                  NgArray<int, PointIndex::BASE> & allowpoint,
+  int ApplyRules (Array<Point3d, PointIndex> & lpoints,
+                  Array<int, PointIndex> & allowpoint,
 		  NgArray<MiniElement2d> & lfaces, INDEX lfacesplit,
 		  INDEX_2_HASHTABLE<int> & connectedpairs,
 		  NgArray<Element> & elements,
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index bb09d193..cb750166 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -151,7 +151,7 @@ namespace netgen
 
 
 
-
+  /*
   class PointIndex
   {
     int i;
@@ -191,8 +191,87 @@ namespace netgen
 #endif  
 
     void DoArchive (Archive & ar) { ar & i; }
-  };
+  }
+    */
 
+  
+  class PointIndex
+  {
+    int i;
+  public:
+    class t_invalid { public: constexpr t_invalid() = default; };
+    static constexpr t_invalid INVALID{};
+    
+    PointIndex () = default;
+    PointIndex (const PointIndex&) = default;
+    PointIndex (PointIndex &&) = default;
+    PointIndex & operator= (const PointIndex&) = default;
+    PointIndex & operator= (PointIndex&&) = default;
+
+    // private:
+    constexpr PointIndex (int ai) : i(ai)
+    {
+#ifdef DEBUG
+      if (ai < PointIndex::BASE)
+        cout << "illegal PointIndex, use PointIndex::INVALID instead" << endl;
+        // throw Exception("illegal PointIndex, use PointIndex::INVALID instead");
+#endif
+    }
+    friend constexpr netgen::PointIndex ngcore::IndexBASE<netgen::PointIndex> ();
+    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);    
+    friend PointIndex operator+ (int, PointIndex);
+    friend PointIndex operator+ (size_t, PointIndex);
+    friend PointIndex operator- (PointIndex, int);
+    friend int operator- (PointIndex, PointIndex);
+    friend bool operator< (PointIndex a, PointIndex b);
+    friend bool operator> (PointIndex a, PointIndex b);
+    friend bool operator>= (PointIndex a, PointIndex b);
+    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; }
+    // private:
+    constexpr operator const int& () const { return i; }
+    explicit constexpr operator int& () { return i; }
+  public:
+    PointIndex operator++ (int) { PointIndex hi(*this); i++; return hi; }
+    PointIndex operator-- (int) { PointIndex hi(*this); i--; return hi; }
+    PointIndex & operator++ () { i++; return *this; }
+    PointIndex operator-- () { i--; return *this; }
+    PointIndex operator+= (int add) { i += add; return *this; }
+    void Invalidate() { i = PointIndex::BASE-1; }
+    bool IsValid() const { return i != PointIndex::BASE-1; }
+#ifdef BASE0
+    static constexpr size_t BASE = 0;
+#else
+    static constexpr size_t BASE = 1;
+#endif  
+
+    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); }  
+  inline PointIndex operator+ (int i, PointIndex pi) { return PointIndex(pi.i+i); }
+  inline PointIndex operator+ (size_t i, PointIndex pi) { return PointIndex(pi.i+i); }
+  inline PointIndex operator- (PointIndex pi, int i) { return PointIndex(pi.i-i); }
+  inline int operator- (PointIndex pa, PointIndex pb) { return PointIndex(pa.i-pb.i); }
+  inline bool operator< (PointIndex a, PointIndex b) { return a.i < b.i; }
+  inline bool operator> (PointIndex a, PointIndex b) { return a.i > b.i; }
+  inline bool operator>= (PointIndex a, PointIndex b) { return a.i >= b.i; }
+  inline bool operator== (PointIndex a, PointIndex b) { return a.i == b.i; }
+  inline bool operator!= (PointIndex a, PointIndex b) { return a.i != b.i; }
+  */
 }
 
 namespace ngcore
@@ -224,6 +303,7 @@ namespace netgen
     PointIndices (PointIndex i1, PointIndex i2) : INDEX_2(i1,i2) { ; } 
     PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); }
     PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_2::operator[](i)); }
+    using INDEX_2::Sort;
     static PointIndices Sort(PointIndex i1, PointIndex i2) { return INDEX_2::Sort(i1, i2); }
     template <size_t J>
     PointIndex get() const { return PointIndex(INDEX_2::operator[](J)); }    
@@ -1627,14 +1707,15 @@ namespace netgen
 
     bool Used (PointIndex pi1, PointIndex pi2)
     {
-      return identifiedpoints.Used (INDEX_2 (pi1, pi2));
+      // return identifiedpoints.Used (INDEX_2 (pi1, pi2));
+      return identifiedpoints.Used (PointIndices<2>(pi1, pi2));
     }
 
     bool UsedSymmetric (PointIndex pi1, PointIndex pi2)
     {
       return 
-	identifiedpoints.Used (INDEX_2 (pi1, pi2)) ||
-	identifiedpoints.Used (INDEX_2 (pi2, pi1));
+	identifiedpoints.Used (PointIndices<2>(pi1, pi2)) ||
+	identifiedpoints.Used (PointIndices<2>(pi2, pi1));
     }
 
     ///
diff --git a/libsrc/meshing/ruler3.cpp b/libsrc/meshing/ruler3.cpp
index adb1fe34..188e6a25 100644
--- a/libsrc/meshing/ruler3.cpp
+++ b/libsrc/meshing/ruler3.cpp
@@ -8,7 +8,7 @@ extern double minother;
 extern double minwithoutother;
 
 
-  static double CalcElementBadness (const NgArray<Point3d, PointIndex::BASE> & points,
+  static double CalcElementBadness (const Array<Point3d, PointIndex> & points,
                                     const Element & elem)
 {
   double vol, l, l4, l5, l6;
@@ -49,8 +49,8 @@ extern double minwithoutother;
 
 int Meshing3 :: ApplyRules 
 (
- NgArray<Point3d, PointIndex::BASE> & lpoints,     // in: local points, out: old+new local points
- NgArray<int, PointIndex::BASE> & allowpoint,      // in: 2 .. it is allowed to use pointi, 1..will be allowed later, 0..no means
+ 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
  INDEX lfacesplit,	       // for local faces in outer radius
  INDEX_2_HASHTABLE<int> & connectedpairs,  // connected pairs for prism-meshing
@@ -659,11 +659,12 @@ int Meshing3 :: ApplyRules
 
 		      // check freezone:
 		      
-		      for (int i = 1; i <= lpoints.Size(); i++)
+		      // for (int i = 1; i <= lpoints.Size(); i++)
+                      for (auto i : lpoints.Range())
 			{
 			  if ( !pused.Get(i) )
 			    {
-			      const Point3d & lp = lpoints.Get(i);
+			      const Point3d & lp = lpoints[i];
 
 			      if (rule->fzbox.IsIn (lp))
 				{
@@ -674,7 +675,7 @@ int Meshing3 :: ApplyRules
 					  (*testout) << "Point " << i 
 						     << " in Freezone" << endl;
 					  snprintf (problems.Elem(ri), 255,
-						   "locpoint %d in Freezone", i);
+                                                    "locpoint %d in Freezone", int(i));
 					}
 				      ok = 0;
 				      break;
@@ -933,11 +934,11 @@ int Meshing3 :: ApplyRules
 
 			  // new points in free-zone ?
 			  for (int i = rule->GetNOldP() + 1; i <= rule->GetNP() && ok; i++)
-			    if (!rule->IsInFreeZone (lpoints.Get(pmap.Get(i))))
+			    if (!rule->IsInFreeZone (lpoints[pmap.Get(i)]))
 			      {
 				if (loktestmode)
 				  {
-				    (*testout) << "Newpoint " << lpoints.Get(pmap.Get(i))
+				    (*testout) << "Newpoint " << lpoints[pmap.Get(i)]
 					       << " outside convex hull" << endl;
 				    snprintf (problems.Elem(ri), 255, "newpoint outside convex hull");
 				  }
@@ -958,9 +959,9 @@ int Meshing3 :: ApplyRules
 			  // Calculate Element badness
 			  
 			  teterr = 0;
-			  for (int i = 1; i <= elements.Size(); i++)
+			  for (auto i : elements.Range())
 			    {
-			      double hf = CalcElementBadness (lpoints, elements.Get(i));
+			      double hf = CalcElementBadness (lpoints, elements[i]);
 			      if (hf > teterr) teterr = hf;
 			    }
 
@@ -1066,25 +1067,29 @@ int Meshing3 :: ApplyRules
 				    {
 				      (*testout) << "P" << i << ": Ref: "
 						 << rule->GetPoint (i) << "  is: "
-						 << lpoints.Get(pmap.Get(i)) << endl;
+						 << lpoints[pmap.Get(i)] << endl;
 				    }
 				}
 			      
 			      tempnewpoints.SetSize (0);
-			      for (int i = noldlp+1; i <= lpoints.Size(); i++)
-				tempnewpoints.Append (lpoints.Get(i));
+			      // for (int i = noldlp+1; i <= lpoints.Size(); i++)
+                              for (auto i : lpoints.Range().Modify(noldlp, 0))
+				tempnewpoints.Append (lpoints[i]);
 			      
 			      tempnewfaces.SetSize (0);
-			      for (int i = noldlf+1; i <= lfaces.Size(); i++)
-				tempnewfaces.Append (lfaces.Get(i));
+			      // for (int i = noldlf+1; i <= lfaces.Size(); i++)
+                              for (auto i : lfaces.Range().Modify(noldlf,0))
+				tempnewfaces.Append (lfaces[i]);
 
 			      tempdelfaces.SetSize (0);
-			      for (int i = 1; i <= delfaces.Size(); i++)
-				tempdelfaces.Append (delfaces.Get(i));
+			      // for (int i = 1; i <= delfaces.Size(); i++)
+                              for (auto i : delfaces.Range())
+				tempdelfaces.Append (delfaces[i]);
 			      
 			      tempelements.SetSize (0);
-			      for (int i = 1; i <= elements.Size(); i++)
-				tempelements.Append (elements.Get(i));
+			      // for (int i = 1; i <= elements.Size(); i++)
+                              for (auto i : elements.Range())
+				tempelements.Append (elements[i]);
 			    }