From 87e472b6fc638cea3663b8b64803af14174c1697 Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Wed, 17 Feb 2021 14:54:14 +0100
Subject: [PATCH] start face-hierarchy in Netgen

---
 libsrc/core/hashtable.hpp   |  7 +++++
 libsrc/meshing/topology.cpp | 59 +++++++++++++++++++++++++++++++++++--
 2 files changed, 64 insertions(+), 2 deletions(-)

diff --git a/libsrc/core/hashtable.hpp b/libsrc/core/hashtable.hpp
index 00739ffb..29f35858 100644
--- a/libsrc/core/hashtable.hpp
+++ b/libsrc/core/hashtable.hpp
@@ -174,6 +174,13 @@ namespace ngcore
     {
       return MakeTupleFromInt<N>()(*this);
     }
+
+    bool Contains (T val)
+    {
+      for (int j = 0; j < N; j++)
+        if (i[j] == val) return true;
+      return false;
+    }
   };
 
   /// sort 2 integers
diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp
index f128c10e..fa81a64f 100644
--- a/libsrc/meshing/topology.cpp
+++ b/libsrc/meshing/topology.cpp
@@ -696,8 +696,8 @@ namespace netgen
 	      {
 		auto verts = edge2vert[i];  // 2 vertices of edge
 
-		if (verts[0] > mesh->mlbetweennodes.Size()+PointIndex::BASE ||
-		    verts[1] > mesh->mlbetweennodes.Size()+PointIndex::BASE)
+		if (verts[0] >= mesh->mlbetweennodes.Size()+PointIndex::BASE ||
+		    verts[1] >= mesh->mlbetweennodes.Size()+PointIndex::BASE)
 		  continue;
                 
 		auto pa0 = mesh->mlbetweennodes[verts[0]]; // two parent vertices of v0
@@ -1436,6 +1436,61 @@ namespace netgen
 	      cout << cnt_err << " elements are not matching !!!" << endl;
 	  }
         // NgProfiler::StopTimer (timer2c);
+
+
+        if (build_hierarchy)
+          {
+            // tets only
+
+            ngcore::ClosedHashTable<INT<3>, int> v2f(nv);
+            for (auto i : Range(face2vert))
+              {
+                auto face = face2vert[i];
+                INT<3> f3(face[0], face[1], face[2]);
+                f3.Sort();
+                v2f[f3] = i;
+              }
+
+            cout << "v2f:" << endl << v2f << endl;
+
+            parent_faces.SetSize (nfa);
+            parent_faces = { -1, { -1, -1, -1, -1 } };
+
+            for (auto i : Range(nfa))
+              {
+                INT<3,PointIndex> f3(face2vert[i][0], face2vert[i][1], face2vert[i][2]);
+                PointIndex vmax = Max(f3);
+                if (vmax >= mesh->mlbetweennodes.Size()+PointIndex::BASE)
+                  continue;
+                
+                auto parents = mesh->mlbetweennodes[vmax];
+
+                // is face part of one parent face (boundary-face) ?
+                for (int j = 0; j < 2; j++)
+                  {
+                    if (f3.Contains(parents[j]))
+                      {
+                        PointIndex v0 = parents[j];
+                        PointIndex v1 = parents[1-j];
+                        
+                        // the third one, on the tip
+                        PointIndex v2 = f3[0]+f3[1]+f3[2] - v0 - v1;
+                        
+                        int classnr = 0;
+                        if (v0 > v1) { Swap (v0, v1); classnr += 1; }
+                        if (v1 > v2) { Swap (v1, v2); classnr += 2; }
+                        if (v0 > v1) { Swap (v0, v1); classnr += 4; }
+                        
+                        INT<3> parentverts(v0, v1, v2);
+                        int pafacenr = v2f[parentverts];
+                        cout << "parent-face = " << pafacenr << endl;
+                        parent_faces[i] = { classnr, { pafacenr, -1, -1, -1 } };
+                      }
+                  }
+              }
+                   
+          }
+        
       }