From 5f93b0584156b8e63cdcdb3574f14fb5ce4fd8d3 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 18 Sep 2009 15:05:54 +0000
Subject: [PATCH] add an option to remove four-triangles vertices

---
 Common/Context.h           |   2 +-
 Common/DefaultOptions.h    |   2 +
 Common/Options.cpp         |  12 ++++
 Common/Options.h           |   1 +
 Fltk/optionWindow.cpp      |   6 ++
 Mesh/meshGFace.cpp         |   2 +
 Mesh/meshGFaceOptimize.cpp | 132 ++++++++++++++++++++++++++++++++++---
 Mesh/meshGFaceOptimize.h   |   1 +
 8 files changed, 147 insertions(+), 11 deletions(-)

diff --git a/Common/Context.h b/Common/Context.h
index 539f3677a1..4a00be4e9a 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -166,7 +166,7 @@ class CTX {
     int points, lines, triangles, quadrangles, tetrahedra, hexahedra, prisms, pyramids;
     int surfacesEdges, surfacesFaces, volumesEdges, volumesFaces, numSubEdges;
     int pointsNum, linesNum, surfacesNum, volumesNum;
-    int optimize, optimizeNetgen, refineSteps, qualityType, labelType;
+    int optimize, optimizeNetgen, refineSteps, qualityType, labelType, remove4triangles;
     double normals, tangents, explode, angleSmoothNormals, allowSwapEdgeAngle;
     double mshFileVersion, labelFrequency, pointSize, lineWidth;
     double qualityInf, qualitySup, radiusInf, radiusSup;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index cadc1fb5c3..39fec36ddf 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1166,6 +1166,8 @@ StringXNumber MeshOptions_Number[] = {
     "RandomFactor * size(triangle)/size(model) approaches machine accuracy)" },
   { F|O, "RefineSteps" , opt_mesh_refine_steps , 10 ,
     "Number of refinement steps in the MeshAdapt-based 2D algorithms" }, 
+  { F|O, "Remove4Triangles" , opt_mesh_remove_4_triangles , 0 ,
+    "Try to remove nodes surrounded by 4 triangles in 2D triangular meshes" }, 
   { F|O, "ReverseAllNormals" , opt_mesh_reverse_all_normals , 0. , 
     "Reverse all the mesh normals (for display)" },
 
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 5a113c3ba8..a9324779f5 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4693,6 +4693,18 @@ double opt_mesh_optimize_netgen(OPT_ARGS_NUM)
 #endif
   return CTX::instance()->mesh.optimizeNetgen;
 }
