1 #include "hxt_tetDelaunayReshape.h"
2 #include "hxt_mesh.h"
3 #include "hxt_message.h"
4 #include "hxt_sort.h"
5 #include "hxt_tetFlag.h"
6 #include "predicates.h"
7 #include <stdint.h>
8
9
getToUndeleteFlag(HXTMesh * mesh,uint64_t tet)10 static inline uint16_t getToUndeleteFlag(HXTMesh* mesh, uint64_t tet) {
11 return getUnusedFlag(mesh, tet, 0);
12 }
13
setToUndeleteFlag(HXTMesh * mesh,uint64_t tet)14 static inline void setToUndeleteFlag(HXTMesh* mesh, uint64_t tet) {
15 setUnusedFlag(mesh, tet, 0);
16 }
17
unsetToUndeleteFlag(HXTMesh * mesh,uint64_t tet)18 static inline void unsetToUndeleteFlag(HXTMesh* mesh, uint64_t tet) {
19 unsetUnusedFlag(mesh, tet, 0);
20 }
21
22 /* check if the cavity is star shaped
23 This isn't useful for pure Delaunay but when we constrain cavity with color,
24 it is useful */
isStarShaped(TetLocal * local,HXTMesh * mesh,const uint32_t vta,uint64_t * blindFaceIndex)25 static HXTStatus isStarShaped(TetLocal* local, HXTMesh* mesh, const uint32_t vta, uint64_t* blindFaceIndex)
26 {
27 HXTVertex* vertices = (HXTVertex*) mesh->vertices.coord;
28 double *vtaCoord = vertices[vta].coord;
29
30 for (uint64_t i=0; i<local->ball.num; i++) {
31 if(local->ball.array[i].node[2]!=HXT_GHOST_VERTEX){
32 double* b = vertices[local->ball.array[i].node[0]].coord;
33 double* c = vertices[local->ball.array[i].node[1]].coord;
34 double* d = vertices[local->ball.array[i].node[2]].coord;
35 if(orient3d(vtaCoord, b, c, d)>=0.0){
36 *blindFaceIndex = i;
37 return HXT_STATUS_INTERNAL;
38 }
39 }
40 }
41 return HXT_STATUS_OK;
42 }
43
44
respectEdgeConstraint(TetLocal * local,HXTMesh * mesh,const uint32_t vta,const uint32_t color,const uint64_t prevDeleted,int * undeleteTet)45 HXTStatus respectEdgeConstraint(TetLocal* local, HXTMesh* mesh, const uint32_t vta, const uint32_t color, const uint64_t prevDeleted, int* undeleteTet)
46 {
47 // HXT_WARNING("a constrained edge was inside the cavity, recovering it");
48
49 *undeleteTet = 0;
50 // all the tetrahedron have the same color 'color', we will use that color to flag them
51 for (uint64_t i=prevDeleted; i<local->deleted.num; i++) {
52 uint64_t delTet = local->deleted.array[i];
53 mesh->tetrahedra.color[delTet] = 0;
54 }
55
56 for (uint64_t i=prevDeleted; i<local->deleted.num; i++) {
57 uint64_t delTet = local->deleted.array[i];
58 for (int edge=0; edge<6; edge++) {
59 if(getEdgeConstraint(mesh, delTet, edge) && (mesh->tetrahedra.color[delTet] & (1U<<edge))==0) {
60 unsigned in_facet;
61 unsigned out_facet;
62
63 getFacetsFromEdge(edge, &in_facet, &out_facet);
64
65 int edgeIsSafe = 0;
66 uint64_t curTet = delTet;
67
68 // first turn
69 do
70 {
71 uint32_t newV = mesh->tetrahedra.node[4*curTet + in_facet];
72
73 // go into the neighbor through out_facet
74 uint64_t neigh = mesh->tetrahedra.neigh[4*curTet + out_facet];
75 curTet = neigh/4;
76 in_facet = neigh%4;
77
78 uint32_t* nodes = mesh->tetrahedra.node + 4*curTet;
79 for (out_facet=0; out_facet<3; out_facet++)
80 if(nodes[out_facet]==newV)
81 break;
82
83 if(getDeletedFlag(mesh, curTet) && !getToUndeleteFlag(mesh, curTet)) {
84 // mark that the edge as been treated
85 #ifdef DEBUG
86 if((mesh->tetrahedra.color[curTet] & (1U<<getEdgeFromFacets(in_facet, out_facet)))!=0)
87 return HXT_ERROR_MSG(HXT_STATUS_ERROR, "the flag says that the tet has already been processed for this edge...");
88 #endif
89 mesh->tetrahedra.color[curTet] |= (1U<<getEdgeFromFacets(in_facet, out_facet));
90 }
91 else {
92 edgeIsSafe=1;
93 }
94
95 } while (curTet!=delTet);
96
97 if(!edgeIsSafe) { // we must find a tetrahedron on the opposite side of vta and mark it as "to undelete".
98 getFacetsFromEdge(edge, &in_facet, &out_facet);
99 curTet = delTet;
100
101 uint64_t tetContainingVta = local->deleted.array[prevDeleted];
102 uint64_t tetToUndelete = HXT_NO_ADJACENT;
103 double distMax = 0.0;
104 double* vtaCoord = mesh->vertices.coord + 4*vta;
105
106 #ifdef DEBUG
107 double* a = mesh->vertices.coord + 4*mesh->tetrahedra.node[4*tetContainingVta];
108 double* b = mesh->vertices.coord + 4*mesh->tetrahedra.node[4*tetContainingVta+1];
109 double* c = mesh->vertices.coord + 4*mesh->tetrahedra.node[4*tetContainingVta+2];
110 double* d = mesh->vertices.coord + 4*mesh->tetrahedra.node[4*tetContainingVta+3];
111
112 if(orient3d(vtaCoord,b,c,d)>0.0 || orient3d(a,vtaCoord,c,d)>0.0 || orient3d(a,b,vtaCoord,d)>0.0 || orient3d(a,b,c,vtaCoord)>0.0) {
113 return HXT_ERROR_MSG(HXT_STATUS_ERROR, "an edge part of a ghost tetrahedron is constrained");
114 }
115 #endif
116
117 // second turn
118 do
119 {
120 uint32_t newV = mesh->tetrahedra.node[4*curTet + in_facet];
121
122 // go into the neighbor through out_facet
123 uint64_t neigh = mesh->tetrahedra.neigh[4*curTet + out_facet];
124 curTet = neigh/4;
125 in_facet = neigh%4;
126
127 uint32_t* nodes = mesh->tetrahedra.node + 4*curTet;
128 for (out_facet=0; out_facet<3; out_facet++)
129 if(nodes[out_facet]==newV)
130 break;
131
132 HXT_ASSERT(newV != HXT_GHOST_VERTEX);
133 HXT_ASSERT(nodes[in_facet] != HXT_GHOST_VERTEX);
134 double* coord1 = mesh->vertices.coord + newV;
135 double* coord2 = mesh->vertices.coord + nodes[in_facet];
136
137 if(curTet!=tetContainingVta) {
138 double dist = 0.0;
139 for (int l=0; l<3; l++) {
140 double meanCoord = (coord1[l]+coord2[l])*0.5;
141 double diff = meanCoord-vtaCoord[l];
142 dist += diff*diff;
143 }
144
145 if(dist>distMax) {
146 dist = distMax;
147 tetToUndelete = curTet;
148 }
149 }
150 } while (curTet!=delTet);
151
152 // printf("marking tet %lu as \"to undelete\" in edgeConstraint\n", tetToUndelete);
153 setToUndeleteFlag(mesh, tetToUndelete);
154 *undeleteTet = 1;
155 }
156 }
157 }
158 }
159
160 for (uint64_t i=prevDeleted; i<local->deleted.num; i++) {
161 uint64_t delTet = local->deleted.array[i];
162 mesh->tetrahedra.color[delTet] = color;
163 }
164
165 return HXT_STATUS_OK;
166 }
167
168
169 // search neighbors in an already initialized hash table. `key` must be unique, and cannot be UINT64_MAX
hashTablePut(HXTGroup2 * pairs,uint64_t mask,uint64_t key,uint64_t value)170 static inline void hashTablePut(HXTGroup2* pairs, uint64_t mask, uint64_t key, uint64_t value) {
171 HXT_ASSERT(key != UINT64_MAX);
172
173 uint64_t hash = hash64(key) & mask;
174 while(pairs[hash].v[0] != UINT64_MAX) {
175 HXT_ASSERT(pairs[hash].v[0] != key); // same key cannot already be there
176 hash = (hash + 1) & mask;
177 }
178
179 pairs[hash].v[0] = key;
180 pairs[hash].v[1] = value;
181 }
182
hashTableGet(HXTGroup2 * pairs,uint64_t mask,uint64_t key)183 static inline uint64_t hashTableGet(HXTGroup2* pairs, uint64_t mask, uint64_t key) {
184 uint64_t hash = hash64(key) & mask;
185 while(pairs[hash].v[0] != key && pairs[hash].v[0] != UINT64_MAX) {
186 hash = (hash + 1) & mask;
187 }
188 if(pairs[hash].v[0] == key) {
189 return pairs[hash].v[1];
190 }
191 return HXT_NO_ADJACENT;
192 }
193
194
bndPushFacet(TetLocal * local,HXTMesh * mesh,uint64_t tet,int facet,const uint64_t neigh)195 static inline void bndPushFacet(TetLocal* local, HXTMesh* mesh, uint64_t tet, int facet, const uint64_t neigh) {
196 uint32_t* nodes = mesh->tetrahedra.node + tet*4;
197 if(facet==0) {
198 bndPush(local, (getFacetConstraint(mesh, tet, 0) ) |
199 (getEdgeConstraint(mesh, tet, 1)>>1) | // constraint on edge 1 (facet 0 2) goes on edge 0
200 (getEdgeConstraint(mesh, tet, 0)<<1) | // constraint on edge 0 (facet 0 1) goes on edge 1
201 (getEdgeConstraint(mesh, tet, 2) ), // constraint on edge 2 (facet 0 3) goes on edge 2
202 nodes[2], nodes[1], nodes[3], neigh);
203 }
204 else if(facet==1) {
205 bndPush(local, (getFacetConstraint(mesh, tet, 1)>>1) |// constraint on facet 1 goes on facet 0
206 (getEdgeConstraint(mesh, tet, 0) ) | // constraint on edge 0 (facet 1 0) goes on edge 0
207 (getEdgeConstraint(mesh, tet, 3)>>2) | // constraint on edge 3 (facet 1 2) goes on edge 1
208 (getEdgeConstraint(mesh, tet, 4)>>2), // constraint on edge 4 (facet 1 3) goes on edge 2
209 nodes[0], nodes[2], nodes[3], neigh);
210 }
211 else if(facet==2) {
212 bndPush(local, (getFacetConstraint(mesh, tet, 2)>>2) |// constraint on facet 2 goes on facet 0
213 (getEdgeConstraint(mesh, tet, 3)>>3) | // constraint on edge 3 (facet 2 1) goes on edge 0
214 (getEdgeConstraint(mesh, tet, 1) ) | // constraint on edge 1 (facet 2 0) goes on edge 1
215 (getEdgeConstraint(mesh, tet, 5)>>3), // constraint on edge 5 (facet 2 3) goes on edge 2
216 nodes[1], nodes[0], nodes[3], neigh);
217 }
218 else {
219 bndPush(local, (getFacetConstraint(mesh, tet, 3)>>3) |// constraint on facet 3 goes on facet 0
220 (getEdgeConstraint(mesh, tet, 2)>>2) | // constraint on edge 2 (facet 3 0) goes on edge 0
221 (getEdgeConstraint(mesh, tet, 4)>>3) | // constraint on edge 4 (facet 3 1) goes on edge 1
222 (getEdgeConstraint(mesh, tet, 5)>>3), // constraint on edge 5 (facet 3 2) goes on edge 2
223 nodes[0], nodes[1], nodes[2], neigh);
224 }
225 }
226
227
228
229
230 // remove the tetrahedra adjacent to the face that does not see the point, progressively, until the cavity is star shaped...
reshapeCavityIfNeeded(TetLocal * local,HXTMesh * mesh,const uint32_t vta,const uint64_t prevNumDeleted,int undeleteTet)231 HXTStatus reshapeCavityIfNeeded(TetLocal* local, HXTMesh* mesh, const uint32_t vta, const uint64_t prevNumDeleted, int undeleteTet)
232 {
233
234 uint64_t blindFace;
235 if(undeleteTet) {
236 blindFace = 0;
237 }
238 else {
239 if(isStarShaped(local, mesh, vta, &blindFace)==HXT_STATUS_OK)
240 return HXT_STATUS_OK;
241 }
242
243
244 // we need an array index for all faces that tells if the face is in the boundary or not, it should correspond to the deleted tet array
245 // In he same array, we can store adjacencies between deleted tet !
246 // => compute those adjacencies with a hash table.
247
248 // 1. put all the deleted tet in a hash table : key = meshIndex, value = deletedIndex
249 // 2. create the faces array
250 // 3. iterate over all deleted tet, and fill the faces array using the hash table. A special bit is set for non deleted neighbors
251 // 4 & 5. reshape the cavity
252
253 // point 1 and 2 can take a lot of memory... can we recycle some ??
254
255 const uint64_t numTet = local->deleted.num - prevNumDeleted;
256 uint64_t* tets = &local->deleted.array[prevNumDeleted];
257
258 uint64_t size = UINT64_C(1) << u64_log2(5*numTet/2);// we want a hash table of at least 2.5*numTet for a good load factor
259 uint64_t mask = size - 1;
260
261 uint64_t* faces;
262 HXTGroup2* pairs;
263
264 HXT_CHECK( hxtMalloc(&pairs, sizeof(HXTGroup2) * size));
265 memset(pairs, -1, size * sizeof(HXTGroup2));
266
267 // 1.
268 for(uint64_t i=0; i<numTet; i++) {
269 hashTablePut(pairs, mask, tets[i], i);
270 }
271
272 // 2.
273 HXT_CHECK( hxtMalloc(&faces, sizeof(uint64_t) * 4 * numTet) );
274
275 // 3.
276 uint64_t curFace = 0;
277 for(uint64_t i=0; i<numTet; i++) {
278 uint64_t curTet = tets[i];
279
280 for(int f=0; f<4; f++) {
281 // search in the table, but without replacing...
282 uint64_t neigh = mesh->tetrahedra.neigh[4 * curTet + f];
283
284 if(getFacetConstraint(mesh, curTet, f) || !getDeletedFlag(mesh, neigh / 4) ) { // it should be the next one in the ball, (we didn't change the order since @diggingACavity)
285 HXT_ASSERT(local->ball.array[curFace].neigh==neigh);
286 local->ball.array[curFace].neigh = 4 * i + f; // recreate the ball array instead
287 faces[4 * i + f] = curFace;
288 curFace++;
289 }
290 else {
291 uint64_t deletedIndex = hashTableGet(pairs, mask, neigh / 4);
292 HXT_ASSERT(deletedIndex != HXT_NO_ADJACENT);
293
294 faces[4 * i + f] = 4 * deletedIndex + neigh % 4;
295 }
296
297 }
298 }
299
300 #ifndef NDEBUG
301 if(curFace != local->ball.num) {
302 return HXT_ERROR_MSG(HXT_STATUS_ERROR, "DEBUG: Faces did not appear in the same order as deleted tet");
303 }
304 #endif
305
306 HXT_CHECK( hxtFree(&pairs) );
307
308 /* 4. Undelete tet that were marked as "tet to undelete" in respectEdgeConstraint.
309 * We have got to update the faces array and the boundary correctly
310 */
311 if(undeleteTet) {
312 for(uint64_t i=0; i<numTet; i++) {
313 uint64_t tetToUndelete = tets[i];
314 if(!getToUndeleteFlag(mesh, tetToUndelete))
315 continue;
316
317 int isBoundary = 0;
318 for(int f=0; f<4; f++) {
319 uint64_t neigh = mesh->tetrahedra.neigh[4 * tetToUndelete + f];
320 if(getFacetConstraint(mesh, tetToUndelete, f) || !getDeletedFlag(mesh, neigh / 4) ) { // we will take car of that after
321 isBoundary = 1;
322 break;
323 }
324 }
325
326 if(isBoundary)
327 continue;
328
329 // 4.1: undelete the tetrahedron
330 // unsigned facet = in % 4;
331 unsetToUndeleteFlag(mesh, tetToUndelete);
332 unsetDeletedFlag(mesh, tetToUndelete);
333 tets[i] = HXT_NO_ADJACENT;// correct deletion from the array is handled at point 6.
334
335 askForBall(local, 3);
336 for(int f=0; f<4; f++) {
337 uint64_t neigh = mesh->tetrahedra.neigh[4 * tetToUndelete + f];
338
339 uint64_t out = faces[4 * i + f];
340 faces[out] = local->ball.num;
341
342 bndPushFacet(local, mesh, tetToUndelete, f, out);
343 }
344 }
345 }
346
347
348 /* 5. iterate over all faces, and undelete (remove from the cavity) the associated tet when the face does not see vta
349 * We have got to update the faces array and the boundary correctly
350 * We also have to update the curFace index such that all faces below it are good
351 * and all faces above or equal are unchecked.
352 */
353 double *vtaCoord = &mesh->vertices.coord[4 * vta];
354 curFace = blindFace;
355 while (curFace<local->ball.num) {
356 uint64_t in = local->ball.array[curFace].neigh;
357 uint64_t tetToUndelete = tets[in / 4];
358
359 if(!getToUndeleteFlag(mesh, tetToUndelete)) {
360 if(local->ball.array[curFace].node[2]==HXT_GHOST_VERTEX){
361 curFace++;
362 continue;
363 }
364 double* b = &mesh->vertices.coord[4 * local->ball.array[curFace].node[0]];
365 double* c = &mesh->vertices.coord[4 * local->ball.array[curFace].node[1]];
366 double* d = &mesh->vertices.coord[4 * local->ball.array[curFace].node[2]];
367 if(orient3d(vtaCoord, b, c, d)<0.0){
368 curFace++;
369 continue;
370 }
371 }
372 else {
373 unsetToUndeleteFlag(mesh, tetToUndelete);
374 }
375
376
377 // 5.1: undelete the tetrahedron
378 // unsigned facet = in % 4;
379 unsetDeletedFlag(mesh, tetToUndelete);
380 tets[in / 4] = HXT_NO_ADJACENT;// correct deletion from the array is handled at point 6.
381
382 askForBall(local, 3);
383 for(int f=0; f<4; f++) {
384 uint64_t neigh = mesh->tetrahedra.neigh[4 * tetToUndelete + f];
385 if(getFacetConstraint(mesh, tetToUndelete, f) || !getDeletedFlag(mesh, neigh / 4) ) { // it is an exterior facet, so we have to remove it from the ball !
386 uint64_t face = faces[4 * (in / 4) + f];
387 if(face < curFace) {
388 curFace--;
389 if(face != curFace) {
390 local->ball.array[face] = local->ball.array[curFace];
391 uint64_t inOther = local->ball.array[face].neigh;
392 HXT_ASSERT(faces[inOther] == curFace);
393 faces[inOther] = face;
394 }
395 face = curFace;
396 }
397
398 // put the last one in here, and change indices in the faces table too !
399 local->ball.num--;
400 if(face != local->ball.num) {
401 local->ball.array[face] = local->ball.array[local->ball.num];
402 uint64_t inOther = local->ball.array[face].neigh;
403 HXT_ASSERT(faces[inOther] == local->ball.num);
404 faces[inOther] = face;
405 }
406 }
407 else { // we have to add a boundary facet at the end
408 uint64_t out = faces[4 * (in / 4) + f];
409 faces[out] = local->ball.num;
410
411 bndPushFacet(local, mesh, tetToUndelete, f, out);
412 }
413 }
414 }
415
416 HXT_CHECK( hxtFree(&faces) );
417
418
419 // 6. put back the correct values in ball neighbors
420 for(uint64_t curFace=0; curFace<local->ball.num; curFace++) {
421 uint64_t in = local->ball.array[curFace].neigh;
422 local->ball.array[curFace].neigh = mesh->tetrahedra.neigh[4*tets[in/4] + in%4];
423 }
424
425
426 // 7. correctly delete the tetrahedra in the deleted array
427 uint64_t shift = 0;
428 for(uint64_t i=0; i<numTet; i++) {
429 if(tets[i]==HXT_NO_ADJACENT)
430 shift++;
431 else if(shift)
432 tets[i-shift] = tets[i];
433 }
434 local->deleted.num -= shift;
435
436 #ifdef DEBUG
437 HXT_ASSERT( isStarShaped(local, mesh, vta, &blindFace)==HXT_STATUS_OK );
438 #endif
439
440 return HXT_STATUS_OK;
441 }
442