1 // Hxt - Copyright (C)
2 // 2016 - 2020 UCLouvain
3 //
4 // See the LICENSE.txt file for license information.
5 //
6 // Contributor(s):
7 // Célestin Marot
8
9 #include "hxt_sort.h"
10 #include "hxt_tetFlag.h"
11
12
13 /* just a function such that:
14 * if hashLess(a,b)==true
15 * then hashLess(b,a)==false
16 *
17 * BUT hashLess(a,b)==true and hashLess(b,c)==true
18 * doesn't imply that hashLess(a,c)==true
19 * and inversely if you replace true with false (no transitivity) */
hashLess(uint64_t a,uint64_t b)20 static inline uint64_t hashLess(uint64_t a, uint64_t b) {
21 return ((a^b)&1)^(a<b);
22 }
23
24 /* good hash function for non-cryptographic purpose
25 * see https://stackoverflow.com/a/12996028 */
hash64(uint64_t x)26 static inline uint64_t hash64(uint64_t x) {
27 x = (x ^ (x >> 30)) * UINT64_C(0xbf58476d1ce4e5b9);
28 x = (x ^ (x >> 27)) * UINT64_C(0x94d049bb133111eb);
29 x = x ^ (x >> 31);
30 return x;
31 }
32
33 /* just a function such that:
34 * if transitiveHashLess(a,b)==true
35 * then transitiveHashLess(b,a)==false
36 *
37 * AND transitiveHashLess(a,b)==true and transitiveHashLess(b,c)==true
38 * imply that transitiveHashLess(a,c)==true
39 * and inversely if you replace true with false (the relation is transitive) */
transitiveHashLess(uint64_t a,uint64_t b)40 static inline uint64_t transitiveHashLess(uint64_t a, uint64_t b) {
41 return hash64(a) < hash64(b);
42 }
43
44
45 #ifndef NDEBUG
countBoundaries(HXTMesh * mesh)46 static inline uint64_t countBoundaries(HXTMesh* mesh) {
47 uint64_t nBoundaries=0;
48
49 #pragma omp parallel for reduction(+:nBoundaries)
50 for (uint64_t i=0; i<4*mesh->tetrahedra.num; i++) {
51 if(mesh->tetrahedra.neigh[i]==HXT_NO_ADJACENT)
52 nBoundaries++;
53 }
54
55 return nBoundaries;
56 }
57 #endif
58
59
60 /**************************************************************************************
61 * Lines --> Triangles MAPPING (for each line, get 3*tri+e)
62 *************************************************************************************/
hxtGetLines2TriMap(HXTMesh * mesh,uint64_t * lines2TriMap,uint64_t * missing)63 HXTStatus hxtGetLines2TriMap(HXTMesh* mesh, uint64_t* lines2TriMap, uint64_t* missing)
64 {
65 HXTGroup2* edgeKey = NULL;
66 uint64_t* numEdges = NULL;
67
68 const int maxThreads = omp_get_max_threads();
69 const uint64_t n = mesh->vertices.num;
70 const uint64_t numEdgesTotal = mesh->triangles.num*3+mesh->lines.num;
71
72 HXT_CHECK( hxtMalloc(&numEdges, maxThreads*sizeof(uint64_t)) );
73 HXT_CHECK( hxtAlignedMalloc(&edgeKey, numEdgesTotal*sizeof(HXTGroup2)) );
74
75 #pragma omp parallel
76 {
77 #pragma omp for nowait
78 for (uint64_t i=0; i<mesh->lines.num; i++) {
79 uint32_t v0 = mesh->lines.node[2*i];
80 uint32_t v1 = mesh->lines.node[2*i+1];
81
82 if(v0<v1) {
83 edgeKey[i].v[0] = v0*n + v1;
84 edgeKey[i].v[1] = 2*i;
85 }
86 else if(v0>v1){
87 edgeKey[i].v[0] = v1*n + v0;
88 edgeKey[i].v[1] = 2*i;
89 }
90 else {
91 edgeKey[i].v[0] = v0*n + v0; // the line begins and ends at the same point...
92 edgeKey[i].v[1] = 1;
93 lines2TriMap[i] = HXT_NO_ADJACENT;
94 }
95 }
96
97 #pragma omp for
98 for (uint64_t i=0; i<mesh->triangles.num; i++) {
99 uint32_t v[3] = {mesh->triangles.node[3*i+0],
100 mesh->triangles.node[3*i+1],
101 mesh->triangles.node[3*i+2]};
102
103 HXTSORT_3_VALUES_INPLACE(uint32_t, v);
104
105 edgeKey[mesh->lines.num+3*i].v[0] = v[0]*n + v[1];
106 edgeKey[mesh->lines.num+3*i].v[1] = 2*(3*i)+1;
107 edgeKey[mesh->lines.num+3*i+1].v[0] = v[0]*n + v[2];
108 edgeKey[mesh->lines.num+3*i+1].v[1] = 2*(3*i+1)+1;
109 edgeKey[mesh->lines.num+3*i+2].v[0] = v[1]*n + v[2];
110 edgeKey[mesh->lines.num+3*i+2].v[1] = 2*(3*i+2)+1;
111 }
112 }
113
114 HXT_CHECK( group2_sort_v0(edgeKey, numEdgesTotal, n*(n-1)-1) );
115
116 #pragma omp parallel
117 {
118 int threadID = omp_get_thread_num();
119 uint64_t localNum = 0;
120
121 #pragma omp for
122 for (uint64_t i=0; i<numEdgesTotal; i++) {
123 if(edgeKey[i].v[1]%2==0) {
124 if(i!=numEdgesTotal-1 && edgeKey[i].v[0]==edgeKey[i+1].v[0]) {
125 #ifndef NDEBUG
126 if(edgeKey[i+1].v[1]%2==0) {
127 HXT_ERROR_MSG(HXT_STATUS_ERROR, "Duplicated line in mesh->lines (" HXTu64 " & " HXTu64 ")\n"
128 "\tThis case is not handled in Release mode, FIX IT !!",
129 edgeKey[i].v[1]/2, edgeKey[i+1].v[1]/2);
130 exit(EXIT_FAILURE);
131 }
132 else
133 #endif
134 {
135 lines2TriMap[edgeKey[i].v[1]/2] = edgeKey[i+1].v[1]/2;
136 }
137 }
138 else /* the edge is not in a triangle */ {
139 localNum++;
140 lines2TriMap[edgeKey[i].v[1]/2] = HXT_NO_ADJACENT;
141 }
142 }
143 }
144
145 numEdges[threadID] = localNum;
146
147 #pragma omp barrier
148 #pragma omp single
149 {
150 int nthreads = omp_get_num_threads();
151 *missing = 0;
152 for (int i=0; i<nthreads; i++) {
153 *missing += numEdges[i];
154 }
155 }
156 }
157
158 HXT_CHECK( hxtFree(&numEdges) );
159 HXT_CHECK( hxtAlignedFree(&edgeKey) );
160
161 return HXT_STATUS_OK;
162 }
163
164
165 /**************************************************************************************
166 * Lines --> Tetrahedra MAPPING (for each line, get 6*tet+e)
167 *************************************************************************************/
hxtGetLines2TetMap(HXTMesh * mesh,uint64_t * lines2TetMap,uint64_t * missing)168 HXTStatus hxtGetLines2TetMap(HXTMesh* mesh, uint64_t* lines2TetMap, uint64_t* missing)
169 {
170 HXT_ASSERT( lines2TetMap!=NULL );
171 HXT_ASSERT( mesh!=NULL );
172
173 #ifndef NDEBUG
174 if(countBoundaries(mesh)!=0)
175 return HXT_ERROR_MSG(HXT_STATUS_ERROR, "a tetrahedron has no neighbor. Add ghosts tetrahedra to use this function");
176 #endif
177
178
179 const int maxThreads = omp_get_max_threads();
180 const uint64_t n = mesh->vertices.num;
181 HXTStatus status = HXT_STATUS_OK;
182 uint64_t numEdgesTotal;
183
184 HXTGroup2* edgeKey = NULL;
185 uint64_t* numEdges;
186 unsigned char* edgeFlag;
187 HXT_CHECK( hxtMalloc(&numEdges, maxThreads*sizeof(uint64_t)) );
188 HXT_CHECK( hxtAlignedMalloc(&edgeFlag, mesh->tetrahedra.num*sizeof(char)) );
189 memset(edgeFlag, 0, mesh->tetrahedra.num*sizeof(char));
190
191
192 #pragma omp parallel
193 {
194 const int threadID = omp_get_thread_num();
195 uint64_t localNum = 0;
196
197 #pragma omp for schedule(static)
198 for (uint64_t tet=0; tet<mesh->tetrahedra.num; tet++) {
199 for (int edge=0; edge<6; edge++) {
200 // get some info on the edge (adjacent facets and nodes)
201
202 unsigned in_facet, out_facet;
203 getFacetsFromEdge(edge, &in_facet, &out_facet);
204
205 uint32_t p0, p1;
206 {
207 unsigned n0, n1;
208 getNodesFromEdge(edge, &n0, &n1);
209 p0 = mesh->tetrahedra.node[4*tet + n0];
210 p1 = mesh->tetrahedra.node[4*tet + n1];
211 }
212
213 if(p0==HXT_GHOST_VERTEX || p1==HXT_GHOST_VERTEX)
214 continue;
215
216 /* visit tetrahedra `tet` by turning around the edge
217 we only want the edge to be registered once,
218 thus we stop as soon as `transitiveHashLess(curTet, tet)==true`
219 */
220 int truth = 1;
221 uint64_t curTet = tet;
222 do
223 {
224 uint32_t newV = mesh->tetrahedra.node[4*curTet + in_facet];
225
226 // go into the neighbor through out_facet
227 uint64_t neigh = mesh->tetrahedra.neigh[4*curTet + out_facet];
228 curTet = neigh/4;
229 in_facet = neigh%4;
230
231 if(transitiveHashLess(curTet, tet)) {
232 truth=0;
233 break;
234 }
235
236 uint32_t* nodes = mesh->tetrahedra.node + 4*curTet;
237 for (out_facet=0; out_facet<3; out_facet++)
238 if(nodes[out_facet]==newV)
239 break;
240
241 } while (curTet!=tet);
242
243 // only the tetrahedron with the biggest hash will register the edge
244 if(truth){
245 edgeFlag[tet] |= 1U<<edge;
246 localNum++;
247 }
248 }
249 }
250
251 numEdges[threadID] = localNum;
252
253 #pragma omp barrier
254 #pragma omp single
255 {
256 // exclusive scan used to compute starting indices
257 int nthreads = omp_get_num_threads();
258 numEdgesTotal = mesh->lines.num;
259 for (int i=0; i<nthreads; i++) {
260 uint32_t tsum = numEdgesTotal + numEdges[i];
261 numEdges[i] = numEdgesTotal;
262 numEdgesTotal = tsum;
263 }
264
265 #ifndef NDEBUG
266 if(numEdgesTotal>2*mesh->tetrahedra.num+mesh->lines.num){
267 HXT_ERROR_MSG(HXT_STATUS_ERROR,
268 "There is less than 2 tetrahedra per edge in average,"
269 "which means the mesh is totally broken !");
270 exit(EXIT_FAILURE);
271 }
272 #endif
273
274 status = hxtAlignedMalloc(&edgeKey, numEdgesTotal*sizeof(HXTGroup2));
275 }
276
277 if(status==HXT_STATUS_OK) {
278 // copy the edges from mesh->lines in the edgeKey struct array
279 #pragma omp for
280 for (uint64_t l=0; l<mesh->lines.num; l++) {
281 uint32_t p0 = mesh->lines.node[2*l+0];
282 uint32_t p1 = mesh->lines.node[2*l+1];
283
284 if(p0<p1) {
285 edgeKey[l].v[0] = p0*n + p1;
286 edgeKey[l].v[1] = 2*l;
287 }
288 else if(p0>p1){
289 edgeKey[l].v[0] = p1*n + p0;
290 edgeKey[l].v[1] = 2*l;
291 }
292 else {
293 edgeKey[l].v[0] = p0*n + p0; // the line begins and ends at the same point...
294 edgeKey[l].v[1] = 1;
295 lines2TetMap[l] = HXT_NO_ADJACENT;
296 }
297
298
299 }
300
301 localNum = numEdges[threadID];
302 #pragma omp for schedule(static)
303 for (uint64_t tet=0; tet<mesh->tetrahedra.num; tet++) {
304 for (unsigned edge=0; edge<6; edge++) {
305 if(edgeFlag[tet] & (1U<<edge)){
306 uint32_t p0, p1;
307 {
308 unsigned n0, n1;
309 getNodesFromEdge(edge, &n0, &n1);
310 p0 = mesh->tetrahedra.node[4*tet + n0];
311 p1 = mesh->tetrahedra.node[4*tet + n1];
312 }
313
314 if(p0<p1){
315 edgeKey[localNum].v[0] = p0*n + p1;
316 }
317 else {
318 edgeKey[localNum].v[0] = p1*n + p0;
319 }
320 edgeKey[localNum].v[1] = 2*(6*tet+edge)+1;
321
322 localNum++;
323 }
324 }
325 }
326 }
327 }
328
329 HXT_CHECK( hxtAlignedFree(&edgeFlag) );
330 HXT_CHECK( status );
331
332 HXT_CHECK( group2_sort_v0(edgeKey, numEdgesTotal, n*(n-1)-1) );
333
334 #pragma omp parallel
335 {
336 const int threadID = omp_get_thread_num();
337 uint64_t localNum = 0;
338
339 #pragma omp for
340 for (uint64_t i=0; i<numEdgesTotal; i++) {
341 if(edgeKey[i].v[1]%2==0) {
342 if(i!=numEdgesTotal-1 && edgeKey[i].v[0]==edgeKey[i+1].v[0]) {
343 #ifndef NDEBUG
344 if(edgeKey[i+1].v[1]%2==0) {
345 HXT_ERROR_MSG(HXT_STATUS_ERROR, "Duplicated line in mesh->lines (" HXTu64 " & " HXTu64 ")\n"
346 "\tThis issue will not produce a direct error if NDEBUG is defined",
347 edgeKey[i].v[1]/2, edgeKey[i+1].v[1]/2);
348 exit(EXIT_FAILURE);
349 }
350 else
351 #endif
352 {
353 lines2TetMap[edgeKey[i].v[1]/2] = edgeKey[i+1].v[1]/2;
354 }
355 }
356 else {
357 lines2TetMap[edgeKey[i].v[1]/2] = HXT_NO_ADJACENT;
358 localNum++;
359 }
360 }
361 }
362
363 numEdges[threadID] = localNum;
364
365 #pragma omp barrier
366 #pragma omp single
367 {
368 int nthreads = omp_get_num_threads();
369 *missing = 0;
370 for (int i=0; i<nthreads; i++) {
371 *missing+=numEdges[i];
372 }
373 }
374 }
375
376 HXT_CHECK( hxtFree(&numEdges) );
377 HXT_CHECK( hxtAlignedFree(&edgeKey) );
378 return HXT_STATUS_OK;
379 }
380
381
382 /**************************************************************************************
383 * Triangles --> Tetrahedra MAPPING
384 *************************************************************************************/
hxtGetTri2TetMap(HXTMesh * mesh,uint64_t * tri2TetMap,uint64_t * missing)385 HXTStatus hxtGetTri2TetMap(HXTMesh* mesh, uint64_t* tri2TetMap, uint64_t* missing)
386 {
387 HXT_ASSERT(tri2TetMap!=NULL);
388
389 if(mesh->triangles.num==0)
390 return HXT_STATUS_OK;
391
392 const uint64_t n = mesh->vertices.num;
393 const int maxThreads = omp_get_max_threads();
394 HXTStatus status = HXT_STATUS_OK;
395 uint64_t numTrianglesTotal;
396
397
398 HXTGroup2 *triKey = NULL;
399 HXTGroup3* pairKey = NULL;
400 uint64_t *numTriangles;
401 HXT_CHECK( hxtMalloc(&numTriangles, maxThreads*sizeof(uint64_t)) );
402
403
404 #ifndef NDEBUG
405 uint64_t nGhosts = 0;
406 uint64_t nBoundaries = countBoundaries(mesh);
407
408 #pragma omp parallel for reduction(+:nGhosts)
409 for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
410 if(mesh->tetrahedra.node[4*i+3]==HXT_GHOST_VERTEX)
411 nGhosts++;
412 }
413
414 if(n <= 2642246){
415 HXT_CHECK( hxtAlignedMalloc(&triKey, (2*(mesh->tetrahedra.num+nBoundaries)-3*(nGhosts+nBoundaries)/2+mesh->triangles.num)*sizeof(HXTGroup2)) );
416 }
417 else{
418 HXT_CHECK( hxtAlignedMalloc(&pairKey, (2*(mesh->tetrahedra.num+nBoundaries)-3*(nGhosts+nBoundaries)/2+mesh->triangles.num)*sizeof(HXTGroup3)) );
419 }
420 #endif
421
422 // create a list of triangles in the mesh:
423 #pragma omp parallel
424 {
425 const int threadID = omp_get_thread_num();
426 uint64_t localNum = 0;
427
428 #pragma omp for schedule(static)
429 for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
430 if(mesh->tetrahedra.node[4*i+3]!=HXT_GHOST_VERTEX) {
431 for (int j=0; j<4; j++) {
432 uint64_t neigh = mesh->tetrahedra.neigh[4*i+j];
433 if(neigh==HXT_NO_ADJACENT || hashLess(i, neigh/4) ||
434 mesh->tetrahedra.node[neigh]==HXT_GHOST_VERTEX)
435 localNum++;
436 }
437 }
438 }
439
440 numTriangles[threadID] = localNum;
441
442 #pragma omp barrier
443 #pragma omp single
444 {
445 // exclusive scan used to compute starting indices
446 int nthreads = omp_get_num_threads();
447 numTrianglesTotal = mesh->triangles.num;
448 for (int i=0; i<nthreads; i++) {
449 uint32_t tsum = numTrianglesTotal + numTriangles[i];
450 numTriangles[i] = numTrianglesTotal;
451 numTrianglesTotal = tsum;
452 }
453
454 #ifndef NDEBUG
455 if(numTrianglesTotal!=2*(mesh->tetrahedra.num+nBoundaries)-3*(nGhosts+nBoundaries)/2+mesh->triangles.num){
456 HXT_ERROR_MSG(HXT_STATUS_ERROR, "you should never go here..."
457 " (" HXTu64 "!=2*" HXTu64 "+(" HXTu64 "-3*" HXTu64 ")/2",
458 numTrianglesTotal-mesh->triangles.num,
459 mesh->tetrahedra.num, nBoundaries, nGhosts);
460 exit(EXIT_FAILURE);
461 }
462 #else
463 if(n <= 2642246){
464 status = hxtAlignedMalloc(&triKey, numTrianglesTotal*sizeof(HXTGroup2));
465 }
466 else{
467 status = hxtAlignedMalloc(&pairKey, numTrianglesTotal*sizeof(HXTGroup3));
468 }
469 #endif
470 }
471
472 if(status==HXT_STATUS_OK) {
473 // copy the triangles from mesh->triangles in the triKey struct array
474 if(n <= 2642246){
475 #pragma omp for
476 for (uint64_t i=0; i<mesh->triangles.num; i++) {
477 uint32_t v[3] = {mesh->triangles.node[3*i+0],
478 mesh->triangles.node[3*i+1],
479 mesh->triangles.node[3*i+2]};
480 HXTSORT_3_VALUES_INPLACE(uint32_t, v);
481 triKey[i].v[1] = 2*i;
482 triKey[i].v[0] = v[0]*(n-1)*n + v[1]*n + v[2];
483 }
484 }
485 else{
486 #pragma omp for
487 for (uint64_t i=0; i<mesh->triangles.num; i++) {
488 uint32_t v[3] = {mesh->triangles.node[3*i+0],
489 mesh->triangles.node[3*i+1],
490 mesh->triangles.node[3*i+2]};
491 HXTSORT_3_VALUES_INPLACE(uint32_t, v);
492 pairKey[i].v[2] = 2*i;
493 pairKey[i].v[1] = v[0]*(n-1) + v[1];
494 pairKey[i].v[0] = v[2];
495 }
496 }
497
498 // add the triangle from the tetrahedral mesh to the triKey struct array
499 localNum = numTriangles[threadID];
500 if(n <= 2642246){
501 #pragma omp for schedule(static)
502 for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
503 if(mesh->tetrahedra.node[4*i+3]!=HXT_GHOST_VERTEX){
504 for (int j=0; j<4; j++) {
505 uint64_t neigh = mesh->tetrahedra.neigh[4*i+j];
506 if(neigh==HXT_NO_ADJACENT || hashLess(i, neigh/4) ||
507 mesh->tetrahedra.node[neigh]==HXT_GHOST_VERTEX) {
508 uint32_t v[3] = {mesh->tetrahedra.node[4*i+((j+1)&3)],
509 mesh->tetrahedra.node[4*i+((j+2)&3)],
510 mesh->tetrahedra.node[4*i+((j+3)&3)]};
511 HXTSORT_3_VALUES_INPLACE(uint32_t, v);
512 triKey[localNum].v[1] = 2*(4*i+j)+1;
513 triKey[localNum].v[0] = v[0]*(n-1)*n + v[1]*n + v[2]; // max: (n-3)(n-1)n + (n-2)n + (n-1)
514 // = (n-2)(n-1)n - 1
515 localNum++;
516 }
517 }
518 }
519 }
520 }
521 else {
522 #pragma omp for schedule(static)
523 for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
524 if(mesh->tetrahedra.node[4*i+3]!=HXT_GHOST_VERTEX){
525 for (int j=0; j<4; j++) {
526 uint64_t neigh = mesh->tetrahedra.neigh[4*i+j];
527 if(neigh==HXT_NO_ADJACENT || hashLess(i, neigh/4) ||
528 mesh->tetrahedra.node[neigh]==HXT_GHOST_VERTEX) {
529 uint32_t v[3] = {mesh->tetrahedra.node[4*i+((j+1)&3)],
530 mesh->tetrahedra.node[4*i+((j+2)&3)],
531 mesh->tetrahedra.node[4*i+((j+3)&3)]};
532 HXTSORT_3_VALUES_INPLACE(uint32_t, v);
533 pairKey[localNum].v[2] = 2*(4*i+j)+1;
534 pairKey[localNum].v[1] = v[0]*(n-1) + v[1]; // max: (n-3)(n-1) + (n-2) = (n-2)(n-1) - 1
535 pairKey[localNum].v[0] = v[2]; // max: n-1
536
537 localNum++;
538 }
539 }
540 }
541 }
542 }
543 }
544 }
545
546 HXT_CHECK(status);
547
548 if(n <= 2642246){
549 // sort triKey
550 HXT_CHECK( group2_sort_v0(triKey, numTrianglesTotal, (n-2)*(n-1)*n-1) );
551 }
552 else{
553 HXT_CHECK( group3_sort_v0(pairKey, numTrianglesTotal, n-1) );
554 HXT_CHECK( group3_sort_v1(pairKey, numTrianglesTotal, (n-2)*(n-1)-1) );
555 }
556
557
558 #pragma omp parallel
559 {
560 const int threadID = omp_get_thread_num();
561 uint64_t localNum = 0;
562
563 if(n <= 2642246){
564 #pragma omp for
565 for (uint64_t i=0; i<numTrianglesTotal; i++) {
566 if(triKey[i].v[1]%2==0) {
567 if(i!=numTrianglesTotal-1 && triKey[i].v[0]==triKey[i+1].v[0]) {
568 #ifndef NDEBUG
569 if(triKey[i+1].v[1]%2==0) {
570 HXT_ERROR_MSG(HXT_STATUS_ERROR, "Duplicated triangle in mesh->triangles (" HXTu64 " & " HXTu64 ")\n"
571 "\tThis case is not handled in Release mode, FIX IT !!",
572 triKey[i].v[1]/2, triKey[i+1].v[1]/2);
573 exit(EXIT_FAILURE);
574 }
575 else
576 #endif
577 {
578 tri2TetMap[triKey[i].v[1]/2] = triKey[i+1].v[1]/2;
579 }
580 }
581 else /* the triangle is missing */ {
582 localNum++;
583 tri2TetMap[triKey[i].v[1]/2] = HXT_NO_ADJACENT;
584 }
585 }
586 }
587 }
588 else{
589 #pragma omp for
590 for (uint64_t i=0; i<numTrianglesTotal; i++) {
591 if(pairKey[i].v[2]%2==0) {
592 if(i!=numTrianglesTotal-1 && pairKey[i].v[0]==pairKey[i+1].v[0] && pairKey[i].v[1]==pairKey[i+1].v[1]) {
593 #ifndef NDEBUG
594 if(pairKey[i+1].v[2]%2==0) {
595 HXT_ERROR_MSG(HXT_STATUS_ERROR, "Duplicated triangle in mesh->triangles (" HXTu64 " & " HXTu64 ")\n"
596 "\tThis case is not handled in Release mode, FIX IT !!",
597 pairKey[i].v[2]/2, pairKey[i+1].v[2]/2);
598 exit(EXIT_FAILURE);
599 }
600 else
601 #endif
602 {
603 tri2TetMap[pairKey[i].v[2]/2] = pairKey[i+1].v[2]/2;
604 }
605 }
606 else /* the triangle is missing */ {
607 localNum++;
608 tri2TetMap[pairKey[i].v[2]/2] = HXT_NO_ADJACENT;
609 }
610 }
611 }
612 }
613
614 numTriangles[threadID] = localNum;
615
616 #pragma omp barrier
617 #pragma omp single
618 {
619 int nthreads = omp_get_num_threads();
620 *missing = 0;
621 for (int i=0; i<nthreads; i++) {
622 *missing += numTriangles[i];
623 }
624 }
625 }
626
627 HXT_CHECK( hxtFree(&numTriangles) );
628 HXT_CHECK( hxtAlignedFree(&triKey) );
629 HXT_CHECK( hxtAlignedFree(&pairKey) );
630
631 return HXT_STATUS_OK;
632 }
633
634
635 /**************************************************************************************
636 * Constrain facets of tetrahedron if it is in tri2TetMap
637 *************************************************************************************/
hxtConstrainTriangles(HXTMesh * mesh,uint64_t * tri2TetMap)638 HXTStatus hxtConstrainTriangles(HXTMesh* mesh, uint64_t* tri2TetMap)
639 {
640 HXT_ASSERT(tri2TetMap!=NULL);
641 HXT_ASSERT(mesh!=NULL);
642
643 char* faceFlag;
644 HXT_CHECK( hxtAlignedMalloc(&faceFlag, 4*mesh->tetrahedra.num*sizeof(char)) );
645 memset(faceFlag, 0, 4*mesh->tetrahedra.num*sizeof(char));
646
647 // fill faceFlag
648 #pragma omp parallel for
649 for (uint64_t i=0; i<mesh->triangles.num; i++) {
650 if(tri2TetMap[i]!=HXT_NO_ADJACENT) {
651 faceFlag[tri2TetMap[i]] = 1;
652 faceFlag[mesh->tetrahedra.neigh[tri2TetMap[i]]] = 1;
653 }
654 }
655
656 // constrain corresponding flag, tetrahedron by tetrahedron to avoid race conditions
657 #pragma omp parallel for
658 for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
659 for (uint64_t j=0; j<4; j++) {
660 if(faceFlag[4*i+j]) {
661 setFacetConstraint(mesh, i, j);
662 }
663 }
664 }
665
666 HXT_CHECK( hxtAlignedFree(&faceFlag) );
667
668 return HXT_STATUS_OK;
669 }
670
671
672 /*********************************************************
673 * Constrain an edge in all tetrahedra surrounding it
674 *********************************************************
675 * In single-thread cases, put edgeFlag==NULL.
676 * This function will directly set the constraint bit corresponding to the edge
677 * '6*firstTet+edge' on all tetrahedra surrounding this edge.
678 *
679 * In multi-threaded case, if multiple thread modify different edges of the
680 * same tetrahedra, modifying the same flag would result in a race condition.
681 * Therefore, in parallel section, you must give an array with a char per edge.
682 * This function will set the edges corresponding to '6*firstTet+edge' to 1.
683 */
hxtConstrainLine(HXTMesh * mesh,uint64_t tet,int edge,char * edgeFlag)684 HXTStatus hxtConstrainLine(HXTMesh* mesh, uint64_t tet, int edge, char* edgeFlag)
685 {
686 uint64_t curTet = tet;
687 unsigned in_facet, out_facet;
688 getFacetsFromEdge(edge, &in_facet, &out_facet);
689
690 // turn around the edge to set edgeFlag of all tetrahedra to 1...
691 do
692 {
693 if(edgeFlag)
694 edgeFlag[6*curTet + getEdgeFromFacets(in_facet, out_facet)] = 1;
695 else
696 setEdgeConstraint(mesh, curTet, getEdgeFromFacets(in_facet, out_facet));
697
698 uint32_t newV = mesh->tetrahedra.node[4*curTet + in_facet];
699
700 // go into the neighbor through out_facet
701 uint64_t neigh = mesh->tetrahedra.neigh[4*curTet + out_facet];
702 curTet = neigh/4;
703 in_facet = neigh%4;
704 uint32_t* nodes = mesh->tetrahedra.node + 4*curTet;
705 for (out_facet=0; out_facet<3; out_facet++)
706 if(nodes[out_facet]==newV)
707 break;
708
709 } while (curTet!=tet);
710
711 return HXT_STATUS_OK;
712 }
713
714
715
716 /**************************************************************************************
717 * Constrain edge of tetrahedron if it is in lines2TetMap but not in lines2TriMap
718 *************************************************************************************/
hxtConstrainLinesNotInTriangles(HXTMesh * mesh,uint64_t * lines2TetMap,uint64_t * lines2TriMap)719 HXTStatus hxtConstrainLinesNotInTriangles(HXTMesh* mesh, uint64_t* lines2TetMap, uint64_t* lines2TriMap)
720 {
721 HXT_ASSERT(lines2TetMap!=NULL);
722 HXT_ASSERT(lines2TriMap!=NULL);
723 HXT_ASSERT(mesh!=NULL);
724
725 char* edgeFlag;
726 HXT_CHECK( hxtAlignedMalloc(&edgeFlag, 6*mesh->tetrahedra.num*sizeof(char)) );
727 memset(edgeFlag, 0, 6*mesh->tetrahedra.num*sizeof(char));
728
729 #pragma omp parallel for
730 for (uint64_t i=0; i<mesh->lines.num; i++) {
731 if(lines2TriMap[i]==HXT_NO_ADJACENT && lines2TetMap[i]!=HXT_NO_ADJACENT) {
732 hxtConstrainLine(mesh, lines2TetMap[i]/6, lines2TetMap[i]%6, edgeFlag);
733 }
734 }
735
736 #pragma omp parallel for
737 for (uint64_t i=0; i<mesh->tetrahedra.num; i++) {
738 for (int j=0; j<6; j++) {
739 if(edgeFlag[6*i+j])
740 setEdgeConstraint(mesh, i, j);
741 }
742 }
743
744
745 HXT_CHECK( hxtAlignedFree(&edgeFlag) );
746
747 return HXT_STATUS_OK;
748 }