+double opt_mesh_remove_4_triangles(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX::instance()->mesh.remove4triangles = (int)val;
+#if defined(HAVE_FLTK)
+  if(FlGui::available() && (action & GMSH_GUI))
+    FlGui::instance()->options->mesh.butt[25]->value
+      (CTX::instance()->mesh.remove4triangles);
+#endif
+  return CTX::instance()->mesh.remove4triangles;
+  
+}
 
 double opt_mesh_refine_steps(OPT_ARGS_NUM)
 {
diff --git a/Common/Options.h b/Common/Options.h
index 5d33eb1d85..1d2e060a11 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -446,6 +446,7 @@ double opt_geometry_match_geom_and_mesh(OPT_ARGS_NUM);
 double opt_mesh_label_frequency(OPT_ARGS_NUM);
 double opt_mesh_optimize(OPT_ARGS_NUM);
 double opt_mesh_optimize_netgen(OPT_ARGS_NUM);
+double opt_mesh_remove_4_triangles(OPT_ARGS_NUM);
 double opt_mesh_refine_steps(OPT_ARGS_NUM);
 double opt_mesh_normals(OPT_ARGS_NUM);
 double opt_mesh_num_sub_edges(OPT_ARGS_NUM);
diff --git a/Fltk/optionWindow.cpp b/Fltk/optionWindow.cpp
index 87fd8894c1..d372f09d25 100644
--- a/Fltk/optionWindow.cpp
+++ b/Fltk/optionWindow.cpp
@@ -400,6 +400,7 @@ static void mesh_options_ok_cb(Fl_Widget *w, void *data)
   opt_mesh_lc_extend_from_boundary(0, GMSH_SET, o->mesh.butt[16]->value());
   opt_mesh_optimize(0, GMSH_SET, o->mesh.butt[2]->value());
   opt_mesh_optimize_netgen(0, GMSH_SET, o->mesh.butt[24]->value());
+  opt_mesh_remove_4_triangles(0,GMSH_SET, o->mesh.butt[25]->value());
   opt_mesh_order(0, GMSH_SET, o->mesh.value[3]->value());
   opt_mesh_smooth_internal_edges(0, GMSH_SET, o->mesh.butt[3]->value());
   opt_mesh_second_order_incomplete(0, GMSH_SET, o->mesh.butt[4]->value());
@@ -1985,6 +1986,11 @@ optionWindow::optionWindow(int deltaFontSize)
         (L + 2 * WB, 2 * WB + 6 * BH, BW, BH, "Optimize high order mesh (2D-plane only)");
       mesh.butt[3]->type(FL_TOGGLE_BUTTON);
       mesh.butt[3]->callback(mesh_options_ok_cb);
+      
+      mesh.butt[25] = new Fl_Check_Button
+        (L + 2 * WB, 2 * WB + 7 * BH, BW, BH, "Try to remove 4 triangles nodes");
+      mesh.butt[25]->type(FL_TOGGLE_BUTTON);
+      mesh.butt[25]->callback(mesh_options_ok_cb);
 
       o->end();
     }
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 94598e5b3b..ed91872a29 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -663,6 +663,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
     sprintf(name, "param%d.pos", gf->tag());
     outputScalarField(m->triangles, name,1);
   }
+  if(CTX::instance()->mesh.remove4triangles)
+    removeFourTrianglesNodes(gf,false);
 
   // delete the mesh
   delete m;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index a7f3b14645..c837123716 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -205,6 +205,128 @@ void parametricCoordinates(MElement *t, GFace *gf, double u[4], double v[4])
   }
 }
 
