From bcbd390f7d632d582b09c43d6c571de98f4a72a4 Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Tue, 31 Dec 2024 21:26:04 +0100
Subject: [PATCH] PointIndex in Identifications

---
 libsrc/csg/zrefine.cpp                 | 23 ++++++--
 libsrc/include/nginterface_v2_impl.hpp |  4 +-
 libsrc/interface/nginterface_v2.cpp    | 10 +++-
 libsrc/interface/writeuser.cpp         | 39 ++++++++++----
 libsrc/meshing/meshclass.cpp           | 34 ++++++------
 libsrc/meshing/meshclass.hpp           |  2 +-
 libsrc/meshing/meshfunc.cpp            |  6 ++-
 libsrc/meshing/meshtype.cpp            | 75 ++++++++++++++++++++------
 libsrc/meshing/meshtype.hpp            | 12 +++--
 libsrc/meshing/python_mesh.cpp         |  5 +-
 libsrc/meshing/secondorder.cpp         |  1 +
 libsrc/meshing/topology.cpp            | 15 ++++--
 libsrc/visualization/vsmesh.cpp        | 45 +++++++++-------
 13 files changed, 188 insertions(+), 83 deletions(-)

diff --git a/libsrc/csg/zrefine.cpp b/libsrc/csg/zrefine.cpp
index 649bb50c..6e154bcb 100644
--- a/libsrc/csg/zrefine.cpp
+++ b/libsrc/csg/zrefine.cpp
@@ -264,14 +264,21 @@ namespace netgen
       {
         auto & identpts =
           mesh.GetIdentifications().GetIdentifiedPoints ();
-        
+
+        /*
 	for (int i = 1; i <= identpts.GetNBags(); i++)
 	  for (int j = 1; j <= identpts.GetBagSize(i); j++)
 	    {
 	      INDEX_3 pair;
 	      int dummy;
 	      identpts.GetData(i, j, pair, dummy);
-              auto idnr = pair[2];
+        */
+        for (auto [hash, val] : identpts)\
+          {
+            auto [hash_pts, idnr] = hash;
+            auto [pi1, pi2] = hash_pts;
+            // auto idnr = pair[2];
+            
 	      const CloseSurfaceIdentification * csid = 
 		dynamic_cast<const CloseSurfaceIdentification*> 
 		(geom->identifications.Get(idnr));
@@ -282,17 +289,25 @@ namespace netgen
 		      if (first_id.Test (idnr))
 			{
 			  first_id.Clear(idnr);
+                          /*
 			  ref_uniform.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels()));
 			  ref_singular.Append (INDEX_3 (pair.I1(), pair.I2(), csid->RefLevels1()));
 			  ref_singular.Append (INDEX_3 (pair.I2(), pair.I1(), csid->RefLevels2()));
+                          */
+			  ref_uniform.Append (INDEX_3 (pi1, pi2, csid->RefLevels()));
+			  ref_singular.Append (INDEX_3 (pi1, pi2, csid->RefLevels1()));
+			  ref_singular.Append (INDEX_3 (pi2, pi1, csid->RefLevels2()));
+                          
 			}
 		    }
 		  else
 		    {   
 		      //const NgArray<double> & slices = csid->GetSlices();
 		      INDEX_4 i4;
-		      i4[0] = pair.I1();
-		      i4[1] = pair.I2();
+		      // i4[0] = pair.I1();
+		      // i4[1] = pair.I2();
+		      i4[0] = pi1; 
+		      i4[1] = pi2; 
 		      i4[2] = idnr;
 		      i4[3] = csid->GetSlices().Size();
 		      ref_slices.Append (i4);
diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp
index 19c41257..d0db9a2b 100644
--- a/libsrc/include/nginterface_v2_impl.hpp
+++ b/libsrc/include/nginterface_v2_impl.hpp
@@ -339,8 +339,8 @@ NGX_INLINE DLL_HEADER Ng_Buffer<int[2]> Ngx_Mesh :: GetPeriodicVertices(int idnr
   mesh->GetIdentifications().GetPairs (idnr+1, apairs);
   for(auto& ind : apairs)
     {
-      ind.I1()--;
-      ind.I2()--;
+      ind.I1() -= IndexBASE<PointIndex>();
+      ind.I2() -= IndexBASE<PointIndex>();
     }
   typedef int ti2[2];
   return { apairs.Size(), (ti2*)(void*)apairs.Release() };
diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp
index 9233c5a7..1e3e7af8 100644
--- a/libsrc/interface/nginterface_v2.cpp
+++ b/libsrc/interface/nginterface_v2.cpp
@@ -1249,9 +1249,15 @@ int Ngx_Mesh::GetElementOrder (int enr) const
 void Ngx_Mesh::GetElementOrders (int enr, int * ox, int * oy, int * oz) const
 {
   if (mesh->GetDimension() == 3)
-    mesh->VolumeElement(enr).GetOrder(*ox, *oy, *oz);
+    {
+      ElementIndex ei = IndexBASE<ElementIndex>() + enr-1;
+      mesh->VolumeElement(ei).GetOrder(*ox, *oy, *oz);
+    }
   else
-    mesh->SurfaceElement(enr).GetOrder(*ox, *oy, *oz);
+    {
+      SurfaceElementIndex sei = IndexBASE<SurfaceElementIndex>() + enr-1;      
+      mesh->SurfaceElement(sei).GetOrder(*ox, *oy, *oz);
+    }
 }
 
 void Ngx_Mesh::SetElementOrder (int enr, int order)
diff --git a/libsrc/interface/writeuser.cpp b/libsrc/interface/writeuser.cpp
index 27370eb0..b4f5f033 100644
--- a/libsrc/interface/writeuser.cpp
+++ b/libsrc/interface/writeuser.cpp
@@ -94,9 +94,13 @@ void WriteNeutralFormat (const Mesh & mesh,
   if (mesh.GetDimension() == 3)
     {
       outfile << ne << "\n";
+      /*
       for (int i = 1; i <= ne; i++)
 	{
 	  Element el = mesh.VolumeElement(i);
+      */
+      for (Element el : mesh.VolumeElements())
+        {
 	  if (inverttets)
 	    el.Invert();
 	  outfile.width(4);
@@ -112,9 +116,13 @@ void WriteNeutralFormat (const Mesh & mesh,
     }
 
   outfile << nse << "\n";
+  /*
   for (int i = 1; i <= nse; i++)
     {
       Element2d el = mesh.SurfaceElement(i);
+  */
+  for (Element2d el : mesh.SurfaceElements())
+    {
       if (invertsurf)
 	el.Invert();
       outfile.width(4);
@@ -583,9 +591,13 @@ void WriteFEPPFormat (const Mesh & mesh,
 
 
       outfile << ne << "\n";
+      /*
       for (i = 1; i <= ne; i++)
 	{
 	  const Element & el = mesh.VolumeElement(i);
+      */
+      for (const Element & el : mesh.VolumeElements())
+        {
 	  outfile.width(4);
 	  outfile << el.GetIndex() << " ";
 	  outfile.width(4);
@@ -677,7 +689,6 @@ void WriteEdgeElementFormat (const Mesh & mesh,
   int nelements = mesh.GetNE();
   int nsurfelem = mesh.GetNSE();
   int nedges = top->GetNEdges();
-  int i, j;
 
   int inverttets = mparam.inverttets;
   int invertsurf = mparam.inverttrigs;
@@ -692,7 +703,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
 
   // vertices with coordinates
   outfile << npoints << "\n";
-  for (i = 1; i <= npoints; i++)
+  for (int i = 1; i <= npoints; i++)
     {
       const Point3d & p = mesh.Point(i);
 
@@ -706,16 +717,24 @@ void WriteEdgeElementFormat (const Mesh & mesh,
 
   // element - edge - list
   outfile << nelements << " " << nedges << "\n";
+  /*
   for (i = 1; i <= nelements; i++)
     {
       Element el = mesh.VolumeElement(i);
+  */
+  for (ElementIndex ei : Range(mesh.VolumeElements()))
+    {
+      int i = ei-IndexBASE(ei)+1;
+      
+      Element el = mesh.VolumeElement(ei);
+
       if (inverttets)
       	el.Invert();
       outfile.width(4);
       outfile << el.GetIndex() << "  ";
       outfile.width(8);
       outfile << el.GetNP();
-      for (j = 1; j <= el.GetNP(); j++)
+      for (int j = 1; j <= el.GetNP(); j++)
 	{
 	  outfile << " ";
 	  outfile.width(8);
@@ -723,11 +742,11 @@ void WriteEdgeElementFormat (const Mesh & mesh,
 	}
 
       // top->GetElementEdges(i,edges);
-      auto eledges = top->GetEdges(ElementIndex(i-1));
+      auto eledges = top->GetEdges(ei);
       outfile << endl << "      ";
       outfile.width(8);
       outfile << eledges.Size();
-      for (j=1; j <= eledges.Size(); j++)
+      for (int j=1; j <= eledges.Size(); j++)
 	{
 	  outfile << " ";
 	  outfile.width(8);
@@ -738,7 +757,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
       // orientation:
       top->GetElementEdgeOrientations(i,edges);
       outfile << "              ";
-      for (j=1; j <= edges.Size(); j++)
+      for (int j=1; j <= edges.Size(); j++)
 	{
 	  outfile << " ";
 	  outfile.width(8);
@@ -749,7 +768,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
 
   // surface element - edge - list (with boundary conditions)
   outfile << nsurfelem << "\n";
-  for (i = 1; i <= nsurfelem; i++)
+  for (int i = 1; i <= nsurfelem; i++)
     {
       Element2d el = mesh.SurfaceElement(i);
       if (invertsurf)
@@ -758,7 +777,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
       outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << "  ";
       outfile.width(8);
       outfile << el.GetNP();
-      for (j = 1; j <= el.GetNP(); j++)
+      for (int j = 1; j <= el.GetNP(); j++)
 	{
 	  outfile << " ";
 	  outfile.width(8);
@@ -769,7 +788,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
       outfile << endl << "      ";
       outfile.width(8);
       outfile << edges.Size();
-      for (j=1; j <= edges.Size(); j++)
+      for (int j=1; j <= edges.Size(); j++)
 	{
 	  outfile << " ";
 	  outfile.width(8);
@@ -782,7 +801,7 @@ void WriteEdgeElementFormat (const Mesh & mesh,
   // int v1, v2;
   // edge - vertex - list
   outfile << nedges << "\n";
-  for (i=1; i <= nedges; i++)
+  for (int i=1; i <= nedges; i++)
     {
       // top->GetEdgeVertices(i,v1,v2);
       auto [v1,v2] = top->GetEdgeVertices(i-1);
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 16ee4a0a..4dae9b84 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -6438,22 +6438,24 @@ namespace netgen
     mapped_points = false;
 
     // Add new points
-    for(auto [p1p2, dummy] : identpts)
+    for(auto [hash, dummy] : identpts)
       {
-        if(p1p2[2] != nr)
+        auto [hash_pts, hash_nr] = hash;
+        if(hash_nr != nr)
           continue;
-        auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}];
-        auto p1 = Point(p1p2.I1());
-        auto p2 = Point(p1p2.I2());
-        ipts.Append(p1p2.I1());
-        mapped_points.SetBit(p1p2.I1());
+        // auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}];
+        auto& ipts = inserted_points[ { hash_pts[0], hash_pts[1] }];
+        auto p1 = Point(hash_pts.I1());
+        auto p2 = Point(hash_pts.I2());
+        ipts.Append(hash_pts.I1());
+        mapped_points.SetBit(hash_pts.I1());
         for(auto slice : slices)
           {
             auto np = p1 + slice * (p2-p1);
             auto npi = AddPoint(np);
             ipts.Append(npi);
           }
-        ipts.Append(p1p2.I2());
+        ipts.Append(hash_pts.I2());
       }
 
     // Split segments
@@ -6644,21 +6646,21 @@ namespace netgen
 
   void Mesh :: CalcMinMaxAngle (double badellimit, double * retvalues) 
   {
-    int i, j;
     int lpi1, lpi2, lpi3, lpi4;
     double phimax = 0, phimin = 10;
     double facephimax = 0, facephimin = 10;
     int illegaltets = 0, negativetets = 0, badtets = 0;
 
-    for (i = 1; i <= GetNE(); i++)
+    // for (int i = 1; i <= GetNE(); i++)
+    for (ElementIndex ei : Range(VolumeElements()))
       {
         int badel = 0;
 
-        Element & el = VolumeElement(i);
+        Element & el = VolumeElement(ei);
 
         if (el.GetType() != TET)
           {
-            VolumeElement(i).Flags().badel = 0;
+            VolumeElement(ei).Flags().badel = 0;
             continue;
           }
 
@@ -6673,8 +6675,8 @@ namespace netgen
           {
             badel = 1;
             illegaltets++;
-            (*testout) << "illegal tet: " << i << " ";
-            for (j = 1; j <= el.GetNP(); j++)
+            (*testout) << "illegal tet: " << ei << " ";
+            for (int j = 1; j <= el.GetNP(); j++)
               (*testout) << el.PNum(j) << " ";
             (*testout) << endl;
           }
@@ -6713,7 +6715,7 @@ namespace netgen
 
 
         // angles in faces
-        for (j = 1; j <= 4; j++)
+        for (int j = 1; j <= 4; j++)
           {
             Element2d face(TRIG);
             el.GetFace (j, face);
@@ -6740,7 +6742,7 @@ namespace netgen
           }
 
 
-        VolumeElement(i).Flags().badel = badel;
+        VolumeElement(ei).Flags().badel = badel;
         if (badel) badtets++;
       }
 
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 90619969..f89dddc1 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -377,7 +377,7 @@ namespace netgen
     const auto & SurfaceElements() const { return surfelements; }
     auto & SurfaceElements() { return surfelements; }
 
-  
+    
     DLL_HEADER void RebuildSurfaceElementLists ();
     DLL_HEADER void GetSurfaceElementsOfFace (int facenr, Array<SurfaceElementIndex> & sei) const;
 
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index 44040fe3..5716f928 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -578,8 +578,10 @@ namespace netgen
              el.SetIndex(m_.domain);
              mesh.AddVolumeElement(el);
          }
-         for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints())
-           mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]);
+         // for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints())
+         // mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]);
+         for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints())         
+           mesh.GetIdentifications().Add( pmap[ get<0>(p1p2)[0] ], pmap[ get<0>(p1p2)[1]] , get<1>(p1p2) );
          for(auto i : Range(m.GetIdentifications().GetMaxNr()))
            {
              mesh.GetIdentifications().SetType(i+1, m.GetIdentifications().GetType(i+1));
diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp
index 24e9f398..2e954f54 100644
--- a/libsrc/meshing/meshtype.cpp
+++ b/libsrc/meshing/meshtype.cpp
@@ -2708,8 +2708,9 @@ namespace netgen
     void Identifications :: DoArchive (Archive & ar)
     {
       ar & maxidentnr;
-      ar & identifiedpoints & identifiedpoints_nr;
-
+      // ar & identifiedpoints & identifiedpoints_nr;
+#pragma message( "Archive CloseHadhTable missing " __FILE__ )
+        
       ar & idpoints_table;
       if (ar.Output())
         {
@@ -2735,11 +2736,11 @@ namespace netgen
   void Identifications :: Add (PointIndex pi1, PointIndex pi2, int identnr)
   {
     //  (*testout) << "Identification::Add, pi1 = " << pi1 << ", pi2 = " << pi2 << ", identnr = " << identnr << endl;
-    INDEX_2 pair (pi1, pi2);
+    PointIndices<2> pair (pi1, pi2);
     identifiedpoints.Set (pair, identnr);
 
-    INDEX_3 tripl (pi1, pi2, identnr);
-    identifiedpoints_nr.Set (tripl, 1);
+    // INDEX_3 tripl (pi1, pi2, identnr);
+    identifiedpoints_nr.Set ( { { pi1, pi2 }, identnr }, 1);
 
     if (identnr > maxidentnr) maxidentnr = identnr;
     names.SetSize(maxidentnr);
@@ -2762,8 +2763,9 @@ namespace netgen
 
   bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const
   {
-    INDEX_3 tripl(pi1, pi2, nr);
-    if (identifiedpoints_nr.Used (tripl))
+    // INDEX_3 tripl(pi1, pi2, nr);
+    // if (identifiedpoints_nr.Used (tripl))
+    if (identifiedpoints_nr.Used ( { { pi1, pi1 }, nr } ) )
       return 1;
     else
       return 0;
@@ -2803,23 +2805,30 @@ namespace netgen
       {
         cout << "getmap, identnr = " << identnr << endl;
 
+        /*
         for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++)
           for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++)
+        */
+        for (auto [hash, val] : identifiedpoints_nr)
             {
+              /*
               INDEX_3 i3;
               int dummy;
               identifiedpoints_nr.GetData (i, j, i3, dummy);
-	    
-              if (i3.I3() == identnr || !identnr)
+              */
+
+              auto [hash_pts, hash_nr] = hash;
+              
+              if (hash_nr == identnr || !identnr)
                 {
                   /*
                   identmap.Elem(i3.I1()) = i3.I2();
                   if(symmetric)
                     identmap.Elem(i3.I2()) = i3.I1();
                   */
-                  identmap[i3.I1()] = i3.I2();
+                  identmap[hash_pts.I1()] = hash_pts.I2();
                   if(symmetric)
-                    identmap[i3.I2()] = i3.I1();
+                    identmap[hash_pts.I2()] = hash_pts.I1();
                 }
             }  
       }
@@ -2830,8 +2839,12 @@ namespace netgen
   Array<INDEX_3> Identifications :: GetPairs () const
   {
     Array<INDEX_3> pairs;
-    for(auto [i3, dummy] : identifiedpoints_nr)
-      pairs.Append(i3);
+    for(auto [hash, dummy] : identifiedpoints_nr)
+      // pairs.Append(i3);
+      {
+        auto [pts,nr] = hash;
+        pairs.Append ( { pts[0], pts[1], nr } );
+      }
     return pairs;
   }
 
@@ -2841,6 +2854,8 @@ namespace netgen
     identpairs.SetSize(0);
   
     if (identnr == 0)
+      {
+        /*
       for (int i = 1; i <= identifiedpoints.GetNBags(); i++)
         for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++)
           {
@@ -2848,8 +2863,14 @@ namespace netgen
             int nr;
             identifiedpoints.GetData (i, j, i2, nr);
             identpairs.Append (i2);
-          }  
+          }
+        */
+        for (auto [hash,val] : identifiedpoints)
+          identpairs.Append (hash);
+      }
     else
+      {
+        /*
       for (int i = 1; i <= identifiedpoints_nr.GetNBags(); i++)
         for (int j = 1; j <= identifiedpoints_nr.GetBagSize(i); j++)
           {
@@ -2859,12 +2880,21 @@ namespace netgen
 	  
             if (i3.I3() == identnr)
               identpairs.Append (INDEX_2(i3.I1(), i3.I2()));
-          }  
+          }
+        */
+        for (auto [hash,val] : identifiedpoints_nr)
+          {
+            auto [hash_pts, hash_nr] = hash;
+            if (hash_nr == identnr)
+              identpairs.Append (hash_pts);
+          }
+      }
   }
 
 
   void Identifications :: SetMaxPointNr (int maxpnum)
   {
+    /*
     for (int i = 1; i <= identifiedpoints.GetNBags(); i++)
       for (int j = 1; j <= identifiedpoints.GetBagSize(i); j++)
         {
@@ -2878,6 +2908,18 @@ namespace netgen
               identifiedpoints.SetData (i, j, i2, -1);	    
             }
         }
+    */
+
+    // can we get data by reference ? 
+    for (auto [hash,data] : identifiedpoints)
+      {
+        if (hash.I1() > IndexBASE<PointIndex>()+maxpnum-1 ||
+            hash.I2() > IndexBASE<PointIndex>()+maxpnum-1)
+          {
+            identifiedpoints[hash] = -1;
+          }
+            
+      }
   }
 
   // Map points in the identifications to new point numbers
@@ -2900,7 +2942,8 @@ namespace netgen
   {
     ost << "Identifications:" << endl;
     ost << "pairs: " << endl << identifiedpoints << endl;
-    ost << "pairs and nr: " << endl << identifiedpoints_nr << endl;
+    // ost << "pairs and nr: " << endl << identifiedpoints_nr << endl;
+#pragma message( "Can't ostream a tuple " __FILE__ )    
     ost << "table: " << endl << idpoints_table << endl;
   }
 
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 0f2e780d..91c5123a 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -244,7 +244,7 @@ namespace netgen
   
 
 
-  class PointIndex : public Index<int,PointIndex,1>
+  class PointIndex : public Index<int,PointIndex,0>
   {
   public:
     using Index::Index;
@@ -1763,13 +1763,15 @@ namespace netgen
     class Mesh & mesh;
 
     /// identify points (thin layers, periodic b.c.)  
-    INDEX_2_HASHTABLE<int> identifiedpoints;
+    // INDEX_2_HASHTABLE<int> identifiedpoints;
+    ClosedHashTable<PointIndices<2>, int> identifiedpoints;
   
     /// the same, with info about the id-nr
-    INDEX_3_HASHTABLE<int> identifiedpoints_nr;
+    // INDEX_3_HASHTABLE<int> identifiedpoints_nr;
+    ClosedHashTable<std::tuple<PointIndices<2>, int>, int> identifiedpoints_nr;
 
     /// sorted by identification nr
-    TABLE<INDEX_2> idpoints_table;
+    TABLE<PointIndices<2>> idpoints_table;
 
     NgArray<ID_TYPE> type;
 
@@ -1805,7 +1807,7 @@ namespace netgen
 
     // bool HasIdentifiedPoints() const { return identifiedpoints != nullptr; } 
     ///
-    INDEX_3_HASHTABLE<int> & GetIdentifiedPoints ()
+    auto & GetIdentifiedPoints ()
     { 
       return identifiedpoints_nr;
     }
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 1a133165..c9c95002 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -1266,7 +1266,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
            py::list points;
            for(const auto& pair : self.GetIdentifications().GetIdentifiedPoints())
              {
-               py::tuple pnts = py::make_tuple(pair.first.I1(), pair.first.I2());
+               // py::tuple pnts = py::make_tuple(pair.first.I1(), pair.first.I2());
+               
+               auto [pi1, pi2] = get<0> (pair.first);
+               py::tuple pnts = py::make_tuple(pi1, pi2);
                points.append(pnts);
              }
            return points;
diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp
index 5527b2a8..56dc7583 100644
--- a/libsrc/meshing/secondorder.cpp
+++ b/libsrc/meshing/secondorder.cpp
@@ -558,6 +558,7 @@ namespace netgen
 		//	      (*testout) << "bad els: " << endl;
 		wrongels = 0;
 		for (int i = 1; i <= ne; i++)
+                  
 		  {
 		    if (!illegalels.Test(i) && 
 			mesh.VolumeElement(i).
diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp
index ce4f9ce7..2af9c76b 100644
--- a/libsrc/meshing/topology.cpp
+++ b/libsrc/meshing/topology.cpp
@@ -1829,8 +1829,8 @@ namespace netgen
   
   void MeshTopology :: GetElementFaces (int elnr, NgArray<int> & elfaces, bool withorientation) const
   {
-    int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
     ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);
+    int nfa = GetNFaces (mesh->VolumeElement(ei).GetType());
     
     elfaces.SetSize (nfa);
 
@@ -1859,7 +1859,8 @@ namespace netgen
 
   void MeshTopology :: GetElementEdgeOrientations (int elnr, NgArray<int> & eorient) const
   {
-    int ned = GetNEdges (mesh->VolumeElement(elnr).GetType());
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);    
+    int ned = GetNEdges (mesh->VolumeElement(ei).GetType());
     eorient.SetSize (ned);
     for (int i = 1; i <= ned; i++)
       // eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1;
@@ -1869,7 +1870,8 @@ namespace netgen
 
   void MeshTopology :: GetElementFaceOrientations (int elnr, NgArray<int> & forient) const
   {
-    int nfa = GetNFaces (mesh->VolumeElement(elnr).GetType());
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);    
+    int nfa = GetNFaces (mesh->VolumeElement(ei).GetType());
     forient.SetSize (nfa);
     for (int i = 1; i <= nfa; i++)
       // forient.Elem(i) = faces.Get(elnr)[i-1].forient;
@@ -2111,7 +2113,9 @@ namespace netgen
 
   int MeshTopology :: GetElementEdgeOrientation (int elnr, int locedgenr) const
   {
-    const Element & el = mesh->VolumeElement (elnr);
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);        
+    
+    const Element & el = mesh->VolumeElement (ei);
     const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType());    
 
     int k = locedgenr;
@@ -2123,7 +2127,8 @@ namespace netgen
   
   int MeshTopology :: GetElementFaceOrientation (int elnr, int locfacenr) const
   {
-    const Element & el = mesh->VolumeElement (elnr);
+    ElementIndex ei = IndexBASE<ElementIndex>() +(elnr-1);        
+    const Element & el = mesh->VolumeElement (ei);
     
     const ELEMENT_FACE * elfaces = MeshTopology::GetFaces0 (el.GetType());
 
diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp
index 0e45aa12..f8f6c5d9 100644
--- a/libsrc/visualization/vsmesh.cpp
+++ b/libsrc/visualization/vsmesh.cpp
@@ -515,20 +515,21 @@ namespace netgen
         if (vispar.drawelementnumbers)
 	  {
 	    NgArray<int> v;
-	    for (int i = 1; i <= mesh->GetNE(); i++)
+	    // for (int i = 1; i <= mesh->GetNE(); i++)
+            for (ElementIndex ei : Range(mesh->VolumeElements()))
 	      {
 		// const ELEMENTTYPE & eltype = mesh->ElementType(i);
 		NgArray<int> pnums;
 
 		Point3d p;
-		const Element & el = mesh->VolumeElement (i);
+		const Element & el = mesh->VolumeElement (ei);
 
 		if ( ! el.PNum(5)) //  eltype == TET )
                   {
 
 		    pnums.SetSize(4);
 		    for( int j = 0; j < pnums.Size(); j++)
-		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);
+		      pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
 
 
 		    const Point3d & p1 = mesh->Point(pnums[0]);
@@ -541,7 +542,7 @@ namespace netgen
                   {
 		    pnums.SetSize(5);
 		    for( int j = 0; j < pnums.Size(); j++)
-		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);
+		      pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
 
 
 		    const Point3d & p1 = mesh->Point(pnums[0]);
@@ -559,7 +560,7 @@ namespace netgen
                   {
 		    pnums.SetSize(6);
 		    for( int j = 0; j < pnums.Size(); j++)
-		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);
+		      pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
 
 		    const Point3d & p1 = mesh->Point(pnums[0]);
 		    const Point3d & p2 = mesh->Point(pnums[1]);
@@ -574,7 +575,7 @@ namespace netgen
                   {
 		    pnums.SetSize(8);
 		    for( int j = 0; j < pnums.Size(); j++)
-		      pnums[j] = mesh->VolumeElement(i).PNum(j+1);
+		      pnums[j] = mesh->VolumeElement(ei).PNum(j+1);
 
 		    const Point3d & p1 = mesh->Point(pnums[0]);
 		    const Point3d & p2 = mesh->Point(pnums[1]);
@@ -589,7 +590,7 @@ namespace netgen
                   }
 
 		glRasterPos3d (p.X(), p.Y(), p.Z());
-		snprintf (buf, size(buf),  "%d", i);
+		snprintf (buf, size(buf),  "%d", ei-IndexBASE(ei));
 		// glCallLists (strlen (buf), GL_UNSIGNED_BYTE, buf);
 		MyOpenGLText (buf);
 
@@ -617,14 +618,15 @@ namespace netgen
 	static float badelcol[] = { 1.0f, 0.0f, 1.0f, 1.0f };
 	glLineWidth (1.0f);
 
-	for (int i = 1; i <= mesh->GetNE(); i++)
+	//for (int i = 1; i <= mesh->GetNE(); i++)
+        for (ElementIndex ei : Range(mesh->VolumeElements()))
 	  {
-            if (mesh->VolumeElement(i).Flags().badel ||
-		mesh->VolumeElement(i).Flags().illegal ||
-		(i == vispar.drawelement))
+            if (mesh->VolumeElement(ei).Flags().badel ||
+		mesh->VolumeElement(ei).Flags().illegal ||
+		(ei-IndexBASE(ei) == vispar.drawelement))
 	      {
 		// copy to be thread-safe
-		Element el = mesh->VolumeElement (i);
+		Element el = mesh->VolumeElement (ei);
 		el.GetSurfaceTriangles (faces);
 
 		glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, badelcol);
@@ -674,9 +676,9 @@ namespace netgen
 	  }
 
 
-	for (int i = 1; i <= mesh->GetNE(); i++)
+        for (ElementIndex ei : Range(mesh->VolumeElements()))
 	  {
-            Element el = mesh->VolumeElement (i);
+            Element el = mesh->VolumeElement (ei);
             int hascp = 0;
             for (int j = 1; j <= el.GetNP(); j++)
 	      if (el.PNum(j) == vispar.centerpoint)
@@ -684,7 +686,7 @@ namespace netgen
 
             if (hascp)
 	      {
-		(*testout) << "draw el " << i << " : ";
+		(*testout) << "draw el " << ei << " : ";
 		for (int j = 1; j <= el.GetNP(); j++)
                   (*testout) << el.PNum(j) << " ";
 		(*testout) << endl;
@@ -860,17 +862,22 @@ namespace netgen
 	      {
                 auto & idpts =
                   mesh->GetIdentifications().GetIdentifiedPoints();
-                
+
+                /*
 		for (int i = 1; i <= idpts.GetNBags(); i++)
                   for (int j = 1; j <= idpts.GetBagSize(i); j++)
 		    {
 		      INDEX_3 pts;
 		      int dummy; // , val;   
-
 		      idpts.GetData (i, j, pts, dummy);
+                */
+                for (auto [hash, val] : idpts)
+                  {
+                    auto [hash_pts, hash_nr] = hash;
+                    auto [pi1, pi2] = hash_pts;
                       // val = pts[2];   
-		      const Point3d & p1 = mesh->Point(pts.I1());
-		      const Point3d & p2 = mesh->Point(pts.I2());
+		      const Point3d & p1 = mesh->Point(pi1);
+		      const Point3d & p2 = mesh->Point(pi2);
 
 		      glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
 				    identifiedcol);