Merge branch 'constexpr_experiments' into 'master'

Constexpr experiments

See merge request ngsolve/netgen!694

libsrc/core/hashtable.hpp

@@ -834,6 +834,12 @@ namespace ngcore
       hash = T_HASH(invalid);
       used = 0;
     }
+
+    void DoArchive (Archive & ar)
+    {
+      ar & hash & cont;
+      ar & size & mask & used;
+    }    
     
     class Iterator
     {

libsrc/csg/edgeflw.cpp

@@ -1714,12 +1714,12 @@ namespace netgen
 	if (oldseg.seginfo == 0)
 	  continue;
 
-	int pi1 = oldseg[0];
+	PointIndex pi1 = oldseg[0];
-	int pi2 = oldseg[1];
+	PointIndex pi2 = oldseg[1];
 
-	int npi1 = geometry.identifications.Get(copyedgeidentification)
+	PointIndex npi1 = geometry.identifications.Get(copyedgeidentification)
 	  -> GetIdentifiedPoint (mesh, pi1);
-	int npi2 = geometry.identifications.Get(copyedgeidentification)
+	PointIndex npi2 = geometry.identifications.Get(copyedgeidentification)
 	  -> GetIdentifiedPoint (mesh, pi2);
 
 	//(*testout) << "copy edge, pts = " << npi1 << " - " << npi2 << endl;
@@ -1885,12 +1885,10 @@ namespace netgen
 
 	    if (seg1.domin != -1 || seg1.domout != -1)
 	      {
-		mesh.AddPoint (p1, layer, EDGEPOINT);
+                seg1[0] = mesh.AddPoint (p1, layer, EDGEPOINT);
-		mesh.AddPoint (p2, layer, EDGEPOINT);
+                seg1[1] = mesh.AddPoint (p2, layer, EDGEPOINT);
-		seg1[0] = mesh.GetNP()-1;
+                seg2[0] = seg1[1];
-		seg1[1] = mesh.GetNP();
+                seg2[1] = seg1[0];
-		seg2[1] = mesh.GetNP()-1;
-		seg2[0] = mesh.GetNP();
 		seg1.geominfo[0].trignum = 1;
 		seg1.geominfo[1].trignum = 1;
 		seg2.geominfo[0].trignum = 1;

libsrc/csg/genmesh.cpp

@@ -35,7 +35,7 @@ namespace netgen
         auto up = geom.GetUserPoint(i);
 	auto pnum = mesh.AddPoint(up);
 	mesh.Points().Last().Singularity (geom.GetUserPointRefFactor(i));
-	mesh.AddLockedPoint (PointIndex (i+1));
+	mesh.AddLockedPoint (pnum);
         int index = up.GetIndex();
         if (index == -1)
           index = mesh.AddCD3Name (up.GetName())+1;

libsrc/csg/identify.cpp

@@ -289,14 +289,14 @@ GetIdentifiedPoint (class Mesh & mesh, PointIndex pi)
   // project to other surface
   snew->Project (hp);
 
-  int newpi = 0;
+  PointIndex newpi(PointIndex::INVALID);
-  for (int i = 1; i <= mesh.GetNP(); i++)
+  for (PointIndex pi : Range(mesh.Points()))
-    if (Dist2 (mesh.Point(i), hp) < 1e-12)
+    if (Dist2 (mesh.Point(pi), hp) < 1e-12)
       {
-	newpi = i;
+	newpi = pi;
 	break;
       }
-  if (!newpi)
+  if (!newpi.IsValid())
     newpi = mesh.AddPoint (hp);
 
   if (snew == s2)
@@ -322,6 +322,7 @@ void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh)
   mesh.GetBox(p1, p2);
   auto eps = 1e-6 * (p2-p1).Length();
 
+  /*
   for (int i = 1; i <= mesh.GetNP(); i++)
     {
       Point<3> p = mesh.Point(i);
@@ -334,13 +335,24 @@ void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh)
 	    if (Dist2(mesh.Point(j), pp) < eps)
 	      {
 		mesh.GetIdentifications().Add (i, j, nr);
-		/*
-		(*testout) << "Identify points(periodic:), nr = " << nr << ": "
-			   << mesh.Point(i) << " - " << mesh.Point(j) << endl;
-		*/
 	      }
 	}
     }
