From e3772bf3af49af3bfa073d17637b06f38bc38550 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Mon, 16 Jan 2023 12:01:17 +0100
Subject: [PATCH] Keep consistent vertex numbering for identified edges

---
 libsrc/meshing/basegeom.cpp | 43 ++++++++++++++++++++++++++-----------
 libsrc/occ/occ_face.cpp     |  6 +++++-
 2 files changed, 36 insertions(+), 13 deletions(-)

diff --git a/libsrc/meshing/basegeom.cpp b/libsrc/meshing/basegeom.cpp
index a4b0b7a4..22c94995 100644
--- a/libsrc/meshing/basegeom.cpp
+++ b/libsrc/meshing/basegeom.cpp
@@ -637,25 +637,44 @@ namespace netgen
                 edge_params[i] = gi.dist;
             }
 
-            // reverse entries if we have decreasing parameters
-            if(edge_params.Size()>=2 && edge_params[0] > edge_params.Last())
-                for(auto i : Range((np-2)/2))
-                {
-                    swap(edge_points[i], edge_points[np-3-i]);
-                    swap(edge_params[i], edge_params[np-3-i]);
-                }
-
             params.SetSize(edge_params.Size()+2);
-            params[0] = 0.;
-            params.Last() = 1.;
 
             for(auto i : Range(edge_params))
                 params[i+1] = edge_params[i];
+
+            if(edge_params.Size()>1)
+            {
+              // Just projecting (code below) does not work for closed edges (startp == endp)
+              // In this case, there are at least 2 inner points which we use to check edge orientation
+              bool reversed = edge_params[1] < edge_params[0];
+              if(reversed)
+              {
+                params[0] = 1.0;
+                params.Last() = 0.0;
+              }
+              else
+              {
+                params.Last() = 1.0;
+                params[0] = 0.0;
+              }
+            }
+            else
+            {
+              for(size_t i : {0UL, pnums_primary.Size()-1})
+              {
+                auto p_mapped = trafo(mesh[pnums_primary[i]]);
+                EdgePointGeomInfo gi;
+                edge->ProjectPoint(p_mapped, &gi);
+                params[i] = gi.dist;
+              }
+            }
         }
 
         pnums.SetSize(edge_points.Size() + 2);
-        pnums[0] = startp;
-        pnums.Last() = endp;
+
+        bool is_reversed = params.Last() < params[0];
+        pnums[0] = is_reversed ? endp : startp;
+        pnums.Last() = is_reversed ? startp : endp;
 
 
         for(auto i : Range(edge_points))
diff --git a/libsrc/occ/occ_face.cpp b/libsrc/occ/occ_face.cpp
index 51af9c58..88fd3fb9 100644
--- a/libsrc/occ/occ_face.cpp
+++ b/libsrc/occ/occ_face.cpp
@@ -116,7 +116,11 @@ namespace netgen
                     seg.epgeominfo[i].v = uv.Y();
                 }
 
-                if(ORIENTATION == REVERSED)
+                bool do_swap = ORIENTATION == REVERSED;
+                if(seg.epgeominfo[1].dist < seg.epgeominfo[0].dist)
+                  do_swap = !do_swap;
+
+                if(do_swap)
                 {
                     swap(seg[0], seg[1]);
                     swap(seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);