From d08e2daa060da23965f68bfc4b083ae8dd271539 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Mon, 8 Jun 2020 10:42:26 +0200
Subject: [PATCH] Do edge swapping for faces individually with tensor product
 meshes

If the mesh contains quads, the edge swapping algorithm switches to
generic improve, which introduces quads everywhere. This is not intended
if one domain contains a tensor product mesh.
Thus, call the optimizer for each face if mesh contains quads but mp.quad
is not set.
---
 libsrc/meshing/meshfunc2d.cpp | 45 ++++++++++++++++++++++++++++++++---
 1 file changed, 42 insertions(+), 3 deletions(-)

diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp
index 88955a36..84f5d399 100644
--- a/libsrc/meshing/meshfunc2d.cpp
+++ b/libsrc/meshing/meshfunc2d.cpp
@@ -20,6 +20,19 @@ namespace netgen
       }
     mesh.Compress();
 
+    bool optimize_swap_separate_faces = false;
+    if(!mp.quad)
+    {
+      bool mixed = false;
+      ParallelFor( Range(mesh.GetNSE()), [&] (auto i) NETGEN_LAMBDA_INLINE
+          {
+            if (mesh[SurfaceElementIndex(i)].GetNP() != 3)
+                mixed = true;
+          });
+      if(mixed)
+        optimize_swap_separate_faces = true;
+    }
+
     const char * optstr = mp.optimize2d.c_str();
     int optsteps = mp.optsteps2d;
 
@@ -31,16 +44,42 @@ namespace netgen
 	    {
 	    case 's': 
 	      {  // topological swap
-		MeshOptimize2d meshopt(mesh);
+                MeshOptimize2d meshopt(mesh);
                 meshopt.SetMetricWeight (mp.elsizeweight);
-		meshopt.EdgeSwapping (0);
+
+                if(optimize_swap_separate_faces)
+                {
+                  for(auto i : Range(1, mesh.GetNFD()+1))
+                  {
+                    meshopt.SetFaceIndex(i);
+                    meshopt.EdgeSwapping (0);
+                  }
+                }
+                else
+                {
+                  meshopt.SetFaceIndex(0);
+                  meshopt.EdgeSwapping (0);
+                }
 		break;
 	      }
 	    case 'S': 
 	      {  // metric swap
 		MeshOptimize2d meshopt(mesh);
                 meshopt.SetMetricWeight (mp.elsizeweight);
-		meshopt.EdgeSwapping (1);
+
+                if(optimize_swap_separate_faces)
+                {
+                  for(auto i : Range(1, mesh.GetNFD()+1))
+                  {
+                    meshopt.SetFaceIndex(i);
+                    meshopt.EdgeSwapping (1);
+                  }
+                }
+                else
+                {
+                  meshopt.SetFaceIndex(0);
+                  meshopt.EdgeSwapping (1);
+                }
 		break;
 	      }
 	    case 'm':