Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
hash_3d.c
Go to the documentation of this file.
1/* =============================================================================
2** This file is part of the mmg software package for the tetrahedral
3** mesh modification.
4** Copyright (c) Bx INP/CNRS/Inria/UBordeaux/UPMC, 2004-
5**
6** mmg is free software: you can redistribute it and/or modify it
7** under the terms of the GNU Lesser General Public License as published
8** by the Free Software Foundation, either version 3 of the License, or
9** (at your option) any later version.
10**
11** mmg is distributed in the hope that it will be useful, but WITHOUT
12** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14** License for more details.
15**
16** You should have received a copy of the GNU Lesser General Public
17** License and of the GNU General Public License along with mmg (in
18** files COPYING.LESSER and COPYING). If not, see
19** <http://www.gnu.org/licenses/>. Please read their terms carefully and
20** use this copy of the mmg distribution only if you accept them.
21** =============================================================================
22*/
23
35#include "libmmg3d.h"
36#include "libmmg3d_private.h"
37
38#define MMG5_KC 13
39
40extern int8_t ddb;
41
51 MMG5_pTetra pt,pt1;
52 MMG5_int k;
53
54 k = 1;
55 do {
56 pt = &mesh->tetra[k];
57 if ( !MG_EOK(pt) ) {
58 pt1 = &mesh->tetra[mesh->ne];
59 assert( pt && pt1 && MG_EOK(pt1) );
60 memcpy(pt,pt1,sizeof(MMG5_Tetra));
61 if ( !MMG3D_delElt(mesh,mesh->ne) ) return 0;
62 }
63 }
64 while ( ++k < mesh->ne );
65
66 /* Recreate nil chain */
67 assert(mesh->ne<=mesh->nemax);
68
69 if ( mesh->ne == mesh->nemax )
70 mesh->nenil = 0;
71 else {
72 mesh->nenil = mesh->ne + 1;
73
74 for(k=mesh->nenil; k<=mesh->nemax-1; k++){
75 mesh->tetra[k].v[3] = k+1;
76 }
77
78 mesh->tetra[mesh->nemax].v[3] = 0;
79 }
80 return 1;
81}
82
84MMG5_int MMG5_hashGetFace(MMG5_Hash *hash,MMG5_int ia,MMG5_int ib,MMG5_int ic) {
85 MMG5_hedge *ph;
86 MMG5_int key;
87 MMG5_int mins,maxs,sum;
88
89 if ( !hash->item ) return 0;
90
91 mins = MG_MIN(ia,MG_MIN(ib,ic));
92 maxs = MG_MAX(ia,MG_MAX(ib,ic));
93
94 /* compute key */
95 sum = ia + ib + ic;
96 key = (MMG5_KA*(int64_t)mins + MMG5_KB*(int64_t)maxs) % hash->siz;
97 ph = &hash->item[key];
98
99 if ( ph->a ) {
100 if ( ph->a == mins && ph->b == maxs && ph->s == sum )
101 return ph->k;
102 else {
103 while ( ph->nxt ) {
104 ph = &hash->item[ph->nxt];
105 if ( ph->a == mins && ph->b == maxs && ph->s == sum ) return ph->k;
106 }
107 }
108 }
109
110 return 0;
111}
112
123 MMG5_pTetra pt,pt1;
124 MMG5_int key;
125 MMG5_int k,kk,pp,l,ll,mins,mins1,maxs,maxs1,sum,sum1,iadr;
126 MMG5_int *hcode,*link,hsize,inival;
127 uint8_t i,ii,i1,i2,i3;
128
129 /* default */
130 if ( mesh->adja ) {
131 return 1;
132 }
133
134 if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug )
135 fprintf(stdout," ** SETTING STRUCTURE\n");
136
137 /* packing : if not hash does not work */
138 if ( pack ) {
139 if ( ! MMG5_paktet(mesh) ) return 0;
140 }
141
142 /* memory alloc */
143 MMG5_ADD_MEM(mesh,(4*mesh->nemax+5)*sizeof(MMG5_int),"adjacency table",
144 fprintf(stderr," Exit program.\n");
145 return 0);
146 MMG5_SAFE_CALLOC(mesh->adja,4*mesh->nemax+5,MMG5_int,return 0);
147 MMG5_SAFE_CALLOC(hcode,mesh->ne+5,MMG5_int,return 0);
148
149 link = mesh->adja;
150 hsize = mesh->ne;
151
152 /* init */
153 if ( mesh->info.ddebug ) fprintf(stdout," h- stage 1: init\n");
154 inival = INT_MAX;
155 iadr = 0;
156 for (k=0; k<=mesh->ne; k++)
157 hcode[k] = -inival;
158
159 /* hash tetras */
160 for (k=1; k<=mesh->ne; k++) {
161 pt = &mesh->tetra[k];
162 if ( !MG_EOK(pt) ) continue;
163 for (i=0; i<4; i++) {
164 i1 = MMG5_idir[i][0];
165 i2 = MMG5_idir[i][1];
166 i3 = MMG5_idir[i][2];
167 mins = MG_MIN(pt->v[i1],MG_MIN(pt->v[i2],pt->v[i3]));
168 maxs = MG_MAX(pt->v[i1],MG_MAX(pt->v[i2],pt->v[i3]));
169
170 /* compute key and insert */
171 sum = pt->v[i1] + pt->v[i2] + pt->v[i3];
172 key = (MMG5_KA*(int64_t)mins+MMG5_KB*(int64_t)maxs+MMG5_KC*(int64_t)sum)%hsize+1;
173 iadr++;
174 link[iadr] = hcode[key];
175 hcode[key] = -iadr;
176 }
177 }
178
179 /* set adjacency */
180 if ( mesh->info.ddebug ) fprintf(stdout," h- stage 2: adjacencies\n");
181 for (l=iadr; l>0; l--) {
182 if ( link[l] >= 0 ) continue;
183
184 /* current element */
185 k = (l-1) / 4 + 1;
186 i = (l-1) % 4;
187 i1 = MMG5_idir[i][0];
188 i2 = MMG5_idir[i][1];
189 i3 = MMG5_idir[i][2];
190 pt = &mesh->tetra[k];
191 mins = MG_MIN(pt->v[i1],MG_MIN(pt->v[i2],pt->v[i3]));
192 maxs = MG_MAX(pt->v[i1],MG_MAX(pt->v[i2],pt->v[i3]));
193 sum = pt->v[i1] + pt->v[i2] + pt->v[i3];
194
195 /* accross link */
196 ll = -link[l];
197 pp = 0;
198 link[l] = 0;
199 while ( ll != inival ) {
200 kk = (ll-1) / 4 + 1;
201 ii = (ll-1) % 4;
202 i1 = MMG5_idir[ii][0];
203 i2 = MMG5_idir[ii][1];
204 i3 = MMG5_idir[ii][2];
205 pt1 = &mesh->tetra[kk];
206 sum1 = pt1->v[i1] + pt1->v[i2] + pt1->v[i3];
207 if ( sum1 == sum ) {
208 mins1 = MG_MIN(pt1->v[i1],MG_MIN(pt1->v[i2],pt1->v[i3]));
209 maxs1 = MG_MAX(pt1->v[i1],MG_MAX(pt1->v[i2],pt1->v[i3]));
210
211 /* adjacent found */
212 if ( mins1 == mins && maxs1 == maxs ) {
213 if ( pp != 0 ) link[pp] = link[ll];
214 link[l] = 4*kk + ii;
215 link[ll] = 4*k + i;
216 break;
217 }
218 }
219 pp = ll;
220 ll = -link[ll];
221 }
222 }
223 MMG5_SAFE_FREE(hcode);
224 return 1;
225}
226
239 MMG5_pPrism pp,pp1;
240 MMG5_int key;
241 MMG5_int k,kk,l,ll,jj;
242 MMG5_int max12,min12,max34,min34,mins,mins1,mins_b, mins_b1,maxs,maxs1;
243 MMG5_int iadr;
244 MMG5_int *hcode,*link,hsize,inival;
245 uint8_t i,ii,i1,i2,i3,i4;
246
247 if ( !mesh->nprism ) return 1;
248
249 /* default */
250 if ( mesh->adjapr ) {
251 if ( abs(mesh->info.imprim) > 3 || mesh->info.ddebug ) {
252 fprintf(stderr,"\n ## Warning: %s: no re-build of adjacencies of prisms. "
253 "mesh->adjapr must be freed to enforce analysis.\n",__func__);
254 }
255 return 1;
256 }
257
258 if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug )
259 fprintf(stdout," ** SETTING PRISMS ADJACENCY\n");
260
261 /* memory alloc */
262 MMG5_ADD_MEM(mesh,(5*mesh->nprism+6)*sizeof(MMG5_int),"prism adjacency table",
263 printf(" Exit program.\n");
264 return 0);
265 MMG5_SAFE_CALLOC(mesh->adjapr,5*mesh->nprism+6,MMG5_int,return 0);
266 MMG5_SAFE_CALLOC(hcode,mesh->nprism+6,MMG5_int,return 0);
267
268 link = mesh->adjapr;
269 hsize = mesh->nprism;
270
271 /* init */
272 if ( mesh->info.ddebug ) fprintf(stdout," h- stage 1: init\n");
273 inival = INT_MAX;
274 iadr = 0;
275 for (k=0; k<=mesh->nprism; k++)
276 hcode[k] = -inival;
277
278 /* hash prism */
279 for (k=1; k<=mesh->nprism; k++) {
280 pp = &mesh->prism[k];
281 assert ( MG_EOK(pp) );
282 for (i=0; i<2; i++) {
283 /* Triangular face */
284 i1 = MMG5_idir_pr[i][0];
285 i2 = MMG5_idir_pr[i][1];
286 i3 = MMG5_idir_pr[i][2];
287
288 min12 = MG_MIN(pp->v[i1],pp->v[i2]);
289 /* mins = minimum index of triangle vertices */
290 mins = MG_MIN(min12,pp->v[i3]);
291
292 max12 = MG_MAX(pp->v[i1],pp->v[i2]);
293 /* maxs = maximum index of triangle vertices */
294 maxs = MG_MAX(max12,pp->v[i3]);
295
296 /* mins_b = second minimum index of triangle vertices */
297 mins_b = pp->v[i1] + pp->v[i2] + pp->v[i3] -mins -maxs;
298
299 /* compute key and insert */
300 key = (MMG5_KA*(int64_t)mins+MMG5_KB*(int64_t)mins_b+MMG5_KC*(int64_t)maxs)%hsize+1;
301 iadr++;
302 link[iadr] = hcode[key];
303 hcode[key] = -iadr;
304 }
305 for ( ; i<5; ++i) {
306 /* Quadrilateral face */
307 i1 = MMG5_idir_pr[i][0];
308 i2 = MMG5_idir_pr[i][1];
309 i3 = MMG5_idir_pr[i][2];
310 i4 = MMG5_idir_pr[i][3];
311
312 min12 = MG_MIN(pp->v[i1],pp->v[i2]);
313 min34 = MG_MIN(pp->v[i3],pp->v[i4]);
314 /* mins = minimum index of quadrilateral vertices */
315 mins = MG_MIN(min12,min34);
316
317 max12 = MG_MAX(pp->v[i1],pp->v[i2]);
318 max34 = MG_MAX(pp->v[i3],pp->v[i4]);
319 /* maxs = maximum index of quadrilateral vertices */
320 maxs = MG_MAX(max12,max34);
321
322 /* mins_b = second minimum index of quadrilateral vertices */
323 mins_b = MG_MIN( MG_MIN(max12,max34),MG_MAX(min12,min34));
324
325 /* compute key and insert */
326 key = (MMG5_KA*(int64_t)mins+MMG5_KB*(int64_t)mins_b+MMG5_KC*(int64_t)maxs)%hsize+1;
327 iadr++;
328 link[iadr] = hcode[key];
329 hcode[key] = -iadr;
330 }
331 }
332
333 /* set adjacency */
334 if ( mesh->info.ddebug ) fprintf(stdout," h- stage 2: adjacencies\n");
335 for (l=iadr; l>0; l--) {
336 if ( link[l] >= 0 ) continue;
337
338 /* current element */
339 k = (l-1) / 5 + 1;
340 i = (l-1) % 5;
341
342 switch (i)
343 {
344 case 0:
345 case 1:
346 i1 = MMG5_idir_pr[i][0];
347 i2 = MMG5_idir_pr[i][1];
348 i3 = MMG5_idir_pr[i][2];
349 pp = &mesh->prism[k];
350
351 min12 = MG_MIN(pp->v[i1],pp->v[i2]);
352 mins = MG_MIN(min12,pp->v[i3]);
353
354 max12 = MG_MAX(pp->v[i1],pp->v[i2]);
355 maxs = MG_MAX(max12,pp->v[i3]);
356
357 mins_b = pp->v[i1] + pp->v[i2] + pp->v[i3] - mins - maxs;
358
359 break;
360
361 default:
362 i1 = MMG5_idir_pr[i][0];
363 i2 = MMG5_idir_pr[i][1];
364 i3 = MMG5_idir_pr[i][2];
365 i4 = MMG5_idir_pr[i][3];
366 pp = &mesh->prism[k];
367
368 min12 = MG_MIN(pp->v[i1],pp->v[i2]);
369 min34 = MG_MIN(pp->v[i3],pp->v[i4]);
370 mins = MG_MIN(min12,min34);
371
372 max12 = MG_MAX(pp->v[i1],pp->v[i2]);
373 max34 = MG_MAX(pp->v[i3],pp->v[i4]);
374 maxs = MG_MAX(max12,max34);
375
376 mins_b = MG_MIN( MG_MIN(max12,max34),MG_MAX(min12,min34));
377 }
378
379 /* accross link */
380 ll = -link[l];
381 jj = 0;
382 link[l] = 0;
383 while ( ll != inival ) {
384 kk = (ll-1) / 5 + 1;
385 ii = (ll-1) % 5;
386
387 switch (ii)
388 {
389 case 0:
390 case 1:
391 i1 = MMG5_idir_pr[ii][0];
392 i2 = MMG5_idir_pr[ii][1];
393 i3 = MMG5_idir_pr[ii][2];
394
395 pp1 = &mesh->prism[kk];
396
397 min12 = MG_MIN(pp1->v[i1],pp1->v[i2]);
398 mins1 = MG_MIN(min12,pp1->v[i3]);
399
400 max12 = MG_MAX(pp1->v[i1],pp1->v[i2]);
401 maxs1 = MG_MAX(max12,pp1->v[i3]);
402
403 mins_b1 = pp1->v[i1] + pp1->v[i2] + pp1->v[i3] - mins1 - maxs1;
404
405 break;
406
407 default:
408 i1 = MMG5_idir_pr[ii][0];
409 i2 = MMG5_idir_pr[ii][1];
410 i3 = MMG5_idir_pr[ii][2];
411 i4 = MMG5_idir_pr[ii][3];
412 pp1 = &mesh->prism[kk];
413
414 min12 = MG_MIN(pp1->v[i1],pp1->v[i2]);
415 min34 = MG_MIN(pp1->v[i3],pp1->v[i4]);
416 mins1 = MG_MIN(min12,min34);
417
418 max12 = MG_MAX(pp1->v[i1],pp1->v[i2]);
419 max34 = MG_MAX(pp1->v[i3],pp1->v[i4]);
420 maxs1 = MG_MAX(max12,max34);
421
422 mins_b1 = MG_MIN( MG_MIN(max12,max34),MG_MAX(min12,min34));
423 }
424
425 /* adjacent found */
426 if ( mins1 == mins && maxs1 == maxs && mins_b1 == mins_b ) {
427 if ( jj != 0 ) link[jj] = link[ll];
428 link[l] = 5*kk + ii;
429 link[ll] = 5*k + i;
430 break;
431 }
432
433 jj = ll;
434 ll = -link[ll];
435 }
436 }
437 MMG5_SAFE_FREE(hcode);
438 return 1;
439}
440
455static inline
457 MMG5_pTetra pt;
458 MMG5_pxTetra pxt;
459 MMG5_pTria ptt;
460 MMG5_hedge *ph;
461 MMG5_int adj,pradj,piv;
462 int64_t list[MMG3D_LMAX+2];
463 MMG5_int key;
464 MMG5_int k,l,i1,i2,na,nb,ia,it1,it2, nr;
465 MMG5_int start;
466 int ilist,nbdy,ipa,ipb;
467 int8_t iface,hasadja,i;
468 static int8_t mmgWarn0=0,mmgWarn1=0;
469
470 nr = 0;
471
472 /* First: seek edges at the interface of two distinct domains and mark it as
473 * required */
474 for (k=1; k<=mesh->nt; k++) {
475 ptt = &mesh->tria[k];
476
477 if ( !MG_EOK(ptt) ) continue;
478
479 for (l=0; l<3; l++) {
480
481 /* Skip parallel edges */
482 if ( (ptt->tag[l] & MG_PARBDY) || (ptt->tag[l] & MG_BDY) ) continue;
483
484 if ( ptt->tag[l] & MG_NOM ) {
485 i1 = MMG5_inxt2[l];
486 i2 = MMG5_iprv2[l];
487
488 /* compute key */
489 na = MG_MIN(ptt->v[i1],ptt->v[i2]);
490 nb = MG_MAX(ptt->v[i1],ptt->v[i2]);
491 key = (MMG5_KA*na + MMG5_KB*nb) % hash->siz;
492 ph = &hash->item[key];
493
494 assert(ph->a);
495 while ( ph->a ) {
496 if ( ph->a == na && ph->b == nb ) break;
497 assert(ph->nxt);
498 ph = &hash->item[ph->nxt];
499 }
500 /* Set edge tag and point tags to MG_REQ if the non-manifold edge shared
501 * separated domains */
502 if ( ph->s > 3 ) {
503 start = ptt->cc/4;
504 assert(start);
505 pt = &mesh->tetra[start];
506
507
508 for (ia=0; ia<6; ++ia) {
509 ipa = MMG5_iare[ia][0];
510 ipb = MMG5_iare[ia][1];
511 if ( (pt->v[ipa] == na && pt->v[ipb] == nb) ||
512 (pt->v[ipa] == nb && pt->v[ipb] == na)) break;
513 }
514 assert(ia<6);
515
516
517 /* Travel throug the shell of the edge until reaching a tetra without adjacent
518 * or until reaching the starting tetra */
519 iface = ptt->cc%4;
520 MMG3D_coquilFaceFirstLoop(mesh,start,na,nb,iface,ia,list,&ilist,&it1,&it2,
521 &piv,&adj,&hasadja,&nbdy,1);
522
523 /* At this point, the first travel, in one direction, of the shell is
524 complete. Now, analyze why the travel ended. */
525 if ( adj == start ) {
526 if ( !it2 ) {
527 if ( !mmgWarn0 ) {
528 mmgWarn0 = 1;
529 fprintf(stderr,"\n ## Warning: %s: at least 1 wrong boundary tag:"
530 " Only 0 or 1 boundary triangles founded in the shell of the edge\n",
531 __func__);
532 }
533 }
534 if ( nbdy < 2 )
535 MMG5_coquilFaceErrorMessage(mesh, it1/4, it2/4);
536 }
537 else {
538 /* A boundary has been detected : slightly different configuration */
539 if ( hasadja ) {
540
541 /* Start back everything from this tetra adj */
542 MMG3D_coquilFaceSecondLoopInit(mesh,piv,&iface,&i,list,&ilist,&it1,
543 &pradj,&adj);
544
545 nbdy = 1;
546 while ( adj ) {
547 pradj = adj;
548
549 if ( MMG5_openCoquilTravel(mesh,na,nb,&adj,&piv,&iface,&i)<0 ) {
550 return 0;
551 }
552
553 /* overflow */
554 if ( ++ilist > MMG3D_LMAX-2 ) {
555 if ( !mmgWarn1 ) {
556 mmgWarn1 = 1;
557 fprintf(stderr,"\n ## Warning: %s: problem in surface remesh"
558 " process. At least 1 shell of edge (%" MMG5_PRId "-%" MMG5_PRId ") contains"
559 " too many elts.\n",__func__,MMG3D_indPt(mesh,na),
560 MMG3D_indPt(mesh,nb));
561 fprintf(stderr,"\n ## Try to modify the hausdorff"
562 " number, or/and the maximum mesh.\n");
563 }
564 return 0;
565 }
566
567 pt = &mesh->tetra[pradj];
568 if ( pt->xt ) {
569 pxt = &mesh->xtetra[pt->xt];
570 if ( pxt->ftag[iface] & MG_BDY ) ++nbdy;
571 }
572 }
573
574 assert(!adj);
575 it2 = 4*pradj + iface;
576
577 if ( (!it1 || !it2) || (it1 == it2) ) {
578 MMG5_coquilFaceErrorMessage(mesh, it1/4, it2/4);
579 return 0;
580 }
581 }
582 }
583
584 /* If ph->s do not match the number of encountred boundaries we have
585 separated domains. */
586 if ( nbdy != ph->s ) {
587 if ( !(ptt->tag[l] & MG_REQ) ) {
588 ptt->tag[l] |= MG_REQ;
589 ptt->tag[l] &= ~MG_NOSURF;
590 ++nr;
591 }
592 mesh->point[ptt->v[MMG5_inxt2[l]]].tag |= MG_REQ;
593 mesh->point[ptt->v[MMG5_iprv2[l]]].tag |= MG_REQ;
594 mesh->point[ptt->v[MMG5_inxt2[l]]].tag &= ~MG_NOSURF;
595 mesh->point[ptt->v[MMG5_iprv2[l]]].tag &= ~MG_NOSURF;
596 }
597
598 /* Work done for this edge: reset ph->s/ */
599 ph->s = 0;
600 }
601 }
602 }
603 }
604 if ( mesh->info.ddebug || abs(mesh->info.imprim) > 3 )
605 fprintf(stdout," %" MMG5_PRId " required edges added\n",nr);
606
607 /* Free the edge hash table */
608 MMG5_DEL_MEM(mesh,hash->item);
609 return 1;
610}
611
623static inline
625 MMG5_pTetra ptet;
626 MMG5_pPoint ppt;
627 MMG5_Hash hash;
628 int i,ier;
629 MMG5_int k,base,np,nc,nm,nre, ng, nrp;
630
631 /* Second: seek the non-required non-manifold points and try to analyse
632 * whether they are corner or required. */
633
634 /* Hash table used by boulernm to store the special edges passing through
635 * a given point */
636 np = 0;
637 for (k=1; k<=mesh->np; ++k) {
638 ppt = &mesh->point[k];
639 if ( (!(ppt->tag & MG_NOM)) || (ppt->tag & MG_REQ) ) continue;
640 ++np;
641 }
642
643 if ( ! MMG5_hashNew(mesh,&hash,np,(MMG5_int)(3.71*np)) ) return 0;
644
645 nc = nre = 0;
646 base = ++mesh->base;
647 for (k=1; k<=mesh->ne; ++k) {
648 ptet = &mesh->tetra[k];
649 if ( !MG_EOK(ptet) ) continue;
650
651 for ( i=0; i<4; ++i ) {
652 ppt = &mesh->point[ptet->v[i]];
653
654 /* Skip parallel points */
655 if ( ppt->tag & MG_PARBDY ) continue;
656
657 if ( (!MG_VOK(ppt)) || (ppt->flag==base) ) continue;
658 ppt->flag = base;
659
660 if ( (!(ppt->tag & MG_NOM)) || (ppt->tag & MG_REQ) ) continue;
661
662 ier = MMG5_boulernm(mesh,&hash, k, i, &ng, &nrp, &nm);
663 if ( ier < 0 ) return 0;
664 else if ( !ier ) continue;
665
666 if ( (ng+nrp+nm) > 2 ) {
667 /* More than 2 feature edges are passing through the point: point is
668 * marked as corner */
669 ppt->tag |= MG_CRN + MG_REQ;
670 ppt->tag &= ~MG_NOSURF;
671 nre++;
672 nc++;
673 }
674 else if ( (ng == 2) || (nrp == 2) || (nm == 2) ) {
675 /* Exactly 2 edges of same type are passing through the point: do
676 * nothing */
677 continue;
678 }
679 else if ( (ng+nrp+nm) == 2 ) {
680 /* 2 edges of different type are passing through the point: point is
681 * marked as required */
682 ppt->tag |= MG_REQ;
683 ppt->tag &= ~MG_NOSURF;
684 nre++;
685 }
686 else if ( ng == 1 && !nrp ){
687 ppt->tag |= MG_CRN + MG_REQ;
688 ppt->tag &= ~MG_NOSURF;
689 nre++;
690 nc++;
691 }
692 else if ( (ng+nrp+nm) == 1 ){
693 /* Only 1 feature edge is passing through the point: point is
694 * marked as corner */
695 assert ( (ng == 1) || (nrp==1) || (nm==1) );
696 ppt->tag |= MG_CRN + MG_REQ;
697 ppt->tag &= ~MG_NOSURF;
698 nre++;
699 nc++;
700 }
701 else {
702 assert ( 0 && "unexpected case");
703 }
704 }
705 }
706
707 /* Free the edge hash table */
708 MMG5_DEL_MEM(mesh,hash.item);
709
710 if ( mesh->info.ddebug || abs(mesh->info.imprim) > 3 )
711 fprintf(stdout," %" MMG5_PRId " corner and %" MMG5_PRId " required vertices added\n",nc,nre);
712
713 return 1;
714}
715
728
729 /* First: seek edges at the interface of two distinct domains and mark it as
730 * required */
731 if ( !MMG5_setEdgeNmTag(mesh,hash) ) return 0;
732
733 /* Second: seek the non-required non-manifold points and try to analyse
734 * whether they are corner or required. */
735 if ( !MMG5_setVertexNmTag(mesh) ) return 0;
736
737 return 1;
738}
739
750
752
753 MMG5_ADD_MEM(mesh,(3*mesh->nt+4)*sizeof(MMG5_int),"surfacic adjacency table",return 0);
754 MMG5_SAFE_CALLOC(mesh->adjt,3*mesh->nt+4,MMG5_int,return 0);
755
756 return MMG5_mmgHashTria(mesh, mesh->adjt, hash, mesh->info.iso) ;
757}
758
759
761int MMG5_hashPop(MMG5_Hash *hash,MMG5_int a,MMG5_int b) {
762 MMG5_hedge *ph,*php;
763 MMG5_int key;
764 MMG5_int ia,ib,iph,iphp;
765
766 ia = MG_MIN(a,b);
767 ib = MG_MAX(a,b);
768 key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
769 ph = &hash->item[key];
770
771 if ( !ph->a ) return 0;
772 else if ( ph->a == ia && ph->b == ib ) {
773 if ( !ph->nxt ) {
774 memset(ph,0,sizeof(MMG5_hedge));
775 return 1;
776 }
777 else {
778 iph = ph->nxt;
779 php = ph;
780 ph = &hash->item[ph->nxt];
781 memcpy(php,ph,sizeof(MMG5_hedge));
782 memset(ph,0,sizeof(MMG5_hedge));
783 ph->nxt = hash->nxt;
784 hash->nxt = iph;
785 return 1;
786 }
787 }
788 while ( ph->nxt ) {
789 php = ph;
790 ph = &hash->item[ph->nxt];
791 if ( ph->a == ia && ph->b == ib ) {
792 if ( !ph->nxt ) {
793 memset(ph,0,sizeof(MMG5_hedge));
794 ph->nxt = hash->nxt;
795 hash->nxt = php->nxt;
796 php->nxt = 0;
797 }
798 else {
799 iph = ph->nxt;
800 iphp = php->nxt;
801 php->nxt = iph;
802 memset(ph,0,sizeof(MMG5_hedge));
803 ph->nxt = hash->nxt;
804 hash->nxt = iphp;
805 }
806 return 1;
807 }
808 }
809 return 0;
810}
811
812
825int MMG5_hTag(MMG5_HGeom *hash,MMG5_int a,MMG5_int b,MMG5_int ref,int16_t tag) {
826 MMG5_hgeom *ph;
827 MMG5_int key;
828 MMG5_int ia,ib;
829
830 ia = MG_MIN(a,b);
831 ib = MG_MAX(a,b);
832 key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
833 ph = &hash->geom[key];
834
835 if ( !ph->a )
836 return 0;
837 else if ( ph->a == ia && ph->b == ib ) {
838 ph->tag |= tag;
839 if ( ref ) {
840 ph->ref = ref;
841 }
842 return 1;
843 }
844 while ( ph->nxt ) {
845 ph = &hash->geom[ph->nxt];
846 if ( ph->a == ia && ph->b == ib ) {
847 ph->tag |= tag;
848 if ( ref ) {
849 ph->ref = ref;
850 }
851 return 1;
852 }
853 }
854 return 0;
855}
856
858int MMG5_hPop(MMG5_HGeom *hash,MMG5_int a,MMG5_int b,MMG5_int *ref,int16_t *tag) {
859 MMG5_hgeom *ph,*php;
860 MMG5_int key;
861 MMG5_int ia,ib,iph,iphp;
862
863 *ref = 0;
864 *tag = 0;
865
866 assert ( hash->siz );
867
868 ia = MG_MIN(a,b);
869 ib = MG_MAX(a,b);
870 key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
871 ph = &hash->geom[key];
872
873 if ( !ph->a ) return 0;
874 else if ( ph->a == ia && ph->b == ib ) {
875 *ref = ph->ref;
876 *tag = ph->tag;
877 if ( !ph->nxt ) {
878 memset(ph,0,sizeof(MMG5_hgeom));
879 }
880 else {
881 iph = ph->nxt;
882 php = ph;
883 ph = &hash->geom[ph->nxt];
884 memcpy(php,ph,sizeof(MMG5_hgeom));
885 memset(ph,0,sizeof(MMG5_hgeom));
886 ph->nxt = hash->nxt;
887 hash->nxt = iph;
888 }
889 return 1;
890 }
891 while ( ph->nxt ) {
892 php = ph;
893 ph = &hash->geom[ph->nxt];
894 if ( ph->a == ia && ph->b == ib ) {
895 *ref = ph->ref;
896 *tag = ph->tag;
897 if ( !ph->nxt ) {
898 memset(ph,0,sizeof(MMG5_hgeom));
899 ph->nxt = hash->nxt;
900 hash->nxt = php->nxt;
901 php->nxt = 0;
902 }
903 else {
904 iph = ph->nxt;
905 iphp = php->nxt;
906 php->nxt = iph;
907 memset(ph,0,sizeof(MMG5_hgeom));
908 ph->nxt = hash->nxt;
909 hash->nxt = iphp;
910 }
911 return 1;
912 }
913 }
914 return 0;
915}
916
918int MMG5_hGet(MMG5_HGeom *hash,MMG5_int a,MMG5_int b,MMG5_int *ref,int16_t *tag) {
919 MMG5_hgeom *ph;
920 MMG5_int key;
921 MMG5_int ia,ib;
922
923 *tag = 0;
924 *ref = 0;
925
926 assert ( hash->siz );
927
928 ia = MG_MIN(a,b);
929 ib = MG_MAX(a,b);
930 key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
931 ph = &hash->geom[key];
932
933 if ( !ph->a ) return 0;
934 else if ( ph->a == ia && ph->b == ib ) {
935 *ref = ph->ref;
936 *tag = ph->tag;
937 return 1;
938 }
939 while ( ph->nxt ) {
940 ph = &hash->geom[ph->nxt];
941 if ( ph->a == ia && ph->b == ib ) {
942 *ref = ph->ref;
943 *tag = ph->tag;
944 return 1;
945 }
946 }
947 return 0;
948}
949
951int MMG5_hEdge(MMG5_pMesh mesh,MMG5_HGeom *hash,MMG5_int a,MMG5_int b,MMG5_int ref,int16_t tag) {
952 MMG5_hgeom *ph;
953 MMG5_int key;
954 MMG5_int ia,ib,j;
955
956 assert ( hash->siz );
957
958 ia = MG_MIN(a,b);
959 ib = MG_MAX(a,b);
960 key = (MMG5_KA*(int64_t)ia + MMG5_KB*(int64_t)ib) % hash->siz;
961 ph = &hash->geom[key];
962
963 if ( ph->a == ia && ph->b == ib )
964 return 1;
965 else if ( ph->a ) {
966 while ( ph->nxt ) {
967 ph = &hash->geom[ph->nxt];
968 if ( ph->a == ia && ph->b == ib ) return 1;
969 }
970 ph->nxt = hash->nxt;
971 ph = &hash->geom[hash->nxt];
972 ph->a = ia; ph->b = ib;
973 ph->ref = ref; ph->tag = tag;
974 hash->nxt = ph->nxt;
975 ph->nxt = 0;
976 if ( hash->nxt >= hash->max ) {
977 if ( mesh->info.ddebug )
978 fprintf(stderr,"\n ## Memory alloc problem (edge): %" MMG5_PRId "\n",hash->max);
980 "larger htab table",
981 fprintf(stderr," Exit program.\n");return 0;);
982 for (j=hash->nxt; j<hash->max; j++) hash->geom[j].nxt = j+1;
983 }
984 return 1;
985 }
986 /* insert new edge */
987 ph->a = ia; ph->b = ib;
988 ph->ref = ref; ph->tag = tag;
989 ph->nxt = 0;
990 return 1;
991}
992
994int MMG5_hNew(MMG5_pMesh mesh,MMG5_HGeom *hash,MMG5_int hsiz,MMG5_int hmax) {
995 MMG5_int k;
996
997 /* adjust hash table params */
998 hash->siz = hsiz + 1;
999 hash->max = hmax + 2;
1000 hash->nxt = hash->siz;
1001
1002 MMG5_ADD_MEM(mesh,(hash->max+1)*sizeof(MMG5_hgeom),"Edge hash table",return 0);
1003 MMG5_SAFE_CALLOC(hash->geom,(hash->max+1),MMG5_hgeom,return 0);
1004
1005 if ( !hash->geom ) {
1006 perror(" ## Memory problem: calloc");
1007 return 0;
1008 }
1009 for (k=hash->siz; k<hash->max; k++)
1010 hash->geom[k].nxt = k+1;
1011
1012 return 1;
1013}
1014
1023 MMG5_pTria pt;
1024 MMG5_pEdge pa;
1025 MMG5_Hash hash;
1026 MMG5_int edg,*adja,k,kk;
1027 int ier;
1028 int16_t tag;
1029 int8_t i,i1,i2;
1030
1031 /* if edges exist in mesh, hash special edges from existing field */
1032 if ( mesh->na ) {
1033 if ( !mesh->htab.geom ) {
1034 mesh->namax = MG_MAX((MMG5_int)(1.5*mesh->na),MMG3D_NAMAX);
1035 if ( !MMG5_hNew(mesh,&mesh->htab,mesh->na,3*mesh->namax) )
1036 return 0;
1037 }
1038 else {
1039 if ( abs(mesh->info.imprim) > 3 || mesh->info.ddebug ) {
1040 fprintf(stderr,"\n ## Warning: %s: no re-hash of edges of mesh. ",
1041 __func__);
1042 fprintf(stderr,"mesh->htab.geom must be freed to enforce analysis.\n");
1043 }
1045 mesh->na = 0;
1046 return 1;
1047 }
1048
1049 /* store initial edges */
1050 for (k=1; k<=mesh->na; k++) {
1051 pa = &mesh->edge[k];
1052 if ( !MMG5_hEdge(mesh,&mesh->htab,pa->a,pa->b,pa->ref,pa->tag) )
1053 return 0;
1054 }
1055
1056 /* now check triangles */
1057 for (k=1; k<=mesh->nt; k++) {
1058 pt = &mesh->tria[k];
1059 for (i=0; i<3; i++) {
1060 if( (pt->tag[i] & MG_PARBDY) && !(pt->tag[i] & MG_PARBDYBDY) ) continue;
1061 i1 = MMG5_inxt2[i];
1062 i2 = MMG5_iprv2[i];
1063 /* transfer non manifold tag to edges */
1064 if ( pt->tag[i] & MG_NOM ) {
1065 ier = MMG5_hTag(&mesh->htab,pt->v[i1],pt->v[i2],pt->edg[i],pt->tag[i]);
1066 if ( !ier ) {
1067 /* The edge is marked as non manifold but doesn't exist in the mesh */
1068 ier = MMG5_hEdge(mesh,&mesh->htab,pt->v[i1],pt->v[i2],pt->edg[i],pt->tag[i]);
1069 if ( !ier )
1070 return 0;
1071 }
1072 }
1073 MMG5_hGet(&mesh->htab,pt->v[i1],pt->v[i2],&edg,&tag);
1074 pt->edg[i] = edg;
1075
1076 /* If we use the nosurf option and the edge is required, we don't want
1077 * to detect it as an edge whose tag has been modified for the option */
1078 if ( mesh->info.nosurf && (tag & MG_REQ) )
1079 pt->tag[i] &= ~MG_NOSURF;
1080
1081 /* Store the edge tag inside the triangle */
1082 pt->tag[i] |= tag;
1083
1084 MMG5_hTag(&mesh->htab,pt->v[i1],pt->v[i2],edg,pt->tag[i]);
1085 }
1086 }
1088 mesh->na = 0;
1089 }
1090 /* else, infer special edges from information carried by triangles */
1091 else {
1092 if ( !mesh->adjt ) {
1093 memset(&hash,0x0,sizeof(MMG5_Hash));
1094 ier = MMG3D_hashTria(mesh,&hash);
1095 MMG5_DEL_MEM(mesh,hash.item);
1096 if ( !ier ) return 0;
1097 }
1098
1099 for (k=1; k<=mesh->nt; k++) {
1100 pt = &mesh->tria[k];
1101 adja = &mesh->adjt[3*(k-1)+1];
1102 for (i=0; i<3; i++) {
1103 if( (pt->tag[i] & MG_PARBDY) && !(pt->tag[i] & MG_PARBDYBDY) ) continue;
1104 kk = adja[i] / 3;
1105 if ( !kk || pt->tag[i] & MG_NOM )
1106 mesh->na++;
1107 else if ( (k < kk) && ( pt->edg[i] || pt->tag[i] ) ) mesh->na++;
1108 }
1109 }
1110
1111 if ( mesh->htab.geom )
1113
1114 mesh->namax = MG_MAX((MMG5_int)(1.5*mesh->na),MMG3D_NAMAX);
1115 if ( !MMG5_hNew(mesh,&mesh->htab,mesh->na,3*mesh->namax) )
1116 return 0;
1117
1118 mesh->na = 0;
1119
1120 /* build hash for edges */
1121 for (k=1; k<=mesh->nt; k++) {
1122 pt = &mesh->tria[k];
1123 adja = &mesh->adjt[3*(k-1)+1];
1124 for (i=0; i<3; i++) {
1125 if( (pt->tag[i] & MG_PARBDY) && !(pt->tag[i] & MG_PARBDYBDY) ) continue;
1126 i1 = MMG5_inxt2[i];
1127 i2 = MMG5_iprv2[i];
1128 kk = adja[i] / 3;
1129 if ( (!kk) || pt->tag[i] & MG_NOM ) {
1130 if ( pt->tag[i] & MG_NOM ) {
1131 if ( mesh->info.iso )
1132 pt->edg[i] = ( pt->edg[i] != 0 ) ? -MMG5_abs(pt->edg[i]) : mesh->info.isoref;
1133 }
1134 if ( !MMG5_hEdge(mesh,&mesh->htab,pt->v[i1],pt->v[i2],pt->edg[i],pt->tag[i]) )
1135 return 0;
1136 }
1137 else if ( k < kk && ( pt->edg[i] || pt->tag[i] ) ) {
1138 if ( !MMG5_hEdge(mesh,&mesh->htab,pt->v[i1],pt->v[i2],pt->edg[i],pt->tag[i]))
1139 return 0;
1140 }
1141 }
1142 }
1143 /* now check triangles */
1144 for (k=1; k<=mesh->nt; k++) {
1145 pt = &mesh->tria[k];
1146 for (i=0; i<3; i++) {
1147 if( (pt->tag[i] & MG_PARBDY) && !(pt->tag[i] & MG_PARBDYBDY) ) continue;
1148 i1 = MMG5_inxt2[i];
1149 i2 = MMG5_iprv2[i];
1150 MMG5_hGet(&mesh->htab,pt->v[i1],pt->v[i2],&edg,&tag);
1151 pt->edg[i] = edg;
1152 pt->tag[i] |= tag;
1153 }
1154 }
1155 }
1156 return 1;
1157}
1158
1176static inline
1177int MMG5_bdryTria(MMG5_pMesh mesh, MMG5_int ntmesh) {
1178 MMG5_pTetra pt,pt1;
1179 MMG5_pPrism pp;
1180 MMG5_pTria ptt;
1181 MMG5_pPoint ppt;
1182 MMG5_pxTetra pxt;
1183 MMG5_pxPrism pxpr;
1184 MMG5_Hash hash;
1185 MMG5_int ref,*adja,adj,k,ia,ib,ic,kt,ntinit;
1186 int tofree=0;
1187 int8_t i,j;
1188
1189 hash.item = NULL;
1190
1191 ntinit = mesh->nt;
1192
1193 if ( mesh->nprism && (ntmesh!=ntinit) ) {
1194 /* If a triangle at the interface between a prism and a tetra is not
1195 * provided, the hashtable is used to recover from the prism a boundary tria
1196 * created by tetra */
1197 if ( ! MMG5_hashNew(mesh,&hash,0.51*ntmesh,1.51*ntmesh) ) return 0;
1198 tofree=1;
1199 }
1200 else if ( mesh->nt ) {
1201 /* Hash given bdry triangles */
1202 if ( ! MMG5_hashNew(mesh,&hash,(MMG5_int)(0.51*mesh->nt),(MMG5_int)(1.51*mesh->nt)) ) return 0;
1203 tofree=1;
1204 }
1205
1206 for (k=1; k<=mesh->nt; k++) {
1207 ptt = &mesh->tria[k];
1208 if ( !MMG5_hashFace(mesh,&hash,ptt->v[0],ptt->v[1],ptt->v[2],k) ) {
1209 MMG5_DEL_MEM(mesh,hash.item);
1210 return 0;
1211 }
1212 for (i=0; i<3; i++) {
1213 ppt = &mesh->point[ptt->v[i]];
1214 if ( !mesh->info.iso ) ppt->tag |= MG_BDY;
1215 }
1216 }
1217
1218 /* Add boundary triangles stored on tetra */
1219 if ( ntmesh != ntinit ) {
1220 for (k=1; k<=mesh->ne; k++) {
1221 pt = &mesh->tetra[k];
1222 if ( !MG_EOK(pt) ) continue;
1223 adja = &mesh->adja[4*(k-1)+1];
1224 pxt = 0;
1225 if ( pt->xt ) pxt = &mesh->xtetra[pt->xt];
1226 for (i=0; i<4; i++) {
1227 adj = adja[i] / 4;
1228 pt1 = &mesh->tetra[adj];
1229
1230 if ( (!mesh->info.opnbdy) || (!mesh->xtetra) ) {
1231 /* Classic mode (no open bdy) or no info stored in xtetra (first
1232 * call) */
1233 if ( adj && ( pt->ref <= pt1->ref) ) continue;
1234 } else {
1235 /* info stored in xtetra and opnbdy mode */
1236 if ( adj && ( (pt->ref<pt1->ref) || (!pt->xt) ||
1237 (!(pxt->ftag[i] & MG_BDY)) || (!MG_GET(pxt->ori,i) ) ) )
1238 continue;
1239 }
1240
1241 ia = pt->v[MMG5_idir[i][0]];
1242 ib = pt->v[MMG5_idir[i][1]];
1243 ic = pt->v[MMG5_idir[i][2]];
1244
1245 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1246 if ( kt ) {
1247 /* Face is already stored */
1248 continue;
1249 }
1250 else if ( mesh->nprism ) {
1251 /* Update the list of boundary trias to be able to recover tria at the
1252 * interface between tet and prisms */
1253 if ( !MMG5_hashFace(mesh,&hash,ia,ib,ic,mesh->nt+1) ) {
1254 MMG5_DEL_MEM(mesh,hash.item);
1255 return 0;
1256 }
1257 }
1258
1259 /* face does not exists: add it in tria array */
1260 mesh->nt++;
1261 ptt = &mesh->tria[mesh->nt];
1262 ptt->v[0] = pt->v[MMG5_idir[i][0]];
1263 mesh->point[ptt->v[0]].tag |= MG_BDY;
1264 ptt->v[1] = pt->v[MMG5_idir[i][1]];
1265 mesh->point[ptt->v[1]].tag |= MG_BDY;
1266 ptt->v[2] = pt->v[MMG5_idir[i][2]];
1267 mesh->point[ptt->v[2]].tag |= MG_BDY;
1268
1269 /* the cc field is used to be able to recover the tetra (and its face)
1270 * from which comes a boundary triangle (when called by packmesh =>
1271 * mesh->nt=0 at the beginning of the function) */
1272 ptt->cc = 4*k + i;
1273
1274 if ( pxt ) {
1275 /* useful only when saving mesh or in ls mode */
1276 for( j = 0; j < 3; j++ ) {
1277 if ( pxt->tag[MMG5_iarf[i][j]] ) {
1278 ptt->tag[j] = pxt->tag[MMG5_iarf[i][j]];
1279 /* Remove redundant boundary tag */
1280 ptt->tag[j] &= ~MG_BDY;
1281 }
1282 if ( pxt->edg[MMG5_iarf[i][j]] )
1283 ptt->edg[j] = pxt->edg[MMG5_iarf[i][j]];
1284 }
1285 }
1286 if ( adj ) {
1287 if ( mesh->info.iso ) {
1288 /* Triangle at the interface between two tets is set to the user-defined ref if any, or else to mesh->info.isoref ref */
1289 if ( pxt && pxt->ftag[i] & MG_BDY )
1290 ptt->ref = pxt->ref[i];
1291 else if( MMG5_isLevelSet(mesh,pt->ref,pt1->ref) )
1292 ptt->ref = mesh->info.isoref;
1293 else
1294 ptt->ref = MG_MIN(pt->ref,pt1->ref);
1295 }
1296 /* useful only when saving mesh or in ls mode */
1297 else {
1298 /* Triangle at the interface between two tet is set its init ref or
1299 * to the min ref of the adja tetra */
1300 ptt->ref = pxt ? pxt->ref[i] : MG_MIN(pt->ref,pt1->ref);
1301 }
1302 }
1303 else {
1304 /* useful only when saving mesh */
1305 ptt->ref = pxt ? pxt->ref[i] : pt->ref;
1306 }
1307 }
1308 }
1309 }
1310
1311 /* Add boundary triangles stored on prisms */
1312 if ( mesh->nprism ) {
1313 for (k=1; k<=mesh->nprism; k++) {
1314 pp = &mesh->prism[k];
1315 if ( !MG_EOK(pp) ) continue;
1316
1317 adja = &mesh->adjapr[5*(k-1)+1];
1318 pxpr = 0;
1319 if ( pp->xpr ) pxpr = &mesh->xprism[pp->xpr];
1320
1321 for (i=0; i<2; i++) {
1322 adj = adja[i]/5;
1323
1324
1325 if ( adj < 0 ) {
1326 if ( !mesh->nt ) continue;
1327
1328 /* Tria at the interface of a prism and a tetra: mark it as required */
1329 ia = pp->v[0+i*3];
1330 ib = pp->v[1+i*3];
1331 ic = pp->v[2+i*3];
1332
1333 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1334 if ( !kt ) continue;
1335
1336 ptt = &mesh->tria[kt];
1337
1338 if ( !(ptt->tag[0] & MG_REQ) ) {
1339 ptt->tag[0] |= MG_REQ;
1340 ptt->tag[0] |= MG_NOSURF;
1341 }
1342 if ( !(ptt->tag[1] & MG_REQ) ) {
1343 ptt->tag[1] |= MG_REQ;
1344 ptt->tag[1] |= MG_NOSURF;
1345 }
1346 if ( !(ptt->tag[2] & MG_REQ) ) {
1347 ptt->tag[2] |= MG_REQ;
1348 ptt->tag[2] |= MG_NOSURF;
1349 }
1350
1351 continue;
1352 }
1353
1354 ref = mesh->prism[adj].ref;
1355 if ( adj && ( pp->ref <= ref) ) continue;
1356
1357 ia = pp->v[0+i*3];
1358 ib = pp->v[1+i*3];
1359 ic = pp->v[2+i*3];
1360
1361 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1362 if ( kt ) {
1363 continue;
1364 }
1365
1366 mesh->nt++;
1367
1368 ptt = &mesh->tria[mesh->nt];
1369 ptt->v[0] = ia;
1370 ptt->v[1] = ib;
1371 ptt->v[2] = ic;
1372 mesh->point[ptt->v[0]].tag |= MG_BDY;
1373 mesh->point[ptt->v[1]].tag |= MG_BDY;
1374 mesh->point[ptt->v[2]].tag |= MG_BDY;
1375
1376 /* the cc field is used to be able to recover the prism (and its face)
1377 * from which comes a boundary triangle (when called by packmesh =>
1378 * mesh->nt=0 at the beginning of the function) */
1379 ptt->cc = 5*k + i;
1380 if ( pxpr ) {
1381 /* useful only when saving mesh */
1382 if ( pxpr->tag[MMG5_iarf_pr[i][0]] ) ptt->tag[0] = pxpr->tag[MMG5_iarf_pr[i][0]];
1383 if ( pxpr->tag[MMG5_iarf_pr[i][1]] ) ptt->tag[1] = pxpr->tag[MMG5_iarf_pr[i][1]];
1384 if ( pxpr->tag[MMG5_iarf_pr[i][2]] ) ptt->tag[2] = pxpr->tag[MMG5_iarf_pr[i][2]];
1385 }
1386 if ( adj ) {
1387 if ( mesh->info.iso ) ptt->ref = mesh->info.isoref;
1388 /* useful only when saving mesh */
1389 else ptt->ref = pxpr ? pxpr->ref[i] : 0;
1390 }
1391 else {
1392 /* useful only when saving mesh */
1393 ptt->ref = pxpr ? pxpr->ref[i] : 0;
1394 }
1395 }
1396 }
1397 }
1398
1399 assert(mesh->nt==ntmesh);
1400
1401 if ( ntmesh != ntinit ) {
1402 /* set point tag */
1403 for (k=1; k<=mesh->nt; k++) {
1404 ptt = &mesh->tria[k];
1405 for (i=0; i<3; i++) {
1406 ppt = &mesh->point[ptt->v[i]];
1407 ppt->tag |= MG_BDY;
1408 }
1409 }
1410 }
1411
1412 if ( tofree ) MMG5_DEL_MEM(mesh,hash.item);
1413
1414 return 1;
1415}
1416
1440 MMG5_pTetra pt,pt1;
1441 MMG5_pPrism pp,pp1;
1442 MMG5_pTria ptt,pttnew;
1443 MMG5_int *adja,adj,k,kk,i,j,ntmesh;
1444 MMG5_int ia,ib,ic, nbl,nt,ntpres;
1445 int iface;
1446 MMG5_Hash hashElt, hashTri;
1447
1449 ntmesh = ntpres = 0;
1450 for (k=1; k<=mesh->ne; k++) {
1451 pt = &mesh->tetra[k];
1452 if ( !MG_EOK(pt) ) continue;
1453 adja = &mesh->adja[4*(k-1)+1];
1454 for (i=0; i<4; i++) {
1455 adj = adja[i];
1456
1457 if ( !adj ) {
1458 ++ntmesh;
1459 continue;
1460 }
1461 adj /= 4;
1462 pt1 = &mesh->tetra[adj];
1463 if ( pt->ref > pt1->ref )
1464 ++ntmesh;
1465 }
1466 }
1467
1468 if ( mesh->info.opnbdy && mesh->xtetra ) {
1469 /* We want to preserve internal triangle and we came from bdryBuild: we need
1470 * to count the preserved boudaries */
1471 for (k=1; k<=mesh->ne; k++) {
1472 pt = &mesh->tetra[k];
1473 if ( !MG_EOK(pt) || !pt->xt ) continue;
1474 adja = &mesh->adja[4*(k-1)+1];
1475 for (i=0; i<4; i++) {
1476 adj = adja[i];
1477
1478 if ( !adj ) continue;
1479
1480 adj /= 4;
1481 pt1 = &mesh->tetra[adj];
1482
1483 if ( pt->ref != pt1->ref ) continue;
1484
1485 if ( (mesh->xtetra[pt->xt].ftag[i] & MG_BDY) &&
1486 (MG_GET(mesh->xtetra[pt->xt].ori,i) ) ) ++ntpres;
1487 }
1488 }
1489 }
1490
1491 if ( mesh->nprism ) {
1492 for (k=1; k<=mesh->nprism; k++) {
1493 pp = &mesh->prism[k];
1494 if ( !MG_EOK(pp) ) continue;
1495
1496 adja = &mesh->adjapr[5*(k-1)+1];
1497 for (i=0; i<2; i++) {
1498 adj = adja[i];
1499
1500 if ( !adj ) {
1501 ++ntmesh;
1502 continue;
1503 }
1504 else if ( adj<0 ) {
1505 continue;
1506 }
1507
1508 adj /= 5;
1509 pp1 = &mesh->prism[adj];
1510 if ( pp->ref > pp1->ref) {
1511 ++ntmesh;
1512 }
1513 }
1514 }
1515
1516 /* Detect the triangles at the interface of the prisms and tetra (they have been
1517 * counted twice) */
1518 if ( ! MMG5_hashNew(mesh,&hashTri,0.51*ntmesh,1.51*ntmesh) ) return 0;
1519 for (k=1; k<=mesh->ne; k++) {
1520 pt = &mesh->tetra[k];
1521 if ( !MG_EOK(pt) ) continue;
1522
1523 adja = &mesh->adja[4*(k-1)+1];
1524 for (i=0; i<4; i++) {
1525 adj = adja[i];
1526 if ( adj ) continue;
1527
1528 ia = pt->v[MMG5_idir[i][0]];
1529 ib = pt->v[MMG5_idir[i][1]];
1530 ic = pt->v[MMG5_idir[i][2]];
1531 if ( !MMG5_hashFace(mesh,&hashTri,ia,ib,ic,5*k+i) ) {
1532 MMG5_DEL_MEM(mesh,hashTri.item);
1533 return 0;
1534 }
1535 }
1536 }
1537
1538 /* Fill the adjacency relationship between prisms and tetra (fill adjapr with
1539 * a negative value to mark this special faces) */
1540 for (k=1; k<=mesh->nprism; k++) {
1541 pp = &mesh->prism[k];
1542 if ( !MG_EOK(pp) ) continue;
1543
1544 adja = &mesh->adjapr[5*(k-1)+1];
1545 for (i=0; i<2; i++) {
1546 adj = adja[i];
1547 if ( adj ) continue;
1548
1549 ia = pp->v[MMG5_idir_pr[i][0]];
1550 ib = pp->v[MMG5_idir_pr[i][1]];
1551 ic = pp->v[MMG5_idir_pr[i][2]];
1552
1553 j = MMG5_hashGetFace(&hashTri,ia,ib,ic);
1554 if ( !j ) continue;
1555
1556 --ntmesh;
1557 adja[i] = -j;
1558 }
1559 }
1560 MMG5_DEL_MEM(mesh,hashTri.item);
1561 }
1562
1565 if ( mesh->nt ) {
1566 if ( ! MMG5_hashNew(mesh,&hashElt,0.51*ntmesh,1.51*ntmesh) ) return 0;
1567 // Hash the boundaries found in the mesh
1568 if ( mesh->info.opnbdy) {
1569 /* We want to keep the internal triangles: we must hash all the tetra faces */
1570 for (k=1; k<=mesh->ne; k++) {
1571 pt = &mesh->tetra[k];
1572 if ( !MG_EOK(pt) ) continue;
1573
1574 for (i=0; i<4; i++) {
1575 ia = pt->v[MMG5_idir[i][0]];
1576 ib = pt->v[MMG5_idir[i][1]];
1577 ic = pt->v[MMG5_idir[i][2]];
1578 if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0;
1579 }
1580 }
1581 } else {
1582 for (k=1; k<=mesh->ne; k++) {
1583 pt = &mesh->tetra[k];
1584 if ( !MG_EOK(pt) ) continue;
1585 adja = &mesh->adja[4*(k-1)+1];
1586 for (i=0; i<4; i++) {
1587 adj = adja[i];
1588 if ( !adj ) {
1589 ia = pt->v[MMG5_idir[i][0]];
1590 ib = pt->v[MMG5_idir[i][1]];
1591 ic = pt->v[MMG5_idir[i][2]];
1592 if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0;
1593 }
1594 adj /= 4;
1595
1596 pt1 = &mesh->tetra[adj];
1597 if ( pt->ref > pt1->ref ) {
1598 ia = pt->v[MMG5_idir[i][0]];
1599 ib = pt->v[MMG5_idir[i][1]];
1600 ic = pt->v[MMG5_idir[i][2]];
1601 if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,4*k+i) ) return 0;
1602 }
1603 }
1604 }
1605 }
1606 for (k=1; k<=mesh->nprism; k++) {
1607 pp = &mesh->prism[k];
1608 if ( !MG_EOK(pp) ) continue;
1609 adja = &mesh->adjapr[5*(k-1)+1];
1610 for (i=0; i<2; i++) {
1611 adj = adja[i];
1612 if ( !adj ) {
1613 ia = pp->v[MMG5_idir_pr[i][0]];
1614 ib = pp->v[MMG5_idir_pr[i][1]];
1615 ic = pp->v[MMG5_idir_pr[i][2]];
1616 if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0;
1617 }
1618 else if ( adj<0 ) continue;
1619
1620 adj /= 5;
1621
1622 pp1 = &mesh->prism[MMG5_abs(adj)];
1623 if ( pp->ref > pp1->ref ) {
1624 ia = pp->v[MMG5_idir_pr[i][0]];
1625 ib = pp->v[MMG5_idir_pr[i][1]];
1626 ic = pp->v[MMG5_idir_pr[i][2]];
1627 if ( !MMG5_hashFace(mesh,&hashElt,ia,ib,ic,5*k+i) ) return 0;
1628 }
1629 }
1630 }
1631
1632
1633 // Travel through the tria, delete those that are not in the hash tab or
1634 // that are stored more that once.
1635 nt=0; nbl=1;
1636
1637 if ( ! MMG5_hashNew(mesh,&hashTri,0.51*mesh->nt,1.51*mesh->nt) ) return 0;
1638
1639 for (k=1; k<=mesh->nt; k++) {
1640 ptt = &mesh->tria[k];
1641
1642 ia = ptt->v[0];
1643 ib = ptt->v[1];
1644 ic = ptt->v[2];
1645
1646 i = MMG5_hashGetFace(&hashElt,ia,ib,ic);
1647 j = MMG5_hashFace(mesh,&hashTri,ia,ib,ic,k);
1648
1649 ptt->cc = i;
1650
1651 if ( !j ) {
1652 MMG5_DEL_MEM(mesh,hashElt.item);
1653 MMG5_DEL_MEM(mesh,hashTri.item);
1654 return 0;
1655 }
1656 else if ( j > 0 ) {
1657 /* the face already exists in the tria table */
1658 continue;
1659 }
1660
1661 if ( !i ) {
1662 /* the triangle is not a boundary tri or a tri at the interface of two
1663 * subdomains with different references and the user don't ask to keep
1664 * it. */
1665 continue;
1666 }
1667
1668 if ( mesh->info.opnbdy ) {
1669 kk = i/4;
1670 iface = i%4;
1671 adj = mesh->adja[4*(kk-1)+1+iface];
1672 /* Check if we have found a triangle at the interface of 2 doms of same
1673 * ref */
1674 if ( adj && mesh->tetra[kk].ref == mesh->tetra[adj/4].ref ) {
1675 ++ntpres;
1676 }
1677 }
1678
1679 ++nt;
1680 if ( k!=nbl ) {
1681 pttnew = &mesh->tria[nbl];
1682 memcpy(pttnew,ptt,sizeof(MMG5_Tria));
1683 }
1684 ++nbl;
1685 }
1686 nbl = mesh->nt-nt;
1687 if ( nbl ) {
1688 fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries provided."
1689 " Ignored\n",__func__,nbl);
1690 MMG5_ADD_MEM(mesh,(-nbl)*sizeof(MMG5_Tria),"triangles",return 0);
1691 MMG5_SAFE_REALLOC(mesh->tria,mesh->nt+1,nt+1,MMG5_Tria,"triangles",return 0);
1692 mesh->nt = nt;
1693 }
1694 MMG5_DEL_MEM(mesh,hashElt.item);
1695 MMG5_DEL_MEM(mesh,hashTri.item);
1696 }
1697 ntmesh +=ntpres;
1698
1701 if ( ntpres && (mesh->info.imprim > 5 || mesh->info.ddebug) )
1702 printf(" %" MMG5_PRId " triangles between 2 tetrahdra with same"
1703 " references\n",ntpres);
1704
1705 if ( mesh->nt==ntmesh && !mesh->nprism ) {
1706 for (k=1; k<=mesh->nt; k++) {
1707 ptt = &mesh->tria[k];
1708 for (i=0; i<3; i++) {
1709 if ( !mesh->info.iso ) mesh->point[ptt->v[i]].tag |= MG_BDY;
1710 }
1711 }
1712 return 1;
1713 }
1714
1715 nbl = 0;
1716 if ( !mesh->nt ) {
1717 MMG5_ADD_MEM(mesh,(ntmesh+1)*sizeof(MMG5_Tria),"triangles",return 0);
1718 MMG5_SAFE_CALLOC(mesh->tria,ntmesh+1,MMG5_Tria,return 0);
1719 }
1720 else {
1721 assert((!mesh->nprism && ntmesh>mesh->nt)||(mesh->nprism && ntmesh>=mesh->nt));
1722 if ( ntmesh > mesh->nt ) {
1723 MMG5_ADD_MEM(mesh,(ntmesh-mesh->nt)*sizeof(MMG5_Tria),"triangles",return 0);
1724 MMG5_SAFE_RECALLOC(mesh->tria,mesh->nt+1,ntmesh+1,MMG5_Tria,"triangles",return 0);
1725 nbl = ntmesh-mesh->nt;
1726 }
1727 }
1728 if ( nbl && (mesh->info.imprim > 5 || mesh->info.ddebug) )
1729 fprintf(stderr,"\n ## Warning: %s: %" MMG5_PRId " extra boundaries founded\n",
1730 __func__,nbl);
1731
1732 /* Fill missing bdy triangles */
1733 return MMG5_bdryTria(mesh,ntmesh);
1734}
1735
1736
1745 MMG5_pTetra pt,pt1;
1746 MMG5_pPrism pp;
1747 MMG5_pTria ptt;
1748 MMG5_pxTetra pxt;
1749 MMG5_pxPrism pxp;
1750 MMG5_Hash hash;
1751 MMG5_int ref,*adja,adj,k,ia,ib,ic,kt,initedg[3];
1752 int j;
1753 int16_t tag,inittag[3];
1754 int8_t i,i1,i2;
1755
1756 if ( !mesh->nt ) return 1;
1757
1758 if ( mesh->xtetra ) {
1759 if ( abs(mesh->info.imprim) > 3 || mesh->info.ddebug ) {
1760 fprintf(stderr,"\n ## Error: %s: mesh->xtetra must be freed.\n",__func__);
1761 }
1762 return 0;
1763 }
1764 if ( mesh->xprism ) {
1765 if ( abs(mesh->info.imprim) > 3 || mesh->info.ddebug ) {
1766 fprintf(stderr,"\n ## Error: %s: mesh->xprism must be freed.\n",__func__);
1767 }
1768 return 0;
1769 }
1770
1771 if ( ! MMG5_hashNew(mesh,&hash,0.51*mesh->nt,1.51*mesh->nt) ) return 0;
1772 for (k=1; k<=mesh->nt; k++) {
1773 ptt = &mesh->tria[k];
1774 if ( !MMG5_hashFace(mesh,&hash,ptt->v[0],ptt->v[1],ptt->v[2],k) ) return 0;
1775 }
1776
1777 mesh->xt = 0;
1778 mesh->xtmax = mesh->ntmax;
1779 assert(mesh->xtmax);
1780
1781 MMG5_ADD_MEM(mesh,(mesh->xtmax+1)*sizeof(MMG5_xTetra),"boundary tetrahedra",
1782 fprintf(stderr," Exit program.\n");
1783 return 0);
1785
1786 /* assign references to tetras faces */
1787 if ( !mesh->info.opnbdy ) {
1788 for (k=1; k<=mesh->ne; k++) {
1789 pt = &mesh->tetra[k];
1790 if ( !MG_EOK(pt) ) continue;
1791 adja = &mesh->adja[4*(k-1)+1];
1792 for (i=0; i<4; i++) {
1793 adj = adja[i] / 4;
1794 pt1 = &mesh->tetra[adj];
1795 if ( !adj || ( pt->ref != pt1->ref) ) {
1796 ia = pt->v[MMG5_idir[i][0]];
1797 ib = pt->v[MMG5_idir[i][1]];
1798 ic = pt->v[MMG5_idir[i][2]];
1799 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1800 assert(kt);
1801 if ( !pt->xt ) {
1802 mesh->xt++;
1803 if ( mesh->xt > mesh->xtmax ) {
1805 "larger xtetra table",
1806 mesh->xt--;
1807 fprintf(stderr," Exit program.\n");return 0;);
1808 }
1809 pt->xt = mesh->xt;
1810 }
1811 ptt = &mesh->tria[kt];
1812 pxt = &mesh->xtetra[pt->xt];
1813 pxt->ref[i] = ptt->ref;
1814 pxt->ftag[i] |= MG_BDY;
1815 pxt->ftag[i] |= (ptt->tag[0] & ptt->tag[1] & ptt->tag[2]);
1816 }
1817 }
1818 }
1819 }
1820 else {
1821 /* Internal triangles preservations */
1822 for (k=1; k<=mesh->ne; k++) {
1823 pt = &mesh->tetra[k];
1824 if ( !MG_EOK(pt) ) continue;
1825
1826 for (i=0; i<4; i++) {
1827 ia = pt->v[MMG5_idir[i][0]];
1828 ib = pt->v[MMG5_idir[i][1]];
1829 ic = pt->v[MMG5_idir[i][2]];
1830 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1831
1832 if ( !kt ) continue;
1833
1834 if ( !pt->xt ) {
1835 mesh->xt++;
1836 if ( mesh->xt > mesh->xtmax ) {
1838 "larger xtetra table",
1839 mesh->xt--;
1840 fprintf(stderr," Exit program.\n");return 0;);
1841 }
1842 pt->xt = mesh->xt;
1843 }
1844 ptt = &mesh->tria[kt];
1845 pxt = &mesh->xtetra[mesh->xt];
1846 pxt->ref[i] = ptt->ref;
1847 pxt->ftag[i] |= MG_BDY;
1848 pxt->ftag[i] |= (ptt->tag[0] & ptt->tag[1] & ptt->tag[2]);
1849 }
1850 }
1851 }
1852
1853 if ( !mesh->info.opnbdy ) {
1854 for (k=1; k<=mesh->ne; k++) {
1855 pt = &mesh->tetra[k];
1856 if ( !MG_EOK(pt) ) continue;
1857 if ( !pt->xt ) continue;
1858 pxt = &mesh->xtetra[pt->xt];
1859 adja = &mesh->adja[4*(k-1)+1];
1860 for (i=0; i<4; i++) {
1861 adj = adja[i] / 4;
1862 pt1 = &mesh->tetra[adj];
1863 /* Set edge tag */
1864 if ( pxt->ftag[i] ) {
1865 if ( adj && (pt->ref == pt1->ref ) ) {
1866 continue;
1867 }
1868 else {
1869 ia = pt->v[MMG5_idir[i][0]];
1870 ib = pt->v[MMG5_idir[i][1]];
1871 ic = pt->v[MMG5_idir[i][2]];
1872 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1873 ptt = &mesh->tria[kt];
1874
1875 /* Set flag to know if tetra has the same orientation than the
1876 * triangle */
1877 if ( ptt->v[0] == ia && ptt->v[1] == ib && ptt->v[2] == ic ) {
1878 MG_SET(pxt->ori,i);
1879 for (j=0; j<3; j++) {
1880 tag = pxt->ftag[i] | ptt->tag[j];
1881 if ( tag ) {
1882 if ( !MMG5_settag(mesh,k,MMG5_iarf[i][j],tag,ptt->edg[j]) )
1883 return 0;
1884 }
1885 }
1886 }
1887 else
1888 MG_CLR(pxt->ori,i);
1889 }
1890 }
1891 }
1892 }
1893 }
1894 else {
1895 for (k=1; k<=mesh->ne; k++) {
1896 pt = &mesh->tetra[k];
1897 if ( !MG_EOK(pt) ) continue;
1898 if ( !pt->xt ) continue;
1899 pxt = &mesh->xtetra[pt->xt];
1900 adja = &mesh->adja[4*(k-1)+1];
1901 for (i=0; i<4; i++) {
1902 adj = adja[i] / 4;
1903 pt1 = &mesh->tetra[adj];
1904
1905 ia = pt->v[MMG5_idir[i][0]];
1906 ib = pt->v[MMG5_idir[i][1]];
1907 ic = pt->v[MMG5_idir[i][2]];
1908 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
1909
1910 if ( !kt ) continue;
1911
1912 ptt = &mesh->tria[kt];
1913
1914 /* Set flag to know if tetra has the same orientation than the triangle
1915 * + force the triangle numbering to match the tetra face numbering */
1916 if ( adj ) {
1917 for ( j=0; j<3; ++j ) {
1918 i1 = MMG5_inxt2[j];
1919 i2 = MMG5_inxt2[i1];
1920 if ( ptt->v[j]==ia && ptt->v[i1]==ib && ptt->v[i2]==ic )
1921 break;
1922 }
1923 if ( j<3 ) {
1924 MG_SET(pxt->ori,i);
1925 if ( j!=0 ) {
1926 /* Triangle vertices+tag/edg reordering */
1927 ptt->v[0] = ia;
1928 ptt->v[1] = ib;
1929 ptt->v[2] = ic;
1930
1931 inittag[0] = ptt->tag[0];
1932 inittag[1] = ptt->tag[1];
1933 inittag[2] = ptt->tag[2];
1934 ptt->tag[0] = inittag[j];
1935 ptt->tag[1] = inittag[i1];
1936 ptt->tag[2] = inittag[i2];
1937
1938 initedg[0] = ptt->edg[0];
1939 initedg[1] = ptt->edg[1];
1940 initedg[2] = ptt->edg[2];
1941 ptt->edg[0] = initedg[j];
1942 ptt->edg[1] = initedg[i1];
1943 ptt->edg[2] = initedg[i2];
1944 }
1945 }
1946 else {
1947 MG_CLR(pxt->ori,i);
1948 }
1949 }
1950 else MG_SET(pxt->ori,i);
1951
1952 /* Set edge tag */
1953 if ( pxt->ftag[i] ) {
1954 if ( adj && ( (pt->ref < pt1->ref) || !MG_GET(pxt->ori,i) ) ) {
1955 continue;
1956 }
1957 else {
1958 for (j=0; j<3; j++) {
1959 tag = pxt->ftag[i] | ptt->tag[j];
1960 if ( tag ) {
1961 if ( !MMG5_settag(mesh,k,MMG5_iarf[i][j],tag,ptt->edg[j]) )
1962 return 0;
1963 }
1964 }
1965 }
1966 }
1967 }
1968 }
1969 }
1970
1971 if ( !mesh->nprism ) {
1972 MMG5_DEL_MEM(mesh,hash.item);
1973 return 1;
1974 }
1975
1976 mesh->xpr = 0;
1977 MMG5_ADD_MEM(mesh,(mesh->nprism+1)*sizeof(MMG5_xPrism),"boundary prisms",
1978 fprintf(stderr," Exit program.\n");
1979 return 0);
1981
1982 /* assign references to prism faces */
1983 for (k=1; k<=mesh->nprism; k++) {
1984 pp = &mesh->prism[k];
1985 if ( !MG_EOK(pp) ) continue;
1986 adja = &mesh->adjapr[5*(k-1)+1];
1987 for (i=0; i<2; i++) {
1988 adj = adja[i] / 5;
1989 if ( adj<0 ) {
1990 ref = mesh->tetra[MMG5_abs(adj)].ref;
1991 } else {
1992 ref = mesh->prism[adj].ref;
1993 }
1994 if ( adj && (pp->ref == ref) ) continue;
1995
1996 ia = pp->v[MMG5_idir_pr[i][0]];
1997 ib = pp->v[MMG5_idir_pr[i][1]];
1998 ic = pp->v[MMG5_idir_pr[i][2]];
1999 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
2000 assert(kt);
2001 if ( !pp->xpr ) {
2002 mesh->xpr++;
2003 pp->xpr = mesh->xpr;
2004 }
2005 ptt = &mesh->tria[kt];
2006 pxp = &mesh->xprism[mesh->xpr];
2007 pxp->ref[i] = ptt->ref;
2008 pxp->ftag[i] |= MG_BDY;
2009 pxp->ftag[i] |= (ptt->tag[0] & ptt->tag[1] & ptt->tag[2]);
2010
2011 for (j=0; j<3; j++) {
2012 pxp->tag[MMG5_iarf[i][j]] |= pxp->ftag[i] | ptt->tag[j];
2013 pxp->edg[MMG5_iarf[i][j]] = ptt->edg[j];
2014 }
2015 }
2016 }
2017 MMG5_ADD_MEM(mesh,(mesh->xpr-mesh->nprism)*sizeof(MMG5_xPrism),"boundary prisms",
2018 fprintf(stderr," Exit program.\n");
2019 return 0);
2021 "boundary prisms",return 0);
2022
2023 MMG5_DEL_MEM(mesh,hash.item);
2024 return 1;
2025}
2026
2045 MMG5_pTetra pt;
2046 MMG5_pTria ptt;
2047 MMG5_pxTetra pxt;
2048 MMG5_Hash hash;
2049 MMG5_int ia,ib,ic,k,kt;
2050 int j;
2051 int16_t tag;
2052 int8_t i;
2053
2054 if ( !mesh->nt ) return 1;
2055 if ( !MMG5_hashNew(mesh,&hash,0.51*mesh->nt,1.51*mesh->nt) ) return 0;
2056 for (k=1; k<=mesh->nt; k++) {
2057 ptt = &mesh->tria[k];
2058 if ( !MMG5_hashFace(mesh,&hash,ptt->v[0],ptt->v[1],ptt->v[2],k) ) {
2059 MMG5_DEL_MEM(mesh,hash.item);
2060 return 0;
2061 }
2062 }
2063
2064 for (k=1; k<=mesh->ne; k++) {
2065 pt = &mesh->tetra[k];
2066 if ( !MG_EOK(pt) ) continue;
2067 if ( pt->tag & MG_REQ ) {
2068 mesh->point[mesh->tetra[k].v[0]].tag |= MG_REQ;
2069 mesh->point[mesh->tetra[k].v[1]].tag |= MG_REQ;
2070 mesh->point[mesh->tetra[k].v[2]].tag |= MG_REQ;
2071 mesh->point[mesh->tetra[k].v[3]].tag |= MG_REQ;
2072 if ( !MMG5_settag(mesh,k,0,MG_REQ,0) ) return 0;
2073 if ( !MMG5_settag(mesh,k,1,MG_REQ,0) ) return 0;
2074 if ( !MMG5_settag(mesh,k,2,MG_REQ,0) ) return 0;
2075 if ( !MMG5_settag(mesh,k,3,MG_REQ,0) ) return 0;
2076 if ( !MMG5_settag(mesh,k,4,MG_REQ,0) ) return 0;
2077 if ( !MMG5_settag(mesh,k,5,MG_REQ,0) ) return 0;
2078 mesh->point[mesh->tetra[k].v[0]].tag &= ~MG_NOSURF;
2079 mesh->point[mesh->tetra[k].v[1]].tag &= ~MG_NOSURF;
2080 mesh->point[mesh->tetra[k].v[2]].tag &= ~MG_NOSURF;
2081 mesh->point[mesh->tetra[k].v[3]].tag &= ~MG_NOSURF;
2082 if ( !MMG5_deltag(mesh,k,0,MG_NOSURF) ) return 0;
2083 if ( !MMG5_deltag(mesh,k,1,MG_NOSURF) ) return 0;
2084 if ( !MMG5_deltag(mesh,k,2,MG_NOSURF) ) return 0;
2085 if ( !MMG5_deltag(mesh,k,3,MG_NOSURF) ) return 0;
2086 if ( !MMG5_deltag(mesh,k,4,MG_NOSURF) ) return 0;
2087 if ( !MMG5_deltag(mesh,k,5,MG_NOSURF) ) return 0;
2088 }
2089
2090 if ( !pt->xt ) continue;
2091 pxt = &mesh->xtetra[pt->xt];
2092
2093 for (i=0; i<4; i++) {
2094 /* Set edge tag */
2095 if ( ! MG_GET(pxt->ori,i) ) continue;
2096 if ( pxt->ftag[i] & MG_BDY ) {
2097 ia = pt->v[MMG5_idir[i][0]];
2098 ib = pt->v[MMG5_idir[i][1]];
2099 ic = pt->v[MMG5_idir[i][2]];
2100 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
2101 assert(kt);
2102 ptt = &mesh->tria[kt];
2103 if ( pt->tag & MG_REQ ) {
2104 pxt->ftag[i] |= MG_REQ;
2105 ptt->tag[0] = MG_REQ;
2106 ptt->tag[1] = MG_REQ;
2107 ptt->tag[2] = MG_REQ;
2108 pxt->ftag[i] &= ~MG_NOSURF;
2109 ptt->tag[0] &= ~MG_NOSURF;
2110 ptt->tag[1] &= ~MG_NOSURF;
2111 ptt->tag[2] &= ~MG_NOSURF;
2112 }
2113 for ( j=0; j<3; j++ ) {
2114 tag = ptt->tag[j];
2115 if ( tag || ptt->edg[j] ) {
2116 if ( !MMG5_settag(mesh,k,MMG5_iarf[i][j],tag,ptt->edg[j]) )
2117 return 0;
2118 }
2119 }
2120 }
2121 }
2122 }
2123 MMG5_DEL_MEM(mesh,hash.item);
2124 return 1;
2125}
2126
2136 MMG5_pTetra pt,pt1;
2137 MMG5_pTria ptt;
2138 MMG5_Hash hash;
2139 MMG5_int *adja,adj,k,kt,ia,ib,ic,nf;
2140 int8_t i;
2141
2142 if ( !mesh->nt ) return 1;
2143
2144 /* store triangles temporarily */
2145 if ( !MMG5_hashNew(mesh,&hash,MG_MAX(0.51*mesh->nt,100),MG_MAX(1.51*mesh->nt,300)) )
2146 return 0;
2147
2148 for (k=1; k<=mesh->nt; k++) {
2149 ptt = &mesh->tria[k];
2150 if ( !MMG5_hashFace(mesh,&hash,ptt->v[0],ptt->v[1],ptt->v[2],k) ) {
2151 MMG5_DEL_MEM(mesh,hash.item);
2152 return 0;
2153 }
2154 }
2155
2156 /* check orientation */
2157 nf = 0;
2158 for (k=1; k<=mesh->ne; k++) {
2159 pt = &mesh->tetra[k];
2160 if ( !MG_EOK(pt) ) continue;
2161 adja = &mesh->adja[4*(k-1)+1];
2162 for (i=0; i<4; i++) {
2163 adj = adja[i] / 4;
2164 pt1 = &mesh->tetra[adj];
2165
2166 if ( adj && (pt->ref <= pt1->ref || pt->ref == MG_PLUS) )
2167 continue;
2168
2169 ia = pt->v[MMG5_idir[i][0]];
2170 ib = pt->v[MMG5_idir[i][1]];
2171 ic = pt->v[MMG5_idir[i][2]];
2172 kt = MMG5_hashGetFace(&hash,ia,ib,ic);
2173 if ( kt ) {
2174 /* check orientation */
2175 ptt = &mesh->tria[kt];
2176 if ( ptt->v[0] == ia && ptt->v[1] == ib && ptt->v[2] == ic )
2177 continue;
2178 else {
2179 ptt->v[0] = ia;
2180 ptt->v[1] = ib;
2181 ptt->v[2] = ic;
2182 nf++;
2183 }
2184 }
2185 }
2186 }
2187 if ( mesh->info.ddebug && nf > 0 )
2188 fprintf(stdout," ## %" MMG5_PRId " faces reoriented\n",nf);
2189
2190 MMG5_DEL_MEM(mesh,hash.item);
2191
2192 return 1;
2193}
int ier
if(!ier) exit(EXIT_FAILURE)
MMG5_pMesh * mesh
int MMG5_settag(MMG5_pMesh mesh, MMG5_int start, int ia, int16_t tag, int edg)
Definition: boulep_3d.c:1228
int MMG5_deltag(MMG5_pMesh mesh, MMG5_int start, int ia, int16_t tag)
Definition: boulep_3d.c:1344
void MMG3D_coquilFaceSecondLoopInit(MMG5_pMesh mesh, MMG5_int piv, int8_t *iface, int8_t *ia, int64_t *list, int *ilist, MMG5_int *it1, MMG5_int *pradj, MMG5_int *adj)
Definition: boulep_3d.c:1787
int MMG3D_coquilFaceFirstLoop(MMG5_pMesh mesh, MMG5_int start, MMG5_int na, MMG5_int nb, int8_t iface, int8_t ia, int64_t *list, int *ilist, MMG5_int *it1, MMG5_int *it2, MMG5_int *piv, MMG5_int *adj, int8_t *hasadja, int *nbdy, int silent)
Definition: boulep_3d.c:1686
void MMG5_coquilFaceErrorMessage(MMG5_pMesh mesh, MMG5_int k1, MMG5_int k2)
Definition: boulep_3d.c:1613
int16_t MMG5_openCoquilTravel(MMG5_pMesh mesh, MMG5_int na, MMG5_int nb, MMG5_int *adj, MMG5_int *piv, int8_t *iface, int8_t *i)
Definition: boulep_3d.c:2011
int MMG5_boulernm(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int start, int ip, MMG5_int *ng, MMG5_int *nr, MMG5_int *nm)
Definition: boulep_3d.c:455
MMG5_int MMG5_hashFace(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int ia, MMG5_int ib, MMG5_int ic, MMG5_int k)
Definition: hash.c:286
int MMG5_mmgHashTria(MMG5_pMesh mesh, MMG5_int *adjt, MMG5_Hash *hash, int chkISO)
Definition: hash.c:58
int MMG5_hashNew(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int hsiz, MMG5_int hmax)
Definition: hash.c:526
int MMG5_paktet(MMG5_pMesh mesh)
Definition: hash_3d.c:50
int MMG3D_hashTetra(MMG5_pMesh mesh, int pack)
Definition: hash_3d.c:122
int MMG5_hGeom(MMG5_pMesh mesh)
Definition: hash_3d.c:1022
int MMG5_hNew(MMG5_pMesh mesh, MMG5_HGeom *hash, MMG5_int hsiz, MMG5_int hmax)
Definition: hash_3d.c:994
int MMG3D_hashPrism(MMG5_pMesh mesh)
Definition: hash_3d.c:238
static int MMG5_bdryTria(MMG5_pMesh mesh, MMG5_int ntmesh)
Definition: hash_3d.c:1177
MMG5_int MMG5_hashGetFace(MMG5_Hash *hash, MMG5_int ia, MMG5_int ib, MMG5_int ic)
Definition: hash_3d.c:84
int MMG5_hGet(MMG5_HGeom *hash, MMG5_int a, MMG5_int b, MMG5_int *ref, int16_t *tag)
Definition: hash_3d.c:918
int MMG5_chkBdryTria(MMG5_pMesh mesh)
Definition: hash_3d.c:1439
int MMG5_bdrySet(MMG5_pMesh mesh)
Definition: hash_3d.c:1744
static int MMG5_setEdgeNmTag(MMG5_pMesh mesh, MMG5_Hash *hash)
Definition: hash_3d.c:456
int8_t ddb
Definition: mmg3d1_delone.c:42
int MMG5_bdryUpdate(MMG5_pMesh mesh)
Definition: hash_3d.c:2044
int MMG5_setNmTag(MMG5_pMesh mesh, MMG5_Hash *hash)
Definition: hash_3d.c:727
int MMG3D_hashTria(MMG5_pMesh mesh, MMG5_Hash *hash)
Definition: hash_3d.c:749
static int MMG5_setVertexNmTag(MMG5_pMesh mesh)
Definition: hash_3d.c:624
int MMG5_hashPop(MMG5_Hash *hash, MMG5_int a, MMG5_int b)
Definition: hash_3d.c:761
int MMG5_hTag(MMG5_HGeom *hash, MMG5_int a, MMG5_int b, MMG5_int ref, int16_t tag)
Definition: hash_3d.c:825
int MMG5_hPop(MMG5_HGeom *hash, MMG5_int a, MMG5_int b, MMG5_int *ref, int16_t *tag)
Definition: hash_3d.c:858
int MMG5_bdryPerm(MMG5_pMesh mesh)
Definition: hash_3d.c:2135
int MMG5_hEdge(MMG5_pMesh mesh, MMG5_HGeom *hash, MMG5_int a, MMG5_int b, MMG5_int ref, int16_t tag)
Definition: hash_3d.c:951
#define MMG5_KC
Definition: hash_3d.c:38
API headers for the mmg3d library.
#define MMG3D_LMAX
Definition: libmmg3d.h:58
MMG5_int MMG3D_indPt(MMG5_pMesh mesh, MMG5_int kp)
Definition: tools_3d.c:916
int MMG3D_delElt(MMG5_pMesh mesh, MMG5_int iel)
Definition: zaldy_3d.c:122
#define MMG3D_NAMAX
static const int8_t MMG5_iarf[4][3]
iarf[i]: edges of face opposite to vertex i
static const uint8_t MMG5_iare[6][2]
vertices of extremities of the edges of the tetra
static const uint8_t MMG5_iarf_pr[5][5]
iarf[i]: edges of face i for a prism
static const uint8_t MMG5_idir_pr[5][4]
idir[i]: vertices of face i for a prism
static const uint8_t MMG5_idir[4][3]
idir[i]: vertices of face opposite to vertex i
#define MG_PLUS
Definition: libmmgtypes.h:71
int MMG5_isLevelSet(MMG5_pMesh mesh, MMG5_int ref0, MMG5_int ref1)
Definition: mmg2.c:471
#define MG_REQ
#define MMG5_SAFE_CALLOC(ptr, size, type, law)
#define MG_EOK(pt)
#define MG_CLR(flag, bit)
#define MG_PARBDY
#define MG_MIN(a, b)
#define MG_MAX(a, b)
#define MMG5_GAP
#define MMG5_ADD_MEM(mesh, size, message, law)
#define MMG5_SAFE_REALLOC(ptr, prevSize, newSize, type, message, law)
static const uint8_t MMG5_iprv2[3]
#define MMG5_KB
#define MMG5_TAB_RECALLOC(mesh, ptr, initSize, wantedGap, type, message, law)
#define MG_GET(flag, bit)
static const uint8_t MMG5_inxt2[6]
#define MG_VOK(ppt)
#define MG_CRN
#define MG_BDY
#define MG_NOSURF
#define MG_NOM
#define MMG5_KA
#define MG_PARBDYBDY
#define MMG5_SAFE_FREE(ptr)
#define MMG5_DEL_MEM(mesh, ptr)
#define MG_SET(flag, bit)
#define MMG5_SAFE_RECALLOC(ptr, prevSize, newSize, type, message, law)
Structure to store edges of a MMG mesh.
Definition: libmmgtypes.h:305
MMG5_int b
Definition: libmmgtypes.h:306
MMG5_int ref
Definition: libmmgtypes.h:307
int16_t tag
Definition: libmmgtypes.h:310
MMG5_int a
Definition: libmmgtypes.h:306
Hash table to store geometric edges.
Definition: libmmgtypes.h:574
MMG5_int max
Definition: libmmgtypes.h:576
MMG5_hgeom * geom
Definition: libmmgtypes.h:575
MMG5_int siz
Definition: libmmgtypes.h:576
MMG5_int nxt
Definition: libmmgtypes.h:576
Identic as MMG5_HGeom but use MMG5_hedge to store edges instead of MMG5_hgeom (memory economy).
Definition: libmmgtypes.h:595
MMG5_int nxt
Definition: libmmgtypes.h:596
MMG5_hedge * item
Definition: libmmgtypes.h:597
MMG5_int siz
Definition: libmmgtypes.h:596
int8_t iso
Definition: libmmgtypes.h:534
int8_t ddebug
Definition: libmmgtypes.h:532
MMG5_int isoref
Definition: libmmgtypes.h:521
uint8_t nosurf
Definition: libmmgtypes.h:546
MMG mesh structure.
Definition: libmmgtypes.h:605
MMG5_int ntmax
Definition: libmmgtypes.h:612
MMG5_Info info
Definition: libmmgtypes.h:651
MMG5_int * adjapr
Definition: libmmgtypes.h:632
MMG5_int xt
Definition: libmmgtypes.h:620
MMG5_int ne
Definition: libmmgtypes.h:612
MMG5_pPoint point
Definition: libmmgtypes.h:641
MMG5_int xpr
Definition: libmmgtypes.h:620
MMG5_int * adja
Definition: libmmgtypes.h:624
MMG5_pPrism prism
Definition: libmmgtypes.h:645
MMG5_int xtmax
Definition: libmmgtypes.h:612
MMG5_int nemax
Definition: libmmgtypes.h:612
MMG5_int nenil
Definition: libmmgtypes.h:622
MMG5_HGeom htab
Definition: libmmgtypes.h:650
MMG5_int namax
Definition: libmmgtypes.h:612
MMG5_int base
Definition: libmmgtypes.h:616
MMG5_pTetra tetra
Definition: libmmgtypes.h:643
MMG5_pxPrism xprism
Definition: libmmgtypes.h:646
MMG5_int nt
Definition: libmmgtypes.h:612
MMG5_pTria tria
Definition: libmmgtypes.h:647
MMG5_int np
Definition: libmmgtypes.h:612
MMG5_pEdge edge
Definition: libmmgtypes.h:649
MMG5_pxTetra xtetra
Definition: libmmgtypes.h:644
MMG5_int nprism
Definition: libmmgtypes.h:613
MMG5_int na
Definition: libmmgtypes.h:612
MMG5_int * adjt
Definition: libmmgtypes.h:628
Structure to store points of a MMG mesh.
Definition: libmmgtypes.h:270
int16_t tag
Definition: libmmgtypes.h:284
MMG5_int flag
Definition: libmmgtypes.h:282
MMG5_int v[6]
Definition: libmmgtypes.h:463
MMG5_int ref
Definition: libmmgtypes.h:464
MMG5_int xpr
Definition: libmmgtypes.h:467
MMG5_int v[4]
Definition: libmmgtypes.h:403
MMG5_int xt
Definition: libmmgtypes.h:407
MMG5_int ref
Definition: libmmgtypes.h:404
int16_t tag
Definition: libmmgtypes.h:410
MMG5_int edg[3]
Definition: libmmgtypes.h:339
int16_t tag[3]
Definition: libmmgtypes.h:342
MMG5_int ref
Definition: libmmgtypes.h:335
MMG5_int v[3]
Definition: libmmgtypes.h:334
MMG5_int cc
Definition: libmmgtypes.h:337
Used to hash edges (memory economy compared to MMG5_hgeom).
Definition: libmmgtypes.h:584
MMG5_int s
Definition: libmmgtypes.h:587
MMG5_int b
Definition: libmmgtypes.h:585
MMG5_int a
Definition: libmmgtypes.h:585
MMG5_int k
Definition: libmmgtypes.h:586
MMG5_int nxt
Definition: libmmgtypes.h:585
Cell of the hash table of geom edges.
Definition: libmmgtypes.h:562
MMG5_int ref
Definition: libmmgtypes.h:565
int16_t tag
Definition: libmmgtypes.h:567
MMG5_int b
Definition: libmmgtypes.h:564
MMG5_int nxt
Definition: libmmgtypes.h:566
MMG5_int a
Definition: libmmgtypes.h:563
Structure to store the surface prism of a MMG mesh.
Definition: libmmgtypes.h:477
int16_t ftag[5]
Definition: libmmgtypes.h:484
int16_t tag[9]
Definition: libmmgtypes.h:486
MMG5_int edg[9]
Definition: libmmgtypes.h:480
MMG5_int ref[5]
Definition: libmmgtypes.h:478
Structure to store the surface tetrahedra of a MMG mesh.
Definition: libmmgtypes.h:418
int16_t ftag[4]
Definition: libmmgtypes.h:423
int16_t tag[6]
Definition: libmmgtypes.h:425
int8_t ori
Definition: libmmgtypes.h:427
MMG5_int edg[6]
Definition: libmmgtypes.h:421
MMG5_int ref[4]
Definition: libmmgtypes.h:419