+  */
+
+  for (auto pi : Range(mesh.Points()))
+    {
+      Point<3> p = mesh[pi];
+      if (s1->PointOnSurface (p))
+	{
+	  Point<3> pp = p;
+          pp = trafo(pp);
+	  s2->Project (pp);
+          for (PointIndex pj : Range(mesh.Points()))
+	    if (Dist2(mesh[pj], pp) < eps)
+              mesh.GetIdentifications().Add (pi, pj, nr);
+	}
+    }
 
   mesh.GetIdentifications().SetType(nr,Identifications::PERIODIC);
 }

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[0] = pair.I1();
-		      i4[1] = pair.I2();
+		      // i4[1] = pair.I2();
+		      i4[0] = pi1; 
+		      i4[1] = pi2; 
 		      i4[2] = idnr;
 		      i4[3] = csid->GetSlices().Size();
 		      ref_slices.Append (i4);

libsrc/general/hashtabl.cpp

@@ -292,7 +292,19 @@ namespace netgen
 
 
 
-
+  BASE_INDEX_3_CLOSED_HASHTABLE ::
+  BASE_INDEX_3_CLOSED_HASHTABLE (size_t size)
+    : hash(RoundUp2(size))
+  {
+    // cout << "orig size = " << size
+    // << ", roundup size = " << hash.Size();
+    size = hash.Size();
+    mask = size-1;
+    // cout << "mask = " << mask << endl;
+    invalid = -1;
+    for (size_t i = 0; i < size; i++)
+      hash[i].I1() = invalid;
+  }
 
 
   void BASE_INDEX_3_CLOSED_HASHTABLE ::

libsrc/general/hashtabl.hpp

@@ -861,9 +861,10 @@ inline ostream & operator<< (ostream & ost, const INDEX_2_CLOSED_HASHTABLE<T> &
   for (int i = 0; i < ht.Size(); i++)
     if (ht.UsedPos(i))
       {
-	INDEX_2 hash;
+	// INDEX_2 hash;
-	T data;
+	// T data;
-	ht.GetData0 (i, hash, data);
+	// ht.GetData0 (i, hash, data);
+        auto [hash,data] = ht.GetBoth(i);
 	ost << "hash = " << hash << ", data = " << data << endl;
       }
   return ost;
@@ -880,7 +881,8 @@ protected:
   size_t mask;
 
 protected: 
-  BASE_INDEX_3_CLOSED_HASHTABLE (size_t size)
+  BASE_INDEX_3_CLOSED_HASHTABLE (size_t size);
+  /*
     : hash(RoundUp2(size))
   {
     // cout << "orig size = " << size
@@ -892,6 +894,7 @@ protected:
     for (size_t i = 0; i < size; i++)
       hash[i].I1() = invalid;
   }
+  */
 
 public:
   int Size() const 
@@ -1073,9 +1076,12 @@ inline ostream & operator<< (ostream & ost, const INDEX_3_CLOSED_HASHTABLE<T> &
   for (int i = 0; i < ht.Size(); i++)
     if (ht.UsedPos(i))
       {
+        /*
 	INDEX_3 hash;
 	T data;
-	ht.GetData (i, hash, data);
+        ht.GetData (i, hash, data);
+        */
+        auto [hash, data] = ht.GetBoth();
 	ost << "hash = " << hash << ", data = " << data << endl;
       }
   return ost;

libsrc/general/template.hpp

@@ -114,8 +114,10 @@ class INDEX_2
 
 public:
   ///
+  // protected:
   INDEX_2 () { }
   INDEX_2 (const INDEX_2&) = default;
+public:
   INDEX_2 (INDEX_2&&) = default;
 
   INDEX_2 & operator= (const INDEX_2&) = default;
@@ -157,7 +159,7 @@ public:
       return INDEX_2 (i1,i2);
   }
 
-
+  operator std::array<INDEX,2>() { return { i[0], i[1] }; }
   ///
   INDEX & I1 () { return i[0]; }
   ///
@@ -213,13 +215,10 @@ public:
   ///
   constexpr INDEX_3 (INDEX ai1, INDEX ai2, INDEX ai3)
     : i{ai1, ai2, ai3} { }
-  // { i[0] = ai1; i[1] = ai2; i[2] = ai3; }
 
-  ///
+  /// 
   constexpr INDEX_3 (const INDEX_3 & in2)
     : i{in2.i[0], in2.i[1], in2.i[2]} { } 
-  // { i[0] = in2.i[0]; i[1] = in2.i[1]; i[2] = in2.i[2]; }
-
 
   static INDEX_3 Sort (INDEX_3 i3)
   {

libsrc/include/nginterface_v2_impl.hpp

@@ -327,7 +327,7 @@ template <> NGX_INLINE DLL_HEADER const Ng_Node<2> Ngx_Mesh :: GetNode<2> (int n
 {
   Ng_Node<2> node;
   node.vertices.ptr = (const int*)mesh->GetTopology().GetFaceVerticesPtr(nr);
-  node.vertices.nv = (node.vertices.ptr[3] == 0) ? 3 : 4;
+  node.vertices.nv = (node.vertices.ptr[3]+1 == PointIndex::BASE) ? 3 : 4;
   node.surface_el = mesh->GetTopology().GetFace2SurfaceElement (nr+1)-1;
   return node;
 }
@@ -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.I1() -= IndexBASE<PointIndex>();
-      ind.I2()--;
+      ind.I2() -= IndexBASE<PointIndex>();
     }
   typedef int ti2[2];
   return { apairs.Size(), (ti2*)(void*)apairs.Release() };

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)

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);

