avoid nasty infinite loop in waling2cavity!

contrib/hxt/README.md

@@ -140,20 +140,18 @@ set_target_properties(my_target hxt_core hxt_predicates hxt_tetBR hxt_tetMesh
 
 ### TODO/future work
 
-core + others:
-
- - extend HXT to 40+ bits nodes, 32 bits colors
-
 tetMesh:
 
- - compute Moore indices on the fly, only when few point to insert or few tetrahedra to optimize
- - try a mesh improvement parallelization based on Geogram point-associated locks
+ - optimize `respectEdgeConstraint` in the same way we did for `reshapeCavityIfNeeded`
+ - optimize the first tet optimization to avoid being stuck if all points are on the same plane
+ - compute Moore indices on the fly for points already in the mesh (find heuristic for when the old tactic should be employed)
+ - try a GSC parallelization based on Geogram point-associated locks
  - add edge-collapse and point insertion
  - modify surface mesh up to a certain distance allowed
  - add a pass of smoothing aimed solely at mesh size control
  - track the order of allocated blocks in reproducible Delaunay
- - simplify the check between Hxt colors and Gmsh colors
  - code the boundary recovery, to replace the tetBR module
+ - extend HXT to 40+ bits nodes
 
 ## References
 

contrib/hxt/tetMesh/src/hxt_tetDelaunay.c

@@ -517,7 +517,6 @@ static inline HXTStatus walking2Cavity(HXTMesh* mesh, HXTPartition* partition, u
 
     // we should pass in orient3d mode here :p
     unsigned index = 4;
-    unsigned outside = 0;
     unsigned randomU = hxtReproducibleLCG(&seed);
     for (unsigned j=0; j<4; j++)
     {
@@ -535,7 +534,6 @@ static inline HXTStatus walking2Cavity(HXTMesh* mesh, HXTPartition* partition, u
           //                        vtaCoord[0], vtaCoord[1], vtaCoord[2]);
           // }
 
-          outside = 1;
           uint32_t node = mesh->tetrahedra.node[curNeigh[i]];
 
           if(node==HXT_GHOST_VERTEX) {
@@ -543,23 +541,21 @@ static inline HXTStatus walking2Cavity(HXTMesh* mesh, HXTPartition* partition, u
             return HXT_STATUS_OK;
           }
 
-          if(!vertexOutOfPartition(vertices, node, rel, partition->startDist)){
+          if(vertexOutOfPartition(vertices, node, rel, partition->startDist))
+            return HXT_STATUS_CONFLICT;
+
           index=i;
           break;
         }
       }
     }
-    }
 
     if(index==4){
       const double* __restrict__ a = vertices[curNode[0]].coord;
       const double* __restrict__ b = vertices[curNode[1]].coord;
       const double* __restrict__ c = vertices[curNode[2]].coord;
       const double* __restrict__ d = vertices[curNode[3]].coord;
-      if(outside) {
-        return HXT_STATUS_CONFLICT;
-      }
-      else if( (orient3d(a,b,c,vtaCoord)>=0) +
+      if( (orient3d(a,b,c,vtaCoord)>=0) +
           (orient3d(a,b,vtaCoord,d)>=0) +
           (orient3d(a,vtaCoord,c,d)>=0) +
           (orient3d(vtaCoord,b,c,d)>=0)>2 ){
@@ -570,6 +566,12 @@ static inline HXTStatus walking2Cavity(HXTMesh* mesh, HXTPartition* partition, u
       return HXT_STATUS_OK;
     }
 
+#ifdef DEBUG
+    nstep++;
+    if(nstep > mesh->tetrahedra.num)
+      return HXT_ERROR_MSG(HXT_STATUS_ERROR, "walk stuck in an infinite loop");
+#endif
+
     enteringFace = curNeigh[index]&3;
     nextTet = curNeigh[index]/4;
   }
@@ -928,6 +930,18 @@ static inline HXTStatus fillingACavity(HXTMesh* mesh, TetLocal* local, unsigned
 
     mesh->tetrahedra.flag[newTet] = local->ball.array[i].flag;
 
+    #ifdef DEBUG
+      if((getFacetConstraint(mesh, neigh / 4, neigh % 4)!=0) ^ (getFacetConstraint(mesh, newTet, 0)!=0))
+        return HXT_ERROR_MSG(HXT_STATUS_ERROR, "constraints are not the same on both side of the new tet");
+      if(getFacetConstraint(mesh, newTet, 1) || getFacetConstraint(mesh, newTet, 2) || getFacetConstraint(mesh, newTet, 3))
+        return HXT_ERROR_MSG(HXT_STATUS_ERROR, "new tetrahedron flag are erroneous (face constraint on new face)");
+      if(getEdgeConstraint(mesh, newTet, getEdgeFromNodes(0, 1)) ||
+         getEdgeConstraint(mesh, newTet, getEdgeFromNodes(0, 2)) ||
+         getEdgeConstraint(mesh, newTet, getEdgeFromNodes(0, 3))
+         )
+        return HXT_ERROR_MSG(HXT_STATUS_ERROR, "new tetrahedron flag are erroneous (edge constraint on new edge)");
+    #endif
+
     // update neighbor's neighbor
     mesh->tetrahedra.neigh[neigh] = 4*newTet;
 
@@ -1450,8 +1464,6 @@ static HXTStatus parallelDelaunay3D(HXTMesh* mesh,
           ******************************************************/
           for (uint32_t i=0; i<localN; i++)
           {
-            // if( i%1000==0)
-            //   printf("%f %%\n", 100.0*i/localN);
             uint32_t passIndex = (localStart+i)%passLength;
             uint32_t vta = nodeInfo[passStart + passIndex].node;
             if(nodeInfo[passStart + passIndex].status==HXT_STATUS_TRYAGAIN){

contrib/hxt/tetMesh/src/hxt_tetDelaunayReshape.c

@@ -387,7 +387,6 @@ HXTStatus reshapeCavityIfNeeded(TetLocal* local, HXTMesh* mesh, const uint32_t v
         uint64_t out = faces[4 * (in / 4) + f];
         faces[out] = local->ball.num;
 
-
         uint64_t* curNeigh = mesh->tetrahedra.neigh + tetToUndelete*4;
         uint32_t* curNode = mesh->tetrahedra.node + tetToUndelete*4;
         if(f==0) {
@@ -442,7 +441,9 @@ HXTStatus reshapeCavityIfNeeded(TetLocal* local, HXTMesh* mesh, const uint32_t v
   }
   local->deleted.num -= shift;
 
-  // HXT_ASSERT( isStarShaped(local, mesh, vta, &blindFace)==HXT_STATUS_OK );
+#ifdef DEBUG
+  HXT_ASSERT( isStarShaped(local, mesh, vta, &blindFace)==HXT_STATUS_OK );
+#endif
 
   return HXT_STATUS_OK;
 }