58 MMG5_int it1,MMG5_int it2,int8_t typchk) {
64 double b0[3],b1[3],n[3],v[3],c[3],ux,uy,uz,ps,disnat,dischg;
65 double cal1,cal2,calnat,calchg,calold,calnew,caltmp,hausd;
66 MMG5_int iel,iel1,iel2,np,nq,na1,na2,k,nminus,nplus,
info,refm;
68 int8_t ifa1,ifa2,ia,ip,iq,ia1,ia2,j,isshell,
ier;
93 for (k=0; k<ilist; k++) {
96 if ( pt->
ref == refm ) {
103 if ( nplus == 1 || nminus == 1 )
return 0;
111 assert( 0<=ifa1 && ifa1<4 &&
"unexpected local face idx");
112 assert( 0<=ifa2 && ifa2<4 &&
"unexpected local face idx");
116 for (ia1=0; ia1<3; ia1++) {
117 if ( (tt1.
v[ia1] != np) && (tt1.
v[ia1] != nq) )
break;
120 if ( ia1==3 )
return 0;
126 for (ia2=0; ia2<3; ia2++) {
127 if ( (tt2.
v[ia2] != np) && (tt2.
v[ia2] != nq) )
break;
131 if ( ia2 ==3 )
return 0;
140 ps = b0[0]*b1[0] + b0[1]*b1[1] + b0[2]*b1[2];
151 if (
ier < 0 )
return -1;
152 else if ( !
ier )
return 0;
153 ps = b0[0]*n[0] + b0[1]*n[1] + b0[2]*n[2];
155 if ( ps < mesh->
info.
dhd )
return 0;
160 if (
ier<0 )
return -1;
161 else if ( !
ier )
return 0;
162 ps = b0[0]*n[0] + b0[1]*n[1] + b0[2]*n[2];
164 if ( ps < mesh->
info.
dhd )
return 0;
169 if (
ier<0 )
return -1;
170 else if ( !
ier )
return 0;
171 ps = b1[0]*n[0] + b1[1]*n[1] + b1[2]*n[2];
173 if ( ps < mesh->
info.
dhd )
return 0;
178 if (
ier<0 )
return -1;
179 else if ( !
ier )
return 0;
180 ps = b1[0]*n[0] + b1[1]*n[1] + b1[2]*n[2];
182 if ( ps < mesh->
info.
dhd )
return 0;
197 if ( tt1.
ref == tt2.
ref ) {
215 if ( tt1.
ref!=par->
ref && tt2.
ref !=par->
ref )
continue;
226 if ( tt1.
ref!=par->
ref && tt2.
ref !=par->
ref )
continue;
244 for ( k=0; k<ilist; ++k ) {
246 if ( par->
ref != pt->
ref )
continue;
257 for ( k=0; k<ilist; ++k ) {
259 if ( par->
ref != pt->
ref )
continue;
267 ux = p1->
c[0] - p0->
c[0];
268 uy = p1->
c[1] - p0->
c[1];
269 uz = p1->
c[2] - p0->
c[2];
276 disnat = c[0]*c[0] + c[1]*c[1] + c[2]*c[2];
282 disnat =
MG_MAX(disnat, c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
286 disnat =
MG_MAX(disnat,hausd * hausd);
290 ux = p1->
c[0] - p0->
c[0];
291 uy = p1->
c[1] - p0->
c[1];
292 uz = p1->
c[2] - p0->
c[2];
301 dischg = c[0]*c[0] + c[1]*c[1] + c[2]*c[2];
307 dischg =
MG_MAX(dischg,c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
308 dischg =
MG_MAX(dischg,hausd * hausd);
310 if ( dischg > disnat )
return 0;
312 if ( typchk==1 && met->
size > 1 && met->
m ) {
317 cal1 = MMG5_caltri(
mesh,met,&tt1);
318 cal2 = MMG5_caltri(
mesh,met,&tt2);
321 calnat =
MG_MIN(cal1,cal2);
322 for (j=0; j<3; j++) {
323 if ( tt1.
v[j] == nq ) tt1.
v[j] = na2;
324 if ( tt2.
v[j] == np ) tt2.
v[j] = na1;
327 if ( typchk==1 && met->
size > 1 && met->
m ) {
332 cal1 = MMG5_caltri(
mesh,met,&tt1);
333 cal2 = MMG5_caltri(
mesh,met,&tt2);
336 calchg =
MG_MIN(cal1,cal2);
337 if ( calchg < 1.01 * calnat )
return 0;
342 ppt0->
c[0] = 0.5*(p0->
c[0] + p1->
c[0]);
343 ppt0->
c[1] = 0.5*(p0->
c[1] + p1->
c[1]);
344 ppt0->
c[2] = 0.5*(p0->
c[2] + p1->
c[2]);
347 if ( typchk == 1 && (met->
size>1) ) {
352 if ( MMG5_intmet(
mesh,met,list[0]/6,list[0]%6,0,0.5) <= 0 )
358 calold = calnew = DBL_MAX;
359 for (k=0; k<ilist; k++) {
364 assert ( isfinite(calold) );
366 ia1 = ia2 = ip = iq = -1;
367 for (j=0; j< 4; j++) {
368 if (pt->
v[j] == np) ip = j;
369 else if (pt->
v[j] == nq) iq = j;
370 else if ( ia1 < 0 ) ia1 = j;
373 assert((ip >= 0) && (iq >= 0) && (ia1 >= 0) && (ia2 >= 0));
374 isshell = (pt->
v[ia1] == na1 || pt->
v[ia2] == na1);
379 if ( typchk==1 && met->
size > 1 && met->
m )
389 MMG5_int adj =
mesh->
adja[4*(iel-1)+1+ip];
402 if ( typchk==1 && met->
size > 1 && met->
m )
407 calnew =
MG_MIN(calnew,caltmp);
412 if ( typchk==1 && met->
size > 1 && met->
m )
421 MMG5_int adj =
mesh->
adja[4*(iel-1)+1+iq];
434 if ( typchk==1 && met->
size > 1 && met->
m )
441 calnew =
MG_MIN(calnew,caltmp);
444 if ( calold <
MMG5_EPSOK && calnew <= calold )
return 0;
446 else if ( calnew < 0.3 * calold )
return 0;
472 MMG5_int iel,np,nq,nm,src,iel1;
474 int8_t ia,iface1,j,ipa,im;
501 if ( (pt1->
v[ipa] != np)&&(pt1->
v[ipa] != nq) ) {
511 c[0] = 0.5*( p0->
c[0] + p1->
c[0]);
512 c[1] = 0.5*( p0->
c[1] + p1->
c[1]);
513 c[2] = 0.5*( p0->
c[2] + p1->
c[2]);
522 fprintf(stderr,
"\n ## Error: %s: unable to allocate a"
523 " new point\n",__func__);
530 if ( typchk == 1 && (met->
size>1) ) {
534 if ( MMG5_intmet(
mesh,met,iel,ia,nm,0.5)<=0 )
return 0;
543 fprintf(stderr,
"\n ## Warning: %s: unable to swap boundary edge.\n",
553 memset(list,0,(
MMG3D_LMAX+2)*
sizeof(MMG5_int));
554 for (j=0; j<3; j++) {
556 if ( pt1->
v[im] == nm )
break;
558 if ( pt1->
v[im] != nm ){
560 fprintf(stderr,
"\n # Warning: %s: pt1->v[im] != nm.\n",__func__);
565 assert(list[0]/4 == iel1);
566 assert(pt1->
v[ipa] == na);
570 fprintf(stderr,
"\n ## Warning: %s: unable to swap boundary edge.\n",
580 assert (
ier &&
"Unable to collapse the point created during the boundary swap");
604 int ifac,
int conf0,MMG5_int adj,
int conf1) {
608 MMG5_int xt1,k1,*adja,iel,np;
609 MMG5_int adj0_2,adj0_3,adj1_1,adj1_2,adj1_3;
611 uint8_t tau0[4],tau1[4];
612 const uint8_t *taued0,*taued1;
617 assert ( ifac>=0 && adj>0 );
643 tau0[0] = 1; tau0[1] = 0; tau0[2] = 3; tau0[3] = 2;
647 tau0[0] = 2; tau0[1] = 0; tau0[2] = 1; tau0[3] = 3;
651 tau0[0] = 3; tau0[1] = 0; tau0[2] = 2; tau0[3] = 1;
657 tau0[0] = 0; tau0[1] = 1; tau0[2] = 2; tau0[3] = 3;
665 assert(pt0->
ref == pt1->
ref);
669 tau1[0] = 0; tau1[1] = 2; tau1[2] = 3; tau1[3] = 1;
673 tau1[0] = 0; tau1[1] = 3; tau1[2] = 1; tau1[3] = 2;
677 tau1[0] = 1; tau1[1] = 0; tau1[2] = 3; tau1[3] = 2;
681 tau1[0] = 1; tau1[1] = 3; tau1[2] = 2; tau1[3] = 0;
685 tau1[0] = 1; tau1[1] = 2; tau1[2] = 0; tau1[3] = 3;
689 tau1[0] = 2; tau1[1] = 0; tau1[2] = 1; tau1[3] = 3;
693 tau1[0] = 2; tau1[1] = 1; tau1[2] = 3; tau1[3] = 0;
697 tau1[0] = 2; tau1[1] = 3; tau1[2] = 0; tau1[3] = 1;
701 tau1[0] = 3; tau1[1] = 0; tau1[2] = 2; tau1[3] = 1;
705 tau1[0] = 3; tau1[1] = 2; tau1[2] = 1; tau1[3] = 0;
709 tau1[0] = 3; tau1[1] = 1; tau1[2] = 0; tau1[3] = 2;
714 tau1[0] = 0; tau1[1] = 1; tau1[2] = 2; tau1[3] = 3;
722 np = pt1->
v[tau1[0]];
728 fprintf(stderr,
"\n ## Error: %s: unable to allocate"
729 " a new element.\n",__func__);
731 fprintf(stderr,
" Exit program.\n");
740 pt0->
v[tau0[1]] = np;
743 pt1->
v[tau0[2]] = np;
746 ptnew->
v[tau0[3]] = np;
756 adj0_2 = adja[tau0[2]];
757 adj0_3 = adja[tau0[3]];
760 adj1_1 = adja[tau1[1]];
761 adj1_2 = adja[tau1[2]];
762 adj1_3 = adja[tau1[3]];
766 adja[tau0[0]] = adj1_1;
767 adja[tau0[2]] = 4*k1 + tau0[1] ;
768 adja[tau0[3]] = 4*iel + tau0[1] ;
770 mesh->
adja[4*(adj1_1/4-1) + 1 + adj1_1%4] = 4*k + tau0[0];
773 adja[tau0[0]] = adj1_3;
774 adja[tau0[1]] = 4*k + tau0[2] ;
775 adja[tau0[2]] = adj0_2;
776 adja[tau0[3]] = 4*iel + tau0[2] ;
778 mesh->
adja[4*(adj1_3/4-1) + 1 + adj1_3%4] = 4*k1 + tau0[0];
780 mesh->
adja[4*(adj0_2/4-1) + 1 + adj0_2%4] = 4*k1 + tau0[2];
783 adja[tau0[0]] = adj1_2;
784 adja[tau0[1]] = 4*k + tau0[3] ;
785 adja[tau0[2]] = 4*k1 + tau0[3] ;
786 adja[tau0[3]] = adj0_3;
788 mesh->
adja[4*(adj1_2/4-1) + 1 + adj1_2%4] = 4*iel + tau0[0];
790 mesh->
adja[4*(adj0_3/4-1) + 1 + adj0_3%4] = 4*iel + tau0[3];
796 xt[0].
tag[taued0[0]] = 0;
797 xt[0].
tag[taued0[3]] = 0;
798 xt[0].
tag[taued0[4]] = 0;
800 xt[0].
edg[taued0[0]] = 0;
801 xt[0].
edg[taued0[3]] = 0;
802 xt[0].
edg[taued0[4]] = 0;
804 xt[0].
ref[ tau0[0]] = 0;
805 xt[0].
ref[ tau0[2]] = 0;
806 xt[0].
ref[ tau0[3]] = 0;
807 xt[0].
ftag[tau0[0]] = 0;
808 xt[0].
ftag[tau0[2]] = 0;
809 xt[0].
ftag[tau0[3]] = 0;
811 MG_SET(xt[0].ori, tau0[0]);
812 MG_SET(xt[0].ori, tau0[2]);
813 MG_SET(xt[0].ori, tau0[3]);
816 xt[1].
tag[taued0[1]] = 0;
817 xt[1].
tag[taued0[3]] = 0;
818 xt[1].
tag[taued0[5]] = 0;
820 xt[1].
edg[taued0[1]] = 0;
821 xt[1].
edg[taued0[3]] = 0;
822 xt[1].
edg[taued0[5]] = 0;
824 xt[1].
ref[ tau0[0]] = 0;
825 xt[1].
ref[ tau0[1]] = 0;
826 xt[1].
ref[ tau0[3]] = 0;
827 xt[1].
ftag[tau0[0]] = 0;
828 xt[1].
ftag[tau0[1]] = 0;
829 xt[1].
ftag[tau0[3]] = 0;
831 MG_SET(xt[1].ori, tau0[0]);
832 MG_SET(xt[1].ori, tau0[1]);
833 MG_SET(xt[1].ori, tau0[3]);
836 xt[1].
tag[taued0[2]] = 0;
837 xt[1].
tag[taued0[4]] = 0;
838 xt[1].
tag[taued0[5]] = 0;
840 xt[1].
edg[taued0[2]] = 0;
841 xt[1].
edg[taued0[4]] = 0;
842 xt[1].
edg[taued0[5]] = 0;
844 xt[1].
ref[ tau0[0]] = 0;
845 xt[1].
ref[ tau0[1]] = 0;
846 xt[1].
ref[ tau0[2]] = 0;
847 xt[1].
ftag[tau0[0]] = 0;
848 xt[1].
ftag[tau0[1]] = 0;
849 xt[1].
ftag[tau0[2]] = 0;
851 MG_SET(xt[1].ori, tau0[0]);
852 MG_SET(xt[1].ori, tau0[1]);
853 MG_SET(xt[1].ori, tau0[2]);
861 xt[0].
tag[taued0[0]] = 0;
862 xt[0].
tag[taued0[3]] = pxt1->
tag[taued1[2]];
863 xt[0].
tag[taued0[4]] = pxt1->
tag[taued1[1]];
865 xt[0].
edg[taued0[0]] = 0;
866 xt[0].
edg[taued0[3]] = pxt1->
edg[taued1[2]];
867 xt[0].
edg[taued0[4]] = pxt1->
edg[taued1[1]];
869 xt[0].
ref[ tau0[0]] = pxt1->
ref[tau1[1]];
870 xt[0].
ref[ tau0[2]] = 0;
871 xt[0].
ref[ tau0[3]] = 0;
872 xt[0].
ftag[tau0[0]] = pxt1->
ftag[tau1[1]];
873 xt[0].
ftag[tau0[2]] = 0;
874 xt[0].
ftag[tau0[3]] = 0;
877 MG_SET(xt[0].ori, tau0[2]);
878 MG_SET(xt[0].ori, tau0[3]);
881 xt[1].
tag[taued0[1]] = 0;
882 xt[1].
tag[taued0[3]] = pxt1->
tag[taued1[0]];
883 xt[1].
tag[taued0[5]] = pxt1->
tag[taued1[2]];
885 xt[1].
edg[taued0[1]] = 0;
886 xt[1].
edg[taued0[3]] = pxt1->
edg[taued1[0]];
887 xt[1].
edg[taued0[5]] = pxt1->
edg[taued1[2]];
889 xt[1].
ref[ tau0[0]] = pxt1->
ref[tau1[3]];
890 xt[1].
ref[ tau0[1]] = 0;
891 xt[1].
ref[ tau0[3]] = 0;
892 xt[1].
ftag[tau0[0]] = pxt1->
ftag[tau1[3]];
893 xt[1].
ftag[tau0[1]] = 0;
894 xt[1].
ftag[tau0[3]] = 0;
897 MG_SET(xt[1].ori, tau0[1]);
898 MG_SET(xt[1].ori, tau0[3]);
901 xt[2].
tag[taued0[2]] = 0;
902 xt[2].
tag[taued0[4]] = pxt1->
tag[taued1[0]];
903 xt[2].
tag[taued0[5]] = pxt1->
tag[taued1[1]];
905 xt[2].
edg[taued0[2]] = 0;
906 xt[2].
edg[taued0[4]] = pxt1->
edg[taued1[0]];
907 xt[2].
edg[taued0[5]] = pxt1->
edg[taued1[1]];
909 xt[2].
ref[ tau0[0]] = pxt1->
ref[tau1[2]];
910 xt[2].
ref[ tau0[1]] = 0;
911 xt[2].
ref[ tau0[2]] = 0;
912 xt[2].
ftag[tau0[0]] = pxt1->
ftag[tau1[2]];
913 xt[2].
ftag[tau0[1]] = 0;
914 xt[2].
ftag[tau0[2]] = 0;
917 MG_SET(xt[2].ori, tau0[1]);
918 MG_SET(xt[2].ori, tau0[2]);
922 isxt[0] = isxt[1] = isxt[2] = 0;
923 for (i=0; i<4; i++) {
924 if ( xt[0].ref[i] || xt[0].ftag[i] ) isxt[0] = 1;
925 if ( xt[1].ref[i] || xt[1].ftag[i] ) isxt[1] = 1;
926 if ( xt[2].ref[i] || xt[2].ftag[i] ) isxt[2] = 1;
929 assert ( (isxt[0] && isxt[1]) || (isxt[1] && isxt[2]) || (isxt[2] && isxt[0]) );
943 "larger xtetra table",
945 fprintf(stderr,
" Exit program.\n");
966 "larger xtetra table",
968 fprintf(stderr,
" Exit program.\n");
982 "larger xtetra table",
984 fprintf(stderr,
" Exit program.\n");
1007 "larger xtetra table",
1009 fprintf(stderr,
" Exit program.\n");
1019 if ( (!metRidTyp) && met->
m && met->
size>1 )
1024 if ( (!metRidTyp) && met->
m && met->
size>1 )
1029 if ( (!metRidTyp) && met->
m && met->
size>1 )
int MMG5_BezierEdge(MMG5_pMesh mesh, MMG5_int ip0, MMG5_int ip1, double b0[3], double b1[3], int8_t ised, double v[3])
int MMG5_boulevolp(MMG5_pMesh mesh, MMG5_int start, int ip, int64_t *list)
MMG5_int MMG5_colver(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ilist, int8_t indq, int8_t typchk)
static double MMG5_orcal(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int iel)
int MMG3D_intmet33_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t i, MMG5_int ip, double s)
API headers for the mmg3d library.
MMG5_int MMG3D_newPt(MMG5_pMesh mesh, double c[3], int16_t tag, MMG5_int src)
int MMG5_split1b(MMG5_pMesh, MMG5_pSol, int64_t *, int, MMG5_int, int, int8_t, int8_t)
MMG5_int MMG3D_newElt(MMG5_pMesh mesh)
int MMG3D_normalAdjaTri(MMG5_pMesh, MMG5_int, int8_t, int, double n[3])
void MMG3D_delPt(MMG5_pMesh mesh, MMG5_int ip)
static const uint8_t MMG5_permedge[12][6]
double MMG5_caltet33_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTetra pt)
static const uint8_t MMG5_iare[6][2]
vertices of extremities of the edges of the tetra
#define MMG3D_POINT_REALLOC(mesh, sol, ip, wantedGap, law, o, tag, src)
void MMG5_tet2tri(MMG5_pMesh mesh, MMG5_int k, int8_t ie, MMG5_Tria *ptt)
int MMG5_norface(MMG5_pMesh mesh, MMG5_int k, int iface, double v[3])
#define MMG3D_TETRA_REALLOC(mesh, jel, wantedGap, law)
static const uint8_t MMG5_idir[4][3]
idir[i]: vertices of face opposite to vertex i
#define MMG5_INCREASE_MEM_MESSAGE()
static const uint8_t MMG5_iprv2[3]
#define MG_EDG_OR_NOM(tag)
#define MMG5_TAB_RECALLOC(mesh, ptr, initSize, wantedGap, type, message, law)
#define MG_GEO_OR_NOM(tag)
#define MG_GET(flag, bit)
static const uint8_t MMG5_inxt2[6]
int MMG5_norpts(MMG5_pMesh, MMG5_int, MMG5_int, MMG5_int, double *)
double MMG5_caltri33_ani(MMG5_pMesh mesh, MMG5_pSol met, MMG5_pTria pt)
#define MG_SET(flag, bit)
Structure to store points of a MMG mesh.
Structure to store the surface tetrahedra of a MMG mesh.
int MMG3D_swap23(MMG5_pMesh mesh, MMG5_pSol met, MMG5_int k, int8_t metRidTyp, int ifac, int conf0, MMG5_int adj, int conf1)
int MMG5_swpbdy(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ret, MMG5_int it1, MMG3D_pPROctree PROctree, int8_t typchk)
int MMG5_chkswpbdy(MMG5_pMesh mesh, MMG5_pSol met, int64_t *list, int ilist, MMG5_int it1, MMG5_int it2, int8_t typchk)