Mmg
Simplicial remeshers (mesh adaptation, isovalue discretization, lagrangian movement)
mmg3d1.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
39#include "libmmg3d.h"
41#include "mmgexterns_private.h"
43
44extern int8_t ddb;
45
60 uint16_t tag,MMG5_int nmref,MMG5_int edgref,
61 double no1[3],double no2[3],double to[3]) {
62
63 if ( MG_EDG_OR_NOM(tag) ) {
64 ppt->ref = nmref;
65 }
66 else {
67 ppt->ref = edgref;
68 }
69
70 MMG5_pxPoint pxp = &mesh->xpoint[ppt->xp];
71 if ( tag & MG_NOM ){
72 memcpy(pxp->n1,no1,3*sizeof(double));
73 memcpy(ppt->n,to,3*sizeof(double));
74 return;
75 }
76 else if ( tag & MG_GEO ) {
77 memcpy(pxp->n1,no1,3*sizeof(double));
78 memcpy(pxp->n2,no2,3*sizeof(double));
79 memcpy(ppt->n,to,3*sizeof(double));
80 return;
81 }
82 else if ( tag & MG_REF ) {
83 memcpy(pxp->n1,no1,3*sizeof(double));
84 memcpy(ppt->n,to,3*sizeof(double));
85 return;
86 }
87 else {
88 memcpy(pxp->n1,no1,3*sizeof(double));
89 return;
90 }
91}
92
102void MMG5_tet2tri(MMG5_pMesh mesh,MMG5_int k,int8_t ie,MMG5_Tria *ptt) {
103 MMG5_pTetra pt;
104 MMG5_pxTetra pxt;
105 int8_t i;
106
107 assert ( 0<=ie && ie<4 && "unexpected local face idx");
108
109 pt = &mesh->tetra[k];
110 memset(ptt,0,sizeof(MMG5_Tria));
111 ptt->v[0] = pt->v[MMG5_idir[ie][0]];
112 ptt->v[1] = pt->v[MMG5_idir[ie][1]];
113 ptt->v[2] = pt->v[MMG5_idir[ie][2]];
114 if ( pt->xt ) {
115 pxt = &mesh->xtetra[pt->xt];
116 ptt->ref = pxt->ref[ie];
117 for (i=0; i<3; i++) {
118 ptt->edg[i] = pxt->edg[MMG5_iarf[ie][i]];
119 ptt->tag[i] = pxt->tag[MMG5_iarf[ie][i]];
120 }
121 }
122 else {
123 for (i=0; i<3; i++) {
124 ptt->edg[i] = 0;
125 ptt->tag[i] = 0;
126 }
127 }
128}
129
140int MMG3D_dichoto(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,MMG5_int *vx) {
141 MMG5_pTetra pt;
142 MMG5_pPoint pa,pb,ps;
143 double o[6][3],p[6][3];
144 float to,tp,t;
145 int ier,it,maxit;
146 MMG5_int ia,ib;
147 int8_t i;
148
149 ier = 1;
150 pt = &mesh->tetra[k];
151
152 /* get point on surface and along segment for edge split */
153 for (i=0; i<6; i++) {
154 memset(p[i],0,3*sizeof(double));
155 memset(o[i],0,3*sizeof(double));
156 if ( vx[i] > 0 ) {
157 ia = pt->v[MMG5_iare[i][0]];
158 ib = pt->v[MMG5_iare[i][1]];
159 pa = &mesh->point[ia];
160 pb = &mesh->point[ib];
161 ps = &mesh->point[vx[i]];
162 o[i][0] = 0.5 * (pa->c[0] + pb->c[0]);
163 o[i][1] = 0.5 * (pa->c[1] + pb->c[1]);
164 o[i][2] = 0.5 * (pa->c[2] + pb->c[2]);
165 p[i][0] = ps->c[0];
166 p[i][1] = ps->c[1];
167 p[i][2] = ps->c[2];
168 }
169 }
170 maxit = 4;
171 it = 0;
172 tp = 1.0;
173 to = 0.0;
174 do {
175 /* compute new position */
176 t = 0.5 * (tp + to);
177 for (i=0; i<6; i++) {
178 if ( vx[i] > 0 ) {
179 ps = &mesh->point[vx[i]];
180 ps->c[0] = o[i][0] + t*(p[i][0] - o[i][0]);
181 ps->c[1] = o[i][1] + t*(p[i][1] - o[i][1]);
182 ps->c[2] = o[i][2] + t*(p[i][2] - o[i][2]);
183 }
184 }
185 switch (pt->flag) {
186 case 1: case 2: case 4: case 8: case 16: case 32:
187 ier = MMG3D_split1_sim(mesh,met,k,vx);
188 break;
189 case 48: case 24: case 40: case 6: case 34: case 36:
190 case 20: case 5: case 17: case 9: case 3: case 10:
191 ier =MMG3D_split2sf_sim(mesh,met,k,vx);
192 break;
193 case 33: case 18: case 12:
194 ier = MMG3D_split2_sim(mesh,met,k,vx);
195 break;
196 case 11: case 21: case 38: case 56:
197 ier = MMG3D_split3_sim(mesh,met,k,vx);
198 break;
199 case 7: case 25: case 42: case 52:
200 ier =MMG3D_split3cone_sim(mesh,met,k,vx);
201 break;
202 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
203 case 14: case 49: case 50: case 44: case 41:
204 ier = MMG3D_split3op_sim(mesh,met,k,vx);
205 break;
206 case 23: case 29: case 53: case 60: case 57: case 58:
207 case 27: case 15: case 43: case 39: case 54: case 46:
208 ier = MMG3D_split4sf_sim(mesh,met,k,vx);
209 break;
210 case 30: case 45: case 51:
211 ier = MMG3D_split4op_sim(mesh,met,k,vx);
212 break;
213 case 62: case 61: case 59: case 55: case 47: case 31:
214 ier = MMG3D_split5_sim(mesh,met,k,vx);
215 break;
216 case 63:
217 ier = MMG3D_split6_sim(mesh,met,k,vx);
218 break;
219 }
220 if ( ier )
221 to = t;
222 else
223 tp = t;
224 }
225 while ( ++it < maxit );
226 /* restore coords of last valid pos. */
227 if ( !ier ) {
228 t = to;
229 for (i=0; i<6; i++) {
230 if ( vx[i] > 0 ) {
231 ps = &mesh->point[vx[i]];
232 ps->c[0] = o[i][0] + t*(p[i][0] - o[i][0]);
233 ps->c[1] = o[i][1] + t*(p[i][1] - o[i][1]);
234 ps->c[2] = o[i][2] + t*(p[i][2] - o[i][2]);
235 }
236 }
237 }
238
239 /* For very ill-shaped elements, we can have no valid position */
240 switch (pt->flag) {
241 case 1: case 2: case 4: case 8: case 16: case 32:
242 ier = MMG3D_split1_sim(mesh,met,k,vx);
243 break;
244 case 48: case 24: case 40: case 6: case 34: case 36:
245 case 20: case 5: case 17: case 9: case 3: case 10:
246 ier =MMG3D_split2sf_sim(mesh,met,k,vx);
247 break;
248 case 33: case 18: case 12:
249 ier = MMG3D_split2_sim(mesh,met,k,vx);
250 break;
251 case 11: case 21: case 38: case 56:
252 ier = MMG3D_split3_sim(mesh,met,k,vx);
253 break;
254 case 7: case 25: case 42: case 52:
255 ier =MMG3D_split3cone_sim(mesh,met,k,vx);
256 break;
257 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
258 case 14: case 49: case 50: case 44: case 41:
259 ier = MMG3D_split3op_sim(mesh,met,k,vx);
260 break;
261 case 23: case 29: case 53: case 60: case 57: case 58:
262 case 27: case 15: case 43: case 39: case 54: case 46:
263 ier = MMG3D_split4sf_sim(mesh,met,k,vx);
264 break;
265 case 30: case 45: case 51:
266 ier = MMG3D_split4op_sim(mesh,met,k,vx);
267 break;
268 case 62: case 61: case 59: case 55: case 47: case 31:
269 ier = MMG3D_split5_sim(mesh,met,k,vx);
270 break;
271 case 63:
272 ier = MMG3D_split6_sim(mesh,met,k,vx);
273 break;
274 }
275 return ier;
276}
277
293int MMG3D_dichoto1b(MMG5_pMesh mesh,MMG5_pSol met,int64_t *list,int ret,MMG5_int ip) {
294 MMG5_pTetra pt;
295 MMG5_pPoint p0,p1,ppt;
296 int it,maxit;
297 MMG5_int iel,np,nq;
298 double m[3],o[3],tp,to,t;
299 int8_t ia,ier;
300
301 iel = list[0] / 6;
302 ia = list[0] % 6;
303 pt = &mesh->tetra[iel];
304
305 np = pt->v[MMG5_iare[ia][0]];
306 nq = pt->v[MMG5_iare[ia][1]];
307 p0 = &mesh->point[np];
308 p1 = &mesh->point[nq];
309
310 /* initial coordinates for new point */
311 ppt = &mesh->point[ip];
312 o[0] = ppt->c[0];
313 o[1] = ppt->c[1];
314 o[2] = ppt->c[2];
315
316 /* midpoint along edge */
317 m[0] = 0.5*(p0->c[0] + p1->c[0]);
318 m[1] = 0.5*(p0->c[1] + p1->c[1]);
319 m[2] = 0.5*(p0->c[2] + p1->c[2]);
320
321 maxit = 4;
322 it = 0;
323 tp = 1.0;
324 to = 0.0;
325 do {
326 t = 0.5*(to + tp);
327 ppt->c[0] = m[0] + t*(o[0]-m[0]);
328 ppt->c[1] = m[1] + t*(o[1]-m[1]);
329 ppt->c[2] = m[2] + t*(o[2]-m[2]);
330
331 ier = MMG3D_simbulgept(mesh,met,list,ret,ip);
332 assert ( (!mesh->info.ddebug) || (mesh->info.ddebug && ier != -1) );
333 if ( ier == 1 )
334 to = t;
335 else
336 tp = t;
337 }
338 while ( ++it < maxit );
339 if ( !ier ) t = to;
340
341 /* For very ill-shaped elements, we can have novalid position */
342 ppt->c[0] = m[0] + t*(o[0]-m[0]);
343 ppt->c[1] = m[1] + t*(o[1]-m[1]);
344 ppt->c[2] = m[2] + t*(o[2]-m[2]);
345
346 return MMG3D_simbulgept(mesh,met,list,ret,ip) ;
347}
348
363int8_t MMG5_chkedg(MMG5_pMesh mesh,MMG5_Tria *pt,int8_t ori, double hmax,
364 double hausd, int locPar) {
365 MMG5_pPoint p[3];
366 MMG5_xPoint *pxp;
367// MMG5_pPar par;
368 double n[3][3],t[3][3],nt[3],*n1,*n2,t1[3],t2[3];
369 double ps,ps2,ux,uy,uz,ll,il,alpha,dis,hma2;
370 MMG5_int ia,ib,ic;//l,info;
371 int8_t i,i1,i2;
372 static int8_t mmgWarn0 = 0, mmgWarn1 = 0;
373
374 ia = pt->v[0];
375 ib = pt->v[1];
376 ic = pt->v[2];
377 p[0] = &mesh->point[ia];
378 p[1] = &mesh->point[ib];
379 p[2] = &mesh->point[ic];
380 pt->flag = 0;
381
382 n1 = n2 = NULL;
383
384 /* normal recovery */
385 for (i=0; i<3; i++) {
386 if ( MG_SIN(p[i]->tag) ) {
387 MMG5_nortri(mesh,pt,n[i]);
388 if(!ori) {
389 n[i][0] *= -1.0;
390 n[i][1] *= -1.0;
391 n[i][2] *= -1.0;
392 }
393 }
394 else if ( (p[i]->tag & MG_NOM) || (p[i]->tag & MG_OPNBDY) ){
395 MMG5_nortri(mesh,pt,n[i]);
396 if(!ori) {
397 n[i][0] *= -1.0;
398 n[i][1] *= -1.0;
399 n[i][2] *= -1.0;
400 }
401 assert(p[i]->xp);
402 memcpy(&t[i],p[i]->n,3*sizeof(double));
403 }
404 else {
405 assert(p[i]->xp);
406 pxp = &mesh->xpoint[p[i]->xp];
407 if ( MG_EDG(p[i]->tag) ) {
408 memcpy(&t[i],p[i]->n,3*sizeof(double));
409 MMG5_nortri(mesh,pt,nt);
410 if(!ori) {
411 nt[0] *= -1.0;
412 nt[1] *= -1.0;
413 nt[2] *= -1.0;
414 }
415 ps = pxp->n1[0]*nt[0] + pxp->n1[1]*nt[1] + pxp->n1[2]*nt[2];
416 ps2 = pxp->n2[0]*nt[0] + pxp->n2[1]*nt[1] + pxp->n2[2]*nt[2];
417 if ( fabs(ps) > fabs(ps2) )
418 memcpy(&n[i],pxp->n1,3*sizeof(double));
419 else
420 memcpy(&n[i],pxp->n2,3*sizeof(double));
421 }
422 else
423 memcpy(&n[i],pxp->n1,3*sizeof(double));
424 }
425 }
426
427 /* analyze edges */
428 for (i=0; i<3; i++) {
429 i1 = MMG5_inxt2[i];
430 i2 = MMG5_iprv2[i];
431
432 /* local parameters at vertices */
433 /* if ( mesh->info.parTyp & MG_Vert ) { */
434
435 /* if ( p[i1]->ref == p[i2]->ref ) { */
436 /* for ( l=0; l<mesh->info.npar; ++l) { */
437 /* par = &mesh->info.par[l]; */
438
439 /* if ( par->elt != MMG5_Vertex || par->ref != p[i1]->ref ) continue; */
440
441 /* if ( !locPar ) { */
442 /* hausd = par->hausd; */
443 /* hmax = par->hmax; */
444 /* } */
445 /* else { */
446 /* hausd = MG_MIN(hausd, par->hausd); */
447 /* hmax = MG_MIN(hmax, par->hmax); */
448 /* } */
449 /* break; */
450 /* } */
451 /* } */
452 /* else { */
453 /* l = 0; */
454 /* info = -1000; */
455 /* do { */
456 /* if ( info >= 0 ) break; */
457
458 /* par = &mesh->info.par[l]; */
459 /* if ( par->elt != MMG5_Vertex ) continue; */
460
461 /* if (p[i1]->ref != par->ref && p[i2]->ref != par->ref ) continue; */
462
463 /* if ( !locPar ) { */
464 /* hausd = par->hausd; */
465 /* hmin = par->hmin; */
466 /* hmax = par->hmax; */
467 /* locPar = 1; */
468 /* } */
469 /* else { */
470 /* hausd = MG_MIN(hausd,par->hausd); */
471 /* hmin = MG_MAX(hmin,par->hmin); */
472 /* hmax = MG_MIN(hmax,par->hmax); */
473 /* } */
474
475 /* info = par->ref; */
476 /* } while ( ++l<mesh->info.npar ); */
477
478 /* for ( ; l<mesh->info.npar; ++l) { */
479 /* par = &mesh->info.par[l]; */
480 /* if ( par->elt != MMG5_Vertex || par->ref == info ) continue; */
481 /* if (p[i1]->ref != par->ref && p[i2]->ref != par->ref ) continue; */
482
483 /* hausd = MG_MIN(hausd,par->hausd); */
484 /* hmax = MG_MIN(hmax,par->hmax); */
485 /* break; */
486 /* } */
487 /* } */
488 /* } */
489
490 hma2 = MMG3D_LLONG*MMG3D_LLONG*hmax*hmax;
491
492 /* check length */
493 ux = p[i2]->c[0] - p[i1]->c[0];
494 uy = p[i2]->c[1] - p[i1]->c[1];
495 uz = p[i2]->c[2] - p[i1]->c[2];
496 ll = ux*ux + uy*uy + uz*uz;
497 if ( ll < MMG5_EPSD ) continue;
498 else if ( ll > hma2 ) {
499 MG_SET(pt->flag,i);
500 continue;
501 }
502 il = 1.0 / sqrt(ll);
503
504 /* Hausdorff w/r tangent direction */
505 if ( MG_EDG_OR_NOM(pt->tag[i]) ) {
506 if ( MG_SIN(p[i1]->tag) ) {
507 t1[0] = il * ux;
508 t1[1] = il * uy;
509 t1[2] = il * uz;
510 }
511 else {
512 if ( !MG_EDG_OR_NOM(p[i1]->tag) ) {
513 if ( !mmgWarn0 ) {
514 fprintf(stderr,"\n ## Warning: %s: a- at least 1 geometrical"
515 " problem: non consistency between point tag (%d) and"
516 " edge tag (%d).\n",__func__,p[i1]->tag,pt->tag[i]);
517 mmgWarn0 = 1;
518 }
519 return -1;
520 }
521 memcpy(t1,t[i1],3*sizeof(double));
522 ps = t1[0]*ux + t1[1]*uy + t1[2]*uz;
523 if ( ps < 0.0 ) {
524 t1[0] *= -1.0;
525 t1[1] *= -1.0;
526 t1[2] *= -1.0;
527 }
528 }
529 if ( MG_SIN(p[i2]->tag) ) {
530 t2[0] = -il * ux;
531 t2[1] = -il * uy;
532 t2[2] = -il * uz;
533 }
534 else {
535 if ( !MG_EDG_OR_NOM(p[i2]->tag) ) {
536 if ( !mmgWarn1 ) {
537 fprintf(stderr,"\n ## Warning: %s: b- at least 1 geometrical"
538 " problem: non consistency between point tag (%d) and"
539 " edge tag (%d).\n",__func__,p[i2]->tag,pt->tag[i]);
540 mmgWarn1 = 1;
541 }
542 return -1;
543 }
544 memcpy(t2,t[i2],3*sizeof(double));
545 ps = - ( t2[0]*ux + t2[1]*uy + t2[2]*uz );
546 if ( ps < 0.0 ) {
547 t2[0] *= -1.0;
548 t2[1] *= -1.0;
549 t2[2] *= -1.0;
550 }
551 }
552 }
553 else {
554 n1 = n[i1];
555 n2 = n[i2];
556 if ( !MMG5_BezierTgt(p[i1]->c,p[i2]->c,n1,n2,t1,t2) ) {
557 t1[0] = ux * il;
558 t1[1] = uy * il;
559 t1[2] = uz * il;
560
561 t2[0] = -ux * il;
562 t2[1] = -uy * il;
563 t2[2] = -uz * il;
564 }
565 }
566 alpha = MMG5_BezierGeod(p[i1]->c,p[i2]->c,t1,t2);
567 ps = t1[0]*ux + t1[1]*uy + t1[2]*uz;
568 ps *= il;
569 dis = alpha*alpha*fabs(1.0 - ps*ps);
570 if ( dis > hausd*hausd ) {
571 MG_SET(pt->flag,i);
572 continue;
573 }
574 ps = -( t2[0]*ux + t2[1]*uy + t2[2]*uz );
575 ps *= il;
576 dis = alpha*alpha*fabs(1.0 - ps*ps);
577
578 if ( dis > hausd*hausd ) {
579 MG_SET(pt->flag,i);
580 continue;
581 }
582 }
583 return pt->flag;
584}
585
598MMG5_int MMG5_swpmsh(MMG5_pMesh mesh,MMG5_pSol met,MMG3D_pPROctree PROctree, int typchk) {
599 MMG5_pTetra pt;
600 MMG5_pxTetra pxt;
601 int it,ilist,ret,maxit;
602 int8_t i,j,ia,ier;
603 MMG5_int k,it1,it2,ns,nns;
604 int64_t list[MMG3D_LMAX+2];
605
606 it = nns = 0;
607 maxit = 2;
608 do {
609 ns = 0;
610 for (k=1; k<=mesh->ne; k++) {
611 pt = &mesh->tetra[k];
612 if ( (!MG_EOK(pt)) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
613 else if ( !pt->xt ) continue;
614 pxt = &mesh->xtetra[pt->xt];
615
616 for (i=0; i<4; i++) {
617 ier = 0;
618 if ( !(pxt->ftag[i] & MG_BDY) ) continue;
619 for (j=0; j<3; j++) {
620 ia = MMG5_iarf[i][j];
621 /* Mark the edge as boundary in case of missing tag */
622 pxt->tag[ia] |= MG_BDY;
623
624 /* No swap of geometric edge */
625 if ( MG_EDG_OR_NOM(pxt->tag[ia]) || (pxt->tag[ia] & MG_REQ) )
626 continue;
627
628 ret = MMG5_coquilface(mesh,k,i,ia,list,&it1,&it2,0);
629 ilist = ret / 2;
630 if ( ret < 0 ) return -1;
631 /* CAUTION: trigger collapse with 2 elements */
632 if ( ilist <= 1 ) continue;
633
634 /* Here, we work on a boundary edge lying along a boundary face */
635 ier = MMG5_chkswpbdy(mesh,met,list,ilist,it1,it2,typchk);
636 if ( ier < 0 )
637 return -1;
638 else if ( ier ) {
639 ier = MMG5_swpbdy(mesh,met,list,ret,it1,PROctree,typchk);
640 if ( ier > 0 ) ns++;
641 else if ( ier < 0 ) return -1;
642 break;
643 }
644 }
645 if ( ier ) break;
646 }
647 }
648 nns += ns;
649 }
650 while ( ++it < maxit && ns > 0 );
651 if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nns > 0 )
652 fprintf(stdout," %8" MMG5_PRId " edge swapped\n",nns);
653
654 return nns;
655}
656
672MMG5_int MMG5_swptet(MMG5_pMesh mesh,MMG5_pSol met,double crit,double declic,
673 MMG3D_pPROctree PROctree,int typchk,MMG5_int testmark) {
674 MMG5_pTetra pt;
675 MMG5_pxTetra pxt;
676 int ilist,it,maxit,ier;
677 int64_t list[MMG3D_LMAX+2];
678 MMG5_int k,ns,nns,nconf;
679 int8_t i;
680
681 maxit = 2;
682 it = nns = 0;
683
684 do {
685 ns = 0;
686 for (k=1; k<=mesh->ne; k++) {
687 pt = &mesh->tetra[k];
688 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
689 else if ( pt->mark < testmark ) continue;
690
691 if ( pt->qual > declic ) continue;
692
693 for (i=0; i<6; i++) {
694 /* Prevent swap of a ref or tagged edge */
695 if ( pt->xt ) {
696 pxt = &mesh->xtetra[pt->xt];
697 if ( pxt->edg[i] || pxt->tag[i] ) continue;
698 }
699
700 nconf = MMG5_chkswpgen(mesh,met,k,i,&ilist,list,crit,typchk);
701 if ( nconf<0 ) return -1;
702 else if ( nconf ) {
703 ier = MMG5_swpgen(mesh,met,nconf,ilist,list,PROctree,typchk);
704 if ( ier > 0 ) ns++;
705 else if ( ier < 0 ) return -1;
706 break;
707 }
708 }
709 }
710 nns += ns;
711 }
712 while ( ++it < maxit && ns > 0 );
713 if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nns > 0 )
714 fprintf(stdout," %8" MMG5_PRId " edge swapped\n",nns);
715
716 return nns;
717}
718
737 double clickSurf,double clickVol,int moveVol, int improveSurf,
738 int improveVolSurf, int improveVol, int maxit,MMG5_int testmark) {
739 MMG5_pTetra pt;
740 MMG5_pPoint ppt;
741 MMG5_pxTetra pxt;
742 MMG5_Tria tt;
743 double *n,caltri;
744 int ier,ilists,ilistv,it,i;
745 MMG5_int k,lists[MMG3D_LMAX+2],nm,nnm,ns,base;
746 int64_t listv[MMG3D_LMAX+2];
747 uint8_t j,i0;
748
749 if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug )
750 fprintf(stdout," ** OPTIMIZING MESH\n");
751
752 base = 1;
753 for (k=1; k<=mesh->np; k++)
754 mesh->point[k].flag = base;
755
756 it = nnm = 0;
757 do {
758 base++;
759 nm = ns = 0;
760 for (k=1; k<=mesh->ne; k++) {
761 pt = &mesh->tetra[k];
762 if ( !MG_EOK(pt) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
763 else if ( pt->mark < testmark ) continue;
764
765 /* point j on face i */
766 for (i=0; i<4; i++) {
767 for (j=0; j<3; j++) {
768 if ( pt->xt ) {
769 pxt = &mesh->xtetra[pt->xt];
770 }
771 else pxt = 0;
772
773 i0 = MMG5_idir[i][j];
774 ppt = &mesh->point[pt->v[i0]];
775 if ( ppt->flag == base ) continue;
776 else if ( MG_SIN(ppt->tag) ) continue;
777
778 if( pt->xt && (pxt->ftag[i] & MG_BDY)) {
779 /* skip required faces */
780 if( pxt->ftag[i] & MG_REQ ) continue;
781
782 assert( 0<=i && i<4 && "unexpected local face idx");
783 MMG5_tet2tri(mesh,k,i,&tt);
784 caltri = MMG5_caltri(mesh,met,&tt);
785
786 if ( caltri >= clickSurf) {
787 j = 3;
788 continue;
789 }
790 }
791 if ( maxit != 1 ) {
792 ppt->flag = base;
793 }
794 ier = 0;
795 if ( ppt->tag & MG_BDY ) {
796 /* Catch a boundary point by an external face, unless point is internal non manifold */
797 if ( (!pt->xt) || !(MG_BDY & pxt->ftag[i]) ) continue;
798 else if ( (ppt->tag & MG_PARBDY) || (ppt->tag & MG_PARBDYBDY) ) continue; /* skip parallel points seen by non-required faces */
799 else if ( ppt->tag & MG_NOM ){
800 if ( ppt->xp && mesh->xpoint[ppt->xp].nnor ) {
801 assert( 0<=i0 && i0<4 && "unexpected local index for vertex");
802 ilistv = MMG5_boulevolp(mesh,k,i0,listv);
803 if ( !ilistv ) continue;
804 /* Iso for now */
805 ier = MMG5_movbdynomintpt_iso(mesh,met,PROctree,listv,ilistv,improveVolSurf);
806 }
807 else {
808 if ( mesh->adja[4*(k-1)+1+i] ) continue;
809
810 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,1);
811 if( !ier ) continue;
812 else if ( ier>0 )
813 ier = MMG5_movbdynompt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
814 else
815 return -1;
816 }
817 }
818 else if ( ppt->tag & MG_GEO ) {
819 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
820 if ( !ier ) continue;
821 else if ( ier>0 )
822 ier = MMG5_movbdyridpt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
823 else
824 return -1;
825 }
826 else if ( ppt->tag & MG_REF ) {
827 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
828 if ( !ier )
829 continue;
830 else if ( ier>0 )
831 ier = MMG5_movbdyrefpt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
832 else
833 return -1;
834 }
835 else {
836 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
837 if ( !ier )
838 continue;
839 else if ( ier<0 )
840 return -1;
841
842 n = &(mesh->xpoint[ppt->xp].n1[0]);
843
844 /* If the orientation of the tetra face is
845 * compatible with the triangle (MG_GET(ori,i)), we know that we
846 * are well orientated. */
847 if ( !MG_GET(pxt->ori,i) ) {
848 if ( !MMG5_directsurfball(mesh,pt->v[i0],lists,ilists,n) )
849 continue;
850 }
851 ier = MMG5_movbdyregpt(mesh,met,PROctree,listv,ilistv,
852 lists,ilists,improveSurf,improveVolSurf);
853 if (ier < 0 ) return -1;
854 else if ( ier ) ns++;
855 }
856 }
857 else if ( moveVol && (pt->qual < clickVol) ) {
858 assert( 0<=i0 && i0<4 && "unexpected local index for vertex");
859 ilistv = MMG5_boulevolp(mesh,k,i0,listv);
860 if ( !ilistv ) continue;
861 ier = MMG5_movintpt(mesh,met,PROctree,listv,ilistv,improveVol);
862 }
863 if ( ier ) {
864 nm++;
865 if(maxit==1){
866 ppt->flag = base;
867 }
868 }
869 }
870 }
871 }
872 nnm += nm;
873 if ( mesh->info.ddebug ) fprintf(stdout," %8" MMG5_PRId " moved, %" MMG5_PRId " geometry\n",nm,ns);
874 }
875 while( ++it < maxit && nm > 0 );
876
877 if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nnm )
878 fprintf(stdout," %8" MMG5_PRId " vertices moved, %d iter.\n",nnm,it);
879
880 return nnm;
881}
882
892static MMG5_int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
893 MMG5_pTetra pt,ptloc;
894 MMG5_pxTetra pxt;
895 MMG5_pPoint p0,p1;
896 MMG5_pPar par;
897 double ll,ux,uy,uz,hmi2;
898 int ilists,ilist;
899 MMG5_int base,k,nc,nnm,lists[MMG3D_LMAX+2],refmin,refplus;
900 int64_t list[MMG3D_LMAX+2];
901 int l,kk,isloc,ifac1;
902
903 uint16_t tag,tag0,tag1,isnm;
904 int16_t isnmint;
905
906 int8_t i,j,ip,iq;
907 int ier, bsret; // function return values/error codes
908
909 nc = nnm = 0;
910
911 /* init point flags */
912 for (k=1; k<=mesh->np; k++) {
913 p0 = &mesh->point[k];
914 p0->flag = 0;
915 }
916
917 for (k=1; k<=mesh->ne; k++) {
918 /* Remark: we can have int32 overflow on large meshes for base field..*/
919 base = ++mesh->base;
920 pt = &mesh->tetra[k];
921 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
922
923 pxt = pt->xt? &mesh->xtetra[pt->xt] : 0;
924
925 for (i=0; i<4; i++) {
926 ier = 0;
927
928 for (j=0; j<3; j++) {
929 if ( pt->xt && (pxt->tag[MMG5_iarf[i][j]] & MG_REQ) ) continue;
930 ip = MMG5_idir[i][MMG5_inxt2[j]];
931 iq = MMG5_idir[i][MMG5_iprv2[j]];
932
933 p0 = &mesh->point[pt->v[ip]];
934 p1 = &mesh->point[pt->v[iq]];
935
936 if ( p0->flag == base ) {
937 /* I think that we can't pass here because we break the loop when base
938 * is setted and just after we increment it */
939 assert(0);
940 continue;
941 }
942 else {
943 /* Ignore OLDPARBDY tag of p0 */
944 uint16_t tag0 = p0->tag, tag1 = p1->tag;
945 tag0 &= ~MG_OLDPARBDY;
946 tag1 &= ~MG_OLDPARBDY;
947 if ( (tag0 > tag1) || (tag0 & MG_REQ) ) {
948
949 /* Unable to merge edge */
950 continue;
951 }
952 }
953
954 /* Ball of point: computed here if needed for the local parameter
955 * evaluation, after length check otherwise (because the ball
956 * computation is time consuming) */
957 ilist = ilists = 0;
958 refmin = refplus = -1;
959 if ( mesh->info.npar ) {
960 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
961 tag = pxt->tag[MMG5_iarf[i][j]];
962 isnm = (tag & MG_NOM);
963 isnmint = ( p0->xp && mesh->xpoint[p0->xp].nnor );
964 /* Ignore MG_OLDPARBDY for tags: commit c51f5f4f86292bdb0e33adb07c21ed327b1e2ec0 (pull request #241)
965 updates edge tags in splits: MG_OLDPARBDY may appear in edge tags and not in corresponding points
966 which could cause a wrong treatment of the next if statement */
967 tag0 = p0->tag & ~MG_OLDPARBDY;
968 tag &= ~MG_OLDPARBDY;
969
970 if ( tag0 > tag ) continue;
971
972 /* Catch an exterior non manifold point by an external face */
973 if ( isnm ) {
974 if ( isnmint ) {
975 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
976 ilist = MMG5_boulevolp(mesh,k,ip,list);
977 }
978 else {
979 if ( mesh->adja[4*(k-1)+1+i] ) continue;
980
981 bsret = MMG5_boulesurfvolpNom(mesh,k,ip,i,
982 list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM);
983 if(bsret==-1 || bsret==-3 || bsret==-4){
984 return -3; // fatal
985 }else if(bsret==-2){
986 continue; // ball computation failed: cannot handle this vertex
987 }
988 assert(bsret==1 && "unexpected return value from MMG5_boulesurfvolpNom");
989 }
990 }
991 else {
992 if (MMG5_boulesurfvolp(mesh,k,ip,i,
993 list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
994 return -2;
995 }
996 }
997 else {
998 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
999 ilist = MMG5_boulevolp(mesh,k,ip,list);
1000 }
1001 }
1002
1003 /* Check length */
1004 if ( typchk == 1 ) {
1005 ux = p1->c[0] - p0->c[0];
1006 uy = p1->c[1] - p0->c[1];
1007 uz = p1->c[2] - p0->c[2];
1008 ll = ux*ux + uy*uy + uz*uz;
1009
1010 /* local parameters*/
1011 hmi2 = mesh->info.hmin;
1012 isloc = 0;
1013
1014 /* Local parameters at vertex */
1015 // take the min of the local hmin at p0 and hmin at p1
1016
1017 /* Local parameters at tetra */
1018 if ( mesh->info.parTyp & MG_Tetra ) {
1019 l = 0;
1020 do
1021 {
1022 if ( isloc ) break;
1023
1024 par = &mesh->info.par[l];
1025 if ( par->elt != MMG5_Tetrahedron ) continue;
1026
1027 for ( kk=0; kk<ilist; ++kk ) {
1028 ptloc = &mesh->tetra[list[kk]/4];
1029 if ( par->ref == ptloc->ref ) {
1030 hmi2 = par->hmin;
1031 isloc = 1;
1032 break;
1033 }
1034 }
1035 } while ( ++l<mesh->info.npar );
1036
1037 for ( ; l<mesh->info.npar; ++l ) {
1038 par = &mesh->info.par[l];
1039 if ( par->elt != MMG5_Tetrahedron ) continue;
1040
1041 for ( kk=0; kk<ilist; ++kk ) {
1042 ptloc = &mesh->tetra[list[kk]/4];
1043 if ( par->ref == ptloc->ref ) {
1044 hmi2 = MG_MAX(hmi2,par->hmin);
1045 break;
1046 }
1047 }
1048 }
1049 }
1050
1051 /* Local parameters at triangle */
1052 if ( mesh->info.parTyp & MG_Tria && ( pt->xt && (pxt->ftag[i] & MG_BDY) )) {
1053 l = 0;
1054 do
1055 {
1056 if ( isloc ) break;
1057
1058 par = &mesh->info.par[l];
1059 if ( par->elt != MMG5_Triangle ) continue;
1060
1061 for ( kk=0; kk<ilists; ++kk ) {
1062 ptloc = &mesh->tetra[lists[kk]/4];
1063 ifac1 = lists[kk] % 4;
1064 assert(ptloc->xt && (mesh->xtetra[ptloc->xt].ftag[ifac1] & MG_BDY) );
1065
1066 if ( par->ref == mesh->xtetra[ptloc->xt].ref[ifac1] ) {
1067 hmi2 = par->hmin;
1068 isloc = 1;
1069 break;
1070 }
1071 }
1072 } while ( ++l<mesh->info.npar );
1073
1074 for ( ; l<mesh->info.npar; ++l ) {
1075 par = &mesh->info.par[l];
1076 if ( par->elt != MMG5_Triangle ) continue;
1077
1078 for ( kk=0; kk<ilists; ++kk ) {
1079 ptloc = &mesh->tetra[lists[kk]/4];
1080 ifac1 = lists[kk] % 4;
1081 assert(ptloc->xt && (mesh->xtetra[ptloc->xt].ftag[ifac1] & MG_BDY) );
1082
1083 if ( par->ref == mesh->xtetra[ptloc->xt].ref[ifac1] ) {
1084 hmi2 = MG_MAX(hmi2,par->hmin);
1085 break;
1086 }
1087 }
1088 }
1089 }
1090
1091 hmi2 = hmi2*hmi2;
1092 if ( ll > hmi2*MMG3D_LSHRT ) continue;
1093 }
1094 else if ( typchk == 2 ) {
1095 ll = MMG5_lenedg(mesh,met,MMG5_iarf[i][j],pt);
1096 // Case of an internal tetra with 4 ridges vertices.
1097 if ( ll == 0 ) continue;
1098 if ( ll > MMG3D_LSHRT ) continue;
1099 }
1100
1101
1102 /* Ball computation if not computed before */
1103 if ( !mesh->info.npar ) {
1104 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
1105 tag = pxt->tag[MMG5_iarf[i][j]];
1106 isnm = (tag & MG_NOM);
1107 isnmint = ( p0->xp && mesh->xpoint[p0->xp].nnor );
1108
1109 tag0 = p0->tag & ~MG_OLDPARBDY;
1110 tag &= ~MG_OLDPARBDY;
1111
1112 if ( tag0 > tag ) continue;
1113 if ( isnm ) {
1114 if ( isnmint ) {
1115 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
1116 ilist = MMG5_boulevolp(mesh,k,ip,list);
1117 }
1118 else {
1119 if ( mesh->adja[4*(k-1)+1+i] ) continue;
1120 bsret = MMG5_boulesurfvolpNom(mesh,k,ip,i,
1121 list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM);
1122 if(bsret==-1 || bsret==-3 || bsret==-4){
1123 return -3; // fatal
1124 }else if(bsret==-2){
1125 continue; // ball computation failed: cannot handle this vertex
1126 }
1127 assert(bsret==1 && "unexpected return value from MMG5_boulesurfvolpNom");
1128 }
1129 }
1130 else {
1131 if (MMG5_boulesurfvolp(mesh,k,ip,i,
1132 list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
1133 return -4;
1134 }
1135 }
1136 else {
1137 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
1138 ilist = MMG5_boulevolp(mesh,k,ip,list);
1139 }
1140 }
1141
1142 /* boundary face: collapse ip on iq */
1143 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
1144 tag = pxt->tag[MMG5_iarf[i][j]];
1145 tag |= MG_BDY;
1146 tag &= ~MG_OLDPARBDY;
1147 tag0 &= ~MG_OLDPARBDY;
1148 if ( tag0 > tag ) continue;
1149
1150 isnm = ( tag & MG_NOM );
1151 isnmint = ( tag & MG_NOM ) ? ( p0->xp && mesh->xpoint[p0->xp].nnor ) : 0;
1152 if ( isnm && (!isnmint) ) {
1153 /* Treat surfacic non manifold points from a surface triangle */
1154 if ( mesh->adja[4*(k-1)+1+i] ) continue;
1155 }
1156 ilist = MMG5_chkcol_bdy(mesh,met,k,i,j,list,ilist,lists,ilists,refmin,refplus,typchk,isnm,isnmint);
1157 }
1158 /* internal face */
1159 else {
1160 isnm = 0;
1161 if ( p0->tag & MG_BDY ) continue;
1162 ilist = MMG5_chkcol_int(mesh,met,k,i,j,list,ilist,typchk);
1163 }
1164
1165 if ( ilist > 0 ) {
1166 ier = MMG5_colver(mesh,met,list,ilist,iq,typchk);
1167 if ( ier < 0 ) return -5;
1168 else if ( ier ) {
1170 break;
1171 }
1172 }
1173 else if (ilist < 0 ) return -6;
1174 }
1175 if ( ier ) {
1176 p1->flag = base;
1177 if ( isnm ) nnm++;
1178 nc++;
1179 break;
1180 }
1181 }
1182 }
1183 if ( nc > 0 && (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) )
1184 fprintf(stdout," %8" MMG5_PRId " vertices removed, %8" MMG5_PRId " non manifold,\n",nc,nnm);
1185
1186 return nc;
1187}
1188
1211 MMG3D_pPROctree *PROctree,MMG5_int k,
1212 int8_t imin,double lmin,MMG5_int* nc) {
1213 MMG5_pTetra pt;
1214 MMG5_pxTetra pxt;
1215 MMG5_pPoint p0,p1;
1216 int64_t list[MMG3D_LMAX+2];
1217 MMG5_int lists[MMG3D_LMAX+2],ip1,ip2;
1218 int ilist,ilists;
1219 int8_t j,i,i1,i2;
1220
1221 if(lmin > MMG3D_LOPTS) {
1222 /* Edge is large enough: nothing to do */
1223 return 3;
1224 }
1225
1226 if ( lmin == 0 ) {
1227 /* Case of an internal tetra with 4 ridges vertices */
1228//#warning is it possible to merge this edge ??
1229 return 0;
1230 }
1231
1232 pt = &mesh->tetra[k];
1233 pxt = pt->xt ? &mesh->xtetra[pt->xt] : 0;
1234
1235 MMG3D_find_bdyface_from_edge(mesh,pt,imin,&i,&j,&i1,&i2,&ip1,&ip2,&p0,&p1);
1236
1237 /* Ignore OLDPARBDY tag of p0 */
1238 uint16_t tag0 = p0->tag;
1239 tag0 &= ~MG_OLDPARBDY;
1240 if ( (tag0 > p1->tag) || (tag0 & MG_REQ) ) {
1241 /* Unable to merge edge: pass to next element */
1242 return 0;
1243 }
1244
1246 ilist = 0;
1247 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
1248 /* Case of a boundary face */
1249 uint16_t tag = pxt->tag[MMG5_iarf[i][j]];
1250 if ( tag & MG_REQ ) {
1251 return 0;
1252 }
1253 tag |= MG_BDY;
1254 tag &= ~MG_OLDPARBDY;
1255 if ( tag0 > tag ) {
1256 return 0;
1257 }
1258 if ( ( tag & MG_NOM ) && (mesh->adja[4*(k-1)+1+i]) ) {
1259 return 0;
1260 }
1261
1262 uint16_t isnm = (p0->tag & MG_NOM);
1263 if (MMG5_boulesurfvolp(mesh,k,i1,i, list,&ilist,lists,&ilists,isnm) < 0 ) {
1264 return -1;
1265 }
1266
1267 ilist = MMG5_chkcol_bdy(mesh,met,k,i,j,list,ilist,lists,ilists,0,0,2,0,0);
1268 }
1269 else {
1270 /* Case of an internal face */
1271 if ( p0->tag & MG_BDY ) {
1272 return 0;
1273 }
1274
1275 ilist = MMG5_boulevolp(mesh,k,i1,list);
1276 ilist = MMG5_chkcol_int(mesh,met,k,i,j,list,ilist,2);
1277 }
1278
1280 if ( ilist > 0 ) {
1281 /* Checks are ok */
1282 int ier = MMG5_colver(mesh,met,list,ilist,i2,2);
1283 if ( ilist < 0 ) {
1284 /* Colver failure */
1285 return 0;
1286 }
1287 if ( ier < 0 ) {
1288 /* Colver failure */
1289 return -1;
1290 }
1291 else if(ier) {
1292 /* Collapse is successful */
1293 if ( PROctree && (*PROctree) ) {
1294 MMG3D_delPROctree(mesh,*PROctree,ier);
1295 }
1297 (*nc)++;
1298 return 2;
1299 }
1300 }
1301 else if (ilist < 0 ) {
1302 /* Checks on chkcol_bdy have failed (error in edge shell computation) */
1303 return -1;
1304 }
1305
1306 return 3;
1307}
1308
1309
1318static inline
1320{
1321 MMG5_pTetra pt;
1322 int ia,i,j;
1323 MMG5_int vx[6],k;
1324
1325 /* Delete the useless added points */
1326 for (k=1; k<=mesh->ne; k++) {
1327 pt = &mesh->tetra[k];
1328 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
1329
1330 memset(vx,0,6*sizeof(MMG5_int));
1331 for (ia=0,i=0; i<3; i++) {
1332 for (j=i+1; j<4; j++,ia++) {
1333 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_REQ) ) continue;
1334 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
1335 if ( vx[ia] > 0 ) {
1336 MMG3D_delPt(mesh,vx[ia]);
1337 if ( !MMG5_hashUpdate(&hash,pt->v[i],pt->v[j],-1) ) {
1338 fprintf(stderr,"\n ## Error: %s: unable to delete point idx"
1339 " along edge %" MMG5_PRId " %" MMG5_PRId ".\n", __func__,
1340 MMG3D_indPt(mesh,pt->v[i]),
1341 MMG3D_indPt(mesh,pt->v[j]));
1342 MMG5_DEL_MEM(mesh,hash.item);
1343 return 0;
1344 }
1345 }
1346 }
1347 }
1348 }
1349 return 1;
1350}
1351
1362static MMG5_int
1364 MMG5_pTetra pt;
1365 MMG5_pPoint p1,p2;
1366 MMG5_xTetra *pxt;
1367 MMG5_Hash hash;
1368 MMG5_pPar par;
1369 double ll,o[3],ux,uy,uz,hma2,mincal;
1370 int l,memlack,ier;
1371 MMG5_int src,vx[6],ip,ip1,ip2,k,ne,ns,nap;
1372 int8_t i,j,ia;
1373
1375 if ( !MMG5_hashNew(mesh,&hash,mesh->np,7*mesh->np) ) return -1;
1376 memlack = ns = nap = 0;
1377
1378 /* Hash all boundary and required edges, and put ip = -1 in hash structure */
1379 for (k=1; k<=mesh->ne; k++) {
1380 pt = &mesh->tetra[k];
1381 if ( !MG_EOK(pt) ) continue;
1382
1383 /* avoid split of edges belonging to a required tet */
1384 if ( pt->tag & MG_REQ ) {
1385 for (i=0; i<6; i++) {
1386 ip1 = pt->v[MMG5_iare[i][0]];
1387 ip2 = pt->v[MMG5_iare[i][1]];
1388 ip = -1;
1389 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
1390 }
1391 continue;
1392 }
1393
1394 if ( !pt->xt ) continue;
1395
1396 pxt = &mesh->xtetra[pt->xt];
1397 for (i=0; i<4; i++) {
1398 if ( pxt->ftag[i] & MG_BDY ) {
1399 for (j=0; j<3; j++) {
1400 ip1 = pt->v[MMG5_idir[i][MMG5_inxt2[j]]];
1401 ip2 = pt->v[MMG5_idir[i][MMG5_iprv2[j]]];
1402 ip = -1;
1403 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
1404 }
1405 }
1406 }
1407 }
1408
1410 mincal = FLT_MAX;
1411 for (k=1; k<=mesh->ne; k++) {
1412 pt = &mesh->tetra[k];
1413 if ( !MG_EOK(pt) ) continue;
1414 pt->flag = 0;
1415 for (i=0; i<6; i++) {
1416 ip = -1;
1417 ip1 = pt->v[MMG5_iare[i][0]];
1418 ip2 = pt->v[MMG5_iare[i][1]];
1419 p1 = &mesh->point[ip1];
1420 p2 = &mesh->point[ip2];
1421 if ( pt->xt ) {
1422 pxt = &mesh->xtetra[pt->xt];
1423 if ( pxt->tag[i] & MG_REQ ) continue;
1424 }
1425 else pxt = 0;
1426
1427 if ( (p1->tag & MG_BDY) && (p2->tag & MG_BDY) ) {
1428 /* Split internal edges connecting bdy points */
1429 ip = MMG5_hashGet(&hash,ip1,ip2);
1430 }
1431 else {
1432 if (typchk == 1) {
1433 ux = p2->c[0] - p1->c[0];
1434 uy = p2->c[1] - p1->c[1];
1435 uz = p2->c[2] - p1->c[2];
1436 ll = ux*ux + uy*uy + uz*uz;
1437
1438 hma2 = mesh->info.hmax;
1439 /* Local parameters */
1440 /* Local param at vertices ip1 and ip2 */
1441 // take the min of the local hmin at ip1 and hmin at ip2
1442
1443 /* Local parameters at tetra */
1444 if ( mesh->info.parTyp & MG_Tetra ) {
1445 for ( l=0; l<mesh->info.npar; ++l ) {
1446 par = &mesh->info.par[l];
1447
1448 if ( par->elt != MMG5_Tetrahedron ) continue;
1449 if ( par->ref != pt->ref ) continue;
1450
1451 hma2 = par->hmax;
1452 break;
1453 }
1454 }
1455 hma2 = MMG3D_LLONG*MMG3D_LLONG*hma2*hma2;
1456 if ( ll > hma2 ) {
1457 ip = MMG5_hashGet(&hash,ip1,ip2);
1458 mincal = MG_MIN(mincal,pt->qual);
1459 }
1460 }
1461 else if ( typchk == 2 ) {
1462 ll = MMG5_lenedg(mesh,met,i,pt);
1463 // Case of an internal tetra with 4 ridges vertices.
1464 if ( ll == 0 ) continue;
1465 if ( ll > MMG3D_LLONG ) {
1466 ip = MMG5_hashGet(&hash,ip1,ip2);
1467 mincal = MG_MIN(mincal,pt->qual);
1468 }
1469 }
1470 }
1471 if ( ip < 0 ) continue;
1472 else if ( !ip ) {
1473 /* new midpoint */
1474 o[0] = 0.5 * (p1->c[0]+p2->c[0]);
1475 o[1] = 0.5 * (p1->c[1]+p2->c[1]);
1476 o[2] = 0.5 * (p1->c[2]+p2->c[2]);
1477
1478#ifdef USE_POINTMAP
1479 src = mesh->point[ip1].src;
1480#else
1481 src = 1;
1482#endif
1483 ip = MMG3D_newPt(mesh,o,0,src);
1484 if ( !ip ) {
1485 /* reallocation of point table */
1486
1488 fprintf(stderr,"\n ## Warning: %s: unable to"
1489 " allocate a new point\n",__func__);
1491 memlack=1;
1492 goto split
1493 ,o,0,src);
1494 }
1495
1496 assert ( met );
1497 if ( met->m ) {
1498 if ( typchk == 1 && (met->size>1) )
1499 ier = MMG3D_intmet33_ani(mesh,met,k,i,ip,0.5);
1500 else
1501 ier = MMG5_intmet(mesh,met,k,i,ip,0.5);
1502
1503 if (!ier) {
1504 // Unable to compute the metric
1505 return -1;
1506 }
1507 else if ( ier < 0 ) {
1508 MMG3D_delPt(mesh,ip);
1509 continue;
1510 }
1511 }
1512
1513 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
1514 MG_SET(pt->flag,i);
1515 nap++;
1516 }
1517 }
1518 }
1519 if ( !nap ) {
1520 MMG5_DEL_MEM(mesh,hash.item);
1521 return 0;
1522 }
1523
1525split:
1526 if ( mincal < MMG5_EPS ) {
1527 /* Delete the useless added points */
1528 if ( !MMG3D_delPatternPts(mesh,hash) ) return -1;
1529
1530 /* Avoid the creation of bad quality elements */
1531 if ( mesh->info.imprim > 5 || mesh->info.ddebug ) {
1532 fprintf(stderr,"\n ## Warning: %s: too bad quality for the worst element."
1533 " Volumic patterns skipped.\n",__func__);
1534 }
1535
1536 MMG5_DEL_MEM(mesh,hash.item);
1537 return 0;
1538 }
1539
1540 ns = 0;
1541 ne = mesh->ne;
1542 for (k=1; k<=ne; k++) {
1543 pt = &mesh->tetra[k];
1544 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
1545 memset(vx,0,6*sizeof(MMG5_int));
1546 pt->flag = 0;
1547 for (ia=0,i=0; i<3; i++) {
1548 for (j=i+1; j<4; j++,ia++) {
1549 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_REQ) ) continue;
1550 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
1551 if ( vx[ia] > 0 ) MG_SET(pt->flag,ia);
1552 }
1553 }
1554
1555 switch (pt->flag) {
1556 case 1: case 2: case 4: case 8: case 16: case 32: /* 1 edge split */
1557 if ( ! MMG5_split1(mesh,met,k,vx,typchk-1) ) return -1;
1558 ns++;
1559 break;
1560 case 48: case 24: case 40: case 6: case 34: case 36:
1561 case 20: case 5: case 17: case 9: case 3: case 10: /* 2 edges (same face) split */
1562 if ( ! MMG5_split2sf(mesh,met,k,vx,typchk-1) ) return -1;
1563 ns++;
1564 break;
1565
1566 case 33: case 18: case 12: /* 2 opposite edges split */
1567 if ( ! MMG5_split2(mesh,met,k,vx,typchk-1) ) return -1;
1568 ns++;
1569 break;
1570
1571 case 11: case 21: case 38: case 56: /* 3 edges on the same faces splitted */
1572 if ( ! MMG5_split3(mesh,met,k,vx,typchk-1) ) return -1;
1573 ns++;
1574 break;
1575
1576 case 7: case 25: case 42: case 52: /* 3 edges on conic configuration splitted */
1577 if ( ! MMG5_split3cone(mesh,met,k,vx,typchk-1) ) return -1;
1578 ns++;
1579 break;
1580
1581 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
1582 case 14: case 49: case 50: case 44: case 41: /* 3 edges on opposite configuration splitted */
1583 if ( ! MMG5_split3op(mesh,met,k,vx,typchk-1) ) return -1;
1584 ns++;
1585 break;
1586
1587 case 23: case 29: case 53: case 60: case 57: case 58:
1588 case 27: case 15: case 43: case 39: case 54: case 46: /* 4 edges with 3 lying on the same face splitted */
1589 if ( ! MMG5_split4sf(mesh,met,k,vx,typchk-1) ) return -1;
1590 ns++;
1591 break;
1592
1593 /* 4 edges with no 3 lying on the same face splitted */
1594 case 30: case 45: case 51:
1595 if ( ! MMG5_split4op(mesh,met,k,vx,typchk-1) ) return -1;
1596 ns++;
1597 break;
1598
1599 case 62: case 61: case 59: case 55: case 47: case 31: /* 5 edges split */
1600 if ( ! MMG5_split5(mesh,met,k,vx,typchk-1) ) return -1;
1601 ns++;
1602 break;
1603
1604 case 63: /* 6 edges split */
1605 if ( ! MMG5_split6(mesh,met,k,vx,typchk-1) ) return -1;
1606 ns++;
1607 break;
1608 }
1609 }
1610
1611 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
1612 fprintf(stdout," %7" MMG5_PRId " splitted\n",nap);
1613
1614 MMG5_DEL_MEM(mesh,hash.item);
1615 if ( memlack ) return -1;
1616 return nap;
1617}
1618
1630static inline int
1632 double dd,to[3];
1633
1634 dd = no[0]*pxp->n1[0]+no[1]*pxp->n1[1]+no[2]*pxp->n1[2];
1635 if ( dd > 1.0-MMG5_EPS ) return 0;
1636
1637 memcpy(pxp->n2,no,3*sizeof(double));
1638
1639 /* a computation of the tangent with respect to these two normals is possible */
1640 to[0] = pxp->n1[1]*pxp->n2[2] - pxp->n1[2]*pxp->n2[1];
1641 to[1] = pxp->n1[2]*pxp->n2[0] - pxp->n1[0]*pxp->n2[2];
1642 to[2] = pxp->n1[0]*pxp->n2[1] - pxp->n1[1]*pxp->n2[0];
1643 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
1644 if ( dd > MMG5_EPSD2 ) {
1645 dd = 1.0 / sqrt(dd);
1646 to[0] *= dd;
1647 to[1] *= dd;
1648 to[2] *= dd;
1649 memcpy(ppt->n,to,3*sizeof(double));
1650 }
1651 else {
1652 /* Detect opposite normals... */
1653 if ( to[0]*to[0]+to[1]*to[1]+to[2]*to[2] > MMG5_EPSD2 ) {
1654 /* Non opposite normals: fail in debug mode */
1655 assert ( dd > MMG5_EPSD2 );
1656 }
1657 else {
1658 /* Opposite normals: unable to compute the tangent, check if we have
1659 * already stored a tangent, fail otherwise */
1660 assert ( ppt->n[0]*ppt->n[0]+ppt->n[1]*ppt->n[1]+ppt->n[2]*ppt->n[2] > MMG5_EPSD2 );
1661 }
1662 }
1663
1664 return 1;
1665}
1666
1685 double no1[3],double no2[3],double to[3] ) {
1686 MMG5_Tria ptt;
1687 double dd;
1688 int ier;
1689 static int8_t warn_n = 0;
1690
1691 assert ( 0<=i && i<4 && "unexpected local idx for face" );
1692 assert ( 0<=j && j<3 && "unexpected local edg odx in face" );
1693
1694 assert( 0<=i && i<4 && "unexpected local face idx");
1695
1696#ifndef NDEBUG
1697 /* Check that face is boundary and has a suitable orientation */
1698 assert ( mesh->tetra[k].xt && "Tetra is not boundary" );
1699
1700 MMG5_pxTetra pxt = &mesh->xtetra[ mesh->tetra[k].xt ];
1701 assert ( (pxt->ftag[i] & MG_BDY) && "Face is not boundary" );
1702 assert ( MG_GET(pxt->ori,i) && "Wrong face orientation" );
1703#endif
1704
1705 MMG5_tet2tri(mesh,k,i,&ptt);
1706 MMG5_nortri(mesh,&ptt,no1);
1707
1708 /* In this case, 'to' orientation depends on the edge processing (so on
1709 * the triangle from which we come) so we can't use it to compute no2. */
1710 ier = MMG3D_normalAdjaTri(mesh,k,i,j,no2);
1711 if ( ier < 0 ) {
1712 return 0;
1713 }
1714 else if ( !ier ) {
1715 if ( !warn_n ) {
1716 warn_n = 1;
1717 fprintf(stderr," ## Warning: %s: %d: error in the computation of normal"
1718 " at triangle.\n",__func__,__LINE__);
1719 }
1720 no2[0] = to[1]*no1[2] - to[2]*no1[1];
1721 no2[1] = to[2]*no1[0] - to[0]*no1[2];
1722 no2[2] = to[0]*no1[1] - to[1]*no1[0];
1723
1724 dd = no2[0]*no2[0] + no2[1]*no2[1] + no2[2]*no2[2];
1725 if ( dd > MMG5_EPSD2 ) {
1726 dd = 1.0 / sqrt(dd);
1727 no2[0] *= dd;
1728 no2[1] *= dd;
1729 no2[2] *= dd;
1730 }
1731 }
1732 else {
1733 assert ( ier==1 );
1734
1735 /* Compute 'to' as intersection of no1 and no2 */
1736 to[0] = no1[1]*no2[2] - no1[2]*no2[1];
1737 to[1] = no1[2]*no2[0] - no1[0]*no2[2];
1738 to[2] = no1[0]*no2[1] - no1[1]*no2[0];
1739 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
1740 if ( dd > MMG5_EPSD2 ) {
1741 dd = 1.0 / sqrt(dd);
1742 to[0] *= dd;
1743 to[1] *= dd;
1744 to[2] *= dd;
1745 }
1746 }
1747 return 1;
1748}
1749
1781 int8_t imax,int8_t i, int8_t j,
1782 MMG5_pxTetra pxt,
1783 MMG5_int ip1,MMG5_int ip2,
1784 MMG5_pPoint p0, MMG5_pPoint p1,
1785 MMG5_int *ref,uint16_t *tag,
1786 double o[3],double to[3],double no1[3],
1787 double no2[3],int64_t *list,int *ilist) {
1788 MMG5_Tria ptt;
1789 double v[3];
1790
1791 if ( (p0->tag & MG_PARBDY) && (p1->tag & MG_PARBDY) ) {
1792 /* Skip edge with extremities on parallel interfaces */
1793 return 0;
1794 }
1795
1796 int8_t ori = MG_GET(pxt->ori,i);
1797 if ( !ori ) {
1798 /* Treat triangles at interface of 2 subdomains from well oriented face */
1799 return 0;
1800 }
1801
1802 *ref = pxt->edg[imax];
1803 *tag = pxt->tag[imax];
1804 if ( (*tag) & MG_REQ ) {
1805 /* No need to split required edges */
1806 return 0;
1807 }
1808
1809 (*tag) |= MG_BDY;
1810 int8_t dummy;
1811 *ilist = MMG5_coquil(mesh,k,imax,list,&dummy);
1812 if ( !(*ilist) ) {
1813 /* On of the tetra of the edge shell is required: we cannot split the edge */
1814 return 0;
1815 }
1816 else if ( (*ilist) < 0 ) {
1817 /* Shell computation has failed */
1818 return -1;
1819 }
1820
1822 if ( (*tag) & MG_NOM ){
1823 /* Edge is non-manifold */
1824 if( !MMG5_BezierNom(mesh,ip1,ip2,0.5,o,no1,to) ) {
1825 /* Unable to compute geom infos at nm edge */
1826 return 0;
1827 }
1828 else if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1829 assert( 0<=i && i<4 && "unexpected local face idx");
1830 MMG5_tet2tri(mesh,k,i,&ptt);
1831 MMG5_nortri(mesh,&ptt,no1);
1832 }
1833 }
1834 else if ( (*tag) & MG_GEO ) {
1835 /* Edge is ridge */
1836 if ( !MMG5_BezierRidge(mesh,ip1,ip2,0.5,o,no1,no2,to) ) {
1837 /* Unable to compute geom infos at ridge */
1838//#warning why a continue here?
1839 return 0;
1840 }
1841 else if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1842 if ( !MMG3D_normalAndTangent_at_sinRidge(mesh,k,i,j,no1,no2,to) ) {
1843 return -1;
1844 }
1845 }
1846 }
1847 else if ( (*tag) & MG_REF ) {
1848 /* Edge is ref */
1849 if ( !MMG5_BezierRef(mesh,ip1,ip2,0.5,o,no1,to) ) {
1850 /* Unable to compute geom infos at ref edge */
1851 return 1;
1852 }
1853 if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1854//#warning creation of sin-sin ref edge to see if it works without normal realloc
1855 assert( 0<=i && i<4 && "unexpected local face idx");
1856 MMG5_tet2tri(mesh,k,i,&ptt);
1857 MMG5_nortri(mesh,&ptt,no1);
1858 }
1859 }
1860 else {
1861 /* Longest edge is regular */
1862 if ( !MMG5_norface(mesh,k,i,v) ) {
1863 /* Unable to treat long edge: try to collapse short one */
1864 return 1;
1865 }
1866 if ( !MMG5_BezierReg(mesh,ip1,ip2,0.5,v,o,no1) ) {
1867 /* Unable to compute geom infos at regular edge */
1868 return 1;
1869 }
1870 else if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1871 assert( 0<=i && i<4 && "unexpected local face idx");
1872 MMG5_tet2tri(mesh,k,i,&ptt);
1873 MMG5_nortri(mesh,&ptt,no1);
1874 }
1875 }
1876 return 2;
1877}
1878
1899 int8_t *i,int8_t *j,int8_t*i1,int8_t*i2,
1900 MMG5_int*ip1,MMG5_int*ip2,MMG5_pPoint *p0,MMG5_pPoint *p1) {
1901
1902 int8_t ifa0 = MMG5_ifar[ied][0];
1903 int8_t ifa1 = MMG5_ifar[ied][1];
1904
1908 /* Default face */
1909 MMG5_pxTetra pxt = pt->xt? &mesh->xtetra[pt->xt] : 0;
1910
1911 (*i) = ifa0;
1912 if ( pt->xt ) {
1913 int16_t is_ifa0_bdy = (pxt->ftag[ifa0] & MG_BDY);
1914 int16_t is_ifa1_bdy = (pxt->ftag[ifa1] & MG_BDY);
1915
1916 if ( is_ifa0_bdy && is_ifa1_bdy ) {
1917 /* Two bdy faces: search if one has a suitable orientation */
1918 int8_t ifa1_ori = MG_GET(pxt->ori,(*i));
1919 (*i) = ifa1_ori ? ifa1 : ifa0;
1920 }
1921 else if ( is_ifa1_bdy ) {
1922 /* only ifa1 is boundary: no need to check for orientation (if it has a
1923 * bad ori we will quit the function later) */
1924 (*i) = ifa1;
1925 }
1926 /* For all other cases (only ifa0 is bdy or no bdy face), we use default
1927 * face (ifa0) */
1928 }
1929
1930 (*j) = MMG5_iarfinv[*i][ied];
1931 (*i1) = MMG5_idir[*i][MMG5_inxt2[*j]];
1932 (*i2) = MMG5_idir[*i][MMG5_iprv2[*j]];
1933 (*ip1) = pt->v[*i1];
1934 (*ip2) = pt->v[*i2];
1935 (*p0) = &mesh->point[*ip1];
1936 (*p1) = &mesh->point[*ip2];
1937
1938}
1939
1958 MMG5_pTetra pt,MMG5_pxTetra pxt,int8_t imax,int8_t typchk,
1959 int8_t chkRidTet,int *warn ) {
1960 MMG5_pPoint p0,p1,ppt;
1961 double o[3],to[3],no1[3],no2[3];
1962 int ilist;
1963 MMG5_int src,ip,ip1,ip2,ref;
1964 int64_t list[MMG3D_LMAX+2];
1965 int ier;
1966 uint16_t tag;
1967 int8_t j,i,i1,i2;
1968
1969 assert ( pxt == &mesh->xtetra[pt->xt] && "suitable xtetra assignation" );
1970
1971 assert ( ((pxt->ftag[MMG5_ifar[imax][0]] & MG_BDY) || (pxt->ftag[MMG5_ifar[imax][1]] & MG_BDY) )
1972 && "Boundary edge has to be splitted from a boundary face" );
1973
1974 /* Mark the edge as MG_BDY to avoid wrong evaluation as an internal edge by
1975 * the interpolation function (see intmet_ani) */
1976 pxt->tag[imax] |= MG_BDY;
1977
1978 /* proceed edges according to lengths */
1979 MMG3D_find_bdyface_from_edge(mesh,pt,imax,&i,&j,&i1,&i2,&ip1,&ip2,&p0,&p1);
1980
1981 ier = MMG3D_build_bezierEdge(mesh,k,imax,i,j,pxt,ip1,ip2,p0,p1,&ref,&tag,
1982 o,to,no1,no2,list,&ilist);
1983 switch ( ier ) {
1984 case -1:
1985 /* Strong failure (due to lack of memory or wrong edge shell) */
1986 return -1;
1987 case 0:
1988 /* We don't want to split the edge (edge shell contains a required tet or
1989 * non manifold edge or ridge reconstruction has failed) */
1990 return 0;
1991 case 1:
1992 /* We don't want to split the edge (ref or regular edge reconstruction has
1993 * failed) */
1994 return 0;
1995 default:
1996 /* We can insert new vertex along edge */
1997 assert ( ier==2 );
1998 }
1999
2000#ifdef USE_POINTMAP
2001 src = mesh->point[ip1].src;
2002#else
2003 src = 1;
2004#endif
2005 ip = MMG3D_newPt(mesh,o,tag,src);
2006 if ( !ip ) {
2007 /* reallocation of point table */
2009 *warn=1;
2010 return 2;
2011 ,o,tag,src);
2012 }
2013 if ( met && met->m ) {
2014 if ( typchk == 1 && (met->size>1) ) {
2015 ier = MMG3D_intmet33_ani(mesh,met,k,imax,ip,0.5);
2016 }
2017 else {
2018 ier = MMG5_intmet(mesh,met,k,imax,ip,0.5);
2019 }
2020 if ( !ier ) {
2021 MMG3D_delPt(mesh,ip);
2022 return -1;
2023 }
2024 else if (ier < 0) {
2025 MMG3D_delPt(mesh,ip);
2026 return 0;
2027 }
2028 }
2029
2030 /* simbulgept needs a valid tangent at ridge point (to build ridge metric in
2031 * order to comute edge lengths). Thus we need to store the geometric info of
2032 * point here. */
2033 ppt = &mesh->point[ip];
2034 MMG3D_set_geom(mesh,ppt,tag,ref,pxt->ref[i],no1,no2,to);
2035
2036 ier = MMG3D_simbulgept(mesh,met,list,ilist,ip);
2037
2038#ifndef NDEBUG
2039 if ( mesh->info.ddebug ) {
2040 assert ( (ier != -1) && "simbulgept failure" );
2041 }
2042#endif
2043
2044 if ( ier < 0 || ier == 2 ) {
2045 MMG3D_delPt(mesh,ip);
2046 return 0;
2047 }
2048 else if ( ier == 0 ) {
2049 ier = MMG3D_dichoto1b(mesh,met,list,ilist,ip);
2050 }
2051 if ( ier == 1 ) { ier = MMG5_split1b(mesh,met,list,ilist,ip,1,typchk-1,chkRidTet); }
2052
2053 /* if we realloc memory in MMG5_split1b pt and pxt pointers are not valid */
2054 if ( ier < 0 ) {
2055 fprintf(stderr,"\n ## Error: %s: unable to split.\n",__func__);
2056 return -1;
2057 }
2058 else if ( ier == 0 || ier == 2 ) {
2059 MMG3D_delPt(mesh,ip);
2060 return 0;
2061 }
2062
2063 return 1;
2064}
2065
2081static
2083 MMG5_pxTetra pxt,int8_t i,MMG5_pTria ptt,int8_t typchk ) {
2084
2085 MMG5_pPar par;
2086 double len,hmax,hausd;
2087 int l;
2088 MMG5_int ip1,ip2;
2089 int8_t isloc,ier;
2090 int8_t j,i1,i2,ia;
2091
2092 if ( typchk == 1 ) {
2093
2094 /* Local parameters for ptt and k */
2095 hmax = mesh->info.hmax;
2096 hausd = mesh->info.hausd;
2097 isloc = 0;
2098
2099 if ( mesh->info.parTyp & MG_Tetra ) {
2100 for ( l=0; l<mesh->info.npar; ++l ) {
2101 par = &mesh->info.par[l];
2102
2103 if ( par->elt != MMG5_Tetrahedron ) continue;
2104 if ( par->ref != pt->ref ) continue;
2105
2106 hmax = par->hmax;
2107 hausd = par->hausd;
2108 isloc = 1;
2109 break;
2110 }
2111 }
2112 if ( mesh->info.parTyp & MG_Tria ) {
2113 if ( isloc ) {
2114 for ( l=0; l<mesh->info.npar; ++l ) {
2115 par = &mesh->info.par[l];
2116
2117 if ( par->elt != MMG5_Triangle ) continue;
2118 if ( par->ref != ptt->ref ) continue;
2119
2120 hmax = MG_MIN(hmax,par->hmax);
2121 hausd = MG_MIN(hausd,par->hausd);
2122 break;
2123 }
2124 }
2125 else {
2126 for ( l=0; l<mesh->info.npar; ++l ) {
2127 par = &mesh->info.par[l];
2128
2129 if ( par->elt != MMG5_Triangle ) continue;
2130 if ( par->ref != ptt->ref ) continue;
2131
2132 hmax = par->hmax;
2133 hausd = par->hausd;
2134 isloc = 1;
2135 break;
2136 }
2137 }
2138 }
2139
2140 ier = MMG5_chkedg(mesh,ptt,MG_GET(pxt->ori,i),hmax,hausd,isloc);
2141 if ( ier < 0 ) {
2142 /* Error */
2143 return 0;
2144 } else if ( ier == 0 ) {
2145 /* Nothing to split */
2146 return 1;
2147 }
2148
2149 /* put back flag on tetra */
2150 for (j=0; j<3; j++){
2151 if ( pxt->tag[MMG5_iarf[i][j]] & MG_REQ ) continue;
2152 if ( MG_GET(ptt->flag,j) ) MG_SET(pt->flag,MMG5_iarf[i][j]);
2153 }
2154 }
2155 else if ( typchk == 2 ) {
2156 for (j=0; j<3; j++) {
2157 ia = MMG5_iarf[i][j];
2158 if ( pxt->tag[ia] & MG_REQ ) continue;
2159 i1 = MMG5_iare[ia][0];
2160 i2 = MMG5_iare[ia][1];
2161 ip1 = pt->v[i1];
2162 ip2 = pt->v[i2];
2163 len = MMG5_lenedg(mesh,met,ia,pt);
2164
2165 assert( isfinite(len) && (len!=-len) );
2166
2167 // Case of an internal tetra with 4 ridges vertices.
2168 if ( len == 0 ) continue;
2169 if ( len > MMG3D_LLONG ) MG_SET(pt->flag,ia);
2170 /* Treat here the ridges coming from a corner (we can not do that after
2171 * because the corner don't have xpoints) */
2172 if ( (mesh->point[ip1].tag & MG_CRN) || (mesh->point[ip2].tag & MG_CRN) ) {
2173 if ( len > MMG3D_LOPTL ) MG_SET(pt->flag,ia);
2174 }
2175 }
2176 }
2177 return 1;
2178}
2179
2190static MMG5_int MMG3D_anatets_ani(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
2191 MMG5_pTetra pt;
2192 MMG5_pxTetra pxt;
2193 MMG5_Tria ptt;
2194 double len,lmax;
2195 double ux,uy,uz;
2196 MMG5_int ns,k,ip1,ip2;
2197 int ier,warn;
2198 int8_t imax,j,i,i1,i2;
2199
2200 assert ( met->m && met->size==6 );
2201
2202 if ( !MMG3D_hashTetra(mesh,1) ) {
2203 fprintf(stderr,"\n ## Hashing problem. Exit program.\n");
2204 return -1;
2205 }
2206
2207 ns = 0;
2208 warn = 0;
2209
2210 for (k=1; k<=mesh->ne; k++) {
2211 pt = &mesh->tetra[k];
2212 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) || !pt->xt ) continue;
2213
2214 /* check boundary face cut w/r Hausdorff or hmax */
2215 pt->flag = 0;
2216 pxt = &mesh->xtetra[pt->xt];
2217
2218 /* Travel well oriented boundary faces to flag edges that have to be cut */
2219 for (i=0; i<4; i++){
2220 if ( pxt->ftag[i] & MG_REQ ) continue;
2221 if ( !(pxt->ftag[i] & MG_BDY) ) continue;
2222
2223 if ( !MG_GET(pxt->ori,i) ) continue;
2224
2225 /* virtual triangle */
2226 assert( 0<=i && i<4 && "unexpected local face idx");
2227 MMG5_tet2tri(mesh,k,i,&ptt);
2228
2229 if ( !MMG3D_chkbdyface(mesh,met,k,pt,pxt,i,&ptt,typchk) ) { continue; }
2230 }
2231
2233 imax = 6;
2234 lmax = 0.;
2235 for ( j=0; j<6; ++j ) {
2236 if ( !MG_GET(pt->flag,j) ) { continue; }
2237
2238 i1 = MMG5_iare[j][0];
2239 i2 = MMG5_iare[j][1];
2240 ip1 = pt->v[i1];
2241 ip2 = pt->v[i2];
2242
2243 ux = mesh->point[ip1].c[0] - mesh->point[ip2].c[0];
2244 uy = mesh->point[ip1].c[1] - mesh->point[ip2].c[1];
2245 uz = mesh->point[ip1].c[2] - mesh->point[ip2].c[2];
2246
2247 len = ux*ux + uy*uy + uz*uz;
2248 if ( len <= lmax) { continue; }
2249 lmax = len;
2250 imax = j;
2251 }
2252
2253 pt->flag = 0;
2254 if ( imax < 6 ) { MG_SET(pt->flag,imax); }
2255
2256 if ( !pt->flag ) continue;
2257
2258 ier = MMG3D_splsurfedge( mesh,met,k,pt,pxt,imax,typchk,1,&warn );
2259
2260 if ( ier==-1 ) { return -1; }
2261 else if ( !ier ) { continue; }
2262 else if ( ier==2 ) {
2263 /* Unable to split due to lack of memory */
2264 return ns;
2265 }
2266
2267 ++ns;
2268 }
2269
2270 return ns;
2271}
2272
2289static MMG5_int
2291 MMG5_pTetra pt;
2292 MMG5_pPoint ppt;
2293 MMG5_Tria ptt,ptt2;
2294 MMG5_xTetra *pxt;
2295 MMG5_xPoint *pxp;
2296 MMG5_Bezier pb,pb2;
2297 MMG5_Hash hash;
2298 double o[3],no[3],to[3],dd;
2299 int ic,it,ier;
2300 MMG5_int ip,vx[6],src,nc,ns,ni,ne,k,ip1,ip2,nap,ixp1,ixp2;
2301 int8_t i,j,j2,ia,i1,i2,ifac,intnom;
2302 static double uv[3][2] = { {0.5,0.5}, {0.,0.5}, {0.5,0.} };
2303 static int8_t mmgWarn = 0, mmgWarn2 = 0;
2304
2306 if ( !MMG5_hashNew(mesh,&hash,mesh->np,7*mesh->np) ) return -1;
2307 ns = nap = 0;
2308
2309 for (k=1; k<=mesh->ne; k++) {
2310 pt = &mesh->tetra[k];
2311 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) || !pt->xt ) continue;
2312
2313 /* check boundary face cut w/r Hausdorff or hmax */
2314 pt->flag = 0;
2315 pxt = &mesh->xtetra[pt->xt];
2316
2317 for (i=0; i<4; i++){
2318 if ( pxt->ftag[i] & MG_REQ ) continue;
2319 if ( !(pxt->ftag[i] & MG_BDY) ) continue;
2320
2321 if ( typchk == 1 ) {
2322 if ( !MG_GET(pxt->ori,i) ) { continue; }
2323 }
2324
2325 /* virtual triangle */
2326 assert( 0<=i && i<4 && "unexpected local face idx");
2327 MMG5_tet2tri(mesh,k,i,&ptt);
2328
2329 if ( !MMG3D_chkbdyface(mesh,met,k,pt,pxt,i,&ptt,typchk) ) { continue; }
2330
2331 if ( !pt->flag ) continue;
2332 ns++;
2333
2334 /* geometric support */
2335 ier = MMG5_bezierCP(mesh,&ptt,&pb,MG_GET(pxt->ori,i));
2336 assert(ier);
2337
2338 /* scan edges in face to split */
2339 for (j=0; j<3; j++) {
2340 ia = MMG5_iarf[i][j];
2341 if ( !MG_GET(pt->flag,ia) ) continue;
2342 if ( pxt->tag[ia] & MG_REQ ) continue;
2343 i1 = MMG5_iare[ia][0];
2344 i2 = MMG5_iare[ia][1];
2345 ip1 = pt->v[i1];
2346 ip2 = pt->v[i2];
2347 ip = MMG5_hashGet(&hash,ip1,ip2);
2348
2349 /* new point along edge */
2350 if ( !ip ) {
2351
2352 if ( ptt.tag[j] & MG_NOM ) {
2353 /* Use BezierNom to ensure that the new nm point has a normal that
2354 * is consistent with the normals at point ip1 and ip2 */
2355 if ( !MMG5_BezierNom(mesh,ip1,ip2,0.5,o,no,to) ) {
2356 continue;
2357 }
2358 }
2359 else {
2360 /* Otherwise we can build the bezier patch */
2361 ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
2362 assert(ier);
2363 }
2364
2365#ifdef USE_POINTMAP
2366 src = mesh->point[ip1].src;
2367#else
2368 src = 1;
2369#endif
2370 ip = MMG3D_newPt(mesh,o,MG_BDY,src);
2371 if ( !ip ) {
2372 /* reallocation of point table */
2374 fprintf(stderr,"\n ## Error: %s: unable to"
2375 " allocate a new point.\n",__func__);
2377 MMG3D_delPatternPts(mesh,hash);return -1;
2378 ,o,MG_BDY,src);
2379 // Now pb->p contain a wrong memory address.
2380 pb.p[0] = &mesh->point[ptt.v[0]];
2381 pb.p[1] = &mesh->point[ptt.v[1]];
2382 pb.p[2] = &mesh->point[ptt.v[2]];
2383 }
2384 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
2385 ppt = &mesh->point[ip];
2386
2387 assert ( met );
2388
2389 // Add MG_BDY tag before interpolation
2390 pxt->tag[ia] |= MG_BDY;
2391
2392 if ( met->m ) {
2393 if ( typchk == 1 && (met->size>1) )
2394 ier = MMG3D_intmet33_ani(mesh,met,k,ia,ip,0.5);
2395 else
2396 ier = MMG5_intmet(mesh,met,k,ia,ip,0.5);
2397
2398 if ( !ier ) {
2399 if ( !mmgWarn ) {
2400 fprintf(stderr,"\n ## Error: %s: unable to interpolate at least"
2401 " 1 metric.\n",__func__);
2402 mmgWarn = 1;
2403 }
2404 return -1;
2405 }
2406 else if ( ier < 0 ) {
2407 MMG3D_delPt(mesh,ip);
2408 continue;
2409 }
2410 }
2411
2412 if ( MG_EDG_OR_NOM(ptt.tag[j]) )
2413 ppt->ref = ptt.edg[j] ? ptt.edg[j] : ptt.ref;
2414 else
2415 ppt->ref = ptt.ref;
2416 ppt->tag |= ptt.tag[j];
2417 pxp = &mesh->xpoint[ppt->xp];
2418
2419 /* Update normal and tangent vectors */
2420 intnom = 0;
2421 if ( ptt.tag[j] & MG_NOM ) {
2422 ixp1 = mesh->point[ip1].xp;
2423 ixp2 = mesh->point[ip2].xp;
2424 intnom = ( mesh->xpoint[ixp1].nnor ) || ( mesh->xpoint[ixp2].nnor );
2425 }
2426 if ( intnom ) pxp->nnor = 1;
2427 else memcpy(pxp->n1,no,3*sizeof(double));
2428
2429 memcpy(ppt->n,to,3*sizeof(double));
2430
2431 if ( mesh->info.fem<typchk ) {
2432 /* A ridge can be at the interface of 2 boundary faces of the same
2433 * tetra: second normal has to be computed */
2434 if ( MG_EDG(ptt.tag[j]) && !(ptt.tag[j] & MG_NOM) ) {
2435 /* Update the second normal and the tangent at point ip if the edge
2436 * is shared by 2 faces (if anatet4 is not called, 1 tetra may have
2437 * 2 bdry faces) */
2438 ifac = (MMG5_ifar[ia][0] == i)? MMG5_ifar[ia][1] : MMG5_ifar[ia][0];
2439 if ( pxt->ftag[ifac] & MG_BDY ) {
2440 j2 = MMG5_iarfinv[ifac][ia];
2441
2442 /* Compute tangent and normal with respect to the face ifac */
2443 /* virtual triangle */
2444 assert( 0<=ifac && ifac<4 && "unexpected local face idx");
2445 MMG5_tet2tri(mesh,k,ifac,&ptt2);
2446
2447 /* geometric support */
2448 ier = MMG5_bezierCP(mesh,&ptt2,&pb2,MG_GET(pxt->ori,ifac));
2449 assert(ier);
2450
2451 ier = MMG3D_bezierInt(&pb2,&uv[j2][0],o,no,to);
2452 assert(ier);
2453
2454 if ( !MMG3D_update_rid_geom(ppt,pxp,no) ) continue;
2455 }
2456 }
2457 }
2458 nap++;
2459 }
2460 else if ( MG_EDG(ptt.tag[j]) && !(ptt.tag[j] & MG_NOM) ) {
2461 /* Point at the interface of 2 boundary faces belonging to different
2462 * tetra : Point has already been created from another tetra so we have
2463 * to store the tangent and the second normal at edge */
2464 ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
2465 assert(ier);
2466
2467 ppt = &mesh->point[ip];
2468 assert(ppt->xp);
2469 pxp = &mesh->xpoint[ppt->xp];
2470
2471 if ( !MMG3D_update_rid_geom(ppt,pxp,no) ) continue;
2472 }
2473 }
2474 }
2475 }
2476
2477 if ( !ns ) {
2478 MMG5_DEL_MEM(mesh,hash.item);
2479 return ns;
2480 }
2481
2484 nc = 0;
2485 for (k=1; k<=mesh->ne; k++) {
2486 pt = &mesh->tetra[k];
2487 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
2488 pxt = pt->xt ? &mesh->xtetra[pt->xt] : 0;
2489
2490 /* update face-edge flag */
2491 for (i=0; i<4; i++) {
2492 /* virtual triangle */
2493 memset(&ptt,0,sizeof(MMG5_Tria));
2494 if ( pt->xt && pxt->ftag[i] && pxt->ftag[i] != MG_OLDPARBDY ) {
2495 assert( 0<=i && i<4 && "unexpected local face idx");
2496 MMG5_tet2tri(mesh,k,i,&ptt);
2497 }
2498
2499 for (j=0; j<3; j++) {
2500 ia = MMG5_iarf[i][j];
2501
2502 /* If edge is already marked, nothing to do */
2503 if ( MG_GET(pt->flag,ia) ) {
2504 continue;
2505 }
2506
2508 /* First: skip edge if required */
2509 if ( pt->xt && (pxt->tag[ia] & MG_REQ) ) continue;
2510
2511 /* Second: if possible treat manifold ridges from a boundary face (to
2512 * ensure the computation of n2) */
2513 if ( pt->xt ) {
2514 if ( pxt->tag[ia] & MG_GEO && (!(pxt->tag[ia] & MG_NOM) ) ) {
2515 int ifac2 = (MMG5_ifar[ia][0]!=i)? MMG5_ifar[ia][1] : MMG5_ifar[ia][0];
2516 if ( (!(pxt->ftag[i] & MG_BDY) ) && (pxt->ftag[ifac2] & MG_BDY) ) {
2517 assert ( ifac2 > i && "lower face is already processed" );
2518 continue;
2519 }
2520 }
2521 }
2522
2523 ip1 = pt->v[MMG5_iare[ia][0]];
2524 ip2 = pt->v[MMG5_iare[ia][1]];
2525 ip = MMG5_hashGet(&hash,ip1,ip2);
2526 if ( ip > 0 ) {
2527
2528 MG_SET(pt->flag,ia);
2529 nc++;
2530
2531 if ( (!(ptt.tag[j] & MG_GEO)) || (ptt.tag[j] & MG_NOM) ) continue;
2532
2533 /* From a boundary face of a boundary tetra we can update normal/tangent
2534 at manifold ridges; */
2535 ppt = &mesh->point[ip];
2536 assert(ppt->xp);
2537 pxp = &mesh->xpoint[ppt->xp];
2538 if ( pt->xt ) ier = MMG5_bezierCP(mesh,&ptt,&pb,MG_GET(pxt->ori,i));
2539 else ier = MMG5_bezierCP(mesh,&ptt,&pb,1);
2540 assert(ier);
2541
2542 ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
2543 assert(ier);
2544
2545 dd = no[0]*pxp->n1[0]+no[1]*pxp->n1[1]+no[2]*pxp->n1[2];
2546 if ( dd > 1.0-MMG5_EPS ) continue;
2547
2548 memcpy(pxp->n2,no,3*sizeof(double));
2549 }
2550 }
2551 }
2552 }
2553 if ( mesh->info.ddebug && nc ) {
2554 fprintf(stdout," %" MMG5_PRId " added\n",nc);
2555 fflush(stdout);
2556 }
2557
2559 for (k=1; k<=mesh->np; k++)
2560 mesh->point[k].flag = 0;
2561
2562 it = 1;
2563 nc = 0;
2564 do {
2565 ni = 0;
2566 for (k=1; k<=mesh->ne; k++) {
2567 pt = &mesh->tetra[k];
2568 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) || !pt->flag ) continue;
2569 memset(vx,0,6*sizeof(MMG5_int));
2570 pt->flag = ic = 0;
2571 for (ia=0,i=0; i<3; i++) {
2572 for (j=i+1; j<4; j++,ia++) {
2573 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_REQ) ) continue;
2574 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
2575 if ( vx[ia] > 0 ) {
2576 MG_SET(pt->flag,ia);
2577 if ( mesh->point[vx[ia]].flag > 2 ) ic = 1;
2578 }
2579 }
2580 }
2581 if ( !pt->flag ) continue;
2582
2583 switch (pt->flag) {
2584 case 1: case 2: case 4: case 8: case 16: case 32:
2585 ier = MMG3D_split1_sim(mesh,met,k,vx);
2586 break;
2587 case 48: case 24: case 40: case 6: case 34: case 36:
2588 case 20: case 5: case 17: case 9: case 3: case 10:
2589 ier =MMG3D_split2sf_sim(mesh,met,k,vx);
2590 break;
2591 case 33: case 18: case 12:
2592 ier = MMG3D_split2_sim(mesh,met,k,vx);
2593 break;
2594 case 11: case 21: case 38: case 56:
2595 ier = MMG3D_split3_sim(mesh,met,k,vx);
2596 break;
2597 case 7: case 25: case 42: case 52:
2598 ier =MMG3D_split3cone_sim(mesh,met,k,vx);
2599 break;
2600 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
2601 case 14: case 49: case 50: case 44: case 41:
2602 ier = MMG3D_split3op_sim(mesh,met,k,vx);
2603 break;
2604 case 23: case 29: case 53: case 60: case 57: case 58:
2605 case 27: case 15: case 43: case 39: case 54: case 46:
2606 ier = MMG3D_split4sf_sim(mesh,met,k,vx);
2607 break;
2608 case 30: case 45: case 51:
2609 ier = MMG3D_split4op_sim(mesh,met,k,vx);
2610 break;
2611 case 62: case 61: case 59: case 55: case 47: case 31:
2612 ier = MMG3D_split5_sim(mesh,met,k,vx);
2613 break;
2614 case 63:
2615 ier = MMG3D_split6_sim(mesh,met,k,vx);
2616 break;
2617 }
2618 if ( ier ) continue;
2619
2620 ni++;
2621 if ( ic == 0 && MMG3D_dichoto(mesh,met,k,vx) ) {
2622 for (ia=0; ia<6; ia++)
2623 if ( vx[ia] > 0 ) {
2624 mesh->point[vx[ia]].flag++;
2625 }
2626 }
2627 else {
2628 for (ia=0; ia<6; ++ia ) {
2629 if ( vx[ia] > 0 ) {
2630 MMG5_hashPop(&hash,pt->v[MMG5_iare[ia][0]],pt->v[MMG5_iare[ia][1]]);
2631 MMG3D_delPt(mesh,vx[ia]);
2632
2633 if ( (mesh->info.ddebug || mesh->info.imprim > 5) && !mmgWarn2 ) {
2634 fprintf(stderr,"\n ## Warning: %s: surfacic pattern: unable to find"
2635 " a valid split for at least 1 point. Point(s) deletion.\n",
2636 __func__ );
2637 mmgWarn2 = 1;
2638 }
2639
2640 }
2641 }
2642 }
2643 }
2644 nc += ni;
2645 }
2646 while( ni > 0 && ++it < 40 );
2647
2648 if ( mesh->info.ddebug && nc ) {
2649 fprintf(stdout," %" MMG5_PRId " corrected, %" MMG5_PRId " invalid\n",nc,ni);
2650 fflush(stdout);
2651 }
2652
2653
2655 ns = 0;
2656 ne = mesh->ne;
2657 for (k=1; k<=ne; k++) {
2658 pt = &mesh->tetra[k];
2659 if ( !MG_EOK(pt) || !pt->flag || (pt->tag & MG_REQ) ) continue;
2660 memset(vx,0,6*sizeof(MMG5_int));
2661 for (ia=0,i=0; i<3; i++) {
2662 for (j=i+1; j<4; j++,ia++) {
2663 if ( MG_GET(pt->flag,ia) ) {
2664 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
2665 assert(vx[ia]);
2666 }
2667 }
2668 }
2669 switch (pt->flag) {
2670 case 1: case 2: case 4: case 8: case 16: case 32: /* 1 edge split */
2671 if ( ! MMG5_split1(mesh,met,k,vx,typchk-1) ) return -1;
2672 ns++;
2673 break;
2674 case 48: case 24: case 40: case 6: case 34: case 36:
2675 case 20: case 5: case 17: case 9: case 3: case 10: /* 2 edges (same face) split */
2676 if ( ! MMG5_split2sf(mesh,met,k,vx,typchk-1) ) return -1;
2677 ns++;
2678 break;
2679
2680 case 33: case 18: case 12: /* 2 opposite edges split */
2681 if ( ! MMG5_split2(mesh,met,k,vx,typchk-1) ) return -1;
2682 ns++;
2683 break;
2684
2685 case 11: case 21: case 38: case 56: /* 3 edges on the same faces splitted */
2686 if ( ! MMG5_split3(mesh,met,k,vx,typchk-1) ) return -1;
2687 ns++;
2688 break;
2689
2690 case 7: case 25: case 42: case 52: /* 3 edges on conic configuration splitted */
2691 if ( ! MMG5_split3cone(mesh,met,k,vx,typchk-1) ) return -1;
2692 ns++;
2693 break;
2694
2695 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
2696 case 14: case 49: case 50: case 44: case 41: /* 3 edges on opposite configuration splitted */
2697 if ( ! MMG5_split3op(mesh,met,k,vx,typchk-1) ) return -1;
2698 ns++;
2699 break;
2700
2701 case 23: case 29: case 53: case 60: case 57: case 58:
2702 case 27: case 15: case 43: case 39: case 54: case 46: /* 4 edges with 3 lying on the same face splitted */
2703 if ( ! MMG5_split4sf(mesh,met,k,vx,typchk-1) ) return -1;
2704 ns++;
2705 break;
2706
2707 /* 4 edges with no 3 lying on the same face splitted */
2708 case 30: case 45: case 51:
2709 if ( ! MMG5_split4op(mesh,met,k,vx,typchk-1) ) return -1;
2710 ns++;
2711 break;
2712
2713 case 62: case 61: case 59: case 55: case 47: case 31: /* 5 edges split */
2714 if ( ! MMG5_split5(mesh,met,k,vx,typchk-1) ) return -1;
2715 ns++;
2716 break;
2717
2718 case 63: /* 6 edges split */
2719 if ( ! MMG5_split6(mesh,met,k,vx,typchk-1) ) return -1;
2720 ns++;
2721 break;
2722 }
2723 }
2724 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
2725 fprintf(stdout," %7" MMG5_PRId " elements splitted\n",nap);
2726
2727 MMG5_DEL_MEM(mesh,hash.item);
2728 return nap;
2729}
2730
2731static MMG5_int (*MMG3D_anatets)(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk);
2732
2755static int MMG3D_anatet4_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,int8_t metRidTyp,
2756 int *ifac,int* conf0,MMG5_int *adj,int *conf1) {
2757 MMG5_pTetra pt,pt1,ptnew;
2758 MMG5_pxTetra pxt0,pxt1;
2759 MMG5_pPoint ppt,ppt0;
2760 double calold0,calold,calnew,calnew0,calnew1,calnew2,calnew3;
2761 double worst_split4bar_cal,worst_swap_cal,cb[4];
2762 int loc_conf0,loc_conf1;
2763 MMG5_int *adja,k1,np;
2764 int nbdy,i,j0,j1;
2765 uint8_t tau0[4],tau1[4];
2766
2767 pt = &mesh->tetra[k];
2768 calold0 = pt->qual;
2769
2770 assert ( pt->xt );
2771
2772 pxt0 = &mesh->xtetra[pt->xt];
2773
2774 ptnew = &mesh->tetra[0];
2775
2777 if ( !mesh->info.noinsert ) {
2778 ppt0 = &mesh->point[0];
2779 memset(ppt0,0,sizeof(MMG5_Point));
2780
2781 for (i=0; i<4; i++) {
2782 ppt = &mesh->point[pt->v[i]];
2783 ppt0->c[0] += ppt->c[0];
2784 ppt0->c[1] += ppt->c[1];
2785 ppt0->c[2] += ppt->c[2];
2786 }
2787 ppt0->c[0] *= 0.25;
2788 ppt0->c[1] *= 0.25;
2789 ppt0->c[2] *= 0.25;
2790
2791 cb[0] = 0.25; cb[1] = 0.25; cb[2] = 0.25; cb[3] = 0.25;
2792
2793 if ( met->m ) {
2794 if ( !metRidTyp && met->size > 1 )
2795 MMG5_interp4bar33_ani(mesh,met,k,0,cb);
2796 else
2797 MMG5_interp4bar(mesh,met,k,0,cb);
2798 }
2799
2800 memcpy(ptnew,pt,sizeof(MMG5_Tetra));
2801 ptnew->v[0] = 0;
2802 if ( (!metRidTyp) && met->m && met->size>1 )
2803 calnew0 = MMG5_caltet33_ani(mesh,met,ptnew);
2804 else
2805 calnew0 = MMG5_orcal(mesh,met,0);
2806
2807 ptnew->v[0] = pt->v[0];
2808 ptnew->v[1] = 0;
2809 if ( (!metRidTyp) && met->m && met->size>1 )
2810 calnew1 = MMG5_caltet33_ani(mesh,met,ptnew);
2811 else
2812 calnew1 = MMG5_orcal(mesh,met,0);
2813
2814 ptnew->v[1] = pt->v[1];
2815 ptnew->v[2] = 0;
2816 if ( (!metRidTyp) && met->m && met->size>1 )
2817 calnew2 = MMG5_caltet33_ani(mesh,met,ptnew);
2818 else
2819 calnew2 = MMG5_orcal(mesh,met,0);
2820
2821 ptnew->v[2] = pt->v[2];
2822 ptnew->v[3] = 0;
2823 if ( (!metRidTyp) && met->m && met->size>1 )
2824 calnew3 = MMG5_caltet33_ani(mesh,met,ptnew);
2825 else
2826 calnew3 = MMG5_orcal(mesh,met,0);
2827
2828 worst_split4bar_cal = MG_MIN(MG_MIN(calnew0,calnew1),MG_MIN(calnew2,calnew3));
2829 }
2830 else
2831 worst_split4bar_cal = 0.;
2832
2834 *ifac = -1;
2835 *adj = -1;
2836 if ( !mesh->info.noswap ) {
2837 worst_swap_cal = 0.;
2838
2839 for (j0=0; j0<4; j0++) {
2840 if ( pxt0->ftag[j0] & MG_BDY ) continue;
2841
2843 adja = &mesh->adja[4*(k-1)+1];
2844 k1 = adja[j0]/4;
2845 j1 = adja[j0]%4;
2846
2847 assert(k1);
2848
2849 /* Search in which configurations are the tetrahedra (default is case 0-0)
2850 *
2851 * 3 2------------- 0
2852 * ,/|`\ |`\ /|
2853 * ,/ | `\ | `\ /.|
2854 * ,/ '. `\ '. `\ / |
2855 * ,/ | `\ | `\ / .|
2856 * ,/ | `\ | /\.|
2857 * 0-----------'.--------2 ' / 3
2858 * `\. | ,/ | / ,/
2859 * `\. | ,/ | /,/
2860 * `\. '. ,/ '. ,/
2861 * `\. |/ |/
2862 * `1 1
2863 */
2864
2865 /* k may be in configuration 0, 3, 6 or 9. Default is case 0 */
2866 loc_conf0 = 3*j0;
2867
2868 switch(loc_conf0) {
2869 case 3:
2870 tau0[0] = 1; tau0[1] = 0; tau0[2] = 3; tau0[3] = 2;
2871 break;
2872 case 6:
2873 tau0[0] = 2; tau0[1] = 0; tau0[2] = 1; tau0[3] = 3;
2874 break;
2875 case 9:
2876 tau0[0] = 3; tau0[1] = 0; tau0[2] = 2; tau0[3] = 1;
2877 break;
2878 default:
2879 tau0[0] = 0; tau0[1] = 1; tau0[2] = 2; tau0[3] = 3;
2880 }
2881
2882 /* k1 may be in configuration j1, j1+1, j1+2 */
2883 pt1 = &mesh->tetra[k1];
2884
2885 if ( pt1->tag & MG_REQ ) continue;
2886
2887 assert(pt->ref == pt1->ref);
2888 for ( i=0; i<3; ++i )
2889 if ( pt->v[MMG5_idir[j0][0]] == pt1->v[MMG5_idir[j1][i]] ) break;
2890
2891 assert(i<3);
2892 loc_conf1 = 3*j1+i;
2893
2894 switch(loc_conf1) {
2895 case 1:
2896 tau1[0] = 0; tau1[1] = 2; tau1[2] = 3; tau1[3] = 1;
2897 break;
2898 case 2:
2899 tau1[0] = 0; tau1[1] = 3; tau1[2] = 1; tau1[3] = 2;
2900 break;
2901 case 3:
2902 tau1[0] = 1; tau1[1] = 0; tau1[2] = 3; tau1[3] = 2;
2903 break;
2904 case 4:
2905 tau1[0] = 1; tau1[1] = 3; tau1[2] = 2; tau1[3] = 0;
2906 break;
2907 case 5:
2908 tau1[0] = 1; tau1[1] = 2; tau1[2] = 0; tau1[3] = 3;
2909 break;
2910 case 6:
2911 tau1[0] = 2; tau1[1] = 0; tau1[2] = 1; tau1[3] = 3;
2912 break;
2913 case 7:
2914 tau1[0] = 2; tau1[1] = 1; tau1[2] = 3; tau1[3] = 0;
2915 break;
2916 case 8:
2917 tau1[0] = 2; tau1[1] = 3; tau1[2] = 0; tau1[3] = 1;
2918 break;
2919 case 9:
2920 tau1[0] = 3; tau1[1] = 0; tau1[2] = 2; tau1[3] = 1;
2921 break;
2922 case 10:
2923 tau1[0] = 3; tau1[1] = 2; tau1[2] = 1; tau1[3] = 0;
2924 break;
2925 case 11:
2926 tau1[0] = 3; tau1[1] = 1; tau1[2] = 0; tau1[3] = 2;
2927 break;
2928 default:
2929 tau1[0] = 0; tau1[1] = 1; tau1[2] = 2; tau1[3] = 3;
2930 }
2931
2933 if ( pt1->xt ) {
2934 pxt1 = &mesh->xtetra[pt1->xt];
2935
2936 nbdy = 0;
2937 if ( pxt1->ftag[tau1[1]] ) ++nbdy;
2938 if ( pxt0->ftag[tau0[1]] ) ++nbdy;
2939 if ( nbdy > 1 ) continue;
2940
2941 nbdy = 0;
2942 if ( pxt1->ftag[tau1[3]] ) ++nbdy;
2943 if ( pxt0->ftag[tau0[2]] ) ++nbdy;
2944 if ( nbdy > 1 ) continue;
2945
2946 nbdy = 0;
2947 if ( pxt1->ftag[tau1[2]] ) ++nbdy;
2948 if ( pxt0->ftag[tau0[3]] ) ++nbdy;
2949 if ( nbdy > 1 ) continue;
2950
2951 }
2952
2954 calold = MG_MIN(calold0,pt1->qual);
2955
2956 ptnew = &mesh->tetra[0];
2957 memcpy(ptnew,pt,sizeof(MMG5_Tetra));
2958 np = pt1->v[tau1[0]];
2959
2960 ptnew->v[tau0[1]] = np;
2961 if ( (!metRidTyp) && met->m && met->size>1 )
2962 calnew0 = MMG5_caltet33_ani(mesh,met,ptnew);
2963 else
2964 calnew0 = MMG5_orcal(mesh,met,0);
2965 if ( calnew0 < MMG5_NULKAL ) continue;
2966
2967 ptnew->v[tau0[1]] = pt->v[tau0[1]];
2968 ptnew->v[tau0[2]] = np;
2969 if ( (!metRidTyp) && met->m && met->size>1 )
2970 calnew1 = MMG5_caltet33_ani(mesh,met,ptnew);
2971 else
2972 calnew1 = MMG5_orcal(mesh,met,0);
2973 if ( calnew1 < MMG5_NULKAL ) continue;
2974
2975 ptnew->v[tau0[2]] = pt->v[tau0[2]];
2976 ptnew->v[tau0[3]] = np;
2977 if ( (!metRidTyp) && met->m && met->size>1 )
2978 calnew2 = MMG5_caltet33_ani(mesh,met,ptnew);
2979 else
2980 calnew2 = MMG5_orcal(mesh,met,0);
2981 if ( calnew2 < MMG5_NULKAL ) continue;
2982
2983 calnew = MG_MIN(calnew0,MG_MIN(calnew1,calnew2));
2984
2985 if ( calold < MMG5_EPSOK ) {
2986 if ( calnew < calold ) continue;
2987 }
2988 else if ( calnew <= MMG5_EPSOK ) continue;
2989
2990 if ( calnew > worst_swap_cal ) {
2991 worst_swap_cal = calnew;
2992 *ifac = j0;
2993 *adj = adja[j0];
2994 *conf0 = loc_conf0;
2995 *conf1 = loc_conf1;
2996 }
2997 }
2998 }
2999 else
3000 worst_swap_cal = 0.;
3001
3002 if ( *ifac < 0 ) {
3003 /* No valid config for the swap23 */
3004 if ( worst_split4bar_cal < MMG5_EPSOK ) return 0;
3005 return 1;
3006 }
3007
3008 assert ( *adj > 0 );
3009 if ( worst_swap_cal < worst_split4bar_cal ) {
3010 if ( worst_split4bar_cal < MMG5_EPSOK ) return 0;
3011 return 1;
3012 }
3013
3014 if ( worst_swap_cal < MMG5_EPSOK ) return 0;
3015
3016 return 2;
3017}
3018
3029static MMG5_int MMG5_anatet4(MMG5_pMesh mesh, MMG5_pSol met,MMG5_int *nf, int8_t typchk) {
3030 MMG5_pTetra pt;
3031 MMG5_pPoint ppt;
3032 MMG5_pxTetra pxt;
3033 int conf0,conf1,ifac,id_op;
3034 MMG5_int ier,ns,k,adj;
3035 int8_t nbdy,j;
3036#ifndef NDEBUG
3037 static int8_t mmgWarn=0;
3038#endif
3039
3040 ns = 0;
3041 for (k=1; k<=mesh->ne; k++) {
3042 pt = &mesh->tetra[k];
3043 if ( !MG_EOK(pt) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
3044 nbdy = 0;
3045 if ( pt->xt ) {
3046 pxt = &mesh->xtetra[pt->xt];
3047 for (j=0; j<4; j++)
3048 if ( ( pxt->ftag[j] & MG_BDY ) && (!(pxt->ftag[j] & MG_PARBDY)) ) nbdy++;
3049 }
3050
3051 /* Check for the number of boundary faces */
3052 if ( nbdy > 1 ) {
3053 id_op = MMG3D_anatet4_sim(mesh,met,k,typchk-1,&ifac,&conf0,&adj,&conf1);
3054 if ( !id_op ) {
3055#ifndef NDEBUG
3056 if ( (!(mesh->info.noswap && mesh->info.noinsert)) && !mmgWarn ) {
3057 mmgWarn=1;
3058 printf("\n ## Warning: %s: unable to swap or split at least 1 tetra"
3059 " with multiple boundary faces.\n",__func__);
3060 }
3061#endif
3062 continue;
3063 }
3064
3065 if ( id_op==1 ) {
3066 /* Split */
3067 ier = MMG5_split4bar(mesh,met,k,typchk-1);
3068 if ( !ier ) return -1;
3069 ns++;
3070 }
3071 else {
3072 assert ( id_op==2 );
3073
3074 /* Swap */
3075 ier = MMG3D_swap23(mesh,met,k,typchk-1,ifac,conf0,adj,conf1);
3076 if ( ier < 0 ) {
3077 return -1;
3078 }
3079 else if ( ier ) ++(*nf);
3080 }
3081 }
3082 /* Check for the number of boundary vertices */
3083 else {
3084 nbdy = 0;
3085 for (j=0; j<4; j++) {
3086 ppt = &mesh->point[pt->v[j]];
3087 if ( (ppt->tag & MG_BDY) && (!(ppt->tag & MG_PARBDY)) ) nbdy++;
3088 }
3089 if ( nbdy == 4 ) {
3090 if ( !mesh->info.noinsert ) {
3091 ier = MMG5_split4bar(mesh,met,k,typchk-1);
3092 if ( !ier ) return -1;
3093 ns++;
3094 }
3095 }
3096 }
3097 }
3098
3099 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
3100 fprintf(stdout," boundary elements: %7" MMG5_PRId " splitted %7" MMG5_PRId " swapped\n",ns,*nf);
3101 return ns;
3102}
3103
3114static MMG5_int MMG5_anatet4rid(MMG5_pMesh mesh, MMG5_pSol met,MMG5_int *nf, int8_t typchk) {
3115 MMG5_pTetra pt;
3116 MMG5_pPoint ppt;
3117 MMG5_int ier;
3118 MMG5_int ns,k;
3119 int8_t nrid,j;
3120
3121 ns = 0;
3122 for (k=1; k<=mesh->ne; k++) {
3123 pt = &mesh->tetra[k];
3124 if ( !MG_EOK(pt) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
3125 nrid = 0;
3126
3127 for (j=0; j<4; j++) {
3128 ppt = &mesh->point[pt->v[j]];
3129 if ( MG_RID(ppt->tag) ) {
3130 /* Non-singular ridge point: metric ridge */
3131 nrid++;
3132 }
3133 }
3134 if ( nrid == 4 ) {
3135 if ( !mesh->info.noinsert ) {
3136 ier = MMG5_split4bar(mesh,met,k,typchk-1);
3137 if ( !ier ) return -1;
3138 ns++;
3139 }
3140 }
3141 }
3142 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
3143 fprintf(stdout," boundary elements: %7" MMG5_PRId " splitted\n",ns);
3144 return ns;
3145}
3146
3147
3159int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
3160 int it,minit,maxit,lastit;
3161 MMG5_int nc,ns,nnc,nns,nnf,ier,nf;
3162
3163 /* pointer to the suitable anatets function */
3164 if ( met->m && met->size==6 ) {
3165 /* if the aniso metric is not compatible with the geometry, the non
3166 * conformal surface operators may create spurious ridges */
3168 }
3169 else {
3171 }
3172
3173 /* analyze tetras : initial splitting */
3174 nns = nnc = nnf = it = 0;
3175 lastit = 0;
3176 minit = 3;
3177 maxit = 6;
3178 mesh->gap = 0.5;
3179 do {
3180 if ( typchk==2 && lastit==1 ) ++mesh->info.fem;
3181
3182
3183 /* split or swap tetra with more than 2 bdry faces */
3184 nf = ier = 0;
3185 if ( mesh->info.fem == typchk ) {
3186 ier = MMG5_anatet4(mesh,met,&nf,typchk);
3187 if ( ier < 0 ) return 0;
3188 }
3189 else if ( (met->size==6) && (typchk == 1) && lastit ) {
3190 ier = MMG5_anatet4rid(mesh,met,&nf,typchk);
3191 if ( ier < 0 ) return 0;
3192 }
3193 else ier = 0;
3194 ns = ier;
3195
3196 /* memory free */
3197 if ( mesh->adja )
3199
3200 if ( !mesh->info.noinsert ) {
3201 /* analyze surface tetras */
3202 ier = MMG3D_anatets(mesh,met,typchk);
3203
3204 if ( ier < 0 ) {
3205 fprintf(stderr,"\n ## Unable to complete surface mesh. Exit program.\n");
3206 return 0;
3207 }
3208 ns += ier;
3209
3210 if ( patternMode ) {
3211 if ( mesh->adja ) { MMG5_DEL_MEM(mesh,mesh->adja); }
3212
3213 /* analyze internal tetras */
3214 ier = MMG5_anatetv(mesh,met,typchk);
3215 if ( ier < 0 ) {
3216 fprintf(stderr,"\n ## Unable to complete volume mesh. Exit program.\n");
3217 return 0;
3218 }
3219 ns += ier;
3220 }
3221 }
3222
3223 if ( (!mesh->adja) && (!MMG3D_hashTetra(mesh,1)) ) {
3224 fprintf(stderr,"\n ## Hashing problem. Exit program.\n");
3225 return 0;
3226 }
3227
3228 /* collapse short edges */
3229 if ( !mesh->info.noinsert ) {
3230 nc = MMG5_coltet(mesh,met,typchk);
3231 if ( nc < 0 ) {
3232 fprintf(stderr,"\n ## Unable to collapse mesh. Exiting.\n");
3233 return 0;
3234 }
3235 }
3236 else nc = 0;
3237
3238 /* attempt to swap */
3239 if ( !mesh->info.noswap ) {
3240 ier = MMG5_swpmsh(mesh,met,NULL,typchk);
3241 if ( ier < 0 ) {
3242 fprintf(stderr,"\n ## Unable to improve mesh. Exiting.\n");
3243 return 0;
3244 }
3245 nf += ier;
3246
3248 if ( ier < 0 ) {
3249 fprintf(stderr,"\n ## Unable to improve mesh. Exiting.\n");
3250 return 0;
3251 }
3252 nf += ier;
3253 }
3254 else nf = 0;
3255
3256 nnc += nc;
3257 nns += ns;
3258 nnf += nf;
3259 if ( (abs(mesh->info.imprim) > 4 || mesh->info.ddebug) && ns+nc+nf > 0 ){
3260#ifndef MMG_PATTERN
3261 fprintf(stdout," ");
3262#endif
3263 fprintf(stdout," %8" MMG5_PRId " splitted, %8" MMG5_PRId " collapsed, %8" MMG5_PRId " swapped\n",ns,nc,nf);
3264 }
3265
3266 if ( it > minit-1 && ( !(ns+nc) || (MMG5_abs(nc-ns) < 0.1 * MG_MAX(nc,ns)) ) ) {
3267 ++lastit;
3268 if ( it > minit && lastit>2 ) break;
3269 }
3270 else if ( it+2 >= maxit ) {
3271 /* Last iteration because we reach the maximal number of iter */
3272 ++lastit;
3273 }
3274 else if ( lastit ) {
3275 /* Avoid the incrementation of mesh->info.fem if we have detected a last
3276 iteration but anatet4 leads to have nc, nf or ns != 0 so we perform a last
3277 iter */
3278 ++lastit;
3279 }
3280 }
3281 while ( ++it < maxit && (ns+nc+nf > 0 || lastit<3) );
3282
3283 if ( mesh->info.imprim > 0 ) {
3284 if ( (abs(mesh->info.imprim) < 5 || mesh->info.ddebug ) && nns+nnc > 0 ) {
3285#ifndef MMG_PATTERN
3286 fprintf(stdout," ");
3287#endif
3288 fprintf(stdout, " %8" MMG5_PRId " splitted, %8" MMG5_PRId " collapsed, %8" MMG5_PRId " swapped, %d iter.\n",nns,nnc,nnf,it);
3289 }
3290 }
3291
3292 return 1;
3293}
int ier
MMG5_pMesh * mesh
int MMG3D_delPROctree(MMG5_pMesh mesh, MMG3D_pPROctree q, const int no)
Definition: PRoctree_3d.c:977
int MMG3D_bezierInt(MMG5_pBezier pb, double uv[2], double o[3], double no[3], double to[3])
Definition: bezier_3d.c:608
int MMG5_BezierTgt(double c1[3], double c2[3], double n1[3], double n2[3], double t1[3], double t2[3])
Definition: bezier_3d.c:53
double MMG5_BezierGeod(double c1[3], double c2[3], double t1[3], double t2[3])
Definition: bezier_3d.c:111
MMG5_Info info
int MMG5_boulevolp(MMG5_pMesh mesh, MMG5_int start, int ip, int64_t *list)
Given a vertex and a tetrahedron, find all tetrahedra in the ball of this vertex.
Definition: boulep_3d.c:57
int MMG5_boulesurfvolp(MMG5_pMesh mesh, MMG5_int start, int ip, int iface, int64_t *listv, int *ilistv, MMG5_int *lists, int *ilists, int isnm)
Definition: boulep_3d.c:610
int MMG5_boulesurfvolpNom(MMG5_pMesh mesh, MMG5_int start, int ip, int iface, int64_t *listv, int *ilistv, MMG5_int *lists, int *ilists, MMG5_int *refmin, MMG5_int *refplus, int isnm)
Definition: boulep_3d.c:780
int MMG5_coquilface(MMG5_pMesh mesh, MMG5_int start, int8_t iface, int ia, int64_t *list, MMG5_int *it1, MMG5_int *it2, int silent)
Definition: boulep_3d.c:1846
int MMG5_coquil(MMG5_pMesh mesh, MMG5_int start, int ia, int64_t *list, int8_t *isbdy)
Definition: boulep_3d.c:1406
MMG5_int MMG5_colver(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ilist, int8_t indq, int8_t typchk)
Definition: colver_3d.c:1154
int MMG5_chkcol_int(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t iface, int8_t iedg, int64_t *list, int ilist, int8_t typchk)
Definition: colver_3d.c:43
int MMG5_chkcol_bdy(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t iface, int8_t iedg, int64_t *listv, int ilistv, MMG5_int *lists, int ilists, MMG5_int refmin, MMG5_int refplus, int8_t typchk, int isnm, int8_t isnmint)
Definition: colver_3d.c:527
#define MMG5_EPSD
#define MMG5_EPS
int MMG5_hashUpdate(MMG5_Hash *hash, MMG5_int a, MMG5_int b, MMG5_int k)
Definition: hash.c:404
int MMG5_hashNew(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int hsiz, MMG5_int hmax)
Definition: hash.c:532
MMG5_int MMG5_hashGet(MMG5_Hash *hash, MMG5_int a, MMG5_int b)
Definition: hash.c:501
int MMG5_hashEdge(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int a, MMG5_int b, MMG5_int k)
Definition: hash.c:346
int MMG3D_hashTetra(MMG5_pMesh mesh, int pack)
Create array of adjacency.
Definition: hash_3d.c:122
int MMG5_hashPop(MMG5_Hash *hash, MMG5_int a, MMG5_int b)
Definition: hash_3d.c:850
static double MMG5_orcal(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int iel)
int MMG5_interp4bar33_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int ip, double cb[4])
Definition: intmet_3d.c:522
int MMG3D_intmet33_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t i, MMG5_int ip, double s)
Definition: intmet_3d.c:126
API headers and documentation for the mmg3d library, for volumetric meshes in 3D.
#define MMG3D_LMAX
Definition: libmmg3d.h:130
int MMG5_split4op(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:3804
int MMG5_split3op(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:2654
int MMG5_split1b(MMG5_pMesh, MMG5_pSol, int64_t *, int, MMG5_int, int, int8_t, int8_t)
Definition: split_3d.c:629
int MMG3D_split6_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:4414
MMG5_int MMG5_chkswpgen(MMG5_pMesh, MMG5_pSol, MMG5_int, int, int *, int64_t *, double, int8_t)
Definition: swapgen_3d.c:57
int MMG5_swpbdy(MMG5_pMesh, MMG5_pSol, int64_t *, int, MMG5_int, MMG3D_pPROctree, int8_t)
Definition: swap_3d.c:489
int MMG3D_split3_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1636
int MMG5_BezierReg(MMG5_pMesh mesh, MMG5_int ip0, MMG5_int ip1, double s, double v[3], double *o, double *no)
Definition: tools_3d.c:696
#define MMG3D_LOPTS
MMG5_int MMG5_split4bar(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t)
Definition: split_3d.c:3050
int MMG5_split5(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:4203
int MMG3D_split3cone_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1921
int MMG3D_normalAdjaTri(MMG5_pMesh, MMG5_int, int8_t, int, double n[3])
Definition: split_3d.c:472
int MMG3D_simbulgept(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ilist, MMG5_int)
Definition: split_3d.c:330
void MMG3D_delPt(MMG5_pMesh mesh, MMG5_int ip)
Definition: zaldy_3d.c:80
int MMG5_split1(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t metRidTyp)
Definition: split_3d.c:134
MMG5_int MMG3D_indPt(MMG5_pMesh mesh, MMG5_int kp)
Definition: tools_3d.c:918
MMG5_int MMG3D_newPt(MMG5_pMesh mesh, double c[3], uint16_t tag, MMG5_int src)
Definition: zaldy_3d.c:39
int MMG5_split4sf(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:3440
static const int8_t MMG5_iarfinv[4][6]
num of the j^th edge in the i^th face
int MMG5_chkswpbdy(MMG5_pMesh, MMG5_pSol, int64_t *, int, MMG5_int, MMG5_int, int8_t)
Definition: swap_3d.c:57
int MMG3D_swap23(MMG5_pMesh, MMG5_pSol, MMG5_int, int8_t, int, int, MMG5_int, int)
Definition: swap_3d.c:625
int MMG5_split6(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:4483
int MMG3D_split2sf_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1103
int MMG5_BezierRidge(MMG5_pMesh mesh, MMG5_int ip0, MMG5_int ip1, double s, double *o, double *no1, double *no2, double *to)
Definition: tools_3d.c:145
int MMG3D_split4op_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:3685
static const int8_t MMG5_iarf[4][3]
iarf[i]: edges of face opposite to vertex i
double MMG5_caltet33_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
Definition: quality_3d.c:109
int MMG5_swpgen(MMG5_pMesh, MMG5_pSol, MMG5_int, int, int64_t *, MMG3D_pPROctree, int8_t)
Definition: swapgen_3d.c:271
static const uint8_t MMG5_iare[6][2]
vertices of extremities of the edges of the tetra
int MMG3D_split1_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:94
int MMG5_split2(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:1489
static const uint8_t MMG5_ifar[6][2]
ifar[i][]: faces sharing the ith edge of the tetra
int MMG5_BezierRef(MMG5_pMesh mesh, MMG5_int ip0, MMG5_int ip1, double s, double *o, double *no, double *to)
Definition: tools_3d.c:352
int MMG3D_split4sf_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:3337
int MMG5_split2sf(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:1254
#define MMG3D_POINT_REALLOC(mesh, sol, ip, wantedGap, law, o, tag, src)
int MMG5_split3cone(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:2059
int MMG5_movbdynomintpt_iso(MMG5_pMesh, MMG5_pSol, MMG3D_pPROctree, int64_t *, int, int)
Definition: movpt_3d.c:1490
#define MMG3D_LOPTL
int MMG3D_split5_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:4111
#define MMG3D_LSWAPIMPROVE
int MMG5_norface(MMG5_pMesh mesh, MMG5_int k, int iface, double v[3])
Definition: tools_3d.c:52
#define MMG3D_SWAP06
int MMG5_split3(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:1707
int MMG5_directsurfball(MMG5_pMesh mesh, MMG5_int ip, MMG5_int *list, int ilist, double n[3])
Definition: tools_3d.c:66
int MMG5_BezierNom(MMG5_pMesh mesh, MMG5_int ip0, MMG5_int ip1, double s, double *o, double *no, double *to)
Definition: tools_3d.c:527
int MMG3D_split3op_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:2522
static const uint8_t MMG5_idir[4][3]
idir[i]: vertices of face opposite to vertex i
#define MMG3D_LLONG
#define MMG3D_LSHRT
int MMG3D_split2_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1427
@ MMG5_Tetrahedron
Definition: libmmgtypes.h:233
@ MMG5_Triangle
Definition: libmmgtypes.h:232
static int MMG3D_delPatternPts(MMG5_pMesh mesh, MMG5_Hash hash)
Definition: mmg3d1.c:1319
int MMG3D_normalAndTangent_at_sinRidge(MMG5_pMesh mesh, MMG5_int k, int i, int j, double no1[3], double no2[3], double to[3])
Definition: mmg3d1.c:1684
static MMG5_int MMG5_anatet4(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int *nf, int8_t typchk)
Definition: mmg3d1.c:3029
MMG5_int MMG5_movtet(MMG5_pMesh mesh, MMG5_pSol met, MMG3D_pPROctree PROctree, double clickSurf, double clickVol, int moveVol, int improveSurf, int improveVolSurf, int improveVol, int maxit, MMG5_int testmark)
Definition: mmg3d1.c:736
static int MMG3D_anatet4_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t metRidTyp, int *ifac, int *conf0, MMG5_int *adj, int *conf1)
Definition: mmg3d1.c:2755
int MMG3D_dichoto1b(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ret, MMG5_int ip)
Definition: mmg3d1.c:293
MMG5_int MMG5_swpmsh(MMG5_pMesh mesh, MMG5_pSol met, MMG3D_pPROctree PROctree, int typchk)
Definition: mmg3d1.c:598
static int MMG3D_update_rid_geom(MMG5_pPoint ppt, MMG5_pxPoint pxp, double no[3])
Definition: mmg3d1.c:1631
int MMG5_anatet(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk, int patternMode)
Definition: mmg3d1.c:3159
static MMG5_int MMG5_anatet4rid(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int *nf, int8_t typchk)
Definition: mmg3d1.c:3114
static MMG5_int MMG3D_anatets_ani(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:2190
int8_t MMG5_chkedg(MMG5_pMesh mesh, MMG5_Tria *pt, int8_t ori, double hmax, double hausd, int locPar)
Definition: mmg3d1.c:363
static MMG5_int(* MMG3D_anatets)(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:2731
static int MMG3D_chkbdyface(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_pTetra pt, MMG5_pxTetra pxt, int8_t i, MMG5_pTria ptt, int8_t typchk)
Definition: mmg3d1.c:2082
int MMG3D_splsurfedge(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_pTetra pt, MMG5_pxTetra pxt, int8_t imax, int8_t typchk, int8_t chkRidTet, int *warn)
Definition: mmg3d1.c:1957
int MMG3D_dichoto(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int *vx)
Definition: mmg3d1.c:140
static MMG5_int MMG3D_anatets_iso(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:2290
void MMG3D_set_geom(MMG5_pMesh mesh, MMG5_pPoint ppt, uint16_t tag, MMG5_int nmref, MMG5_int edgref, double no1[3], double no2[3], double to[3])
Definition: mmg3d1.c:59
int8_t ddb
Definition: mmg3d1_delone.c:42
int MMG3D_adpcoledg(MMG5_pMesh mesh, MMG5_pSol met, MMG3D_pPROctree *PROctree, MMG5_int k, int8_t imin, double lmin, MMG5_int *nc)
Definition: mmg3d1.c:1210
void MMG5_tet2tri(MMG5_pMesh mesh, MMG5_int k, int8_t ie, MMG5_Tria *ptt)
Definition: mmg3d1.c:102
static MMG5_int MMG5_anatetv(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:1363
MMG5_int MMG5_swptet(MMG5_pMesh mesh, MMG5_pSol met, double crit, double declic, MMG3D_pPROctree PROctree, int typchk, MMG5_int testmark)
Definition: mmg3d1.c:672
int8_t MMG3D_build_bezierEdge(MMG5_pMesh mesh, MMG5_int k, int8_t imax, int8_t i, int8_t j, MMG5_pxTetra pxt, MMG5_int ip1, MMG5_int ip2, MMG5_pPoint p0, MMG5_pPoint p1, MMG5_int *ref, uint16_t *tag, double o[3], double to[3], double no1[3], double no2[3], int64_t *list, int *ilist)
Definition: mmg3d1.c:1780
static MMG5_int MMG5_coltet(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:892
void MMG3D_find_bdyface_from_edge(MMG5_pMesh mesh, MMG5_pTetra pt, int8_t ied, int8_t *i, int8_t *j, int8_t *i1, int8_t *i2, MMG5_int *ip1, MMG5_int *ip2, MMG5_pPoint *p0, MMG5_pPoint *p1)
Definition: mmg3d1.c:1898
#define MG_Tria
#define MG_REQ
#define MG_GEO
#define MMG5_EPSOK
#define MG_EOK(pt)
#define MG_OPNBDY
#define MMG5_INCREASE_MEM_MESSAGE()
#define MG_PARBDY
#define MG_MIN(a, b)
#define MG_MAX(a, b)
#define MG_EDG(tag)
static const uint8_t MMG5_iprv2[3]
#define MG_OLDPARBDY
#define MG_EDG_OR_NOM(tag)
#define MG_SIN(tag)
#define MMG5_NULKAL
#define MG_Tetra
#define MG_GET(flag, bit)
#define MG_RID(tag)
static const uint8_t MMG5_inxt2[6]
int MMG5_nortri(MMG5_pMesh mesh, MMG5_pTria pt, double *n)
Definition: tools.c:209
#define MG_CRN
#define MG_BDY
#define MG_NOM
#define MMG5_EPSD2
#define MG_REF
#define MG_PARBDYBDY
#define MMG5_DEL_MEM(mesh, ptr)
#define MG_SET(flag, bit)
MMG5_pPoint p[3]
Identic as MMG5_HGeom but use MMG5_hedge to store edges instead of MMG5_hgeom (memory economy).
Definition: libmmgtypes.h:603
MMG5_hedge * item
Definition: libmmgtypes.h:605
int8_t parTyp
Definition: libmmgtypes.h:548
int8_t ddebug
Definition: libmmgtypes.h:539
uint8_t noswap
Definition: libmmgtypes.h:553
double hmin
Definition: libmmgtypes.h:525
uint8_t noinsert
Definition: libmmgtypes.h:553
double hmax
Definition: libmmgtypes.h:525
MMG5_pPar par
Definition: libmmgtypes.h:524
double hausd
Definition: libmmgtypes.h:525
int8_t fem
Definition: libmmgtypes.h:546
MMG mesh structure.
Definition: libmmgtypes.h:613
MMG5_Info info
Definition: libmmgtypes.h:659
MMG5_int ne
Definition: libmmgtypes.h:620
MMG5_pPoint point
Definition: libmmgtypes.h:649
MMG5_int mark
Definition: libmmgtypes.h:626
MMG5_int * adja
Definition: libmmgtypes.h:632
MMG5_pxPoint xpoint
Definition: libmmgtypes.h:650
double gap
Definition: libmmgtypes.h:616
MMG5_int base
Definition: libmmgtypes.h:624
MMG5_pTetra tetra
Definition: libmmgtypes.h:651
MMG5_int np
Definition: libmmgtypes.h:620
MMG5_pxTetra xtetra
Definition: libmmgtypes.h:652
Local parameters for a specific entity and reference.
Definition: libmmgtypes.h:263
double hmin
Definition: libmmgtypes.h:264
double hmax
Definition: libmmgtypes.h:265
double hausd
Definition: libmmgtypes.h:266
MMG5_int ref
Definition: libmmgtypes.h:267
int8_t elt
Definition: libmmgtypes.h:268
Structure to store vertices of an MMG mesh.
Definition: libmmgtypes.h:276
double n[3]
Definition: libmmgtypes.h:278
double c[3]
Definition: libmmgtypes.h:277
uint16_t tag
Definition: libmmgtypes.h:290
MMG5_int xp
Definition: libmmgtypes.h:285
MMG5_int ref
Definition: libmmgtypes.h:284
MMG5_int flag
Definition: libmmgtypes.h:288
double * m
Definition: libmmgtypes.h:680
Structure to store tetrahedra of an MMG mesh.
Definition: libmmgtypes.h:407
MMG5_int v[4]
Definition: libmmgtypes.h:409
MMG5_int xt
Definition: libmmgtypes.h:413
MMG5_int flag
Definition: libmmgtypes.h:416
MMG5_int ref
Definition: libmmgtypes.h:410
MMG5_int mark
Definition: libmmgtypes.h:412
uint16_t tag
Definition: libmmgtypes.h:417
double qual
Definition: libmmgtypes.h:408
Structure to store triangles of a MMG mesh.
Definition: libmmgtypes.h:338
MMG5_int edg[3]
Definition: libmmgtypes.h:345
MMG5_int ref
Definition: libmmgtypes.h:341
MMG5_int flag
Definition: libmmgtypes.h:347
uint16_t tag[3]
Definition: libmmgtypes.h:348
MMG5_int v[3]
Definition: libmmgtypes.h:340
Structure to store surface vertices of an MMG mesh.
Definition: libmmgtypes.h:300
int8_t nnor
Definition: libmmgtypes.h:303
double n2[3]
Definition: libmmgtypes.h:301
double n1[3]
Definition: libmmgtypes.h:301
Structure to store additional information for the surface tetrahedra of an MMG mesh.
Definition: libmmgtypes.h:425
uint16_t tag[6]
Definition: libmmgtypes.h:432
int8_t ori
Definition: libmmgtypes.h:434
MMG5_int edg[6]
Definition: libmmgtypes.h:428
MMG5_int ref[4]
Definition: libmmgtypes.h:426
uint16_t ftag[4]
Definition: libmmgtypes.h:430