libsrc/meshing/adfront2.hpp

@@ -95,7 +95,7 @@ namespace netgen
   {
   private:
     /// Point Indizes
-    INDEX_2 l;            
+    INDEX_2 l;  // want to replace by std::array<int,2> l;
     /// quality class 
     int lineclass;      
     /// geometry specific data
@@ -109,23 +109,12 @@ namespace netgen
 
     ///
     FrontLine (const INDEX_2 & al)
-    {
+      : l(al), lineclass(1) { } 
-      l = al;
-      lineclass = 1;
-    }
-
 
     ///
-    const INDEX_2 & L () const
+    const auto & L () const { return l; }
-    {
-      return l;
-    }
-
     ///
-    int LineClass() const
+    int LineClass() const { return lineclass; }
-    {
-      return lineclass;
-    }
 
     ///
     void IncrementClass ()
@@ -141,13 +130,13 @@ namespace netgen
     ///
     bool Valid () const
     {
-      return l.I1() != -1;
+      return l[0] != -1;
     }
     ///
     void Invalidate ()
     {
-      l.I1() = -1;
+      l[0] = -1;
-      l.I2() = -1;
+      l[1] = -1;
       lineclass = 1000;
     }
 

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& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}];
-        auto p1 = Point(p1p2.I1());
+        auto& ipts = inserted_points[ { hash_pts[0], hash_pts[1] }];
-        auto p2 = Point(p1p2.I2());
+        auto p1 = Point(hash_pts.I1());
-        ipts.Append(p1p2.I1());
+        auto p2 = Point(hash_pts.I2());
-        mapped_points.SetBit(p1p2.I1());
+        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 << " ";
+            (*testout) << "illegal tet: " << ei << " ";
-            for (j = 1; j <= el.GetNP(); j++)
+            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++;
       }
 

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;
 
@@ -507,7 +507,7 @@ namespace netgen
     ///
     LocalH & LocalHFunction (int layer=1) { return * lochfunc[layer-1]; }
 
