From b1182e48b8a6e4ae43dcd7d7e86dee8461d6bdc2 Mon Sep 17 00:00:00 2001
From: Matthias Hochsteger <mhochsteger@cerbsim.com>
Date: Tue, 1 Oct 2024 13:28:14 +0200
Subject: [PATCH] Fix RemoveIllegalElements

- Only search in relevant domain
- Break if number of illegal elements increases (avoids infinite loop)
  -> This shouldn't actually happen and is just a workaround until the
  optimization routines are fixed
---
 libsrc/meshing/meshclass.cpp |  4 ++--
 libsrc/meshing/meshclass.hpp |  2 +-
 libsrc/meshing/meshfunc.cpp  | 16 ++++++++++------
 libsrc/meshing/meshfunc.hpp  |  2 +-
 4 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index 0a291e91..adad7339 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -6739,7 +6739,7 @@ namespace netgen
   }
 
 
-  int Mesh :: MarkIllegalElements ()
+  int Mesh :: MarkIllegalElements (int domain)
   {
     if(!boundaryedges)
       BuildBoundaryEdges();
@@ -6749,7 +6749,7 @@ namespace netgen
     {
       int cnt_local = 0;
       for(auto & el : volelements.Range(myrange))
-        if (!LegalTet (el))
+        if ((domain==0 || el.GetIndex() == domain) && !LegalTet (el))
           cnt_local++;
       cnt += cnt_local;
     });
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 41b4832b..5d7c7e7a 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -649,7 +649,7 @@ namespace netgen
       Marks elements which are dangerous to refine
       return: number of illegal elements
     */
-	DLL_HEADER int MarkIllegalElements ();
+	DLL_HEADER int MarkIllegalElements (int domain=0);
 
     /// orient surface mesh, for one sub-domain only
 	DLL_HEADER void SurfaceMeshOrientation ();
diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp
index c551301b..fc9ca918 100644
--- a/libsrc/meshing/meshfunc.cpp
+++ b/libsrc/meshing/meshfunc.cpp
@@ -521,7 +521,7 @@ namespace netgen
            throw NgException ("Stop meshing since surface mesh not consistent");
         }
      }
-     RemoveIllegalElements (mesh);
+     RemoveIllegalElements (mesh, domain);
   }
 
   void MergeMeshes( Mesh & mesh, Array<MeshingData> & md )
@@ -728,13 +728,10 @@ namespace netgen
 
 
 
-  void RemoveIllegalElements (Mesh & mesh3d)
+  void RemoveIllegalElements (Mesh & mesh3d, int domain)
   {
     static Timer t("RemoveIllegalElements"); RegionTimer reg(t);
     
-    int it = 10;
-    int nillegal, oldn;
-
     // return, if non-pure tet-mesh
     /*
       if (!mesh3d.PureTetMesh())
@@ -742,12 +739,16 @@ namespace netgen
     */
     mesh3d.CalcSurfacesOfNode();
 
-    nillegal = mesh3d.MarkIllegalElements();
+    int nillegal = mesh3d.MarkIllegalElements(domain);
     if(nillegal)
       PrintMessage (1, "Remove Illegal Elements");
 
+    int oldn = nillegal;
+    int nillegal_min = nillegal;
+
     MeshingParameters dummymp;
     MeshOptimize3d optmesh(mesh3d, dummymp, OPT_LEGAL);
+    int it = 10;
     while (nillegal && (it--) > 0)
       {
 	if (multithread.terminate)
@@ -763,6 +764,9 @@ namespace netgen
 
 	oldn = nillegal;
 	nillegal = mesh3d.MarkIllegalElements();
+        nillegal_min = min(nillegal_min, nillegal);
+        if(nillegal > nillegal_min)
+          break;
 
 	if (oldn != nillegal)
 	  it = 10;
diff --git a/libsrc/meshing/meshfunc.hpp b/libsrc/meshing/meshfunc.hpp
index 81157eee..854cdd8f 100644
--- a/libsrc/meshing/meshfunc.hpp
+++ b/libsrc/meshing/meshfunc.hpp
@@ -30,7 +30,7 @@ DLL_HEADER MESHING3_RESULT MeshVolume (const MeshingParameters & mp, Mesh& mesh3
 DLL_HEADER MESHING3_RESULT OptimizeVolume (const MeshingParameters & mp, Mesh& mesh3d);
 //			       const CSGeometry * geometry = NULL);
 
-DLL_HEADER void RemoveIllegalElements (Mesh & mesh3d);
+DLL_HEADER void RemoveIllegalElements (Mesh & mesh3d, int domain = 0);
 
 
 enum MESHING_STEP {