+double surfaceTriangleUV(MElement *t,GFace *gf){
+  double u[4],v[4];
+  parametricCoordinates(t,gf,u,v);
+  return 0.5*fabs((u[1]-u[0])*(v[2]-v[0])-(u[2]-u[0])*(v[1]-v[0]));
+}
+
+double surfaceTriangleUV(MVertex *v1, MVertex *v2, MVertex *v3,           
+                         const std::vector<double> &Us,
+                         const std::vector<double> &Vs)
+{
+  const double v12[2] = {Us[v2->getIndex()] - Us[v1->getIndex()],
+                         Vs[v2->getIndex()] - Vs[v1->getIndex()]};
+  const double v13[2] = {Us[v3->getIndex()] - Us[v1->getIndex()],
+                         Vs[v3->getIndex()] - Vs[v1->getIndex()]};
+  return 0.5 * fabs (v12[0] * v13[1] - v12[1] * v13[0]);
+}
+
+int _removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
+{
+
+  v2t_cont adj;
+  buildVertexToElement(gf->triangles,adj);
+  v2t_cont :: iterator it = adj.begin();
+  int n=0;
+  std::set<MElement*> touched;
+  while (it != adj.end()) {
+    bool skip=false;
+    double surfaceRef=0;
+    if(it->second.size()==4) {
+      const std::vector<MElement*> &lt = it->second;
+      MVertex* edges[4][2];
+      for(int i=0;i<4;i++) {
+        if(touched.find(lt[i])!=touched.end() || lt[i]->getNumVertices()!=3){
+          skip=true;
+          break;
+        }
+        int j;
+
+        surfaceRef+=surfaceTriangleUV(lt[i],gf);
+        for(j=0;j<3;j++) {
+          if(lt[i]->getVertex(j)==it->first) {
+            edges[i][0] = lt[i]->getVertex((j+1)%3);
+            edges[i][1] = lt[i]->getVertex((j+2)%3);
+            break;
+          }
+        }
+        if(j==3)
+          throw;
+      }
+      if(skip){
+        it++;
+        continue;
+      }
+
+      for(int i=1;i<3;i++) {
+        for(int j=i+1;j<4;j++) {
+          if(edges[j][0]==edges[i-1][1]){
+            MVertex *buf[2]={edges[i][0],edges[i][1]};
+            edges[i][0]=edges[j][0];
+            edges[i][1]=edges[j][1];
+            edges[j][0]=buf[0];
+            edges[j][1]=buf[1];
+            break;
+          }
+        }
+      }
+      if(edges[0][1]==edges[1][0] && edges[1][1]==edges[2][0] && edges[2][1] == edges[3][0] && edges[3][1]==edges[0][0]) {
+        if(replace_by_quads){
+          gf->quadrangles.push_back(new MQuadrangle(edges[0][0],edges[1][0],edges[2][0],edges[3][0]));
+        }else{
+          MTriangle *newt[4];
+          double surf[4],qual[4];
+          for(int i=0;i<4;i++){
+            newt[i] = new MTriangle(edges[i][0],edges[(i+1)%4][0],edges[(i+2)%4][0]);
+            surf[i]=surfaceTriangleUV(newt[i],gf);
+            qual[i]=qmTriangle(newt[i],QMTRI_RHO);
+          }
+          double q02=(fabs((surf[0]+surf[2]-surfaceRef)/surfaceRef)<1e-8) ? std::min(qual[0],qual[2]) : -1;
+          double q13=(fabs((surf[1]+surf[3]-surfaceRef)/surfaceRef)<1e-8) ? std::min(qual[1],qual[3]) : -1;
+          if(q02>q13 && q02 >0) {
+            delete newt[1];
+            delete newt[3];
+            gf->triangles.push_back(newt[0]);
+            gf->triangles.push_back(newt[2]);
+          } else if (q13 >0) {
+            delete newt[0];
+            delete newt[2];
+            gf->triangles.push_back(newt[1]);
+            gf->triangles.push_back(newt[3]);
+          } else {
+            it++;
+            throw;
+            continue;
+          }
+        }
+        for(int i=0;i<4;i++) {
+          touched.insert(lt[i]);
+        }
+        n++;
+      }
+    }
+    it++;
+  }
+  std::vector<MTriangle*> triangles2;
+  for(unsigned int i = 0; i < gf->triangles.size(); i++){
+    if(touched.find(gf->triangles[i]) == touched.end()){
+      triangles2.push_back(gf->triangles[i]);
+    }
+    else {
+      delete gf->triangles[i];
+    }    
+  } 
+  gf->triangles = triangles2;
+  Msg::Debug("%i four-triangles vertices removed",n);
+  return n;
+}
+
+void removeFourTrianglesNodes(GFace *gf,bool replace_by_quads)
+{
+  while(_removeFourTrianglesNodes(gf,replace_by_quads));
+}
+
 void laplaceSmoothing(GFace *gf)
 {
   v2t_cont adj;
@@ -252,16 +374,6 @@ void laplaceSmoothing(GFace *gf)
   }
 }
 
-double surfaceTriangleUV(MVertex *v1, MVertex *v2, MVertex *v3,           
-                         const std::vector<double> &Us,
-                         const std::vector<double> &Vs)
-{
-  const double v12[2] = {Us[v2->getIndex()] - Us[v1->getIndex()],
-                         Vs[v2->getIndex()] - Vs[v1->getIndex()]};
-  const double v13[2] = {Us[v3->getIndex()] - Us[v1->getIndex()],
-                         Vs[v3->getIndex()] - Vs[v1->getIndex()]};
-  return 0.5 * fabs (v12[0] * v13[1] - v12[1] * v13[0]);
-}
 
 bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
                   std::vector<MTri3*> &newTris, const gmshSwapCriterion &cr,               
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 67f7fa2f00..6c65013e6c 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -41,6 +41,7 @@ int edgeSplitPass(double maxLC, GFace *gf,
                   std::vector<double> &Vs,
                   std::vector<double> &vSizes ,
                   std::vector<double> &vSizesBGM);
+void removeFourTrianglesNodes(GFace *gf, bool replace_by_quads);
 int edgeCollapsePass(double minLC, GFace *gf, 
                      std::set<MTri3*, compareTri3Ptr> &allTris,
                      std::vector<double> &Us,
-- 
GitLab