Merge branch 'fixDivideAndConquerAlgorithm' into 'master'

Fix divide and conquer algorithm

See merge request gmsh/gmsh!41

Mesh/DivideAndConquer.cpp

@@ -333,7 +333,7 @@ int DocRecord::BuildDelaunay()
 
 // This routine insert the point 'newPoint' in the list dlist,
 // respecting the clock-wise orientation
-int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPoint)
+int DocRecord::DListInsert(PointNumero centerPoint, PointNumero newPoint)
 {
   DListRecord *p, *newp;
   double alpha1, alpha, beta, xx, yy;
@@ -342,6 +342,7 @@ int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPo
   newp = new DListRecord;
   newp->point_num = newPoint;
 
+  DListRecord **dlist = &points[centerPoint].adjacent;
   if(*dlist == NULL) {
     *dlist = newp;
     Pred(*dlist) = newp;
@@ -362,6 +363,8 @@ int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPo
   p = *dlist;
   first = p->point_num;
 
+  DPoint center = points[centerPoint].where;
+
   // first, compute polar coord. of the first point
   yy = (double)(points[first].where.v - center.v);
   xx = (double)(points[first].where.h - center.h);
@@ -378,15 +381,39 @@ int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPo
     yy = (double)(points[Succ(p)->point_num].where.v - center.v);
     xx = (double)(points[Succ(p)->point_num].where.h - center.h);
     alpha = atan2(yy, xx) - alpha1;
-    if(alpha <= 1.e-15)
+    if(alpha <= -1.e-15 || Succ(p)->point_num == first)
+      alpha += 2. * M_PI;
+    else if(abs(alpha) <= 1e-15 && IsRightOf(centerPoint, first, Succ(p)->point_num))
       alpha += 2. * M_PI;
-    if(alpha >= beta) {
+    if(alpha >= beta+1e-15) {
       Succ(newp) = Succ(p);
       Succ(p) = newp;
       Pred(newp) = p;
       Pred(Succ(newp)) = newp;
       return 1;
     }
+    else if (alpha >= beta - 1e-15)
+    {
+      // Angles more or less the same. To keep consistency with rest of algorithm
+      // check with IsLeftOf to see which point should be first
+      if (IsRightOf(centerPoint, Succ(p)->point_num, newPoint))
+      {
+        // New point should become before Succ(p)
+        Succ(newp) = Succ(p);
+        Succ(p) = newp;
+        Pred(newp) = p;
+        Pred(Succ(newp)) = newp;
+      }
+      else
+      {
+        // New point should become after Succ(p)
+        Succ(newp) = Succ(Succ(p));
+        Succ(Succ(p)) = newp;
+        Pred(newp) = Succ(p);
+        Pred(Succ(newp)) = newp;
+      }
+      return 1;
+    }
     p = Succ(p);
   } while(p != *dlist);
 
@@ -398,8 +425,8 @@ int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPo
 // the point 'b' in the adjency list of 'a'
 int DocRecord::Insert(PointNumero a, PointNumero b)
 {
-  int rslt = DListInsert(&points[a].adjacent, points[a].where, b);
-  rslt &= DListInsert(&points[b].adjacent, points[b].where, a);
+  int rslt = DListInsert(a, b);
+  rslt &= DListInsert(b, a);
   return rslt;
 }
 
@@ -605,6 +632,36 @@ void DocRecord::makePosView(std::string fileName, GFace *gf)
     }
     fprintf(f,"};\n");
   }
+  else
+  {
+    fprintf(f, "View \"scalar\" {\n");
+    for (PointNumero iPoint = 0; iPoint < numPoints; iPoint++){
+      DListPeek p = points[iPoint].adjacent;
+      if (!p){
+        continue;
+      }
+      std::vector<PointNumero> adjacentPoints;
+      do {
+        adjacentPoints.push_back(p->point_num);
+        p = Pred(p);
+      } while (p != points[iPoint].adjacent);
+      adjacentPoints.push_back(p->point_num);
+
+      for (size_t iTriangle = 0; iTriangle < adjacentPoints.size() - 1; iTriangle++) {
+        const PointNumero jPoint = adjacentPoints[iTriangle];
+        const PointNumero kPoint = adjacentPoints[iTriangle + 1];
+        if ((jPoint > iPoint) && (kPoint > iPoint) &&
+          (IsRightOf(iPoint, jPoint, kPoint))) {
+          fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+            points[iPoint].where.h, points[iPoint].where.v, 0.0,
+            points[jPoint].where.h, points[jPoint].where.v, 0.0,
+            points[kPoint].where.h, points[kPoint].where.v, 0.0,
+            (double)iPoint, (double)jPoint, (double)kPoint);
+        }
+      }
+    }
+    fprintf(f, "};\n");
+  }
   fclose(f);
 }
 

Mesh/DivideAndConquer.h

@@ -75,7 +75,7 @@ class DocRecord{
   int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k);
   int Merge(DT vl, DT vr);
   DT RecurTrig(PointNumero left, PointNumero right);
-  int DListInsert(DListRecord **dlist, DPoint center, PointNumero newPoint);
+  int DListInsert(PointNumero centerPoint, PointNumero newPoint);
   int Insert(PointNumero a, PointNumero b);
   int DListDelete(DListPeek *dlist, PointNumero oldPoint);
   int Delete(PointNumero a, PointNumero b);