From f9d7d3a4fdc23f713809188d11573083a185daf4 Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Mon, 24 Jun 2024 09:56:44 +0200
Subject: [PATCH] store newest vertex from bisection

---
 libsrc/include/nginterface_v2.hpp      |  1 +
 libsrc/include/nginterface_v2_impl.hpp |  2 ++
 libsrc/meshing/bisect.cpp              | 10 +++++++++-
 libsrc/meshing/meshtype.hpp            | 10 +++++++++-
 4 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp
index 5c48fe68..39aa3107 100644
--- a/libsrc/include/nginterface_v2.hpp
+++ b/libsrc/include/nginterface_v2.hpp
@@ -160,6 +160,7 @@ namespace netgen
     FlatArray<T_FACE2> faces;    
     Ng_Facets facets;
     bool is_curved;
+    int8_t newest_vertex;
   };
 
   
diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp
index a6f56df8..facd6c81 100644
--- a/libsrc/include/nginterface_v2_impl.hpp
+++ b/libsrc/include/nginterface_v2_impl.hpp
@@ -191,6 +191,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<2> (size_t nr) const
       ret.facets.ptr = ret.edges.Data();
     }
   ret.is_curved = el.IsCurved();
+  ret.newest_vertex = el.NewestVertex();
   return ret;
 }
 
@@ -226,6 +227,7 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<3> (size_t nr) const
   ret.facets.ptr = ret.faces.Data();
 
   ret.is_curved = el.IsCurved();
+  ret.newest_vertex = el.NewestVertex();
   return ret;
 }
 
diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index 02fceb22..409750f3 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -34,6 +34,7 @@ namespace netgen
     // unsigned char faceedges[4];
     bool incorder;
     unsigned int order:6;
+    int8_t newest_vertex;
 
     MarkedTet() = default;
     /*
@@ -195,6 +196,7 @@ namespace netgen
 
     bool incorder;
     unsigned int order:6;
+    int8_t newest_vertex;
   };
   
   ostream & operator<< (ostream & ost, const MarkedTri & mt)
@@ -1258,6 +1260,8 @@ namespace netgen
     newtet1.marked = nm;
     newtet2.marked = nm;
 
+    newtet1.newest_vertex = oldtet.newest_vertex;
+
 #ifdef DEBUG
     *testout << "newtet1,before = " << newtet1 << endl;
     *testout << "newtet2,before = " << newtet2 << endl;
@@ -1267,6 +1271,7 @@ namespace netgen
       {
 	if (i == oldtet.tetedge1)
 	  {
+            newtet2.newest_vertex = i;
 	    newtet2.pnums[i] = newp;
 	    newtet2.faceedges[i] = oldtet.faceedges[i];  // inherited face
 	    newtet2.faceedges[vis1] = i;        // cut faces
@@ -1463,11 +1468,12 @@ namespace netgen
     newtri1.pnums[pe2] = newp;
     newtri1.pgeominfo[pe2] = newpgi;
     newtri1.markededge = pe2;
+    newtri1.newest_vertex = oldtri.newest_vertex;
 
     newtri2.pnums[pe1] = newp;
     newtri2.pgeominfo[pe1] = newpgi;
     newtri2.markededge = pe1;
-
+    newtri2.newest_vertex = pe1;
 
     newtri1.surfid = oldtri.surfid;
     newtri2.surfid = oldtri.surfid;
@@ -3715,6 +3721,7 @@ namespace netgen
             el.SetOrder (tet.order);
             for (int j = 0; j < 4; j++)
               el[j] = tet.pnums[j];
+            el.NewestVertex() = tet.newest_vertex;
             mesh.SetVolumeElement (ElementIndex(i), el);
           }
        });
@@ -3811,6 +3818,7 @@ namespace netgen
                 el[j] = trig.pnums[j];
                 el.GeomInfoPi(j+1) = trig.pgeominfo[j];
               }
+            el.NewestVertex() = trig.newest_vertex;
             mesh.SetSurfaceElement (SurfaceElementIndex(i), el);
           }
        });
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index d765980f..d814b850 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -425,6 +425,7 @@ namespace netgen
     // control whether it is visible or not
     bool visible:1;  // element visible
     bool is_curved;   // element is (high order) curved
+    int8_t newest_vertex = -1; // from refinement via bisection
     /// order for hp-FEM
     unsigned int orderx:6;
     unsigned int ordery:6;
@@ -562,6 +563,9 @@ namespace netgen
     ///
     const PointGeomInfo & GeomInfoPiMod (int i) const { return geominfo[(i-1) % np]; }
 
+    auto & NewestVertex() { return newest_vertex; }
+    auto NewestVertex() const { return newest_vertex; }
+
     void DoArchive (Archive & ar)
     {
       short _np, _typ;
@@ -731,7 +735,8 @@ namespace netgen
     ELEMENT_TYPE typ;
     /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism)
     int8_t np;
-
+    int8_t newest_vertex = -1; // from refinement via bisection
+    
     /// sub-domain index
     int index;
     /// order for hp-FEM
@@ -856,6 +861,9 @@ namespace netgen
     ///
     const PointIndex & PNumMod (int i) const { return pnum[(i-1) % np]; }
 
+    auto & NewestVertex() { return newest_vertex; }
+    auto NewestVertex() const { return newest_vertex; }
+
     void DoArchive (Archive & ar)
     {
       short _np, _typ;