diff --git a/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c b/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
index c31e3eee3f8013efc3eeb0e626b083183874f77c..30067c1902599722d710cd512290359bcb14241d 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetDelaunay.c
@@ -206,7 +206,7 @@ unsigned computePasses(uint32_t passes[12],
* initialisation of the TetLocal structure
******************************************/
static inline HXTStatus localInit(TetLocal* local){
- local->ball.size = 1020;
+ local->ball.size = 1024;
local->ball.num = 0;
HXT_CHECK( hxtAlignedMalloc(&local->ball.array, sizeof(cavityBnd_t)*local->ball.size));
@@ -770,56 +770,132 @@ HXT_ASSERT(((size_t) verticesID)%SIMD_ALIGN==0);
/**************************************************************
* compute adjacencies by sorting => O(n log n)
**************************************************************/
+
+static inline uint64_t hash64(uint64_t x) {
+ x = (x ^ (x >> 30)) * UINT64_C(0xbf58476d1ce4e5b9);
+ x = (x ^ (x >> 27)) * UINT64_C(0x94d049bb133111eb);
+ return x ^ (x >> 31);
+}
+
+// key must be unique, and it cannot be 0
+static inline uint64_t hashTableSearch(HXTGroup2* pairs, uint64_t mask, uint64_t key, uint64_t neigh) {
+ uint64_t hash = hash64(key) & mask;
+ while(pairs[hash].v[0] != key && pairs[hash].v[0] != 0) {
+ hash = (hash + 1) & mask;
+ }
+ if(pairs[hash].v[0] == key) {
+ return pairs[hash].v[1];
+ }
+
+ pairs[hash].v[0] = key;
+ pairs[hash].v[1] = neigh;
+ return HXT_NO_ADJACENT;
+}
+
+
static inline HXTStatus computeAdjacenciesSlow(HXTMesh* mesh, TetLocal* local, const uint64_t start, const uint64_t blength)
{
- HXTGroup2* edges;
-
- HXT_CHECK( askForBall(local, 2 * blength) );
- edges = (HXTGroup2*) local->ball.array; // recycle array
+ uint64_t size = UINT64_C(1) << u64_log2(5*blength);// 3*blength/2 entries => we want a hash table of at least 3*blength for a good load factor
+ uint64_t mask = size - 1;
+
+ // recycle memory
+ if(sizeof(cavityBnd_t) * local->ball.size < sizeof(HXTGroup2) * size) {
+ uint64_t capacity = (sizeof(HXTGroup2) * size + sizeof(cavityBnd_t) - 1) / sizeof(cavityBnd_t);
+ capacity = MAX(2 * local->ball.size, capacity);
+ HXT_CHECK( hxtAlignedFree(&local->ball.array) );
+ HXT_CHECK( hxtAlignedMalloc(&local->ball.array, sizeof(cavityBnd_t) * capacity) );
+ local->ball.size = capacity;
+ }
+ HXTGroup2* pairs = (HXTGroup2*) local->ball.array;
+ memset(pairs, 0, size * sizeof(HXTGroup2));
const uint8_t index[4] = {2,3,1,2};
- const uint32_t n = mesh->vertices.num + 1;
- uint64_t maxKey = 0;
+#ifdef DEBUG
+ int found = 0;
+#endif
for (uint64_t i=0; i<blength; i++)
{
- uint64_t curTet = local->deleted.array[start+ i];
+ uint64_t curTet = local->deleted.array[start + i];
const uint32_t* __restrict__ Node = mesh->tetrahedra.node + 4*curTet;
for (unsigned j=0; j<3; j++)
{
- uint32_t n0 = Node[index[j]] + 1;
- uint32_t n1 = Node[index[j+1]] + 1;
-
- if(n0 < n1)
- edges[3*i + j].v[0] = ((uint64_t) n0 * n) + n1;
- else
- edges[3*i + j].v[0] = ((uint64_t) n1 * n) + n0;
-
- if(edges[3*i + j].v[0] > maxKey)
- maxKey = edges[3*i + j].v[0];
-
- edges[3*i + j].v[1] = 4*curTet+j+1;
+ uint32_t n0 = Node[index[j]];
+ uint32_t n1 = Node[index[j+1]];
+ uint64_t key = (uint64_t) (n0 ^ n1) << 32 | MAX(n0, n1);
+
+ uint64_t neigh = hashTableSearch(pairs, mask, key, 4 * curTet + j + 1);
+ if(neigh != HXT_NO_ADJACENT) {
+#ifdef DEBUG
+ found++;
+#endif
+ mesh->tetrahedra.neigh[4 * curTet + j + 1] = neigh;
+ mesh->tetrahedra.neigh[neigh] = 4 * curTet + j + 1;
+ }
}
}
- HXT_CHECK( group2_sort_v0(edges, 3*blength, maxKey) );
-
- for(int i=0; i<3*blength; i+=2) {
-#ifndef NDEBUG
- if(edges[i].v[0] != edges[i+1].v[0])
- return HXT_ERROR_MSG(HXT_STATUS_ERROR, "Failed to compute adjacencies (b)");
+#ifdef DEBUG
+ if(2 * found != 3 * blength)
+ return HXT_ERROR_MSG(HXT_STATUS_ERROR, "Failed to compute adjacencies (b)");
#endif
- uint64_t f0 = edges[i].v[1];
- uint64_t f1 = edges[i+1].v[1];
- mesh->tetrahedra.neigh[f0] = f1;
- mesh->tetrahedra.neigh[f1] = f0;
- }
return HXT_STATUS_OK;
}
+
+// static inline HXTStatus computeAdjacenciesSlow(HXTMesh* mesh, TetLocal* local, const uint64_t start, const uint64_t blength)
+// {
+// HXTGroup2* edges;
+
+// HXT_CHECK( askForBall(local, 2 * blength) );
+// edges = (HXTGroup2*) local->ball.array; // recycle array
+
+// const uint8_t index[4] = {2,3,1,2};
+
+// const uint32_t n = mesh->vertices.num + 1;
+// uint64_t maxKey = 0;
+
+// for (uint64_t i=0; i<blength; i++)
+// {
+// uint64_t curTet = local->deleted.array[start+ i];
+// const uint32_t* __restrict__ Node = mesh->tetrahedra.node + 4*curTet;
+// for (unsigned j=0; j<3; j++)
+// {
+// uint32_t n0 = Node[index[j]] + 1;
+// uint32_t n1 = Node[index[j+1]] + 1;
+
+// if(n0 < n1)
+// edges[3*i + j].v[0] = ((uint64_t) n0 * n) + n1;
+// else
+// edges[3*i + j].v[0] = ((uint64_t) n1 * n) + n0;
+
+// if(edges[3*i + j].v[0] > maxKey)
+// maxKey = edges[3*i + j].v[0];
+
+// edges[3*i + j].v[1] = 4*curTet+j+1;
+// }
+// }
+
+// HXT_CHECK( group2_sort_v0(edges, 3*blength, maxKey) );
+
+// for(int i=0; i<3*blength; i+=2) {
+// #ifndef NDEBUG
+// if(edges[i].v[0] != edges[i+1].v[0])
+// return HXT_ERROR_MSG(HXT_STATUS_ERROR, "Failed to compute adjacencies (b)");
+// #endif
+// uint64_t f0 = edges[i].v[1];
+// uint64_t f1 = edges[i+1].v[1];
+
+// mesh->tetrahedra.neigh[f0] = f1;
+// mesh->tetrahedra.neigh[f1] = f0;
+// }
+// return HXT_STATUS_OK;
+// }
+
+
/****************************************
* filling back the cavity (DelaunayBall)
****************************************/
diff --git a/contrib/hxt/tetMesh/src/hxt_tetDelaunayReshape.c b/contrib/hxt/tetMesh/src/hxt_tetDelaunayReshape.c
index bd43205e44bd1d75ba51299877b4d11950e6d4f6..87d13c029d0dbbed2d61a938bcfb8f5ec3b30a0a 100644
--- a/contrib/hxt/tetMesh/src/hxt_tetDelaunayReshape.c
+++ b/contrib/hxt/tetMesh/src/hxt_tetDelaunayReshape.c
@@ -287,118 +287,118 @@ HXTStatus hashTableDestroy(HXTHashTable* map) {
// remove the tetrahedra adjacent to the face that does not see the point, progressively, until the cavity is star shaped...
HXTStatus reshapeCavityIfNeeded(TetLocal* local, HXTMesh* mesh, const uint32_t vta, uint64_t prevNumDeleted)
{
- uint64_t count = 0;
+// uint64_t count = 0;
- // STEP 1: for each ball faces, we need the index of the tet in the deleted array, instead of the index of the tet in the mesh
- // for each deleted tet, we need the index of the 4 faces in the ball, or the index of the neighbor in the deleted array
- // if the facet is not a facet of the cavity. To differentiate both, BND_FACE_FLAG is set for indices of the ball.
- uint64_t numTet = local->deleted.num - prevNumDeleted;
- uint64_t* tets = &local->deleted.array[prevNumDeleted];
+// // STEP 1: for each ball faces, we need the index of the tet in the deleted array, instead of the index of the tet in the mesh
+// // for each deleted tet, we need the index of the 4 faces in the ball, or the index of the neighbor in the deleted array
+// // if the facet is not a facet of the cavity. To differentiate both, BND_FACE_FLAG is set for indices of the ball.
+// uint64_t numTet = local->deleted.num - prevNumDeleted;
+// uint64_t* tets = &local->deleted.array[prevNumDeleted];
- HXTHashTable map;
- HXT_CHECK( hashTableInit(&map, numTet) );
+// HXTHashTable map;
+// HXT_CHECK( hashTableInit(&map, numTet) );
- for(uint64_t i=0; i<numTet; i++) {
- hashTablePut(&map, tets[i], i);
- }
+// for(uint64_t i=0; i<numTet; i++) {
+// hashTablePut(&map, tets[i], i);
+// }
- uint64_t* faces;
- HXT_CHECK( hxtMalloc(&faces, sizeof(uint64_t) * 4 * numTet) );
+// uint64_t* faces;
+// HXT_CHECK( hxtMalloc(&faces, sizeof(uint64_t) * 4 * numTet) );
- // just copy the adjacency into faces
- uint64_t curFace = 0;
- for(uint64_t i=0; i<numTet; i++) {
- uint64_t curTet = tets[i];
- for(uint64_t j=0; j<4; j++) {
+// // just copy the adjacency into faces
+// uint64_t curFace = 0;
+// for(uint64_t i=0; i<numTet; i++) {
+// uint64_t curTet = tets[i];
+// for(uint64_t j=0; j<4; j++) {
- uint64_t neigh = mesh->tetrahedra.neigh[4 * curTet + j];
- if(neigh != HXT_NO_ADJACENT && getDeletedFlag(mesh, neigh / 4)) {
- faces[4 * i + j] = 4 * hashTableGet(&map, neigh / 4) + neigh % 4;
- }
- else { // it should be the next one in the ball, (we didn't change the order since @diggingACavity)
- HXT_ASSERT(local->ball.array[curFace].neigh==neigh);
- faces[4 * i + j] = curFace;
- curFace++;
- }
- }
- }
+// uint64_t neigh = mesh->tetrahedra.neigh[4 * curTet + j];
+// if(neigh != HXT_NO_ADJACENT && getDeletedFlag(mesh, neigh / 4)) {
+// faces[4 * i + j] = 4 * hashTableGet(&map, neigh / 4) + neigh % 4;
+// }
+// else { // it should be the next one in the ball, (we didn't change the order since @diggingACavity)
+// HXT_ASSERT(local->ball.array[curFace].neigh==neigh);
+// faces[4 * i + j] = curFace;
+// curFace++;
+// }
+// }
+// }
- HXT_CHECK( hashTableDestroy(&map) );
+// HXT_CHECK( hashTableDestroy(&map) );
-#ifndef NDEBUG
- if(curFace != local->ball.num) {
- return HXT_ERROR_MSG(HXT_STATUS_ERROR, "DEBUG: Faces did not appear in the same order as deleted tet");
- }
-#endif
+// #ifndef NDEBUG
+// if(curFace != local->ball.num) {
+// return HXT_ERROR_MSG(HXT_STATUS_ERROR, "DEBUG: Faces did not appear in the same order as deleted tet");
+// }
+// #endif
- // double *vtaCoord = &mesh->vertices.coord[4 * vta];
+// // double *vtaCoord = &mesh->vertices.coord[4 * vta];
- // // STEP 2: for each face that does not see vta, delete the associated tet, and update the boundary of the cavity
- // for (uint64_t curFace=0; curFace<local->ball.num; curFace++) {
- // if(local->ball.array[curFace].node[2]==HXT_GHOST_VERTEX){
- // continue;
- // }
+// // // STEP 2: for each face that does not see vta, delete the associated tet, and update the boundary of the cavity
+// // for (uint64_t curFace=0; curFace<local->ball.num; curFace++) {
+// // if(local->ball.array[curFace].node[2]==HXT_GHOST_VERTEX){
+// // continue;
+// // }
- // double* b = &mesh->vertices.coord[4 * local->ball.array[curFace].node[0]];
- // double* c = &mesh->vertices.coord[4 * local->ball.array[curFace].node[1]];
- // double* d = &mesh->vertices.coord[4 * local->ball.array[curFace].node[2]];
- // if(orient3d(vtaCoord, b, c, d)<0.0){
- // continue;
- // }
+// // double* b = &mesh->vertices.coord[4 * local->ball.array[curFace].node[0]];
+// // double* c = &mesh->vertices.coord[4 * local->ball.array[curFace].node[1]];
+// // double* d = &mesh->vertices.coord[4 * local->ball.array[curFace].node[2]];
+// // if(orient3d(vtaCoord, b, c, d)<0.0){
+// // continue;
+// // }
- // // 2.1 update the facets so that no facet points to the current tet... does not work with our structure. SHITTT !!
- // // => we are forced to use an O(n^2) algo if we do not want to use a super large table
+// // // 2.1 update the facets so that no facet points to the current tet... does not work with our structure. SHITTT !!
+// // // => we are forced to use an O(n^2) algo if we do not want to use a super large table
- // // 2.2: undelete the tetrahedron
- // uint64_t index = local->ball.array[curFace].neigh;
- // uint64_t tetToUndelete = local->deleted.array[index / 4];
- // unsigned facet = index % 4;
- // unsetDeletedFlag(mesh, tetToUndelete);
+// // // 2.2: undelete the tetrahedron
+// // uint64_t index = local->ball.array[curFace].neigh;
+// // uint64_t tetToUndelete = local->deleted.array[index / 4];
+// // unsigned facet = index % 4;
+// // unsetDeletedFlag(mesh, tetToUndelete);
- // // second part: remove the faces of the tet that were on the boundary
- // // do not forget to update curFace to point just before the next face to process
- // uint64_t bndFaces[4] = {HXT_NO_ADJACENT, HXT_NO_ADJACENT, HXT_NO_ADJACENT, HXT_NO_ADJACENT};
- // int nbndFace = 0;
+// // // second part: remove the faces of the tet that were on the boundary
+// // // do not forget to update curFace to point just before the next face to process
+// // uint64_t bndFaces[4] = {HXT_NO_ADJACENT, HXT_NO_ADJACENT, HXT_NO_ADJACENT, HXT_NO_ADJACENT};
+// // int nbndFace = 0;
- // // the way that the boundary faces are contructed,
- // // TODO: use a bidirectional structure to be able to find the faces of deleted tet.
+// // // the way that the boundary faces are contructed,
+// // // TODO: use a bidirectional structure to be able to find the faces of deleted tet.
- // // face to delete: curface
+// // // face to delete: curface
- // for (uint64_t i=local->ball.num; nbndFace<4 && i>0; i--) {
- // if(mesh->tetrahedra.neigh[local->ball.array[i-1].neigh]/4==tetToUndelete) {
- // bndFaces[nbndFace++] = local->ball.array[i-1].neigh;
- // local->ball.num--;
- // local->ball.array[i-1] = local->ball.array[local->ball.num];
- // }
- // }
+// // for (uint64_t i=local->ball.num; nbndFace<4 && i>0; i--) {
+// // if(mesh->tetrahedra.neigh[local->ball.array[i-1].neigh]/4==tetToUndelete) {
+// // bndFaces[nbndFace++] = local->ball.array[i-1].neigh;
+// // local->ball.num--;
+// // local->ball.array[i-1] = local->ball.array[local->ball.num];
+// // }
+// // }
- // // third part: add new facets to the boundary
+// // // third part: add new facets to the boundary
- // count++;
- // }
+// // count++;
+// // }
- // // if(count)
- // // printf("nbr of tet to undelete = %lu, nbr of deleted tet = %lu\n", count, local->deleted.num);
+// // // if(count)
+// // // printf("nbr of tet to undelete = %lu, nbr of deleted tet = %lu\n", count, local->deleted.num);
- // // // STEP 3: compute all the flags at once ! ??
+// // // // STEP 3: compute all the flags at once ! ??
- // STEP 4: put back the right neighbor value inside the ball struct
- for(uint64_t i=0; i<local->ball.num; i++) {
- uint64_t index = local->ball.array[i].neigh;
- uint64_t interiorTet = tets[index / 4];
- unsigned facet = index % 4;
- local->ball.array[i].neigh = mesh->tetrahedra.neigh[4 * interiorTet + facet];
- }
+// // STEP 4: put back the right neighbor value inside the ball struct
+// for(uint64_t i=0; i<local->ball.num; i++) {
+// uint64_t index = local->ball.array[i].neigh;
+// uint64_t interiorTet = tets[index / 4];
+// unsigned facet = index % 4;
+// local->ball.array[i].neigh = mesh->tetrahedra.neigh[4 * interiorTet + facet];
+// }
- HXT_CHECK( hxtFree(&faces) );
+// HXT_CHECK( hxtFree(&faces) );
uint64_t blindFace = 0;