76 ps = nt[0]*n[0] + nt[1]*n[1] + nt[2]*n[2];
77 if ( ps > 0.0 )
return 1;
79 for (j=1; j<=(ilist-1)/2; j++) {
81 list[j] = list[ilist -j];
102 if(pt->
v[ip] == nump)
break;
110 if(pt->
v[ipt] == numq)
return 1;
114 assert(pt->
v[ipt] == numq);
117 for(l=0;l<ilist-1;l++){
146 double *no1,
double *no2,
double *to ) {
148 double ux,uy,uz,n01[3],n02[3],n11[3],n12[3],t0[3],t1[3];
149 double ps,ps2,b0[3],b1[3],bn[3],ll,il,ntemp[3],dd,alpha;
157 ux = p1->
c[0] - p0->
c[0];
158 uy = p1->
c[1] - p0->
c[1];
159 uz = p1->
c[2] - p0->
c[2];
160 ll = ux*ux + uy*uy + uz*uz;
173 memcpy(t0,&(p0->
n[0]),3*
sizeof(
double));
174 ps = t0[0]*ux + t0[1]*uy + t0[2]*uz;
189 memcpy(t1,&(p1->
n[0]),3*
sizeof(
double));
190 ps = - ( t1[0]*ux + t1[1]*uy + t1[2]*uz );
199 b0[0] = p0->
c[0] + alpha * t0[0];
200 b0[1] = p0->
c[1] + alpha * t0[1];
201 b0[2] = p0->
c[2] + alpha * t0[2];
203 b1[0] = p1->
c[0] + alpha * t1[0];
204 b1[1] = p1->
c[1] + alpha * t1[1];
205 b1[2] = p1->
c[2] + alpha * t1[2];
207 o[0] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[0] + 3.0*s*(1.0-s)*(1.0-s)*b0[0] + \
208 3.0*s*s*(1.0-s)*b1[0] + s*s*s*p1->
c[0];
210 o[1] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[1] + 3.0*s*(1.0-s)*(1.0-s)*b0[1] + \
211 3.0*s*s*(1.0-s)*b1[1] + s*s*s*p1->
c[1];
213 o[2] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[2] + 3.0*s*(1.0-s)*(1.0-s)*b0[2] + \
214 3.0*s*s*(1.0-s)*b1[2] + s*s*s*p1->
c[2];
219 memcpy(to,t0,3*
sizeof(
double));
241 ps = n01[0] * n11[0] + n01[1] * n11[1] + n01[2] * n11[2];
242 ps2 = n01[0] * n12[0] + n01[1] * n12[1] + n01[2] * n12[2];
244 memcpy(ntemp,n11,3*
sizeof(
double));
245 memcpy(n11,n12,3*
sizeof(
double));
246 memcpy(n12,ntemp,3*
sizeof(
double));
251 ps = ux*(n01[0] + n11[0]) + uy*(n01[1] + n11[1]) + uz*(n01[2] + n11[2]);
254 bn[0] = n01[0] + n11[0] -ps*ux;
255 bn[1] = n01[1] + n11[1] -ps*uy;
256 bn[2] = n01[2] + n11[2] -ps*uz;
258 dd = bn[0]*bn[0] + bn[1]*bn[1] + bn[2]*bn[2];
265 no1[0] = (1.0-s)*(1.0-s)*n01[0] + 2.0*s*(1.0-s)*bn[0] + s*s*n11[0];
266 no1[1] = (1.0-s)*(1.0-s)*n01[1] + 2.0*s*(1.0-s)*bn[1] + s*s*n11[1];
267 no1[2] = (1.0-s)*(1.0-s)*n01[2] + 2.0*s*(1.0-s)*bn[2] + s*s*n11[2];
268 dd = no1[0]*no1[0] + no1[1]*no1[1] + no1[2]*no1[2];
277 ps = ux*(n02[0] + n12[0]) + uy*(n02[1] + n12[1]) + uz*(n02[2] + n12[2]);
280 bn[0] = n02[0] + n12[0] -ps*ux;
281 bn[1] = n02[1] + n12[1] -ps*uy;
282 bn[2] = n02[2] + n12[2] -ps*uz;
284 dd = bn[0]*bn[0] + bn[1]*bn[1] + bn[2]*bn[2];
291 no2[0] = (1.0-s)*(1.0-s)*n02[0] + 2.0*s*(1.0-s)*bn[0] + s*s*n12[0];
292 no2[1] = (1.0-s)*(1.0-s)*n02[1] + 2.0*s*(1.0-s)*bn[1] + s*s*n12[1];
293 no2[2] = (1.0-s)*(1.0-s)*n02[2] + 2.0*s*(1.0-s)*bn[2] + s*s*n12[2];
294 dd = no2[0]*no2[0] + no2[1]*no2[1] + no2[2]*no2[2];
301 to[0] = no1[1]*no2[2] - no1[2]*no2[1];
302 to[1] = no1[2]*no2[0] - no1[0]*no2[2];
303 to[2] = no1[0]*no2[1] - no1[1]*no2[0];
307 ps = t0[0]*t1[0] + t0[1]*t1[1] + t0[2]*t1[2];
313 to[0] = (1.0-s)*t0[0] + s*t1[0];
314 to[1] = (1.0-s)*t0[1] + s*t1[1];
315 to[2] = (1.0-s)*t0[2] + s*t1[2];
318 ps = to[0]*no1[0] + to[1]*no1[1] + to[2]*no1[2];
319 to[0] = to[0] -ps*no1[0];
320 to[1] = to[1] -ps*no1[1];
321 to[2] = to[2] -ps*no1[2];
324 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
353 double *no,
double *to) {
355 double ux,uy,uz,n01[3],n02[3],n11[3],n12[3],ntemp[3],t0[3],t1[3];
356 double ps,b0[3],b1[3],bn[3],ll,il,dd,alpha,ps2;
361 ux = p1->
c[0] - p0->
c[0];
362 uy = p1->
c[1] - p0->
c[1];
363 uz = p1->
c[2] - p0->
c[2];
364 ll = ux*ux + uy*uy + uz*uz;
379 memcpy(t0,&(p0->
n[0]),3*
sizeof(
double));
380 ps = t0[0]*ux + t0[1]*uy + t0[2]*uz;
394 memcpy(t1,&(p1->
n[0]),3*
sizeof(
double));
395 ps = - ( t1[0]*ux + t1[1]*uy + t1[2]*uz );
405 b0[0] = p0->
c[0] + alpha * t0[0];
406 b0[1] = p0->
c[1] + alpha * t0[1];
407 b0[2] = p0->
c[2] + alpha * t0[2];
409 b1[0] = p1->
c[0] + alpha * t1[0];
410 b1[1] = p1->
c[1] + alpha * t1[1];
411 b1[2] = p1->
c[2] + alpha * t1[2];
413 o[0] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[0] + 3.0*s*(1.0-s)*(1.0-s)*b0[0] + \
414 3.0*s*s*(1.0-s)*b1[0] + s*s*s*p1->
c[0];
416 o[1] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[1] + 3.0*s*(1.0-s)*(1.0-s)*b0[1] + \
417 3.0*s*s*(1.0-s)*b1[1] + s*s*s*p1->
c[1];
419 o[2] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[2] + 3.0*s*(1.0-s)*(1.0-s)*b0[2] + \
420 3.0*s*s*(1.0-s)*b1[2] + s*s*s*p1->
c[2];
424 memcpy(to,t0,3*
sizeof(
double));
447 ps = n01[0] * n11[0] + n01[1] * n11[1] + n01[2] * n11[2];
448 ps2 = n01[0] * n12[0] + n01[1] * n12[1] + n01[2] * n12[2];
450 memcpy(ntemp,n11,3*
sizeof(
double));
451 memcpy(n11,n12,3*
sizeof(
double));
452 memcpy(n12,ntemp,3*
sizeof(
double));
456 ps = ux*(n01[0] + n11[0]) + uy*(n01[1] + n11[1]) + uz*(n01[2] + n11[2]);
459 bn[0] = n01[0] + n11[0] -ps*ux;
460 bn[1] = n01[1] + n11[1] -ps*uy;
461 bn[2] = n01[2] + n11[2] -ps*uz;
463 dd = bn[0]*bn[0] + bn[1]*bn[1] + bn[2]*bn[2];
470 no[0] = (1.0-s)*(1.0-s)*n01[0] + 2.0*s*(1.0-s)*bn[0] + s*s*n11[0];
471 no[1] = (1.0-s)*(1.0-s)*n01[1] + 2.0*s*(1.0-s)*bn[1] + s*s*n11[1];
472 no[2] = (1.0-s)*(1.0-s)*n01[2] + 2.0*s*(1.0-s)*bn[2] + s*s*n11[2];
474 dd = no[0]*no[0] + no[1]*no[1] + no[2]*no[2];
483 ps = t0[0]*t1[0] + t0[1]*t1[1] + t0[2]*t1[2];
489 to[0] = (1.0-s)*t0[0] + s*t1[0];
490 to[1] = (1.0-s)*t0[1] + s*t1[1];
491 to[2] = (1.0-s)*t0[2] + s*t1[2];
494 ps = to[0]*no[0] + to[1]*no[1] + to[2]*no[2];
495 to[0] = to[0] -ps*no[0];
496 to[1] = to[1] -ps*no[1];
497 to[2] = to[2] -ps*no[2];
499 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
529 double ux,uy,uz,il,ll,ps,alpha,dd;
530 double t0[3],t1[3],b0[3],b1[3],n0[3],n1[3],bn[3];
536 ux = p1->
c[0] - p0->
c[0];
537 uy = p1->
c[1] - p0->
c[1];
538 uz = p1->
c[2] - p0->
c[2];
539 ll = ux*ux + uy*uy + uz*uz;
553 memcpy(t0,&(p0->
n[0]),3*
sizeof(
double));
554 ps = t0[0]*ux + t0[1]*uy + t0[2]*uz;
567 memcpy(t1,&(p1->
n[0]),3*
sizeof(
double));
568 ps = - ( t1[0]*ux + t1[1]*uy + t1[2]*uz );
578 b0[0] = p0->
c[0] + alpha * t0[0];
579 b0[1] = p0->
c[1] + alpha * t0[1];
580 b0[2] = p0->
c[2] + alpha * t0[2];
582 b1[0] = p1->
c[0] + alpha * t1[0];
583 b1[1] = p1->
c[1] + alpha * t1[1];
584 b1[2] = p1->
c[2] + alpha * t1[2];
586 o[0] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[0] + 3.0*s*(1.0-s)*(1.0-s)*b0[0] + \
587 3.0*s*s*(1.0-s)*b1[0] + s*s*s*p1->
c[0];
589 o[1] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[1] + 3.0*s*(1.0-s)*(1.0-s)*b0[1] + \
590 3.0*s*s*(1.0-s)*b1[1] + s*s*s*p1->
c[1];
592 o[2] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[2] + 3.0*s*(1.0-s)*(1.0-s)*b0[2] + \
593 3.0*s*s*(1.0-s)*b1[2] + s*s*s*p1->
c[2];
597 memcpy(to,t0,3*
sizeof(
double));
618 ps = ux*(n0[0] + n1[0]) + uy*(n0[1] + n1[1]) + uz*(n0[2] + n1[2]);
621 bn[0] = n0[0] + n1[0] -ps*ux;
622 bn[1] = n0[1] + n1[1] -ps*uy;
623 bn[2] = n0[2] + n1[2] -ps*uz;
625 dd = bn[0]*bn[0] + bn[1]*bn[1] + bn[2]*bn[2];
632 no[0] = (1.0-s)*(1.0-s)*n0[0] + 2.0*s*(1.0-s)*bn[0] + s*s*n1[0];
633 no[1] = (1.0-s)*(1.0-s)*n0[1] + 2.0*s*(1.0-s)*bn[1] + s*s*n1[1];
634 no[2] = (1.0-s)*(1.0-s)*n0[2] + 2.0*s*(1.0-s)*bn[2] + s*s*n1[2];
636 dd = no[0]*no[0] + no[1]*no[1] + no[2]*no[2];
646 ps = t0[0]*t1[0] + t0[1]*t1[1] + t0[2]*t1[2];
652 to[0] = (1.0-s)*t0[0] + s*t1[0];
653 to[1] = (1.0-s)*t0[1] + s*t1[1];
654 to[2] = (1.0-s)*t0[2] + s*t1[2];
658 ps = to[0]*no[0] + to[1]*no[1] + to[2]*no[2];
659 to[0] = to[0] -ps*no[0];
660 to[1] = to[1] -ps*no[1];
661 to[2] = to[2] -ps*no[2];
664 dd = to[0]*to[0] + to[1]*to[1] + to[2]*to[2];
696 double b0[3],b1[3],bn[3],t0[3],t1[3],np0[3],np1[3],alpha,ux,uy,uz,ps1,ps2,ll;
697 double il,dd,*n1,*n2;
702 ux = p1->
c[0] - p0->
c[0];
703 uy = p1->
c[1] - p0->
c[1];
704 uz = p1->
c[2] - p0->
c[2];
706 ll = ux*ux + uy*uy + uz*uz;
708 np0[0] = np0[1] = np0[2] = 0;
709 np1[0] = np1[1] = np1[2] = 0;
713 o[0] = 0.5*( p0->
c[0] + p1->
c[0] );
714 o[1] = 0.5*( p0->
c[1] + p1->
c[1] );
715 o[2] = 0.5*( p0->
c[2] + p1->
c[2] );
717 memcpy(no,v,3*
sizeof(
double));
728 ps1 = n1[0]*v[0] + n1[1]*v[1] + n1[2]*v[2];
729 ps2 = n2[0]*v[0] + n2[1]*v[1] + n2[2]*v[2];
731 if(fabs(ps1) < fabs(ps2)){
741 memcpy(np0,np1,3*
sizeof(
double));
747 ps1 = n1[0]*v[0] + n1[1]*v[1] + n1[2]*v[2];
748 ps2 = n2[0]*v[0] + n2[1]*v[1] + n2[2]*v[2];
750 if(fabs(ps1) < fabs(ps2)){
760 memcpy(np1,np0,3*
sizeof(
double));
766 ps1 = n1[0]*v[0] + n1[1]*v[1] + n1[2]*v[2];
767 ps2 = n2[0]*v[0] + n2[1]*v[1] + n2[2]*v[2];
769 if(fabs(ps1) < fabs(ps2)){
783 ps1 = n1[0]*v[0] + n1[1]*v[1] + n1[2]*v[2];
784 ps2 = n2[0]*v[0] + n2[1]*v[1] + n2[2]*v[2];
786 if(fabs(ps1) < fabs(ps2)){
799 ps1 = ux*(np0[0] + np1[0]) + uy*(np0[1] + np1[1]) + uz*(np0[2] + np1[2]);
802 bn[0] = np0[0] + np1[0] -ps1*ux;
803 bn[1] = np0[1] + np1[1] -ps1*uy;
804 bn[2] = np0[2] + np1[2] -ps1*uz;
806 dd = bn[0]*bn[0] + bn[1]*bn[1] + bn[2]*bn[2];
814 no[0] = (1.0-s)*(1.0-s)*np0[0] + 2.0*s*(1.0-s)*bn[0] + s*s*np1[0];
815 no[1] = (1.0-s)*(1.0-s)*np0[1] + 2.0*s*(1.0-s)*bn[1] + s*s*np1[1];
816 no[2] = (1.0-s)*(1.0-s)*np0[2] + 2.0*s*(1.0-s)*bn[2] + s*s*np1[2];
818 dd = no[0]*no[0] + no[1]*no[1] + no[2]*no[2];
839 b0[0] = p0->
c[0] + alpha * t0[0];
840 b0[1] = p0->
c[1] + alpha * t0[1];
841 b0[2] = p0->
c[2] + alpha * t0[2];
843 b1[0] = p1->
c[0] + alpha * t1[0];
844 b1[1] = p1->
c[1] + alpha * t1[1];
845 b1[2] = p1->
c[2] + alpha * t1[2];
847 o[0] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[0] + 3.0*s*(1.0-s)*(1.0-s)*b0[0] + \
848 3.0*s*s*(1.0-s)*b1[0] + s*s*s*p1->
c[0];
850 o[1] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[1] + 3.0*s*(1.0-s)*(1.0-s)*b0[1] + \
851 3.0*s*s*(1.0-s)*b1[1] + s*s*s*p1->
c[1];
853 o[2] = (1.0-s)*(1.0-s)*(1.0-s)*p0->
c[2] + 3.0*s*(1.0-s)*(1.0-s)*b0[2] + \
854 3.0*s*s*(1.0-s)*b1[2] + s*s*s*p1->
c[2];
879 while ( !
MG_EOK(pt) && k < ne ) {
904 for (k=1; k<=
mesh->
ne; k++) {
908 if ( k == kel )
return ne;
934 while ( !
MG_VOK(ppt) && k < np ) {
960 for (k=1; k<=
mesh->
np; k++) {
964 if ( k == kp )
return np;
979 inm = fopen(fileName,
"w");
981 fprintf(inm,
"----------> %" MMG5_PRId
" MMG5_TETRAHEDRAS <----------\n",
mesh->
ne);
982 for(k=1; k<=
mesh->
ne; k++) {
984 fprintf(inm,
"num %" MMG5_PRId
" -> %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
"\n",k,pt->
v[0],pt->
v[1],
986 fprintf(inm,
"ref,tag,xt -> %" MMG5_PRId
" %d %" MMG5_PRId
"\n",pt->
ref,pt->
tag,pt->
xt);
989 fprintf(inm,
"tag -> %d %d %d %d %d %d\n",pxt->
tag[0],pxt->
tag[1],
991 fprintf(inm,
"edg -> %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
"\n",pxt->
edg[0],pxt->
edg[1],
993 fprintf(inm,
"ftag -> %d %d %d %d\n",pxt->
ftag[0],pxt->
ftag[1],
995 fprintf(inm,
"ref -> %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
" %" MMG5_PRId
"\n",pxt->
ref[0],pxt->
ref[1],
997 fprintf(inm,
"ori -> %d \n",pxt->
ori);
1001 fprintf(inm,
"---------> END MMG5_TETRAHEDRAS <--------\n");
1024 MMG5_int *lists,
int ilists,
1025 double* hausd_ip,
double *hmin_ip,
double *hmax_ip) {
1029 double hausd, hmin, hmax;
1030 int l,k,isloc,ifac1;
1048 for ( k=0; k<ilistv; ++k ) {
1050 if ( par->
ref == pt->
ref ) {
1064 for ( k=0; k<ilistv; ++k ) {
1066 if ( par->
ref == pt->
ref ) {
1086 for ( k=0; k<ilists; ++k ) {
1088 ifac1 = lists[k] % 4;
1105 for ( k=0; k<ilists; ++k ) {
1107 ifac1 = lists[k] % 4;
1121 if ( hausd_ip ) *hausd_ip = hausd;
1122 if ( hmin_ip ) *hmin_ip = hmin;
1123 if ( hmax_ip ) *hmax_ip = hmax;
1144 double* hausd_ip,
double *hmin_ip,
double *hmax_ip) {
1149 double hausd, hmin, hmax;
1153 MMG5_int ifac1,ifac2;
1154 static int8_t mmgWarn0;
1172 ifac1 = ifac2 = 4*iel + iface;
1182 fprintf(stderr,
" ## Warning: %s: unable to take into account local"
1183 " parameters at at least 1 vertex.\n",__func__ );
1192 if ( pxt->
ref[iface]!=par->
ref )
continue;
1247 for ( k=0; k<ilistv; ++k ) {
1249 if ( par->
ref != pt->
ref )
continue;
1262 for ( k=0; k<ilistv; ++k ) {
1264 if ( par->
ref != pt->
ref )
continue;
1274 if ( hausd_ip ) *hausd_ip = hausd;
1275 if ( hmin_ip ) *hmin_ip = hmin;
1276 if ( hmax_ip ) *hmax_ip = hmax;
1297 for ( k=1; k<=
mesh->
np; k++ ) {
1302 ppt->
tag &= ~MG_NUL;
1306 for ( k=1; k<=
mesh->
ne; k++ ) {
1308 if ( !
MG_EOK(pt) )
continue;
1310 for (i=0; i<4; i++) {
1312 ppt->
tag &= ~MG_NUL;
1318 if ( !
MG_EOK(pq) )
continue;
1320 for (i=0; i<6; i++) {
1322 ppt->
tag &= ~MG_NUL;
1345 MMG5_int k,*adja,iadr,iadrv;
1348 for ( k=1 ; k <=
mesh->
ne ; k++) {
1351 if ( !
MG_EOK(pt) )
continue;
1359 if ( pt->
ref == nsd )
continue;
1364 iadr = nfac*(k-1) + 1;
1366 for ( i=0; i<nfac; ++i ) {
1373 mesh->
adja[nfac*(iadrv-1)+1+iv] = 0;
1398 fprintf(stdout,
"\n -- ONLY KEEP DOMAIN OF REF %d\n",nsd );
int MMG5_BezierTgt(double c1[3], double c2[3], double n1[3], double n2[3], double t1[3], double t2[3])
double MMG5_BezierGeod(double c1[3], double c2[3], double t1[3], double t2[3])
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)
API headers for the mmg3d library.
static const int8_t MMG5_idirinv[4][4]
void MMG3D_delPt(MMG5_pMesh mesh, MMG5_int ip)
int MMG3D_delElt(MMG5_pMesh mesh, MMG5_int iel)
static const uint8_t MMG5_idir[4][3]
idir[i]: vertices of face opposite to vertex i
static const uint8_t MMG5_iprv2[3]
static const uint8_t MMG5_inxt2[6]
void MMG5_mark_verticesAsUnused(MMG5_pMesh mesh)
#define MG_SIN_OR_NOM(tag)
int MMG5_norpts(MMG5_pMesh, MMG5_int, MMG5_int, MMG5_int, double *)
Structure to store points of a MMG mesh.
Structure to store the surface tetrahedra of a MMG mesh.