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 int16_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
622 /* No swap of geometric edge */
623 if ( MG_EDG_OR_NOM(pxt->tag[ia]) || (pxt->tag[ia] & MG_REQ) )
624 continue;
625
626 ret = MMG5_coquilface(mesh,k,i,ia,list,&it1,&it2,0);
627 ilist = ret / 2;
628 if ( ret < 0 ) return -1;
629 /* CAUTION: trigger collapse with 2 elements */
630 if ( ilist <= 1 ) continue;
631 ier = MMG5_chkswpbdy(mesh,met,list,ilist,it1,it2,typchk);
632 if ( ier < 0 )
633 return -1;
634 else if ( ier ) {
635 ier = MMG5_swpbdy(mesh,met,list,ret,it1,PROctree,typchk);
636 if ( ier > 0 ) ns++;
637 else if ( ier < 0 ) return -1;
638 break;
639 }
640 }
641 if ( ier ) break;
642 }
643 }
644 nns += ns;
645 }
646 while ( ++it < maxit && ns > 0 );
647 if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nns > 0 )
648 fprintf(stdout," %8" MMG5_PRId " edge swapped\n",nns);
649
650 return nns;
651}
652
668MMG5_int MMG5_swptet(MMG5_pMesh mesh,MMG5_pSol met,double crit,double declic,
669 MMG3D_pPROctree PROctree,int typchk,MMG5_int testmark) {
670 MMG5_pTetra pt;
671 MMG5_pxTetra pxt;
672 int ilist,it,maxit,ier;
673 int64_t list[MMG3D_LMAX+2];
674 MMG5_int k,ns,nns,nconf;
675 int8_t i;
676
677 maxit = 2;
678 it = nns = 0;
679
680 do {
681 ns = 0;
682 for (k=1; k<=mesh->ne; k++) {
683 pt = &mesh->tetra[k];
684 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
685 else if ( pt->mark < testmark ) continue;
686
687 if ( pt->qual > declic ) continue;
688
689 for (i=0; i<6; i++) {
690 /* Prevent swap of a ref or tagged edge */
691 if ( pt->xt ) {
692 pxt = &mesh->xtetra[pt->xt];
693 if ( pxt->edg[i] || pxt->tag[i] ) continue;
694 }
695
696 nconf = MMG5_chkswpgen(mesh,met,k,i,&ilist,list,crit,typchk);
697 if ( nconf<0 ) return -1;
698 else if ( nconf ) {
699 ier = MMG5_swpgen(mesh,met,nconf,ilist,list,PROctree,typchk);
700 if ( ier > 0 ) ns++;
701 else if ( ier < 0 ) return -1;
702 break;
703 }
704 }
705 }
706 nns += ns;
707 }
708 while ( ++it < maxit && ns > 0 );
709 if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nns > 0 )
710 fprintf(stdout," %8" MMG5_PRId " edge swapped\n",nns);
711
712 return nns;
713}
714
733 double clickSurf,double clickVol,int moveVol, int improveSurf,
734 int improveVolSurf, int improveVol, int maxit,MMG5_int testmark) {
735 MMG5_pTetra pt;
736 MMG5_pPoint ppt;
737 MMG5_pxTetra pxt;
738 MMG5_Tria tt;
739 double *n,caltri;
740 int ier,ilists,ilistv,it,i;
741 MMG5_int k,lists[MMG3D_LMAX+2],nm,nnm,ns,base;
742 int64_t listv[MMG3D_LMAX+2];
743 uint8_t j,i0;
744
745 if ( abs(mesh->info.imprim) > 5 || mesh->info.ddebug )
746 fprintf(stdout," ** OPTIMIZING MESH\n");
747
748 base = 1;
749 for (k=1; k<=mesh->np; k++)
750 mesh->point[k].flag = base;
751
752 it = nnm = 0;
753 do {
754 base++;
755 nm = ns = 0;
756 for (k=1; k<=mesh->ne; k++) {
757 pt = &mesh->tetra[k];
758 if ( !MG_EOK(pt) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
759 else if ( pt->mark < testmark ) continue;
760
761 /* point j on face i */
762 for (i=0; i<4; i++) {
763 for (j=0; j<3; j++) {
764 if ( pt->xt ) {
765 pxt = &mesh->xtetra[pt->xt];
766 }
767 else pxt = 0;
768
769 i0 = MMG5_idir[i][j];
770 ppt = &mesh->point[pt->v[i0]];
771 if ( ppt->flag == base ) continue;
772 else if ( MG_SIN(ppt->tag) ) continue;
773
774 if( pt->xt && (pxt->ftag[i] & MG_BDY)) {
775 /* skip required faces */
776 if( pxt->ftag[i] & MG_REQ ) continue;
777
778 assert( 0<=i && i<4 && "unexpected local face idx");
779 MMG5_tet2tri(mesh,k,i,&tt);
780 caltri = MMG5_caltri(mesh,met,&tt);
781
782 if ( caltri >= clickSurf) {
783 j = 3;
784 continue;
785 }
786 }
787 if ( maxit != 1 ) {
788 ppt->flag = base;
789 }
790 ier = 0;
791 if ( ppt->tag & MG_BDY ) {
792 /* Catch a boundary point by an external face, unless point is internal non manifold */
793 if ( (!pt->xt) || !(MG_BDY & pxt->ftag[i]) ) continue;
794 else if ( (ppt->tag & MG_PARBDY) || (ppt->tag & MG_PARBDYBDY) ) continue; /* skip parallel points seen by non-required faces */
795 else if ( ppt->tag & MG_NOM ){
796 if ( ppt->xp && mesh->xpoint[ppt->xp].nnor ) {
797 assert( 0<=i0 && i0<4 && "unexpected local index for vertex");
798 ilistv = MMG5_boulevolp(mesh,k,i0,listv);
799 if ( !ilistv ) continue;
800 /* Iso for now */
801 ier = MMG5_movbdynomintpt_iso(mesh,met,PROctree,listv,ilistv,improveVolSurf);
802 }
803 else {
804 if ( mesh->adja[4*(k-1)+1+i] ) continue;
805
806 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,1);
807 if( !ier ) continue;
808 else if ( ier>0 )
809 ier = MMG5_movbdynompt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
810 else
811 return -1;
812 }
813 }
814 else if ( ppt->tag & MG_GEO ) {
815 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
816 if ( !ier ) continue;
817 else if ( ier>0 )
818 ier = MMG5_movbdyridpt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
819 else
820 return -1;
821 }
822 else if ( ppt->tag & MG_REF ) {
823 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
824 if ( !ier )
825 continue;
826 else if ( ier>0 )
827 ier = MMG5_movbdyrefpt(mesh,met,PROctree,listv,ilistv,lists,ilists,improveVolSurf);
828 else
829 return -1;
830 }
831 else {
832 ier=MMG5_boulesurfvolp(mesh,k,i0,i,listv,&ilistv,lists,&ilists,0);
833 if ( !ier )
834 continue;
835 else if ( ier<0 )
836 return -1;
837
838 n = &(mesh->xpoint[ppt->xp].n1[0]);
839
840 /* If the orientation of the tetra face is
841 * compatible with the triangle (MG_GET(ori,i)), we know that we
842 * are well orientated. */
843 if ( !MG_GET(pxt->ori,i) ) {
844 if ( !MMG5_directsurfball(mesh,pt->v[i0],lists,ilists,n) )
845 continue;
846 }
847 ier = MMG5_movbdyregpt(mesh,met,PROctree,listv,ilistv,
848 lists,ilists,improveSurf,improveVolSurf);
849 if (ier < 0 ) return -1;
850 else if ( ier ) ns++;
851 }
852 }
853 else if ( moveVol && (pt->qual < clickVol) ) {
854 assert( 0<=i0 && i0<4 && "unexpected local index for vertex");
855 ilistv = MMG5_boulevolp(mesh,k,i0,listv);
856 if ( !ilistv ) continue;
857 ier = MMG5_movintpt(mesh,met,PROctree,listv,ilistv,improveVol);
858 }
859 if ( ier ) {
860 nm++;
861 if(maxit==1){
862 ppt->flag = base;
863 }
864 }
865 }
866 }
867 }
868 nnm += nm;
869 if ( mesh->info.ddebug ) fprintf(stdout," %8" MMG5_PRId " moved, %" MMG5_PRId " geometry\n",nm,ns);
870 }
871 while( ++it < maxit && nm > 0 );
872
873 if ( (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) && nnm )
874 fprintf(stdout," %8" MMG5_PRId " vertices moved, %d iter.\n",nnm,it);
875
876 return nnm;
877}
878
888static int MMG5_coltet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
889 MMG5_pTetra pt,ptloc;
890 MMG5_pxTetra pxt;
891 MMG5_pPoint p0,p1;
892 MMG5_pPar par;
893 double ll,ux,uy,uz,hmi2;
894 int ilists,ilist;
895 MMG5_int base,k,nc,nnm,lists[MMG3D_LMAX+2],refmin,refplus;
896 int64_t list[MMG3D_LMAX+2];
897 int l,kk,isloc,ifac1;
898 int16_t tag,isnm,isnmint;
899 int8_t i,j,ip,iq;
900 int ier;
901
902 nc = nnm = 0;
903
904 /* init point flags */
905 for (k=1; k<=mesh->np; k++) {
906 p0 = &mesh->point[k];
907 p0->flag = 0;
908 }
909
910 for (k=1; k<=mesh->ne; k++) {
911 /* Remark: we can have int32 overflow on large meshes for base field..*/
912 base = ++mesh->base;
913 pt = &mesh->tetra[k];
914 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
915
916 pxt = pt->xt? &mesh->xtetra[pt->xt] : 0;
917
918 for (i=0; i<4; i++) {
919 ier = 0;
920
921 for (j=0; j<3; j++) {
922 if ( pt->xt && (pxt->tag[MMG5_iarf[i][j]] & MG_REQ) ) continue;
923 ip = MMG5_idir[i][MMG5_inxt2[j]];
924 iq = MMG5_idir[i][MMG5_iprv2[j]];
925
926 p0 = &mesh->point[pt->v[ip]];
927 p1 = &mesh->point[pt->v[iq]];
928
929 if ( p0->flag == base ) {
930 /* I think that we can't pass here because we break the loop when base
931 * is setted and just after we increment it */
932 assert(0);
933 continue;
934 }
935 else {
936 /* Ignore OLDPARBDY tag of p0 */
937 int16_t tag = p0->tag;
938 tag &= ~MG_OLDPARBDY;
939 if ( (tag > p1->tag) || (tag & MG_REQ) ) {
940 /* Unable to merge edge */
941 continue;
942 }
943 }
944
945 /* Ball of point: computed here if needed for the local parameter
946 * evaluation, after length check otherwise (because the ball
947 * computation is time consuming) */
948 ilist = ilists = 0;
949 refmin = refplus = -1;
950 if ( mesh->info.npar ) {
951 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
952 tag = pxt->tag[MMG5_iarf[i][j]];
953 isnm = (tag & MG_NOM);
954 isnmint = ( p0->xp && mesh->xpoint[p0->xp].nnor );
955
956 if ( p0->tag > tag ) continue;
957
958 /* Catch an exterior non manifold point by an external face */
959 if ( isnm ) {
960 if ( isnmint ) {
961 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
962 ilist = MMG5_boulevolp(mesh,k,ip,list);
963 }
964 else {
965 if ( mesh->adja[4*(k-1)+1+i] ) continue;
966 if (MMG5_boulesurfvolpNom(mesh,k,ip,i,
967 list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 )
968 return -1;
969 }
970 }
971 else {
972 if (MMG5_boulesurfvolp(mesh,k,ip,i,
973 list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
974 return -1;
975 }
976 }
977 else {
978 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
979 ilist = MMG5_boulevolp(mesh,k,ip,list);
980 }
981 }
982
983 /* Check length */
984 if ( typchk == 1 ) {
985 ux = p1->c[0] - p0->c[0];
986 uy = p1->c[1] - p0->c[1];
987 uz = p1->c[2] - p0->c[2];
988 ll = ux*ux + uy*uy + uz*uz;
989
990 /* local parameters*/
991 hmi2 = mesh->info.hmin;
992 isloc = 0;
993
994 /* Local parameters at vertex */
995 // take the min of the local hmin at p0 and hmin at p1
996
997 /* Local parameters at tetra */
998 if ( mesh->info.parTyp & MG_Tetra ) {
999 l = 0;
1000 do
1001 {
1002 if ( isloc ) break;
1003
1004 par = &mesh->info.par[l];
1005 if ( par->elt != MMG5_Tetrahedron ) continue;
1006
1007 for ( kk=0; kk<ilist; ++kk ) {
1008 ptloc = &mesh->tetra[list[kk]/4];
1009 if ( par->ref == ptloc->ref ) {
1010 hmi2 = par->hmin;
1011 isloc = 1;
1012 break;
1013 }
1014 }
1015 } while ( ++l<mesh->info.npar );
1016
1017 for ( ; l<mesh->info.npar; ++l ) {
1018 par = &mesh->info.par[l];
1019 if ( par->elt != MMG5_Tetrahedron ) continue;
1020
1021 for ( kk=0; kk<ilist; ++kk ) {
1022 ptloc = &mesh->tetra[list[kk]/4];
1023 if ( par->ref == ptloc->ref ) {
1024 hmi2 = MG_MAX(hmi2,par->hmin);
1025 break;
1026 }
1027 }
1028 }
1029 }
1030
1031 /* Local parameters at triangle */
1032 if ( mesh->info.parTyp & MG_Tria && ( pt->xt && (pxt->ftag[i] & MG_BDY) )) {
1033 l = 0;
1034 do
1035 {
1036 if ( isloc ) break;
1037
1038 par = &mesh->info.par[l];
1039 if ( par->elt != MMG5_Triangle ) continue;
1040
1041 for ( kk=0; kk<ilists; ++kk ) {
1042 ptloc = &mesh->tetra[lists[kk]/4];
1043 ifac1 = lists[kk] % 4;
1044 assert(ptloc->xt && (mesh->xtetra[ptloc->xt].ftag[ifac1] & MG_BDY) );
1045
1046 if ( par->ref == mesh->xtetra[ptloc->xt].ref[ifac1] ) {
1047 hmi2 = par->hmin;
1048 isloc = 1;
1049 break;
1050 }
1051 }
1052 } while ( ++l<mesh->info.npar );
1053
1054 for ( ; l<mesh->info.npar; ++l ) {
1055 par = &mesh->info.par[l];
1056 if ( par->elt != MMG5_Triangle ) continue;
1057
1058 for ( kk=0; kk<ilists; ++kk ) {
1059 ptloc = &mesh->tetra[lists[kk]/4];
1060 ifac1 = lists[kk] % 4;
1061 assert(ptloc->xt && (mesh->xtetra[ptloc->xt].ftag[ifac1] & MG_BDY) );
1062
1063 if ( par->ref == mesh->xtetra[ptloc->xt].ref[ifac1] ) {
1064 hmi2 = MG_MAX(hmi2,par->hmin);
1065 break;
1066 }
1067 }
1068 }
1069 }
1070
1071 hmi2 = hmi2*hmi2;
1072 if ( ll > hmi2*MMG3D_LSHRT ) continue;
1073 }
1074 else if ( typchk == 2 ) {
1075 ll = MMG5_lenedg(mesh,met,MMG5_iarf[i][j],pt);
1076 // Case of an internal tetra with 4 ridges vertices.
1077 if ( ll == 0 ) continue;
1078 if ( ll > MMG3D_LSHRT ) continue;
1079 }
1080
1081
1082 /* Ball computation if not computed before */
1083 if ( !mesh->info.npar ) {
1084 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
1085 tag = pxt->tag[MMG5_iarf[i][j]];
1086 isnm = (tag & MG_NOM);
1087 isnmint = ( p0->xp && mesh->xpoint[p0->xp].nnor );
1088
1089 if ( p0->tag > tag ) continue;
1090 if ( isnm ) {
1091 if ( isnmint ) {
1092 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
1093 ilist = MMG5_boulevolp(mesh,k,ip,list);
1094 }
1095 else {
1096 if ( mesh->adja[4*(k-1)+1+i] ) continue;
1097 if (MMG5_boulesurfvolpNom(mesh,k,ip,i,
1098 list,&ilist,lists,&ilists,&refmin,&refplus,p0->tag & MG_NOM) < 0 )
1099 return -1;
1100 }
1101 }
1102 else {
1103 if (MMG5_boulesurfvolp(mesh,k,ip,i,
1104 list,&ilist,lists,&ilists,p0->tag & MG_NOM) < 0 )
1105 return -1;
1106 }
1107 }
1108 else {
1109 assert( 0<=ip && ip<4 && "unexpected local index for vertex");
1110 ilist = MMG5_boulevolp(mesh,k,ip,list);
1111 }
1112 }
1113
1114 /* boundary face: collapse ip on iq */
1115 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
1116 tag = pxt->tag[MMG5_iarf[i][j]];
1117 tag |= MG_BDY;
1118 if ( p0->tag > tag ) continue;
1119
1120 isnm = ( tag & MG_NOM );
1121 isnmint = ( tag & MG_NOM ) ? ( p0->xp && mesh->xpoint[p0->xp].nnor ) : 0;
1122 if ( isnm && (!isnmint) ) {
1123 /* Treat surfacic non manifold points from a surface triangle */
1124 if ( mesh->adja[4*(k-1)+1+i] ) continue;
1125 }
1126 ilist = MMG5_chkcol_bdy(mesh,met,k,i,j,list,ilist,lists,ilists,refmin,refplus,typchk,isnm,isnmint);
1127 }
1128 /* internal face */
1129 else {
1130 isnm = 0;
1131 if ( p0->tag & MG_BDY ) continue;
1132 ilist = MMG5_chkcol_int(mesh,met,k,i,j,list,ilist,typchk);
1133 }
1134
1135 if ( ilist > 0 ) {
1136 ier = MMG5_colver(mesh,met,list,ilist,iq,typchk);
1137 if ( ier < 0 ) return -1;
1138 else if ( ier ) {
1140 break;
1141 }
1142 }
1143 else if (ilist < 0 ) return -1;
1144 }
1145 if ( ier ) {
1146 p1->flag = base;
1147 if ( isnm ) nnm++;
1148 nc++;
1149 break;
1150 }
1151 }
1152 }
1153 if ( nc > 0 && (abs(mesh->info.imprim) > 5 || mesh->info.ddebug) )
1154 fprintf(stdout," %8" MMG5_PRId " vertices removed, %8" MMG5_PRId " non manifold,\n",nc,nnm);
1155
1156 return nc;
1157}
1158
1181 MMG3D_pPROctree *PROctree,MMG5_int k,
1182 int8_t imin,double lmin,MMG5_int* nc) {
1183 MMG5_pTetra pt;
1184 MMG5_pxTetra pxt;
1185 MMG5_pPoint p0,p1;
1186 int64_t list[MMG3D_LMAX+2];
1187 MMG5_int lists[MMG3D_LMAX+2],ip1,ip2;
1188 int ilist,ilists;
1189 int8_t j,i,i1,i2;
1190
1191 if(lmin > MMG3D_LOPTS) {
1192 /* Edge is large enough: nothing to do */
1193 return 3;
1194 }
1195
1196 if ( lmin == 0 ) {
1197 /* Case of an internal tetra with 4 ridges vertices */
1198//#warning is it possible to merge this edge ??
1199 return 0;
1200 }
1201
1202 pt = &mesh->tetra[k];
1203 pxt = pt->xt ? &mesh->xtetra[pt->xt] : 0;
1204
1205 MMG3D_find_bdyface_from_edge(mesh,pt,imin,&i,&j,&i1,&i2,&ip1,&ip2,&p0,&p1);
1206
1207 /* Ignore OLDPARBDY tag of p0 */
1208 int16_t tag0 = p0->tag;
1209 tag0 &= ~MG_OLDPARBDY;
1210 if ( (tag0 > p1->tag) || (tag0 & MG_REQ) ) {
1211 /* Unable to merge edge: pass to next element */
1212 return 0;
1213 }
1214
1216 ilist = 0;
1217 if ( pt->xt && (pxt->ftag[i] & MG_BDY) ) {
1218 /* Case of a boundary face */
1219 int16_t tag = pxt->tag[MMG5_iarf[i][j]];
1220 if ( tag & MG_REQ ) {
1221 return 0;
1222 }
1223 tag |= MG_BDY;
1224 tag &= ~MG_OLDPARBDY;
1225 if ( tag0 > tag ) {
1226 return 0;
1227 }
1228 if ( ( tag & MG_NOM ) && (mesh->adja[4*(k-1)+1+i]) ) {
1229 return 0;
1230 }
1231
1232 int16_t isnm = (p0->tag & MG_NOM);
1233 if (MMG5_boulesurfvolp(mesh,k,i1,i, list,&ilist,lists,&ilists,isnm) < 0 ) {
1234 return -1;
1235 }
1236
1237 ilist = MMG5_chkcol_bdy(mesh,met,k,i,j,list,ilist,lists,ilists,0,0,2,0,0);
1238 }
1239 else {
1240 /* Case of an internal face */
1241 if ( p0->tag & MG_BDY ) {
1242 return 0;
1243 }
1244
1245 ilist = MMG5_boulevolp(mesh,k,i1,list);
1246 ilist = MMG5_chkcol_int(mesh,met,k,i,j,list,ilist,2);
1247 }
1248
1250 if ( ilist > 0 ) {
1251 /* Checks are ok */
1252 int ier = MMG5_colver(mesh,met,list,ilist,i2,2);
1253 if ( ilist < 0 ) {
1254 /* Colver failure */
1255 return 0;
1256 }
1257 if ( ier < 0 ) {
1258 /* Colver failure */
1259 return -1;
1260 }
1261 else if(ier) {
1262 /* Collapse is successful */
1263 if ( PROctree && (*PROctree) ) {
1264 MMG3D_delPROctree(mesh,*PROctree,ier);
1265 }
1267 (*nc)++;
1268 return 2;
1269 }
1270 }
1271 else if (ilist < 0 ) {
1272 /* Checks on chkcol_bdy have failed (error in edge shell computation) */
1273 return -1;
1274 }
1275
1276 return 3;
1277}
1278
1279
1288static inline
1290{
1291 MMG5_pTetra pt;
1292 int ia,i,j;
1293 MMG5_int vx[6],k;
1294
1295 /* Delete the useless added points */
1296 for (k=1; k<=mesh->ne; k++) {
1297 pt = &mesh->tetra[k];
1298 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
1299
1300 memset(vx,0,6*sizeof(MMG5_int));
1301 for (ia=0,i=0; i<3; i++) {
1302 for (j=i+1; j<4; j++,ia++) {
1303 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_REQ) ) continue;
1304 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
1305 if ( vx[ia] > 0 ) {
1306 MMG3D_delPt(mesh,vx[ia]);
1307 if ( !MMG5_hashUpdate(&hash,pt->v[i],pt->v[j],-1) ) {
1308 fprintf(stderr,"\n ## Error: %s: unable to delete point idx"
1309 " along edge %" MMG5_PRId " %" MMG5_PRId ".\n", __func__,
1310 MMG3D_indPt(mesh,pt->v[i]),
1311 MMG3D_indPt(mesh,pt->v[j]));
1312 MMG5_DEL_MEM(mesh,hash.item);
1313 return 0;
1314 }
1315 }
1316 }
1317 }
1318 }
1319 return 1;
1320}
1321
1332static MMG5_int
1334 MMG5_pTetra pt;
1335 MMG5_pPoint p1,p2;
1336 MMG5_xTetra *pxt;
1337 MMG5_Hash hash;
1338 MMG5_pPar par;
1339 double ll,o[3],ux,uy,uz,hma2,mincal;
1340 int l,memlack,ier;
1341 MMG5_int src,vx[6],ip,ip1,ip2,k,ne,ns,nap;
1342 int8_t i,j,ia;
1343
1345 if ( !MMG5_hashNew(mesh,&hash,mesh->np,7*mesh->np) ) return -1;
1346 memlack = ns = nap = 0;
1347
1348 /* Hash all boundary and required edges, and put ip = -1 in hash structure */
1349 for (k=1; k<=mesh->ne; k++) {
1350 pt = &mesh->tetra[k];
1351 if ( !MG_EOK(pt) ) continue;
1352
1353 /* avoid split of edges belonging to a required tet */
1354 if ( pt->tag & MG_REQ ) {
1355 for (i=0; i<6; i++) {
1356 ip1 = pt->v[MMG5_iare[i][0]];
1357 ip2 = pt->v[MMG5_iare[i][1]];
1358 ip = -1;
1359 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
1360 }
1361 continue;
1362 }
1363
1364 if ( !pt->xt ) continue;
1365
1366 pxt = &mesh->xtetra[pt->xt];
1367 for (i=0; i<4; i++) {
1368 if ( pxt->ftag[i] & MG_BDY ) {
1369 for (j=0; j<3; j++) {
1370 ip1 = pt->v[MMG5_idir[i][MMG5_inxt2[j]]];
1371 ip2 = pt->v[MMG5_idir[i][MMG5_iprv2[j]]];
1372 ip = -1;
1373 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
1374 }
1375 }
1376 }
1377 }
1378
1380 mincal = FLT_MAX;
1381 for (k=1; k<=mesh->ne; k++) {
1382 pt = &mesh->tetra[k];
1383 if ( !MG_EOK(pt) ) continue;
1384 pt->flag = 0;
1385 for (i=0; i<6; i++) {
1386 ip = -1;
1387 ip1 = pt->v[MMG5_iare[i][0]];
1388 ip2 = pt->v[MMG5_iare[i][1]];
1389 p1 = &mesh->point[ip1];
1390 p2 = &mesh->point[ip2];
1391 if ( pt->xt ) {
1392 pxt = &mesh->xtetra[pt->xt];
1393 if ( pxt->tag[i] & MG_REQ ) continue;
1394 }
1395 else pxt = 0;
1396
1397 if ( (p1->tag & MG_BDY) && (p2->tag & MG_BDY) ) {
1398 /* Split internal edges connecting bdy points */
1399 ip = MMG5_hashGet(&hash,ip1,ip2);
1400 }
1401 else {
1402 if (typchk == 1) {
1403 ux = p2->c[0] - p1->c[0];
1404 uy = p2->c[1] - p1->c[1];
1405 uz = p2->c[2] - p1->c[2];
1406 ll = ux*ux + uy*uy + uz*uz;
1407
1408 hma2 = mesh->info.hmax;
1409 /* Local parameters */
1410 /* Local param at vertices ip1 and ip2 */
1411 // take the min of the local hmin at ip1 and hmin at ip2
1412
1413 /* Local parameters at tetra */
1414 if ( mesh->info.parTyp & MG_Tetra ) {
1415 for ( l=0; l<mesh->info.npar; ++l ) {
1416 par = &mesh->info.par[l];
1417
1418 if ( par->elt != MMG5_Tetrahedron ) continue;
1419 if ( par->ref != pt->ref ) continue;
1420
1421 hma2 = par->hmax;
1422 break;
1423 }
1424 }
1425 hma2 = MMG3D_LLONG*MMG3D_LLONG*hma2*hma2;
1426 if ( ll > hma2 ) {
1427 ip = MMG5_hashGet(&hash,ip1,ip2);
1428 mincal = MG_MIN(mincal,pt->qual);
1429 }
1430 }
1431 else if ( typchk == 2 ) {
1432 ll = MMG5_lenedg(mesh,met,i,pt);
1433 // Case of an internal tetra with 4 ridges vertices.
1434 if ( ll == 0 ) continue;
1435 if ( ll > MMG3D_LLONG ) {
1436 ip = MMG5_hashGet(&hash,ip1,ip2);
1437 mincal = MG_MIN(mincal,pt->qual);
1438 }
1439 }
1440 }
1441 if ( ip < 0 ) continue;
1442 else if ( !ip ) {
1443 /* new midpoint */
1444 o[0] = 0.5 * (p1->c[0]+p2->c[0]);
1445 o[1] = 0.5 * (p1->c[1]+p2->c[1]);
1446 o[2] = 0.5 * (p1->c[2]+p2->c[2]);
1447
1448#ifdef USE_POINTMAP
1449 src = mesh->point[ip1].src;
1450#else
1451 src = 1;
1452#endif
1453 ip = MMG3D_newPt(mesh,o,0,src);
1454 if ( !ip ) {
1455 /* reallocation of point table */
1456
1458 fprintf(stderr,"\n ## Warning: %s: unable to"
1459 " allocate a new point\n",__func__);
1461 memlack=1;
1462 goto split
1463 ,o,0,src);
1464 }
1465
1466 assert ( met );
1467 if ( met->m ) {
1468 if ( typchk == 1 && (met->size>1) )
1469 ier = MMG3D_intmet33_ani(mesh,met,k,i,ip,0.5);
1470 else
1471 ier = MMG5_intmet(mesh,met,k,i,ip,0.5);
1472
1473 if (!ier) {
1474 // Unable to compute the metric
1475 return -1;
1476 }
1477 else if ( ier < 0 ) {
1478 MMG3D_delPt(mesh,ip);
1479 continue;
1480 }
1481 }
1482
1483 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
1484 MG_SET(pt->flag,i);
1485 nap++;
1486 }
1487 }
1488 }
1489 if ( !nap ) {
1490 MMG5_DEL_MEM(mesh,hash.item);
1491 return 0;
1492 }
1493
1495split:
1496 if ( mincal < MMG5_EPS ) {
1497 /* Delete the useless added points */
1498 if ( !MMG3D_delPatternPts(mesh,hash) ) return -1;
1499
1500 /* Avoid the creation of bad quality elements */
1501 if ( mesh->info.imprim > 5 || mesh->info.ddebug ) {
1502 fprintf(stderr,"\n ## Warning: %s: too bad quality for the worst element."
1503 " Volumic patterns skipped.\n",__func__);
1504 }
1505
1506 MMG5_DEL_MEM(mesh,hash.item);
1507 return 0;
1508 }
1509
1510 ns = 0;
1511 ne = mesh->ne;
1512 for (k=1; k<=ne; k++) {
1513 pt = &mesh->tetra[k];
1514 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
1515 memset(vx,0,6*sizeof(MMG5_int));
1516 pt->flag = 0;
1517 for (ia=0,i=0; i<3; i++) {
1518 for (j=i+1; j<4; j++,ia++) {
1519 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_REQ) ) continue;
1520 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
1521 if ( vx[ia] > 0 ) MG_SET(pt->flag,ia);
1522 }
1523 }
1524
1525 switch (pt->flag) {
1526 case 1: case 2: case 4: case 8: case 16: case 32: /* 1 edge split */
1527 if ( ! MMG5_split1(mesh,met,k,vx,typchk-1) ) return -1;
1528 ns++;
1529 break;
1530 case 48: case 24: case 40: case 6: case 34: case 36:
1531 case 20: case 5: case 17: case 9: case 3: case 10: /* 2 edges (same face) split */
1532 if ( ! MMG5_split2sf(mesh,met,k,vx,typchk-1) ) return -1;
1533 ns++;
1534 break;
1535
1536 case 33: case 18: case 12: /* 2 opposite edges split */
1537 if ( ! MMG5_split2(mesh,met,k,vx,typchk-1) ) return -1;
1538 ns++;
1539 break;
1540
1541 case 11: case 21: case 38: case 56: /* 3 edges on the same faces splitted */
1542 if ( ! MMG5_split3(mesh,met,k,vx,typchk-1) ) return -1;
1543 ns++;
1544 break;
1545
1546 case 7: case 25: case 42: case 52: /* 3 edges on conic configuration splitted */
1547 if ( ! MMG5_split3cone(mesh,met,k,vx,typchk-1) ) return -1;
1548 ns++;
1549 break;
1550
1551 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
1552 case 14: case 49: case 50: case 44: case 41: /* 3 edges on opposite configuration splitted */
1553 if ( ! MMG5_split3op(mesh,met,k,vx,typchk-1) ) return -1;
1554 ns++;
1555 break;
1556
1557 case 23: case 29: case 53: case 60: case 57: case 58:
1558 case 27: case 15: case 43: case 39: case 54: case 46: /* 4 edges with 3 lying on the same face splitted */
1559 if ( ! MMG5_split4sf(mesh,met,k,vx,typchk-1) ) return -1;
1560 ns++;
1561 break;
1562
1563 /* 4 edges with no 3 lying on the same face splitted */
1564 case 30: case 45: case 51:
1565 if ( ! MMG5_split4op(mesh,met,k,vx,typchk-1) ) return -1;
1566 ns++;
1567 break;
1568
1569 case 62: case 61: case 59: case 55: case 47: case 31: /* 5 edges split */
1570 if ( ! MMG5_split5(mesh,met,k,vx,typchk-1) ) return -1;
1571 ns++;
1572 break;
1573
1574 case 63: /* 6 edges split */
1575 if ( ! MMG5_split6(mesh,met,k,vx,typchk-1) ) return -1;
1576 ns++;
1577 break;
1578 }
1579 }
1580
1581 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
1582 fprintf(stdout," %7" MMG5_PRId " splitted\n",nap);
1583
1584 MMG5_DEL_MEM(mesh,hash.item);
1585 if ( memlack ) return -1;
1586 return nap;
1587}
1588
1600static inline int
1602 double dd,to[3];
1603
1604 dd = no[0]*pxp->n1[0]+no[1]*pxp->n1[1]+no[2]*pxp->n1[2];
1605 if ( dd > 1.0-MMG5_EPS ) return 0;
1606
1607 memcpy(pxp->n2,no,3*sizeof(double));
1608
1609 /* a computation of the tangent with respect to these two normals is possible */
1610 to[0] = pxp->n1[1]*pxp->n2[2] - pxp->n1[2]*pxp->n2[1];
1611 to[1] = pxp->n1[2]*pxp->n2[0] - pxp->n1[0]*pxp->n2[2];
1612 to[2] = pxp->n1[0]*pxp->n2[1] - pxp->n1[1]*pxp->n2[0];
1613 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
1614 if ( dd > MMG5_EPSD2 ) {
1615 dd = 1.0 / sqrt(dd);
1616 to[0] *= dd;
1617 to[1] *= dd;
1618 to[2] *= dd;
1619 memcpy(ppt->n,to,3*sizeof(double));
1620 }
1621 else {
1622 /* Detect opposite normals... */
1623 if ( to[0]*to[0]+to[1]*to[1]+to[2]*to[2] > MMG5_EPSD2 ) {
1624 /* Non opposite normals: fail in debug mode */
1625 assert ( dd > MMG5_EPSD2 );
1626 }
1627 else {
1628 /* Opposite normals: unable to compute the tangent, check if we have
1629 * already stored a tangent, fail otherwise */
1630 assert ( ppt->n[0]*ppt->n[0]+ppt->n[1]*ppt->n[1]+ppt->n[2]*ppt->n[2] > MMG5_EPSD2 );
1631 }
1632 }
1633
1634 return 1;
1635}
1636
1655 double no1[3],double no2[3],double to[3] ) {
1656 MMG5_Tria ptt;
1657 double dd;
1658 int ier;
1659 static int8_t warn_n = 0;
1660
1661 assert ( 0<=i && i<4 && "unexpected local idx for face" );
1662 assert ( 0<=j && j<3 && "unexpected local edg odx in face" );
1663
1664 assert( 0<=i && i<4 && "unexpected local face idx");
1665
1666#ifndef NDEBUG
1667 /* Check that face is boundary and has a suitable orientation */
1668 assert ( mesh->tetra[k].xt && "Tetra is not boundary" );
1669
1670 MMG5_pxTetra pxt = &mesh->xtetra[ mesh->tetra[k].xt ];
1671 assert ( (pxt->ftag[i] & MG_BDY) && "Face is not boundary" );
1672 assert ( MG_GET(pxt->ori,i) && "Wrong face orientation" );
1673#endif
1674
1675 MMG5_tet2tri(mesh,k,i,&ptt);
1676 MMG5_nortri(mesh,&ptt,no1);
1677
1678 /* In this case, 'to' orientation depends on the edge processing (so on
1679 * the triangle from which we come) so we can't use it to compute no2. */
1680 ier = MMG3D_normalAdjaTri(mesh,k,i,j,no2);
1681 if ( ier < 0 ) {
1682 return 0;
1683 }
1684 else if ( !ier ) {
1685 if ( !warn_n ) {
1686 warn_n = 1;
1687 fprintf(stderr," ## Warning: %s: %d: error in the computation of normal"
1688 " at triangle.\n",__func__,__LINE__);
1689 }
1690 no2[0] = to[1]*no1[2] - to[2]*no1[1];
1691 no2[1] = to[2]*no1[0] - to[0]*no1[2];
1692 no2[2] = to[0]*no1[1] - to[1]*no1[0];
1693
1694 dd = no2[0]*no2[0] + no2[1]*no2[1] + no2[2]*no2[2];
1695 if ( dd > MMG5_EPSD2 ) {
1696 dd = 1.0 / sqrt(dd);
1697 no2[0] *= dd;
1698 no2[1] *= dd;
1699 no2[2] *= dd;
1700 }
1701 }
1702 else {
1703 assert ( ier==1 );
1704
1705 /* Compute 'to' as intersection of no1 and no2 */
1706 to[0] = no1[1]*no2[2] - no1[2]*no2[1];
1707 to[1] = no1[2]*no2[0] - no1[0]*no2[2];
1708 to[2] = no1[0]*no2[1] - no1[1]*no2[0];
1709 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
1710 if ( dd > MMG5_EPSD2 ) {
1711 dd = 1.0 / sqrt(dd);
1712 to[0] *= dd;
1713 to[1] *= dd;
1714 to[2] *= dd;
1715 }
1716 }
1717 return 1;
1718}
1719
1751 int8_t imax,int8_t i, int8_t j,
1752 MMG5_pxTetra pxt,
1753 MMG5_int ip1,MMG5_int ip2,
1754 MMG5_pPoint p0, MMG5_pPoint p1,
1755 MMG5_int *ref,int16_t *tag,
1756 double o[3],double to[3],double no1[3],
1757 double no2[3],int64_t *list,int *ilist) {
1758 MMG5_Tria ptt;
1759 double v[3];
1760
1761 if ( (p0->tag & MG_PARBDY) && (p1->tag & MG_PARBDY) ) {
1762 /* Skip edge with extremities on parallel interfaces */
1763 return 0;
1764 }
1765
1766 int8_t ori = MG_GET(pxt->ori,i);
1767 if ( !ori ) {
1768 /* Treat triangles at interface of 2 subdomains from well oriented face */
1769 return 0;
1770 }
1771
1772 *ref = pxt->edg[imax];
1773 *tag = pxt->tag[imax];
1774 if ( (*tag) & MG_REQ ) {
1775 /* No need to split required edges */
1776 return 0;
1777 }
1778
1779 (*tag) |= MG_BDY;
1780 int8_t dummy;
1781 *ilist = MMG5_coquil(mesh,k,imax,list,&dummy);
1782 if ( !(*ilist) ) {
1783 /* On of the tetra of the edge shell is required: we cannot split the edge */
1784 return 0;
1785 }
1786 else if ( (*ilist) < 0 ) {
1787 /* Shell computation has failed */
1788 return -1;
1789 }
1790
1792 if ( (*tag) & MG_NOM ){
1793 /* Edge is non-manifold */
1794 if( !MMG5_BezierNom(mesh,ip1,ip2,0.5,o,no1,to) ) {
1795 /* Unable to compute geom infos at nm edge */
1796 return 0;
1797 }
1798 else if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1799 assert( 0<=i && i<4 && "unexpected local face idx");
1800 MMG5_tet2tri(mesh,k,i,&ptt);
1801 MMG5_nortri(mesh,&ptt,no1);
1802 }
1803 }
1804 else if ( (*tag) & MG_GEO ) {
1805 /* Edge is ridge */
1806 if ( !MMG5_BezierRidge(mesh,ip1,ip2,0.5,o,no1,no2,to) ) {
1807 /* Unable to compute geom infos at ridge */
1808//#warning why a continue here?
1809 return 0;
1810 }
1811 else if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1812 if ( !MMG3D_normalAndTangent_at_sinRidge(mesh,k,i,j,no1,no2,to) ) {
1813 return -1;
1814 }
1815 }
1816 }
1817 else if ( (*tag) & MG_REF ) {
1818 /* Edge is ref */
1819 if ( !MMG5_BezierRef(mesh,ip1,ip2,0.5,o,no1,to) ) {
1820 /* Unable to compute geom infos at ref edge */
1821 return 1;
1822 }
1823 if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1824//#warning creation of sin-sin ref edge to see if it works without normal realloc
1825 assert( 0<=i && i<4 && "unexpected local face idx");
1826 MMG5_tet2tri(mesh,k,i,&ptt);
1827 MMG5_nortri(mesh,&ptt,no1);
1828 }
1829 }
1830 else {
1831 /* Longest edge is regular */
1832 if ( !MMG5_norface(mesh,k,i,v) ) {
1833 /* Unable to treat long edge: try to collapse short one */
1834 return 1;
1835 }
1836 if ( !MMG5_BezierReg(mesh,ip1,ip2,0.5,v,o,no1) ) {
1837 /* Unable to compute geom infos at regular edge */
1838 return 1;
1839 }
1840 else if ( MG_SIN(p0->tag) && MG_SIN(p1->tag) ) {
1841 assert( 0<=i && i<4 && "unexpected local face idx");
1842 MMG5_tet2tri(mesh,k,i,&ptt);
1843 MMG5_nortri(mesh,&ptt,no1);
1844 }
1845 }
1846 return 2;
1847}
1848
1869 int8_t *i,int8_t *j,int8_t*i1,int8_t*i2,
1870 MMG5_int*ip1,MMG5_int*ip2,MMG5_pPoint *p0,MMG5_pPoint *p1) {
1871
1872 int8_t ifa0 = MMG5_ifar[ied][0];
1873 int8_t ifa1 = MMG5_ifar[ied][1];
1874
1878 /* Default face */
1879 MMG5_pxTetra pxt = pt->xt? &mesh->xtetra[pt->xt] : 0;
1880
1881 (*i) = ifa0;
1882 if ( pt->xt ) {
1883 int16_t is_ifa0_bdy = (pxt->ftag[ifa0] & MG_BDY);
1884 int16_t is_ifa1_bdy = (pxt->ftag[ifa1] & MG_BDY);
1885
1886 if ( is_ifa0_bdy && is_ifa1_bdy ) {
1887 /* Two bdy faces: search if one has a suitable orientation */
1888 int8_t ifa1_ori = MG_GET(pxt->ori,(*i));
1889 (*i) = ifa1_ori ? ifa1 : ifa0;
1890 }
1891 else if ( is_ifa1_bdy ) {
1892 /* only ifa1 is boundary: no need to check for orientation (if it has a
1893 * bad ori we will quit the function later) */
1894 (*i) = ifa1;
1895 }
1896 /* For all other cases (only ifa0 is bdy or no bdy face), we use default
1897 * face (ifa0) */
1898 }
1899
1900 (*j) = MMG5_iarfinv[*i][ied];
1901 (*i1) = MMG5_idir[*i][MMG5_inxt2[*j]];
1902 (*i2) = MMG5_idir[*i][MMG5_iprv2[*j]];
1903 (*ip1) = pt->v[*i1];
1904 (*ip2) = pt->v[*i2];
1905 (*p0) = &mesh->point[*ip1];
1906 (*p1) = &mesh->point[*ip2];
1907
1908}
1909
1928 MMG5_pTetra pt,MMG5_pxTetra pxt,int8_t imax,int8_t typchk,
1929 int8_t chkRidTet,int *warn ) {
1930 MMG5_pPoint p0,p1,ppt;
1931 double o[3],to[3],no1[3],no2[3];
1932 int ilist;
1933 MMG5_int src,ip,ip1,ip2,ref;
1934 int64_t list[MMG3D_LMAX+2];
1935 int ier;
1936 int16_t tag;
1937 int8_t j,i,i1,i2;
1938
1939 assert ( pxt == &mesh->xtetra[pt->xt] );
1940
1941 /* proceed edges according to lengths */
1942 MMG3D_find_bdyface_from_edge(mesh,pt,imax,&i,&j,&i1,&i2,&ip1,&ip2,&p0,&p1);
1943
1944 ier = MMG3D_build_bezierEdge(mesh,k,imax,i,j,pxt,ip1,ip2,p0,p1,&ref,&tag,
1945 o,to,no1,no2,list,&ilist);
1946 switch ( ier ) {
1947 case -1:
1948 /* Strong failure (due to lack of memory or wrong edge shell) */
1949 return -1;
1950 case 0:
1951 /* We don't want to split the edge (edge shell contains a required tet or
1952 * non manifold edge or ridge reconstruction has failed) */
1953 return 0;
1954 case 1:
1955 /* We don't want to split the edge (ref or regular edge reconstruction has
1956 * failed) */
1957 return 0;
1958 default:
1959 /* We can insert new vertex along edge */
1960 assert ( ier==2 );
1961 }
1962
1963#ifdef USE_POINTMAP
1964 src = mesh->point[ip1].src;
1965#else
1966 src = 1;
1967#endif
1968 ip = MMG3D_newPt(mesh,o,tag,src);
1969 if ( !ip ) {
1970 /* reallocation of point table */
1972 *warn=1;
1973 return 2;
1974 ,o,tag,src);
1975 }
1976 if ( met && met->m ) {
1977 if ( typchk == 1 && (met->size>1) ) {
1978 ier = MMG3D_intmet33_ani(mesh,met,k,imax,ip,0.5);
1979 }
1980 else {
1981 ier = MMG5_intmet(mesh,met,k,imax,ip,0.5);
1982 }
1983 if ( !ier ) {
1984 MMG3D_delPt(mesh,ip);
1985 return -1;
1986 }
1987 else if (ier < 0) {
1988 MMG3D_delPt(mesh,ip);
1989 return 0;
1990 }
1991 }
1992
1993 /* simbulgept needs a valid tangent at ridge point (to build ridge metric in
1994 * order to comute edge lengths). Thus we need to store the geometric info of
1995 * point here. */
1996 ppt = &mesh->point[ip];
1997 MMG3D_set_geom(mesh,ppt,tag,ref,pxt->ref[i],no1,no2,to);
1998
1999 ier = MMG3D_simbulgept(mesh,met,list,ilist,ip);
2000
2001#ifndef NDEBUG
2002 if ( mesh->info.ddebug ) {
2003 assert ( (ier != -1) && "simbulgept failure" );
2004 }
2005#endif
2006
2007 if ( ier < 0 || ier == 2 ) {
2008 MMG3D_delPt(mesh,ip);
2009 return 0;
2010 }
2011 else if ( ier == 0 ) {
2012 ier = MMG3D_dichoto1b(mesh,met,list,ilist,ip);
2013 }
2014 if ( ier == 1 ) { ier = MMG5_split1b(mesh,met,list,ilist,ip,1,typchk-1,chkRidTet); }
2015
2016 /* if we realloc memory in MMG5_split1b pt and pxt pointers are not valid */
2017 if ( ier < 0 ) {
2018 fprintf(stderr,"\n ## Error: %s: unable to split.\n",__func__);
2019 return -1;
2020 }
2021 else if ( ier == 0 || ier == 2 ) {
2022 MMG3D_delPt(mesh,ip);
2023 return 0;
2024 }
2025
2026 return 1;
2027}
2028
2044static
2046 MMG5_pxTetra pxt,int8_t i,MMG5_pTria ptt,int8_t typchk ) {
2047
2048 MMG5_pPar par;
2049 double len,hmax,hausd;
2050 int l;
2051 MMG5_int ip1,ip2;
2052 int8_t isloc,ier;
2053 int8_t j,i1,i2,ia;
2054
2055 if ( typchk == 1 ) {
2056
2057 /* Local parameters for ptt and k */
2058 hmax = mesh->info.hmax;
2059 hausd = mesh->info.hausd;
2060 isloc = 0;
2061
2062 if ( mesh->info.parTyp & MG_Tetra ) {
2063 for ( l=0; l<mesh->info.npar; ++l ) {
2064 par = &mesh->info.par[l];
2065
2066 if ( par->elt != MMG5_Tetrahedron ) continue;
2067 if ( par->ref != pt->ref ) continue;
2068
2069 hmax = par->hmax;
2070 hausd = par->hausd;
2071 isloc = 1;
2072 break;
2073 }
2074 }
2075 if ( mesh->info.parTyp & MG_Tria ) {
2076 if ( isloc ) {
2077 for ( l=0; l<mesh->info.npar; ++l ) {
2078 par = &mesh->info.par[l];
2079
2080 if ( par->elt != MMG5_Triangle ) continue;
2081 if ( par->ref != ptt->ref ) continue;
2082
2083 hmax = MG_MIN(hmax,par->hmax);
2084 hausd = MG_MIN(hausd,par->hausd);
2085 break;
2086 }
2087 }
2088 else {
2089 for ( l=0; l<mesh->info.npar; ++l ) {
2090 par = &mesh->info.par[l];
2091
2092 if ( par->elt != MMG5_Triangle ) continue;
2093 if ( par->ref != ptt->ref ) continue;
2094
2095 hmax = par->hmax;
2096 hausd = par->hausd;
2097 isloc = 1;
2098 break;
2099 }
2100 }
2101 }
2102
2103 ier = MMG5_chkedg(mesh,ptt,MG_GET(pxt->ori,i),hmax,hausd,isloc);
2104 if ( ier < 0 ) {
2105 /* Error */
2106 return 0;
2107 } else if ( ier == 0 ) {
2108 /* Nothing to split */
2109 return 1;
2110 }
2111
2112 /* put back flag on tetra */
2113 for (j=0; j<3; j++){
2114 if ( pxt->tag[MMG5_iarf[i][j]] & MG_REQ ) continue;
2115 if ( MG_GET(ptt->flag,j) ) MG_SET(pt->flag,MMG5_iarf[i][j]);
2116 }
2117 }
2118 else if ( typchk == 2 ) {
2119 for (j=0; j<3; j++) {
2120 ia = MMG5_iarf[i][j];
2121 if ( pxt->tag[ia] & MG_REQ ) continue;
2122 i1 = MMG5_iare[ia][0];
2123 i2 = MMG5_iare[ia][1];
2124 ip1 = pt->v[i1];
2125 ip2 = pt->v[i2];
2126 len = MMG5_lenedg(mesh,met,ia,pt);
2127
2128 assert( isfinite(len) && (len!=-len) );
2129
2130 // Case of an internal tetra with 4 ridges vertices.
2131 if ( len == 0 ) continue;
2132 if ( len > MMG3D_LLONG ) MG_SET(pt->flag,ia);
2133 /* Treat here the ridges coming from a corner (we can not do that after
2134 * because the corner don't have xpoints) */
2135 if ( (mesh->point[ip1].tag & MG_CRN) || (mesh->point[ip2].tag & MG_CRN) ) {
2136 if ( len > MMG3D_LOPTL ) MG_SET(pt->flag,ia);
2137 }
2138 }
2139 }
2140 return 1;
2141}
2142
2153static MMG5_int MMG3D_anatets_ani(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk) {
2154 MMG5_pTetra pt;
2155 MMG5_pxTetra pxt;
2156 MMG5_Tria ptt;
2157 double len,lmax;
2158 double ux,uy,uz;
2159 MMG5_int ns,k,ip1,ip2;
2160 int ier,warn;
2161 int8_t imax,j,i,i1,i2;
2162
2163 assert ( met->m && met->size==6 );
2164
2165 if ( !MMG3D_hashTetra(mesh,1) ) {
2166 fprintf(stderr,"\n ## Hashing problem. Exit program.\n");
2167 return -1;
2168 }
2169
2170 ns = 0;
2171 warn = 0;
2172
2173 for (k=1; k<=mesh->ne; k++) {
2174 pt = &mesh->tetra[k];
2175 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) || !pt->xt ) continue;
2176
2177 /* check boundary face cut w/r Hausdorff or hmax */
2178 pt->flag = 0;
2179 pxt = &mesh->xtetra[pt->xt];
2180
2181 /* Travel well oriented boundary faces to flag edges that have to be cut */
2182 for (i=0; i<4; i++){
2183 if ( pxt->ftag[i] & MG_REQ ) continue;
2184 if ( !(pxt->ftag[i] & MG_BDY) ) continue;
2185
2186 if ( !MG_GET(pxt->ori,i) ) continue;
2187
2188 /* virtual triangle */
2189 assert( 0<=i && i<4 && "unexpected local face idx");
2190 MMG5_tet2tri(mesh,k,i,&ptt);
2191
2192 if ( !MMG3D_chkbdyface(mesh,met,k,pt,pxt,i,&ptt,typchk) ) { continue; }
2193 }
2194
2196 imax = 6;
2197 lmax = 0.;
2198 for ( j=0; j<6; ++j ) {
2199 if ( !MG_GET(pt->flag,j) ) { continue; }
2200
2201 i1 = MMG5_iare[j][0];
2202 i2 = MMG5_iare[j][1];
2203 ip1 = pt->v[i1];
2204 ip2 = pt->v[i2];
2205
2206 ux = mesh->point[ip1].c[0] - mesh->point[ip2].c[0];
2207 uy = mesh->point[ip1].c[1] - mesh->point[ip2].c[1];
2208 uz = mesh->point[ip1].c[2] - mesh->point[ip2].c[2];
2209
2210 len = ux*ux + uy*uy + uz*uz;
2211 if ( len <= lmax) { continue; }
2212 lmax = len;
2213 imax = j;
2214 }
2215
2216 pt->flag = 0;
2217 if ( imax < 6 ) { MG_SET(pt->flag,imax); }
2218
2219 if ( !pt->flag ) continue;
2220
2221 ier = MMG3D_splsurfedge( mesh,met,k,pt,pxt,imax,typchk,1,&warn );
2222
2223 if ( ier==-1 ) { return -1; }
2224 else if ( !ier ) { continue; }
2225 else if ( ier==2 ) {
2226 /* Unable to split due to lack of memory */
2227 return ns;
2228 }
2229
2230 ++ns;
2231 }
2232
2233 return ns;
2234}
2235
2252static MMG5_int
2254 MMG5_pTetra pt;
2255 MMG5_pPoint ppt;
2256 MMG5_Tria ptt,ptt2;
2257 MMG5_xTetra *pxt;
2258 MMG5_xPoint *pxp;
2259 MMG5_Bezier pb,pb2;
2260 MMG5_Hash hash;
2261 double o[3],no[3],to[3],dd;
2262 int ic,it,ier;
2263 MMG5_int ip,vx[6],src,nc,ns,ni,ne,k,ip1,ip2,nap,ixp1,ixp2;
2264 int8_t i,j,j2,ia,i1,i2,ifac,intnom;
2265 static double uv[3][2] = { {0.5,0.5}, {0.,0.5}, {0.5,0.} };
2266 static int8_t mmgWarn = 0, mmgWarn2 = 0;
2267
2269 if ( !MMG5_hashNew(mesh,&hash,mesh->np,7*mesh->np) ) return -1;
2270 ns = nap = 0;
2271
2272 for (k=1; k<=mesh->ne; k++) {
2273 pt = &mesh->tetra[k];
2274 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) || !pt->xt ) continue;
2275
2276 /* check boundary face cut w/r Hausdorff or hmax */
2277 pt->flag = 0;
2278 pxt = &mesh->xtetra[pt->xt];
2279
2280 for (i=0; i<4; i++){
2281 if ( pxt->ftag[i] & MG_REQ ) continue;
2282 if ( !(pxt->ftag[i] & MG_BDY) ) continue;
2283
2284 if ( typchk == 1 ) {
2285 if ( !MG_GET(pxt->ori,i) ) { continue; }
2286 }
2287
2288 /* virtual triangle */
2289 assert( 0<=i && i<4 && "unexpected local face idx");
2290 MMG5_tet2tri(mesh,k,i,&ptt);
2291
2292 if ( !MMG3D_chkbdyface(mesh,met,k,pt,pxt,i,&ptt,typchk) ) { continue; }
2293
2294 if ( !pt->flag ) continue;
2295 ns++;
2296
2297 /* geometric support */
2298 ier = MMG5_bezierCP(mesh,&ptt,&pb,MG_GET(pxt->ori,i));
2299 assert(ier);
2300
2301 /* scan edges in face to split */
2302 for (j=0; j<3; j++) {
2303 ia = MMG5_iarf[i][j];
2304 if ( !MG_GET(pt->flag,ia) ) continue;
2305 if ( pxt->tag[ia] & MG_REQ ) continue;
2306 i1 = MMG5_iare[ia][0];
2307 i2 = MMG5_iare[ia][1];
2308 ip1 = pt->v[i1];
2309 ip2 = pt->v[i2];
2310 ip = MMG5_hashGet(&hash,ip1,ip2);
2311
2312 /* new point along edge */
2313 if ( !ip ) {
2314
2315 if ( ptt.tag[j] & MG_NOM ) {
2316 /* Use BezierNom to ensure that the new nm point has a normal that
2317 * is consistent with the normals at point ip1 and ip2 */
2318 if ( !MMG5_BezierNom(mesh,ip1,ip2,0.5,o,no,to) ) {
2319 continue;
2320 }
2321 }
2322 else {
2323 /* Otherwise we can build the bezier patch */
2324 ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
2325 assert(ier);
2326 }
2327
2328#ifdef USE_POINTMAP
2329 src = mesh->point[ip1].src;
2330#else
2331 src = 1;
2332#endif
2333 ip = MMG3D_newPt(mesh,o,MG_BDY,src);
2334 if ( !ip ) {
2335 /* reallocation of point table */
2337 fprintf(stderr,"\n ## Error: %s: unable to"
2338 " allocate a new point.\n",__func__);
2340 MMG3D_delPatternPts(mesh,hash);return -1;
2341 ,o,MG_BDY,src);
2342 // Now pb->p contain a wrong memory address.
2343 pb.p[0] = &mesh->point[ptt.v[0]];
2344 pb.p[1] = &mesh->point[ptt.v[1]];
2345 pb.p[2] = &mesh->point[ptt.v[2]];
2346 }
2347 if ( !MMG5_hashEdge(mesh,&hash,ip1,ip2,ip) ) return -1;
2348 ppt = &mesh->point[ip];
2349
2350 assert ( met );
2351 if ( met->m ) {
2352 if ( typchk == 1 && (met->size>1) )
2353 ier = MMG3D_intmet33_ani(mesh,met,k,ia,ip,0.5);
2354 else
2355 ier = MMG5_intmet(mesh,met,k,ia,ip,0.5);
2356
2357 if ( !ier ) {
2358 if ( !mmgWarn ) {
2359 fprintf(stderr,"\n ## Error: %s: unable to interpolate at least"
2360 " 1 metric.\n",__func__);
2361 mmgWarn = 1;
2362 }
2363 return -1;
2364 }
2365 else if ( ier < 0 ) {
2366 MMG3D_delPt(mesh,ip);
2367 continue;
2368 }
2369 }
2370
2371 if ( MG_EDG_OR_NOM(ptt.tag[j]) )
2372 ppt->ref = ptt.edg[j] ? ptt.edg[j] : ptt.ref;
2373 else
2374 ppt->ref = ptt.ref;
2375 ppt->tag |= ptt.tag[j];
2376 pxp = &mesh->xpoint[ppt->xp];
2377
2378 /* Update normal and tangent vectors */
2379 intnom = 0;
2380 if ( ptt.tag[j] & MG_NOM ) {
2381 ixp1 = mesh->point[ip1].xp;
2382 ixp2 = mesh->point[ip2].xp;
2383 intnom = ( mesh->xpoint[ixp1].nnor ) || ( mesh->xpoint[ixp2].nnor );
2384 }
2385 if ( intnom ) pxp->nnor = 1;
2386 else memcpy(pxp->n1,no,3*sizeof(double));
2387
2388 memcpy(ppt->n,to,3*sizeof(double));
2389
2390 if ( mesh->info.fem<typchk ) {
2391 /* A ridge can be at the interface of 2 boundary faces of the same
2392 * tetra: second normal has to be computed */
2393 if ( MG_EDG(ptt.tag[j]) && !(ptt.tag[j] & MG_NOM) ) {
2394 /* Update the second normal and the tangent at point ip if the edge
2395 * is shared by 2 faces (if anatet4 is not called, 1 tetra may have
2396 * 2 bdry faces) */
2397 ifac = (MMG5_ifar[ia][0] == i)? MMG5_ifar[ia][1] : MMG5_ifar[ia][0];
2398 if ( pxt->ftag[ifac] & MG_BDY ) {
2399 j2 = MMG5_iarfinv[ifac][ia];
2400
2401 /* Compute tangent and normal with respect to the face ifac */
2402 /* virtual triangle */
2403 assert( 0<=ifac && ifac<4 && "unexpected local face idx");
2404 MMG5_tet2tri(mesh,k,ifac,&ptt2);
2405
2406 /* geometric support */
2407 ier = MMG5_bezierCP(mesh,&ptt2,&pb2,MG_GET(pxt->ori,ifac));
2408 assert(ier);
2409
2410 ier = MMG3D_bezierInt(&pb2,&uv[j2][0],o,no,to);
2411 assert(ier);
2412
2413 if ( !MMG3D_update_rid_geom(ppt,pxp,no) ) continue;
2414 }
2415 }
2416 }
2417 nap++;
2418 }
2419 else if ( MG_EDG(ptt.tag[j]) && !(ptt.tag[j] & MG_NOM) ) {
2420 /* Point at the interface of 2 boundary faces belonging to different
2421 * tetra : Point has alredy been created from another tetra so we have
2422 * to store the tangent and the second normal at edge */
2423 ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
2424 assert(ier);
2425
2426 ppt = &mesh->point[ip];
2427 assert(ppt->xp);
2428 pxp = &mesh->xpoint[ppt->xp];
2429
2430 if ( !MMG3D_update_rid_geom(ppt,pxp,no) ) continue;
2431 }
2432 }
2433 }
2434 }
2435
2436 if ( !ns ) {
2437 MMG5_DEL_MEM(mesh,hash.item);
2438 return ns;
2439 }
2440
2443 nc = 0;
2444 for (k=1; k<=mesh->ne; k++) {
2445 pt = &mesh->tetra[k];
2446 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) ) continue;
2447 pxt = pt->xt ? &mesh->xtetra[pt->xt] : 0;
2448
2449 /* update face-edge flag */
2450 for (i=0; i<4; i++) {
2451 /* virtual triangle */
2452 memset(&ptt,0,sizeof(MMG5_Tria));
2453 if ( pt->xt && pxt->ftag[i] && pxt->ftag[i] != MG_OLDPARBDY ) {
2454 assert( 0<=i && i<4 && "unexpected local face idx");
2455 MMG5_tet2tri(mesh,k,i,&ptt);
2456 }
2457
2458 for (j=0; j<3; j++) {
2459 ia = MMG5_iarf[i][j];
2460
2461 /* If edge is already marked, nothing to do */
2462 if ( MG_GET(pt->flag,ia) ) {
2463 continue;
2464 }
2465
2467 /* First: skip edge if required */
2468 if ( pt->xt && (pxt->tag[ia] & MG_REQ) ) continue;
2469
2470 /* Second: if possible treat manifold ridges from a boundary face (to
2471 * ensure the computation of n2) */
2472 if ( pt->xt ) {
2473 if ( pxt->tag[ia] & MG_GEO && (!(pxt->tag[ia] & MG_NOM) ) ) {
2474 int ifac2 = (MMG5_ifar[ia][0]!=i)? MMG5_ifar[ia][1] : MMG5_ifar[ia][0];
2475 if ( (!(pxt->ftag[i] & MG_BDY) ) && (pxt->ftag[ifac2] & MG_BDY) ) {
2476 assert ( ifac2 > i && "lower face is already processed" );
2477 continue;
2478 }
2479 }
2480 }
2481
2482 ip1 = pt->v[MMG5_iare[ia][0]];
2483 ip2 = pt->v[MMG5_iare[ia][1]];
2484 ip = MMG5_hashGet(&hash,ip1,ip2);
2485 if ( ip > 0 ) {
2486
2487 MG_SET(pt->flag,ia);
2488 nc++;
2489
2490 if ( (!(ptt.tag[j] & MG_GEO)) || (ptt.tag[j] & MG_NOM) ) continue;
2491
2492 /* From a boundary face of a boundary tetra we can update normal/tangent
2493 at manifold ridges; */
2494 ppt = &mesh->point[ip];
2495 assert(ppt->xp);
2496 pxp = &mesh->xpoint[ppt->xp];
2497 if ( pt->xt ) ier = MMG5_bezierCP(mesh,&ptt,&pb,MG_GET(pxt->ori,i));
2498 else ier = MMG5_bezierCP(mesh,&ptt,&pb,1);
2499 assert(ier);
2500
2501 ier = MMG3D_bezierInt(&pb,&uv[j][0],o,no,to);
2502 assert(ier);
2503
2504 dd = no[0]*pxp->n1[0]+no[1]*pxp->n1[1]+no[2]*pxp->n1[2];
2505 if ( dd > 1.0-MMG5_EPS ) continue;
2506
2507 memcpy(pxp->n2,no,3*sizeof(double));
2508 }
2509 }
2510 }
2511 }
2512 if ( mesh->info.ddebug && nc ) {
2513 fprintf(stdout," %" MMG5_PRId " added\n",nc);
2514 fflush(stdout);
2515 }
2516
2518 for (k=1; k<=mesh->np; k++)
2519 mesh->point[k].flag = 0;
2520
2521 it = 1;
2522 nc = 0;
2523 do {
2524 ni = 0;
2525 for (k=1; k<=mesh->ne; k++) {
2526 pt = &mesh->tetra[k];
2527 if ( !MG_EOK(pt) || (pt->tag & MG_REQ) || !pt->flag ) continue;
2528 memset(vx,0,6*sizeof(MMG5_int));
2529 pt->flag = ic = 0;
2530 for (ia=0,i=0; i<3; i++) {
2531 for (j=i+1; j<4; j++,ia++) {
2532 if ( pt->xt && (mesh->xtetra[pt->xt].tag[ia] & MG_REQ) ) continue;
2533 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
2534 if ( vx[ia] > 0 ) {
2535 MG_SET(pt->flag,ia);
2536 if ( mesh->point[vx[ia]].flag > 2 ) ic = 1;
2537 }
2538 }
2539 }
2540 if ( !pt->flag ) continue;
2541
2542 switch (pt->flag) {
2543 case 1: case 2: case 4: case 8: case 16: case 32:
2544 ier = MMG3D_split1_sim(mesh,met,k,vx);
2545 break;
2546 case 48: case 24: case 40: case 6: case 34: case 36:
2547 case 20: case 5: case 17: case 9: case 3: case 10:
2548 ier =MMG3D_split2sf_sim(mesh,met,k,vx);
2549 break;
2550 case 33: case 18: case 12:
2551 ier = MMG3D_split2_sim(mesh,met,k,vx);
2552 break;
2553 case 11: case 21: case 38: case 56:
2554 ier = MMG3D_split3_sim(mesh,met,k,vx);
2555 break;
2556 case 7: case 25: case 42: case 52:
2557 ier =MMG3D_split3cone_sim(mesh,met,k,vx);
2558 break;
2559 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
2560 case 14: case 49: case 50: case 44: case 41:
2561 ier = MMG3D_split3op_sim(mesh,met,k,vx);
2562 break;
2563 case 23: case 29: case 53: case 60: case 57: case 58:
2564 case 27: case 15: case 43: case 39: case 54: case 46:
2565 ier = MMG3D_split4sf_sim(mesh,met,k,vx);
2566 break;
2567 case 30: case 45: case 51:
2568 ier = MMG3D_split4op_sim(mesh,met,k,vx);
2569 break;
2570 case 62: case 61: case 59: case 55: case 47: case 31:
2571 ier = MMG3D_split5_sim(mesh,met,k,vx);
2572 break;
2573 case 63:
2574 ier = MMG3D_split6_sim(mesh,met,k,vx);
2575 break;
2576 }
2577 if ( ier ) continue;
2578
2579 ni++;
2580 if ( ic == 0 && MMG3D_dichoto(mesh,met,k,vx) ) {
2581 for (ia=0; ia<6; ia++)
2582 if ( vx[ia] > 0 ) mesh->point[vx[ia]].flag++;
2583 }
2584 else {
2585 for (ia=0; ia<6; ++ia ) {
2586 if ( vx[ia] > 0 ) {
2587 MMG5_hashPop(&hash,pt->v[MMG5_iare[ia][0]],pt->v[MMG5_iare[ia][1]]);
2588 MMG3D_delPt(mesh,vx[ia]);
2589
2590 if ( (mesh->info.ddebug || mesh->info.imprim > 5) && !mmgWarn2 ) {
2591 fprintf(stderr,"\n ## Warning: %s: surfacic pattern: unable to find"
2592 " a valid split for at least 1 point. Point(s) deletion.\n",
2593 __func__ );
2594 mmgWarn2 = 1;
2595 }
2596
2597 }
2598 }
2599 }
2600 }
2601 nc += ni;
2602 }
2603 while( ni > 0 && ++it < 40 );
2604
2605 if ( mesh->info.ddebug && nc ) {
2606 fprintf(stdout," %" MMG5_PRId " corrected, %" MMG5_PRId " invalid\n",nc,ni);
2607 fflush(stdout);
2608 }
2609
2610
2612 ns = 0;
2613 ne = mesh->ne;
2614 for (k=1; k<=ne; k++) {
2615 pt = &mesh->tetra[k];
2616 if ( !MG_EOK(pt) || !pt->flag || (pt->tag & MG_REQ) ) continue;
2617 memset(vx,0,6*sizeof(MMG5_int));
2618 for (ia=0,i=0; i<3; i++) {
2619 for (j=i+1; j<4; j++,ia++) {
2620 if ( MG_GET(pt->flag,ia) ) {
2621 vx[ia] = MMG5_hashGet(&hash,pt->v[i],pt->v[j]);
2622 assert(vx[ia]);
2623 }
2624 }
2625 }
2626 switch (pt->flag) {
2627 case 1: case 2: case 4: case 8: case 16: case 32: /* 1 edge split */
2628 if ( ! MMG5_split1(mesh,met,k,vx,typchk-1) ) return -1;
2629 ns++;
2630 break;
2631 case 48: case 24: case 40: case 6: case 34: case 36:
2632 case 20: case 5: case 17: case 9: case 3: case 10: /* 2 edges (same face) split */
2633 if ( ! MMG5_split2sf(mesh,met,k,vx,typchk-1) ) return -1;
2634 ns++;
2635 break;
2636
2637 case 33: case 18: case 12: /* 2 opposite edges split */
2638 if ( ! MMG5_split2(mesh,met,k,vx,typchk-1) ) return -1;
2639 ns++;
2640 break;
2641
2642 case 11: case 21: case 38: case 56: /* 3 edges on the same faces splitted */
2643 if ( ! MMG5_split3(mesh,met,k,vx,typchk-1) ) return -1;
2644 ns++;
2645 break;
2646
2647 case 7: case 25: case 42: case 52: /* 3 edges on conic configuration splitted */
2648 if ( ! MMG5_split3cone(mesh,met,k,vx,typchk-1) ) return -1;
2649 ns++;
2650 break;
2651
2652 case 35: case 19: case 13: case 37: case 22: case 28: case 26:
2653 case 14: case 49: case 50: case 44: case 41: /* 3 edges on opposite configuration splitted */
2654 if ( ! MMG5_split3op(mesh,met,k,vx,typchk-1) ) return -1;
2655 ns++;
2656 break;
2657
2658 case 23: case 29: case 53: case 60: case 57: case 58:
2659 case 27: case 15: case 43: case 39: case 54: case 46: /* 4 edges with 3 lying on the same face splitted */
2660 if ( ! MMG5_split4sf(mesh,met,k,vx,typchk-1) ) return -1;
2661 ns++;
2662 break;
2663
2664 /* 4 edges with no 3 lying on the same face splitted */
2665 case 30: case 45: case 51:
2666 if ( ! MMG5_split4op(mesh,met,k,vx,typchk-1) ) return -1;
2667 ns++;
2668 break;
2669
2670 case 62: case 61: case 59: case 55: case 47: case 31: /* 5 edges split */
2671 if ( ! MMG5_split5(mesh,met,k,vx,typchk-1) ) return -1;
2672 ns++;
2673 break;
2674
2675 case 63: /* 6 edges split */
2676 if ( ! MMG5_split6(mesh,met,k,vx,typchk-1) ) return -1;
2677 ns++;
2678 break;
2679 }
2680 }
2681 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
2682 fprintf(stdout," %7" MMG5_PRId " elements splitted\n",nap);
2683
2684 MMG5_DEL_MEM(mesh,hash.item);
2685 return nap;
2686}
2687
2688static MMG5_int (*MMG3D_anatets)(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk);
2689
2712static int MMG3D_anatet4_sim(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int k,int8_t metRidTyp,
2713 int *ifac,int* conf0,MMG5_int *adj,int *conf1) {
2714 MMG5_pTetra pt,pt1,ptnew;
2715 MMG5_pxTetra pxt0,pxt1;
2716 MMG5_pPoint ppt,ppt0;
2717 double calold0,calold,calnew,calnew0,calnew1,calnew2,calnew3;
2718 double worst_split4bar_cal,worst_swap_cal,cb[4];
2719 int loc_conf0,loc_conf1;
2720 MMG5_int *adja,k1,np;
2721 int nbdy,i,j0,j1;
2722 uint8_t tau0[4],tau1[4];
2723
2724 pt = &mesh->tetra[k];
2725 calold0 = pt->qual;
2726
2727 assert ( pt->xt );
2728
2729 pxt0 = &mesh->xtetra[pt->xt];
2730
2731 ptnew = &mesh->tetra[0];
2732
2734 if ( !mesh->info.noinsert ) {
2735 ppt0 = &mesh->point[0];
2736 memset(ppt0,0,sizeof(MMG5_Point));
2737
2738 for (i=0; i<4; i++) {
2739 ppt = &mesh->point[pt->v[i]];
2740 ppt0->c[0] += ppt->c[0];
2741 ppt0->c[1] += ppt->c[1];
2742 ppt0->c[2] += ppt->c[2];
2743 }
2744 ppt0->c[0] *= 0.25;
2745 ppt0->c[1] *= 0.25;
2746 ppt0->c[2] *= 0.25;
2747
2748 cb[0] = 0.25; cb[1] = 0.25; cb[2] = 0.25; cb[3] = 0.25;
2749
2750 if ( met->m ) {
2751 if ( !metRidTyp && met->size > 1 )
2752 MMG5_interp4bar33_ani(mesh,met,k,0,cb);
2753 else
2754 MMG5_interp4bar(mesh,met,k,0,cb);
2755 }
2756
2757 memcpy(ptnew,pt,sizeof(MMG5_Tetra));
2758 ptnew->v[0] = 0;
2759 if ( (!metRidTyp) && met->m && met->size>1 )
2760 calnew0 = MMG5_caltet33_ani(mesh,met,ptnew);
2761 else
2762 calnew0 = MMG5_orcal(mesh,met,0);
2763
2764 ptnew->v[0] = pt->v[0];
2765 ptnew->v[1] = 0;
2766 if ( (!metRidTyp) && met->m && met->size>1 )
2767 calnew1 = MMG5_caltet33_ani(mesh,met,ptnew);
2768 else
2769 calnew1 = MMG5_orcal(mesh,met,0);
2770
2771 ptnew->v[1] = pt->v[1];
2772 ptnew->v[2] = 0;
2773 if ( (!metRidTyp) && met->m && met->size>1 )
2774 calnew2 = MMG5_caltet33_ani(mesh,met,ptnew);
2775 else
2776 calnew2 = MMG5_orcal(mesh,met,0);
2777
2778 ptnew->v[2] = pt->v[2];
2779 ptnew->v[3] = 0;
2780 if ( (!metRidTyp) && met->m && met->size>1 )
2781 calnew3 = MMG5_caltet33_ani(mesh,met,ptnew);
2782 else
2783 calnew3 = MMG5_orcal(mesh,met,0);
2784
2785 worst_split4bar_cal = MG_MIN(MG_MIN(calnew0,calnew1),MG_MIN(calnew2,calnew3));
2786 }
2787 else
2788 worst_split4bar_cal = 0.;
2789
2791 *ifac = -1;
2792 *adj = -1;
2793 if ( !mesh->info.noswap ) {
2794 worst_swap_cal = 0.;
2795
2796 for (j0=0; j0<4; j0++) {
2797 if ( pxt0->ftag[j0] & MG_BDY ) continue;
2798
2800 adja = &mesh->adja[4*(k-1)+1];
2801 k1 = adja[j0]/4;
2802 j1 = adja[j0]%4;
2803
2804 assert(k1);
2805
2806 /* Search in which configurations are the tetrahedra (default is case 0-0)
2807 *
2808 * 3 2------------- 0
2809 * ,/|`\ |`\ /|
2810 * ,/ | `\ | `\ /.|
2811 * ,/ '. `\ '. `\ / |
2812 * ,/ | `\ | `\ / .|
2813 * ,/ | `\ | /\.|
2814 * 0-----------'.--------2 ' / 3
2815 * `\. | ,/ | / ,/
2816 * `\. | ,/ | /,/
2817 * `\. '. ,/ '. ,/
2818 * `\. |/ |/
2819 * `1 1
2820 */
2821
2822 /* k may be in configuration 0, 3, 6 or 9. Default is case 0 */
2823 loc_conf0 = 3*j0;
2824
2825 switch(loc_conf0) {
2826 case 3:
2827 tau0[0] = 1; tau0[1] = 0; tau0[2] = 3; tau0[3] = 2;
2828 break;
2829 case 6:
2830 tau0[0] = 2; tau0[1] = 0; tau0[2] = 1; tau0[3] = 3;
2831 break;
2832 case 9:
2833 tau0[0] = 3; tau0[1] = 0; tau0[2] = 2; tau0[3] = 1;
2834 break;
2835 default:
2836 tau0[0] = 0; tau0[1] = 1; tau0[2] = 2; tau0[3] = 3;
2837 }
2838
2839 /* k1 may be in configuration j1, j1+1, j1+2 */
2840 pt1 = &mesh->tetra[k1];
2841
2842 if ( pt1->tag & MG_REQ ) continue;
2843
2844 assert(pt->ref == pt1->ref);
2845 for ( i=0; i<3; ++i )
2846 if ( pt->v[MMG5_idir[j0][0]] == pt1->v[MMG5_idir[j1][i]] ) break;
2847
2848 assert(i<3);
2849 loc_conf1 = 3*j1+i;
2850
2851 switch(loc_conf1) {
2852 case 1:
2853 tau1[0] = 0; tau1[1] = 2; tau1[2] = 3; tau1[3] = 1;
2854 break;
2855 case 2:
2856 tau1[0] = 0; tau1[1] = 3; tau1[2] = 1; tau1[3] = 2;
2857 break;
2858 case 3:
2859 tau1[0] = 1; tau1[1] = 0; tau1[2] = 3; tau1[3] = 2;
2860 break;
2861 case 4:
2862 tau1[0] = 1; tau1[1] = 3; tau1[2] = 2; tau1[3] = 0;
2863 break;
2864 case 5:
2865 tau1[0] = 1; tau1[1] = 2; tau1[2] = 0; tau1[3] = 3;
2866 break;
2867 case 6:
2868 tau1[0] = 2; tau1[1] = 0; tau1[2] = 1; tau1[3] = 3;
2869 break;
2870 case 7:
2871 tau1[0] = 2; tau1[1] = 1; tau1[2] = 3; tau1[3] = 0;
2872 break;
2873 case 8:
2874 tau1[0] = 2; tau1[1] = 3; tau1[2] = 0; tau1[3] = 1;
2875 break;
2876 case 9:
2877 tau1[0] = 3; tau1[1] = 0; tau1[2] = 2; tau1[3] = 1;
2878 break;
2879 case 10:
2880 tau1[0] = 3; tau1[1] = 2; tau1[2] = 1; tau1[3] = 0;
2881 break;
2882 case 11:
2883 tau1[0] = 3; tau1[1] = 1; tau1[2] = 0; tau1[3] = 2;
2884 break;
2885 default:
2886 tau1[0] = 0; tau1[1] = 1; tau1[2] = 2; tau1[3] = 3;
2887 }
2888
2890 if ( pt1->xt ) {
2891 pxt1 = &mesh->xtetra[pt1->xt];
2892
2893 nbdy = 0;
2894 if ( pxt1->ftag[tau1[1]] ) ++nbdy;
2895 if ( pxt0->ftag[tau0[1]] ) ++nbdy;
2896 if ( nbdy > 1 ) continue;
2897
2898 nbdy = 0;
2899 if ( pxt1->ftag[tau1[3]] ) ++nbdy;
2900 if ( pxt0->ftag[tau0[2]] ) ++nbdy;
2901 if ( nbdy > 1 ) continue;
2902
2903 nbdy = 0;
2904 if ( pxt1->ftag[tau1[2]] ) ++nbdy;
2905 if ( pxt0->ftag[tau0[3]] ) ++nbdy;
2906 if ( nbdy > 1 ) continue;
2907
2908 }
2909
2911 calold = MG_MIN(calold0,pt1->qual);
2912
2913 ptnew = &mesh->tetra[0];
2914 memcpy(ptnew,pt,sizeof(MMG5_Tetra));
2915 np = pt1->v[tau1[0]];
2916
2917 ptnew->v[tau0[1]] = np;
2918 if ( (!metRidTyp) && met->m && met->size>1 )
2919 calnew0 = MMG5_caltet33_ani(mesh,met,ptnew);
2920 else
2921 calnew0 = MMG5_orcal(mesh,met,0);
2922 if ( calnew0 < MMG5_NULKAL ) continue;
2923
2924 ptnew->v[tau0[1]] = pt->v[tau0[1]];
2925 ptnew->v[tau0[2]] = np;
2926 if ( (!metRidTyp) && met->m && met->size>1 )
2927 calnew1 = MMG5_caltet33_ani(mesh,met,ptnew);
2928 else
2929 calnew1 = MMG5_orcal(mesh,met,0);
2930 if ( calnew1 < MMG5_NULKAL ) continue;
2931
2932 ptnew->v[tau0[2]] = pt->v[tau0[2]];
2933 ptnew->v[tau0[3]] = np;
2934 if ( (!metRidTyp) && met->m && met->size>1 )
2935 calnew2 = MMG5_caltet33_ani(mesh,met,ptnew);
2936 else
2937 calnew2 = MMG5_orcal(mesh,met,0);
2938 if ( calnew2 < MMG5_NULKAL ) continue;
2939
2940 calnew = MG_MIN(calnew0,MG_MIN(calnew1,calnew2));
2941
2942 if ( calold < MMG5_EPSOK ) {
2943 if ( calnew < calold ) continue;
2944 }
2945 else if ( calnew <= MMG5_EPSOK ) continue;
2946
2947 if ( calnew > worst_swap_cal ) {
2948 worst_swap_cal = calnew;
2949 *ifac = j0;
2950 *adj = adja[j0];
2951 *conf0 = loc_conf0;
2952 *conf1 = loc_conf1;
2953 }
2954 }
2955 }
2956 else
2957 worst_swap_cal = 0.;
2958
2959 if ( *ifac < 0 ) {
2960 /* No valid config for the swap23 */
2961 if ( worst_split4bar_cal < MMG5_EPSOK ) return 0;
2962 return 1;
2963 }
2964
2965 assert ( *adj > 0 );
2966 if ( worst_swap_cal < worst_split4bar_cal ) {
2967 if ( worst_split4bar_cal < MMG5_EPSOK ) return 0;
2968 return 1;
2969 }
2970
2971 if ( worst_swap_cal < MMG5_EPSOK ) return 0;
2972
2973 return 2;
2974}
2975
2986static MMG5_int MMG5_anatet4(MMG5_pMesh mesh, MMG5_pSol met,MMG5_int *nf, int8_t typchk) {
2987 MMG5_pTetra pt;
2988 MMG5_pPoint ppt;
2989 MMG5_pxTetra pxt;
2990 int conf0,conf1,ifac,id_op;
2991 MMG5_int ier,ns,k,adj;
2992 int8_t nbdy,j;
2993#ifndef NDEBUG
2994 static int8_t mmgWarn=0;
2995#endif
2996
2997 ns = 0;
2998 for (k=1; k<=mesh->ne; k++) {
2999 pt = &mesh->tetra[k];
3000 if ( !MG_EOK(pt) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
3001 nbdy = 0;
3002 if ( pt->xt ) {
3003 pxt = &mesh->xtetra[pt->xt];
3004 for (j=0; j<4; j++)
3005 if ( ( pxt->ftag[j] & MG_BDY ) && (!(pxt->ftag[j] & MG_PARBDY)) ) nbdy++;
3006 }
3007
3008 /* Check for the number of boundary faces */
3009 if ( nbdy > 1 ) {
3010 id_op = MMG3D_anatet4_sim(mesh,met,k,typchk-1,&ifac,&conf0,&adj,&conf1);
3011 if ( !id_op ) {
3012#ifndef NDEBUG
3013 if ( (!(mesh->info.noswap && mesh->info.noinsert)) && !mmgWarn ) {
3014 mmgWarn=1;
3015 printf("\n ## Warning: %s: unable to swap or split at least 1 tetra"
3016 " with multiple boundary faces.\n",__func__);
3017 }
3018#endif
3019 continue;
3020 }
3021
3022 if ( id_op==1 ) {
3023 /* Split */
3024 ier = MMG5_split4bar(mesh,met,k,typchk-1);
3025 if ( !ier ) return -1;
3026 ns++;
3027 }
3028 else {
3029 assert ( id_op==2 );
3030
3031 /* Swap */
3032 ier = MMG3D_swap23(mesh,met,k,typchk-1,ifac,conf0,adj,conf1);
3033 if ( ier < 0 ) {
3034 return -1;
3035 }
3036 else if ( ier ) ++(*nf);
3037 }
3038 }
3039 /* Check for the number of boundary vertices */
3040 else {
3041 nbdy = 0;
3042 for (j=0; j<4; j++) {
3043 ppt = &mesh->point[pt->v[j]];
3044 if ( (ppt->tag & MG_BDY) && (!(ppt->tag & MG_PARBDY)) ) nbdy++;
3045 }
3046 if ( nbdy == 4 ) {
3047 if ( !mesh->info.noinsert ) {
3048 ier = MMG5_split4bar(mesh,met,k,typchk-1);
3049 if ( !ier ) return -1;
3050 ns++;
3051 }
3052 }
3053 }
3054 }
3055
3056 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
3057 fprintf(stdout," boundary elements: %7" MMG5_PRId " splitted %7" MMG5_PRId " swapped\n",ns,*nf);
3058 return ns;
3059}
3060
3071static MMG5_int MMG5_anatet4rid(MMG5_pMesh mesh, MMG5_pSol met,MMG5_int *nf, int8_t typchk) {
3072 MMG5_pTetra pt;
3073 MMG5_pPoint ppt;
3074 MMG5_int ier;
3075 MMG5_int ns,k;
3076 int8_t nrid,j;
3077
3078 ns = 0;
3079 for (k=1; k<=mesh->ne; k++) {
3080 pt = &mesh->tetra[k];
3081 if ( !MG_EOK(pt) || pt->ref < 0 || (pt->tag & MG_REQ) ) continue;
3082 nrid = 0;
3083
3084 for (j=0; j<4; j++) {
3085 ppt = &mesh->point[pt->v[j]];
3086 if ( MG_RID(ppt->tag) ) {
3087 /* Non-singular ridge point: metric ridge */
3088 nrid++;
3089 }
3090 }
3091 if ( nrid == 4 ) {
3092 if ( !mesh->info.noinsert ) {
3093 ier = MMG5_split4bar(mesh,met,k,typchk-1);
3094 if ( !ier ) return -1;
3095 ns++;
3096 }
3097 }
3098 }
3099 if ( (mesh->info.ddebug || abs(mesh->info.imprim) > 5) && ns > 0 )
3100 fprintf(stdout," boundary elements: %7" MMG5_PRId " splitted\n",ns);
3101 return ns;
3102}
3103
3104
3116int MMG5_anatet(MMG5_pMesh mesh,MMG5_pSol met,int8_t typchk, int patternMode) {
3117 int it,minit,maxit,lastit;
3118 MMG5_int nc,ns,nnc,nns,nnf,ier,nf;
3119
3120 /* pointer toward the suitable anatets function */
3121 if ( met->m && met->size==6 ) {
3122 /* if the aniso metric is not compatible with the geometry, the non
3123 * conformal surface operators may create spurious ridges */
3125 }
3126 else {
3128 }
3129
3130 /* analyze tetras : initial splitting */
3131 nns = nnc = nnf = it = 0;
3132 lastit = 0;
3133 minit = 3;
3134 maxit = 6;
3135 mesh->gap = 0.5;
3136 do {
3137 if ( typchk==2 && lastit==1 ) ++mesh->info.fem;
3138
3139
3140 /* split or swap tetra with more than 2 bdry faces */
3141 nf = ier = 0;
3142 if ( mesh->info.fem == typchk ) {
3143 ier = MMG5_anatet4(mesh,met,&nf,typchk);
3144 if ( ier < 0 ) return 0;
3145 }
3146 else if ( (met->size==6) && (typchk == 1) && lastit ) {
3147 ier = MMG5_anatet4rid(mesh,met,&nf,typchk);
3148 if ( ier < 0 ) return 0;
3149 }
3150 else ier = 0;
3151 ns = ier;
3152
3153 /* memory free */
3154 if ( mesh->adja )
3156
3157 if ( !mesh->info.noinsert ) {
3158 /* analyze surface tetras */
3159 ier = MMG3D_anatets(mesh,met,typchk);
3160
3161 if ( ier < 0 ) {
3162 fprintf(stderr,"\n ## Unable to complete surface mesh. Exit program.\n");
3163 return 0;
3164 }
3165 ns += ier;
3166
3167 if ( patternMode ) {
3168 if ( mesh->adja ) { MMG5_DEL_MEM(mesh,mesh->adja); }
3169
3170 /* analyze internal tetras */
3171 ier = MMG5_anatetv(mesh,met,typchk);
3172 if ( ier < 0 ) {
3173 fprintf(stderr,"\n ## Unable to complete volume mesh. Exit program.\n");
3174 return 0;
3175 }
3176 ns += ier;
3177 }
3178 }
3179
3180 if ( (!mesh->adja) && (!MMG3D_hashTetra(mesh,1)) ) {
3181 fprintf(stderr,"\n ## Hashing problem. Exit program.\n");
3182 return 0;
3183 }
3184
3185 /* collapse short edges */
3186 if ( !mesh->info.noinsert ) {
3187 nc = MMG5_coltet(mesh,met,typchk);
3188 if ( nc < 0 ) {
3189 fprintf(stderr,"\n ## Unable to collapse mesh. Exiting.\n");
3190 return 0;
3191 }
3192 }
3193 else nc = 0;
3194
3195 /* attempt to swap */
3196 if ( !mesh->info.noswap ) {
3197 ier = MMG5_swpmsh(mesh,met,NULL,typchk);
3198 if ( ier < 0 ) {
3199 fprintf(stderr,"\n ## Unable to improve mesh. Exiting.\n");
3200 return 0;
3201 }
3202 nf += ier;
3203
3205 if ( ier < 0 ) {
3206 fprintf(stderr,"\n ## Unable to improve mesh. Exiting.\n");
3207 return 0;
3208 }
3209 nf += ier;
3210 }
3211 else nf = 0;
3212
3213 nnc += nc;
3214 nns += ns;
3215 nnf += nf;
3216 if ( (abs(mesh->info.imprim) > 4 || mesh->info.ddebug) && ns+nc+nf > 0 ){
3217#ifndef MMG_PATTERN
3218 fprintf(stdout," ");
3219#endif
3220 fprintf(stdout," %8" MMG5_PRId " splitted, %8" MMG5_PRId " collapsed, %8" MMG5_PRId " swapped\n",ns,nc,nf);
3221 }
3222
3223 if ( it > minit-1 && ( !(ns+nc) || (MMG5_abs(nc-ns) < 0.1 * MG_MAX(nc,ns)) ) ) {
3224 ++lastit;
3225 if ( it > minit && lastit>2 ) break;
3226 }
3227 else if ( it+2 >= maxit ) {
3228 /* Last iteration because we reach the maximal number of iter */
3229 ++lastit;
3230 }
3231 else if ( lastit ) {
3232 /* Avoid the incrementation of mesh->info.fem if we have detected a last
3233 iteration but anatet4 leads to have nc, nf or ns != 0 so we perform a last
3234 iter */
3235 ++lastit;
3236 }
3237 }
3238 while ( ++it < maxit && (ns+nc+nf > 0 || lastit<3) );
3239
3240 if ( mesh->info.imprim > 0 ) {
3241 if ( (abs(mesh->info.imprim) < 5 || mesh->info.ddebug ) && nns+nnc > 0 ) {
3242#ifndef MMG_PATTERN
3243 fprintf(stdout," ");
3244#endif
3245 fprintf(stdout, " %8" MMG5_PRId " splitted, %8" MMG5_PRId " collapsed, %8" MMG5_PRId " swapped, %d iter.\n",nns,nnc,nnf,it);
3246 }
3247 }
3248
3249 return 1;
3250}
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)
Definition: boulep_3d.c:54
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:607
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:771
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:1843
int MMG5_coquil(MMG5_pMesh mesh, MMG5_int start, int ia, int64_t *list, int8_t *isbdy)
Definition: boulep_3d.c:1403
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:1144
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:515
#define MMG5_EPSD
#define MMG5_EPS
int MMG5_hashUpdate(MMG5_Hash *hash, MMG5_int a, MMG5_int b, MMG5_int k)
Definition: hash.c:400
int MMG5_hashNew(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int hsiz, MMG5_int hmax)
Definition: hash.c:526
MMG5_int MMG5_hashGet(MMG5_Hash *hash, MMG5_int a, MMG5_int b)
Definition: hash.c:495
int MMG5_hashEdge(MMG5_pMesh mesh, MMG5_Hash *hash, MMG5_int a, MMG5_int b, MMG5_int k)
Definition: hash.c:345
int MMG3D_hashTetra(MMG5_pMesh mesh, int pack)
Definition: hash_3d.c:122
int MMG5_hashPop(MMG5_Hash *hash, MMG5_int a, MMG5_int b)
Definition: hash_3d.c:761
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:462
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:101
API headers for the mmg3d library.
#define MMG3D_LMAX
Definition: libmmg3d.h:58
MMG5_int MMG3D_newPt(MMG5_pMesh mesh, double c[3], int16_t tag, MMG5_int src)
Definition: zaldy_3d.c:39
int MMG5_split4op(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:3676
int MMG5_split3op(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:2574
int MMG5_split1b(MMG5_pMesh, MMG5_pSol, int64_t *, int, MMG5_int, int, int8_t, int8_t)
Definition: split_3d.c:617
int MMG3D_split6_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:4258
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:467
int MMG3D_split3_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1571
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:694
#define MMG3D_LOPTS
MMG5_int MMG5_split4bar(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t)
Definition: split_3d.c:2964
int MMG5_split5(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:4053
int MMG3D_split3cone_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1790
int MMG3D_normalAdjaTri(MMG5_pMesh, MMG5_int, int8_t, int, double n[3])
Definition: split_3d.c:466
int MMG3D_simbulgept(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ilist, MMG5_int)
Definition: split_3d.c:324
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:136
MMG5_int MMG3D_indPt(MMG5_pMesh mesh, MMG5_int kp)
Definition: tools_3d.c:916
int MMG5_split4sf(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:3354
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:603
int MMG5_split6(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:4327
int MMG3D_split2sf_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:1075
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:3555
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:1430
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:3251
int MMG5_split2sf(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6], int8_t)
Definition: split_3d.c:1228
#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:1971
int MMG5_movbdynomintpt_iso(MMG5_pMesh, MMG5_pSol, MMG3D_pPROctree, int64_t *, int, int)
Definition: movpt_3d.c:1482
#define MMG3D_LOPTL
int MMG3D_split5_sim(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int vx[6])
Definition: split_3d.c:3962
#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:1642
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:2442
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:1368
@ MMG5_Tetrahedron
Definition: libmmgtypes.h:227
@ MMG5_Triangle
Definition: libmmgtypes.h:226
static int MMG3D_delPatternPts(MMG5_pMesh mesh, MMG5_Hash hash)
Definition: mmg3d1.c:1289
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:1654
static MMG5_int MMG5_anatet4(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int *nf, int8_t typchk)
Definition: mmg3d1.c:2986
void MMG3D_set_geom(MMG5_pMesh mesh, MMG5_pPoint ppt, int16_t tag, MMG5_int nmref, MMG5_int edgref, double no1[3], double no2[3], double to[3])
Definition: mmg3d1.c:59
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:732
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:2712
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:1601
int MMG5_anatet(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk, int patternMode)
Definition: mmg3d1.c:3116
static MMG5_int MMG5_anatet4rid(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int *nf, int8_t typchk)
Definition: mmg3d1.c:3071
static MMG5_int MMG3D_anatets_ani(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:2153
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:2688
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:2045
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:1927
int MMG3D_dichoto(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, MMG5_int *vx)
Definition: mmg3d1.c:140
static int MMG5_coltet(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:888
static MMG5_int MMG3D_anatets_iso(MMG5_pMesh mesh, MMG5_pSol met, int8_t typchk)
Definition: mmg3d1.c:2253
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:1180
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:1333
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:668
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, int16_t *tag, double o[3], double to[3], double no1[3], double no2[3], int64_t *list, int *ilist)
Definition: mmg3d1.c:1750
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:1868
#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:595
MMG5_hedge * item
Definition: libmmgtypes.h:597
int8_t parTyp
Definition: libmmgtypes.h:541
int8_t ddebug
Definition: libmmgtypes.h:532
uint8_t noswap
Definition: libmmgtypes.h:546
double hmin
Definition: libmmgtypes.h:518
uint8_t noinsert
Definition: libmmgtypes.h:546
double hmax
Definition: libmmgtypes.h:518
MMG5_pPar par
Definition: libmmgtypes.h:517
double hausd
Definition: libmmgtypes.h:518
int8_t fem
Definition: libmmgtypes.h:539
MMG mesh structure.
Definition: libmmgtypes.h:605
MMG5_Info info
Definition: libmmgtypes.h:651
MMG5_int ne
Definition: libmmgtypes.h:612
MMG5_pPoint point
Definition: libmmgtypes.h:641
MMG5_int mark
Definition: libmmgtypes.h:618
MMG5_int * adja
Definition: libmmgtypes.h:624
MMG5_pxPoint xpoint
Definition: libmmgtypes.h:642
double gap
Definition: libmmgtypes.h:608
MMG5_int base
Definition: libmmgtypes.h:616
MMG5_pTetra tetra
Definition: libmmgtypes.h:643
MMG5_int np
Definition: libmmgtypes.h:612
MMG5_pxTetra xtetra
Definition: libmmgtypes.h:644
double hmin
Definition: libmmgtypes.h:258
double hmax
Definition: libmmgtypes.h:259
double hausd
Definition: libmmgtypes.h:260
MMG5_int ref
Definition: libmmgtypes.h:261
int8_t elt
Definition: libmmgtypes.h:262
Structure to store points of a MMG mesh.
Definition: libmmgtypes.h:270
double n[3]
Definition: libmmgtypes.h:272
int16_t tag
Definition: libmmgtypes.h:284
double c[3]
Definition: libmmgtypes.h:271
MMG5_int xp
Definition: libmmgtypes.h:279
MMG5_int ref
Definition: libmmgtypes.h:278
MMG5_int flag
Definition: libmmgtypes.h:282
double * m
Definition: libmmgtypes.h:671
MMG5_int v[4]
Definition: libmmgtypes.h:403
MMG5_int xt
Definition: libmmgtypes.h:407
MMG5_int flag
Definition: libmmgtypes.h:409
MMG5_int ref
Definition: libmmgtypes.h:404
MMG5_int mark
Definition: libmmgtypes.h:406
double qual
Definition: libmmgtypes.h:402
int16_t tag
Definition: libmmgtypes.h:410
MMG5_int edg[3]
Definition: libmmgtypes.h:339
int16_t tag[3]
Definition: libmmgtypes.h:342
MMG5_int ref
Definition: libmmgtypes.h:335
MMG5_int flag
Definition: libmmgtypes.h:341
MMG5_int v[3]
Definition: libmmgtypes.h:334
Structure to store surface points of a MMG mesh.
Definition: libmmgtypes.h:294
int8_t nnor
Definition: libmmgtypes.h:297
double n2[3]
Definition: libmmgtypes.h:295
double n1[3]
Definition: libmmgtypes.h:295
Structure to store the surface tetrahedra of a MMG mesh.
Definition: libmmgtypes.h:418
int16_t ftag[4]
Definition: libmmgtypes.h:423
int16_t tag[6]
Definition: libmmgtypes.h:425
int8_t ori
Definition: libmmgtypes.h:427
MMG5_int edg[6]
Definition: libmmgtypes.h:421
MMG5_int ref[4]
Definition: libmmgtypes.h:419