-    shared_ptr<LocalH> GetLocalH(int layer=1) const
+    shared_ptr<LocalH> & GetLocalH(int layer=1) const
     {
       if(lochfunc.Size() == 1)
         return lochfunc[0];

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())
+         // for(const auto& [p1p2, dummy] : m.GetIdentifications().GetIdentifiedPoints())
-           mesh.GetIdentifications().Add(pmap[p1p2[0]], pmap[p1p2[1]], p1p2[2]);
+         // 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));

libsrc/meshing/meshtype.cpp

@@ -2709,8 +2709,8 @@ namespace netgen
     {
       ar & maxidentnr;
       ar & identifiedpoints & identifiedpoints_nr;
-
       ar & idpoints_table;
+      
       if (ar.Output())
         {
           size_t s = type.Size();
@@ -2735,11 +2735,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);
+    // INDEX_3 tripl (pi1, pi2, identnr);
-    identifiedpoints_nr.Set (tripl, 1);
+    identifiedpoints_nr.Set ( { { pi1, pi2 }, identnr }, 1);
 
     if (identnr > maxidentnr) maxidentnr = identnr;
     names.SetSize(maxidentnr);
@@ -2762,8 +2762,9 @@ namespace netgen
 
   bool Identifications :: Get (PointIndex pi1, PointIndex pi2, int nr) const
   {
-    INDEX_3 tripl(pi1, pi2, nr);
+    // INDEX_3 tripl(pi1, pi2, nr);
-    if (identifiedpoints_nr.Used (tripl))
+    // if (identifiedpoints_nr.Used (tripl))
+    if (identifiedpoints_nr.Used ( { { pi1, pi1 }, nr } ) )
       return 1;
     else
       return 0;
@@ -2803,23 +2804,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 +2838,12 @@ namespace netgen
   Array<INDEX_3> Identifications :: GetPairs () const
   {
     Array<INDEX_3> pairs;
-    for(auto [i3, dummy] : identifiedpoints_nr)
+    for(auto [hash, dummy] : identifiedpoints_nr)
-      pairs.Append(i3);
+      // pairs.Append(i3);
+      {
+        auto [pts,nr] = hash;
+        pairs.Append ( { pts[0], pts[1], nr } );
+      }
     return pairs;
   }
 
@@ -2841,6 +2853,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 +2862,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 +2879,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 +2907,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 +2941,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;
   }
 

libsrc/meshing/meshtype.hpp

@@ -292,6 +292,9 @@ namespace netgen
     constexpr PointIndex operator[] (int i) const { return PointIndex(INDEX_2::operator[](i)); }
     PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_2::operator[](i)); }
 
+    template <typename ARCHIVE>
+    void DoArchive(ARCHIVE& ar) { ar.Do(&I1(), 2); }
+    
     PointIndex & I1 () { return (*this)[0]; }
     PointIndex & I2 () { return (*this)[1]; }
     PointIndex I1 () const { return (*this)[0]; }
@@ -314,6 +317,10 @@ namespace netgen
     constexpr 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)); }
+
+    template <typename ARCHIVE>
+    void DoArchive(ARCHIVE& ar) { ar.Do(&I1(), 3); }
+    
     PointIndex & I1 () { return (*this)[0]; }
     PointIndex & I2 () { return (*this)[1]; }
     PointIndex & I3 () { return (*this)[2]; }
@@ -336,6 +343,9 @@ namespace netgen
     PointIndex operator[] (int i) const { return PointIndex(INDEX_4::operator[](i)); }
     PointIndex & operator[] (int i) { return reinterpret_cast<PointIndex&>(INDEX_4::operator[](i)); }
 
+    template <typename ARCHIVE>
+    void DoArchive(ARCHIVE& ar) { ar.Do(&I1(), 4); }
+    
     PointIndex & I1 () { return (*this)[0]; }
     PointIndex & I2 () { return (*this)[1]; }
     PointIndex & I3 () { return (*this)[2]; }
@@ -1763,13 +1773,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 +1817,7 @@ namespace netgen
 
     // bool HasIdentifiedPoints() const { return identifiedpoints != nullptr; } 
     ///
-    INDEX_3_HASHTABLE<int> & GetIdentifiedPoints ()
+    auto & GetIdentifiedPoints ()
     { 
       return identifiedpoints_nr;
     }

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;

