bug solved, adjacencies more efficient

src/hxt_tetrahedra.c

@@ -24,6 +24,7 @@ Author: Célestin Marot (celestin.marot@uclouvain.be)                        */
 #include <float.h>
 
 #define ABS(x) ((x) >= 0 ? (x) : -(x))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
 
 static const uint64_t NEIGH_H = UINT64_C(0xFFFFFFFFFFFFFFFC);
 
@@ -39,7 +40,7 @@ typedef struct HXT_Delaunay_struct {
   uint64_t num_deleted;
   uint64_t size_deleted;
 
-  uint64_t* buffer;
+  void* buffer;
 } Delaunay_t;
 
 static inline uint64_t walking2cavity(mesh_t* mesh, uint64_t curTet, const uint32_t vta);
@@ -226,7 +227,7 @@ status_t HXT_tetrahedra_compute(mesh_t* mesh){
 
   HXT_CHECK( HXT_malloc(&Del.tmp, Del.size_tmp*sizeof(Del.tmp[0])) );
   HXT_CHECK( HXT_malloc(&Del.deleted, Del.size_deleted*sizeof(uint64_t)) );
-  HXT_CHECK( HXT_malloc(&Del.buffer, 32*(mesh->num_vertices+1)*sizeof(uint64_t)) );
+  HXT_CHECK( HXT_malloc(&Del.buffer, MAX(sizeof(uint64_t)*1024, sizeof(uint32_t)*(mesh->num_vertices+1))) );
 
 
   uint64_t curTet = 0; // we always begin with the first tet. (index 0)
@@ -569,52 +570,53 @@ status_t diggingacavity(mesh_t* mesh, Delaunay_t* Del, uint64_t curTet, const ui
 status_t computeAdjacencies(mesh_t* mesh, Delaunay_t* Del, const uint64_t start, const uint64_t blength){
   if(blength<=60){ // N+2<=32 => N<=30 => 2N<=60
     // for each point p on the surface of the cavity, compute p-vti
-
-    uint64_t* buffer = Del->buffer;
-
-
-    uint64_t i;
-    for (i=0; i<blength; i++){
-      uint32_t* Node = Del->tmp[i].node;
-      buffer[Node[1]+1] = UINT32_MAX;
-      buffer[Node[2]+1] = UINT32_MAX;
-      buffer[Node[3]+1] = UINT32_MAX;
-    }
-
-    uint64_t npts = 0;
-    for (i=0; i<blength; i++)
     {
-      uint32_t* Node = Del->tmp[i].node;
-      if(buffer[Node[1]+1]>npts){
-        buffer[Node[1]+1] = npts++;
-      }
-      Node[1] = buffer[Node[1]+1];
-      if(buffer[Node[2]+1]>npts){
-        buffer[Node[2]+1] = npts++;
-      }
-      Node[2] = buffer[Node[2]+1];
-      if(buffer[Node[3]+1]>npts){
-        buffer[Node[3]+1] = npts++;
-      }
-      Node[3] = buffer[Node[3]+1];
+        uint32_t* bufferu32 = Del->buffer;
+        uint64_t i;
+        for (i=0; i<blength; i++){
+          uint32_t* Node = Del->tmp[i].node;
+          bufferu32[Node[1]+1] = UINT32_MAX;
+          bufferu32[Node[2]+1] = UINT32_MAX;
+          bufferu32[Node[3]+1] = UINT32_MAX;
+        }
 
-      HXT_ASSERT(Node[1]<=mesh->num_vertices);
-      HXT_ASSERT(Node[2]<=mesh->num_vertices);
-      HXT_ASSERT(Node[3]<=mesh->num_vertices);
+        uint64_t npts = 0;
+        for (i=0; i<blength; i++)
+        {
+          uint32_t* Node = Del->tmp[i].node;
+          if(bufferu32[Node[1]+1]>npts){
+            bufferu32[Node[1]+1] = npts++;
+          }
+          Node[1] = bufferu32[Node[1]+1];
+          if(bufferu32[Node[2]+1]>npts){
+            bufferu32[Node[2]+1] = npts++;
+          }
+          Node[2] = bufferu32[Node[2]+1];
+          if(bufferu32[Node[3]+1]>npts){
+            bufferu32[Node[3]+1] = npts++;
+          }
+          Node[3] = bufferu32[Node[3]+1];
+
+          HXT_ASSERT(Node[1]<=mesh->num_vertices);
+          HXT_ASSERT(Node[2]<=mesh->num_vertices);
+          HXT_ASSERT(Node[3]<=mesh->num_vertices);
+        }
     }
 
+    uint64_t* bufferu64 = Del->buffer;
+    uint64_t i;
     for (i=0; i<blength; i++)
     {
-      buffer[Del->tmp[i].node[1]*32 + Del->tmp[i].node[2]] = Del->tmp[i].bnd + 3;
-      buffer[Del->tmp[i].node[2]*32 + Del->tmp[i].node[3]] = Del->tmp[i].bnd + 1;
-      buffer[Del->tmp[i].node[3]*32 + Del->tmp[i].node[1]] = Del->tmp[i].bnd + 2;
+      bufferu64[Del->tmp[i].node[1]*32 + Del->tmp[i].node[2]] = Del->tmp[i].bnd + 3;
+      bufferu64[Del->tmp[i].node[2]*32 + Del->tmp[i].node[3]] = Del->tmp[i].bnd + 1;
+      bufferu64[Del->tmp[i].node[3]*32 + Del->tmp[i].node[1]] = Del->tmp[i].bnd + 2;
     }
 
     for (i=0; i<blength; i++)
     {
-      mesh->tetrahedra.neigh[Del->tmp[i].bnd+ 1] = buffer[Del->tmp[i].node[3]*32 + Del->tmp[i].node[2]];
-      mesh->tetrahedra.neigh[Del->tmp[i].bnd+ 2] = buffer[Del->tmp[i].node[1]*32 + Del->tmp[i].node[3]];
-      mesh->tetrahedra.neigh[Del->tmp[i].bnd+ 3] = buffer[Del->tmp[i].node[2]*32 + Del->tmp[i].node[1]];
+      mesh->tetrahedra.neigh[Del->tmp[i].bnd+ 1] = bufferu64[Del->tmp[i].node[3]*32 + Del->tmp[i].node[2]];
+      mesh->tetrahedra.neigh[Del->tmp[i].bnd+ 2] = bufferu64[Del->tmp[i].node[1]*32 + Del->tmp[i].node[3]];
+      mesh->tetrahedra.neigh[Del->tmp[i].bnd+ 3] = bufferu64[Del->tmp[i].node[2]*32 + Del->tmp[i].node[1]];
     }
   }
   else{