From 960fe8a342e605b7fc72a7e8140981c1198b2ef0 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Mon, 18 Sep 2023 13:58:43 +0200
Subject: [PATCH] Reset MeshTopology before 2d optimize

---
 libsrc/geom2d/genmesh2d.cpp    |  2 +-
 libsrc/meshing/meshfunc2d.cpp  | 21 +++++++++------------
 libsrc/meshing/python_mesh.cpp |  8 ++++----
 3 files changed, 14 insertions(+), 17 deletions(-)

diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp
index 7ab995fc..6d40bd0d 100644
--- a/libsrc/geom2d/genmesh2d.cpp
+++ b/libsrc/geom2d/genmesh2d.cpp
@@ -5,7 +5,7 @@ namespace netgen
 {
   // extern DLL_HEADER MeshingParameters mparam;
 
-  extern void Optimize2d (Mesh & mesh, MeshingParameters & mp);
+  extern void Optimize2d (Mesh & mesh, MeshingParameters & mp, int faceindex=0);
 
 
 
diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp
index 84f5d399..d3a2542b 100644
--- a/libsrc/meshing/meshfunc2d.cpp
+++ b/libsrc/meshing/meshfunc2d.cpp
@@ -4,7 +4,7 @@
 namespace netgen
 {
 
-  DLL_HEADER void Optimize2d (Mesh & mesh, MeshingParameters & mp)
+  DLL_HEADER void Optimize2d (Mesh & mesh, MeshingParameters & mp, int faceindex)
   {
     static Timer timer("optimize2d"); RegionTimer reg(timer);
 
@@ -33,19 +33,25 @@ namespace netgen
         optimize_swap_separate_faces = true;
     }
 
+    if(faceindex)
+      optimize_swap_separate_faces = false;
+
     const char * optstr = mp.optimize2d.c_str();
     int optsteps = mp.optsteps2d;
 
+    // reset topology
+    mesh.GetTopology() = MeshTopology(mesh);
     for (int i = 1; i <= optsteps; i++)
       for (size_t j = 1; j <= strlen(optstr); j++)
 	{
 	  if (multithread.terminate) break;
+          MeshOptimize2d meshopt(mesh);
+          meshopt.SetMetricWeight (mp.elsizeweight);
+          meshopt.SetFaceIndex(faceindex);
 	  switch (optstr[j-1])
 	    {
 	    case 's': 
 	      {  // topological swap
-                MeshOptimize2d meshopt(mesh);
-                meshopt.SetMetricWeight (mp.elsizeweight);
 
                 if(optimize_swap_separate_faces)
                 {
@@ -57,16 +63,12 @@ namespace netgen
                 }
                 else
                 {
-                  meshopt.SetFaceIndex(0);
                   meshopt.EdgeSwapping (0);
                 }
 		break;
 	      }
 	    case 'S': 
 	      {  // metric swap
-		MeshOptimize2d meshopt(mesh);
-                meshopt.SetMetricWeight (mp.elsizeweight);
-
                 if(optimize_swap_separate_faces)
                 {
                   for(auto i : Range(1, mesh.GetNFD()+1))
@@ -77,22 +79,17 @@ namespace netgen
                 }
                 else
                 {
-                  meshopt.SetFaceIndex(0);
                   meshopt.EdgeSwapping (1);
                 }
 		break;
 	      }
 	    case 'm': 
 	      {
-		MeshOptimize2d meshopt(mesh);
-                meshopt.SetMetricWeight (mp.elsizeweight);
 		meshopt.ImproveMesh(mp);
 		break;
 	      }
 	    case 'c': 
 	      {
-		MeshOptimize2d meshopt(mesh);
-                meshopt.SetMetricWeight (mp.elsizeweight);
 		meshopt.CombineImprove();
 		break;
 	      }
diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp
index 1dbef786..34379844 100644
--- a/libsrc/meshing/python_mesh.cpp
+++ b/libsrc/meshing/python_mesh.cpp
@@ -78,7 +78,7 @@ namespace netgen
 {
   extern bool netgen_executable_started;
   extern shared_ptr<NetgenGeometry> ng_geometry;
-  extern void Optimize2d (Mesh & mesh, MeshingParameters & mp);
+  extern void Optimize2d (Mesh & mesh, MeshingParameters & mp, int faceindex=0);
 }
 
 
@@ -1283,7 +1283,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
             OptimizeVolume (mp, self);
           }, py::arg("mp"), py::call_guard<py::gil_scoped_release>())
 
-    .def ("OptimizeMesh2d", [](Mesh & self, MeshingParameters* pars)
+    .def ("OptimizeMesh2d", [](Mesh & self, MeshingParameters* pars, int faceindex)
           {
             self.CalcLocalH(0.5);
             MeshingParameters mp;
@@ -1291,8 +1291,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
             else mp.optsteps2d = 5;
             if(!self.GetGeometry())
               throw Exception("Cannot optimize surface mesh without geometry!");
-            Optimize2d (self, mp);
-          }, py::arg("mp")=nullptr, py::call_guard<py::gil_scoped_release>())
+            Optimize2d (self, mp, faceindex);
+          }, py::arg("mp")=nullptr, py::arg("faceindex")=0, py::call_guard<py::gil_scoped_release>())
     
     .def ("Refine", FunctionPointer
           ([](Mesh & self, bool adaptive)