libsrc/meshing/secondorder.cpp

@@ -174,8 +174,8 @@ namespace netgen
 	int nnp = newel.GetNP();
 	for (int j = 0; j < nnp-onp; j++)
 	  {
-	    int pi1 = newel[betw[j][0]];
+	    PointIndex pi1 = newel[betw[j][0]];
-	    int pi2 = newel[betw[j][1]];
+	    PointIndex pi2 = newel[betw[j][1]];
 	  
 	    INDEX_2 i2 = INDEX_2::Sort (pi1, pi2);
 	  
@@ -424,7 +424,6 @@ namespace netgen
   {
     PrintMessage (3, "Validate mesh");
     int np = mesh.GetNP();
-    int ne = mesh.GetNE();
     // int i, j;
     NgArray<INDEX_2> parents(np);
   
@@ -558,6 +557,7 @@ namespace netgen
 		//	      (*testout) << "bad els: " << endl;
 		wrongels = 0;
 		for (int i = 1; i <= ne; i++)
+                  
 		  {
 		    if (!illegalels.Test(i) && 
 			mesh.VolumeElement(i).

libsrc/meshing/topology.cpp

@@ -160,7 +160,7 @@ namespace netgen
             
             { // triangle
               PointIndices<4> face(el[elfaces[j][0]], el[elfaces[j][1]], 
-                                   el[elfaces[j][2]], 0);
+                                   el[elfaces[j][2]], PointIndex(PointIndex::INVALID));
 
 
               [[maybe_unused]] int facedir = 0;
@@ -266,7 +266,7 @@ namespace netgen
             
             PointIndices<4> face(el.PNum(elfaces[0][0]),
                                  el.PNum(elfaces[0][1]),
-                                 el.PNum(elfaces[0][2]),0);
+                                 el.PNum(elfaces[0][2]), PointIndex(PointIndex::INVALID));
             
             // facedir = 0;
             if (face[0] > face[1])
@@ -1161,7 +1161,7 @@ namespace netgen
                       size_t pos;
                       if (vert2face.PositionCreate(face, pos))
                         {
-                          face2vert[nfa] = { face[0], face[1], face[2], 0 }; // i4;
+                          face2vert[nfa] = { face[0], face[1], face[2], PointIndex::BASE-1 }; // i4;
                           vert2face.SetData (pos, face, nfa);
                           nfa++;
                         }
@@ -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());
 
@@ -2131,7 +2136,7 @@ namespace netgen
     if (elfaces[j][3] < 0)
       { // triangle
         INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], 
-                     el[elfaces[j][2]], 0);
+                     el[elfaces[j][2]], PointIndex::BASE-1 );
         
         int facedir = 0;
         if (face.I1() > face.I2())
@@ -2198,7 +2203,7 @@ namespace netgen
     if (elfaces[j][3] < 0)
       { // triangle
         INDEX_4 face(el[elfaces[j][0]], el[elfaces[j][1]], 
-                     el[elfaces[j][2]], 0);
+                     el[elfaces[j][2]], PointIndex(PointIndex::INVALID));
         
         int facedir = 0;
         if (face.I1() > face.I2())
@@ -2267,7 +2272,7 @@ namespace netgen
     vertices.SetSize(4);
     for (int i = 0; i < 4; i++)
       vertices[i] = face2vert[fnr-1][i];
-    if (vertices[3] == 0)
+    if (vertices[3]+1==PointIndex::BASE)
       vertices.SetSize(3);
   }
 

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 ||
+            if (mesh->VolumeElement(ei).Flags().badel ||
-		mesh->VolumeElement(i).Flags().illegal ||
+		mesh->VolumeElement(ei).Flags().illegal ||
-		(i == vispar.drawelement))
+		(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 & p1 = mesh->Point(pi1);
-		      const Point3d & p2 = mesh->Point(pts.I2());
+		      const Point3d & p2 = mesh->Point(pi2);
 
 		      glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE,
 				    identifiedcol);