📄 ibg2.cxx
字号:
ibg2CellType t=cct(c);
int s1,s2,sr,sl,so,c1,co,cu,ce,cn,n1,n2,nh,e,i,rc;
int *cn0=ccn(c);
wzFloat *cxx;
wzFloat x,y,z,xh,yh,zh,x2,y2,z2,d;
if(t!=wzCellType3Tetrahedron) return 0;
for(e=0;e<wzCellTypeLines(t);e++){
if(ibg2Incorrect != ibg2iStatusOfLine(ibg2Geom,cnd(cn0[wzCellTypeEPoint1(t,e)])
,cnd(cn0[wzCellTypeEPoint2(t,e)]))) continue;
n1=cn0[s1=wzCellTypeEPoint1(t,e)];
n2=cn0[s2=wzCellTypeEPoint2(t,e)];
if(cnt(n1)!=ibg2SRegion) continue;
if(cnt(n2)!=ibg2SRegion) continue;
sr=wzCellTypeESideR(t,e);
sl=wzCellTypeESideL(t,e);
goto newpoint;
}
for(e=0;e<wzCellTypeLines(t);e++){
if(ibg2Incorrect != ibg2iStatusOfLine(ibg2Geom,cnd(cn0[wzCellTypeEPoint1(t,e)])
,cnd(cn0[wzCellTypeEPoint2(t,e)]))) continue;
n1=cn0[s1=wzCellTypeEPoint1(t,e)];
n2=cn0[s2=wzCellTypeEPoint2(t,e)];
if(cnt(n1)!=ibg2SRegion) {
if(cnt(n2)!=ibg2SRegion) continue;
nh=n1; n1=n2; n2=nh;
sh=s1; s1=s2; s2=sh;
sr=wzCellTypeESideL(t,e);
sl=wzCellTypeESideR(t,e);
} else {sr=wzCellTypeESideR(t,e);
sl=wzCellTypeESideL(t,e);
}
goto newpoint;
}
for(e=0;e<wzCellTypeLines(t);e++){
if(ibg2Incorrect != ibg2iStatusOfLine(ibg2Geom,cnd(cn0[wzCellTypeEPoint1(t,e)])
,cnd(cn0[wzCellTypeEPoint2(t,e)]))) continue;
n1=cn0[s1=wzCellTypeEPoint1(t,e)];
n2=cn0[s2=wzCellTypeEPoint2(t,e)];
wzAssert(0);
sr=wzCellTypeESideR(t,e);
sl=wzCellTypeESideL(t,e);
goto newpoint;
}
wzAssert(0);
return 0;
newpoint:
wzGridNGet(*this,nh);
rc = ibg2iFaceWithLine(ibg2Geom,cnd(nh),cnd(n1),cnd(n2));
if(rc != ibg2FaceFound) {
wzAssert(0);
return 0;
}
cxx = cnx(n1); x = cxx[0]; y = cxx[1]; z = cxx[2];
cxx = cnx(nh); xh = cxx[0]-x; yh = cxx[1]-y; zh = cxx[2]-z;
cxx = cnx(n2); x2 = cxx[0]-x; y2 = cxx[1]-y; z2 = cxx[2]-z;
if((d=(xh*x2+yh*y2+zh*z2)/(x2*x2+y2*y2+z2*z2))<g3sdelmin){
cxx = cnx(n1);
x = cxx[0]; y = cxx[1]; z = cxx[2];
*cnd(n1) = *cnd(nh);
cxx[0] = x; cxx[1] = y; cxx[2] = z;
wzGridNFree(*this,nh);
if(cRegion(c)==ibg2NotFound) return 0;
return 1;
}else if(d>g3sdelmax){
cxx = cnx(n2);
x = cxx[0]; y = cxx[1]; z = cxx[2];
*cnd(n2) = *cnd(nh);
cxx[0] = x; cxx[1] = y; cxx[2] = z;
wzGridNFree(*this,nh);
if(cRegion(c)==ibg2NotFound) return 0;
return 1;
}
c1 = c; co = 0;
// cr = ccs(c1)[sr];
do{
g0CGet(&cn,wzCellType3Tetrahedron);
if(cn<=0) {wzAssert(0); return 0;}
ccn(cn)[s1]=nh;
ccn(cn)[s2]=n2;
ccn(cn)[sr]=ccn(c1)[sr];
ccn(cn)[sl]=ccn(c1)[sl];
ccn(c1)[s2]=nh;
ccs(cn)[s1]=ce=ccs(c1)[s1];
for(i=0;i<wzCellTypeSides(cct(ce));i++)
if(ccs(ce)[i]==c1) {ccs(ce)[i]=cn; break;}
ccs(cn)[s2]=c1;
if(co>0){
ccs(cn)[sl]=co;
ccs(co)[so]=cn;
}
ccs(c1)[s1]=cn;
if(cRegion(c1)==ibg2NotFound){
if(st[c1]==0){st[c1] = nosegment0; nosegment0 = c1;}
}
if(cRegion(cn)==ibg2NotFound){
if(st[cn]==0){ st[cn] = nosegment0; nosegment0 = cn;}
}
g3cvol(c1); g3cvol(cn);
co = cn; so = sr; cu = c1; c1 = ccs(c1)[sr];
s1=s2=sr=sl= -1;
for(i=0;i<4;i++){
if (ccn(c1)[i]==n1) s1=i;
else if(ccn(c1)[i]==n2) s2=i;
else if(ccs(c1)[i]==cu) sl=i;
else sr=i;
}
wzAssert(s1>=0); if(c1!=c) wzAssert(s2>=0); wzAssert(sr>=0); wzAssert(sl>=0);
}while(c1!=c);
cn = ccs(c1)[s1];
ccs(cn)[sl]=co;
ccs(co)[so]=cn;
return 1;
}
int wzDelaunayGenerator::g3SDelete(int cb, int sb)
{ibg2CellType tb=cct(cb),ts;
int cs,cn,co,ss,ss0,ss1,ss2,sb0,sb1,sb2,sh,nb,n0,n1,n2,ns,nh,i;
wzFloat *cxx;
wzFloat x,y,z,x0,y0,z0,x1,y1,z1,x2,y2,z2,xs,ys,zs,det,det0,vol0;
cs=ccs(cb)[sb]; ts=cct(cs);
if(ccout(cs)) return 0;
for(ss=0;ss<wzCellTypeSides(ts);ss++) if(ccs(cs)[ss]==cb) break;
nb = ccn(cb)[sb];
ns = ccn(cs)[ss];
n0 = ccn(cb)[sb0=wzCellTypeSPoint1(tb,sb,0)];
n1 = ccn(cb)[sb1=wzCellTypeSPoint1(tb,sb,1)];
n2 = ccn(cb)[sb2=wzCellTypeSPoint1(tb,sb,2)];
if(ibg2Incorrect == ibg2iStatusOfLine(ibg2Geom,cnd(nb),cnd(ns))) return 0;
cxx = cnx(nb); x = cxx[0]; y = cxx[1]; z = cxx[2];
if(delaunay_on){ /* Umkugeltest:
if((ccv(cs).x1[0]-x)*(ccv(cs).x2[0]-x)
+ (ccv(cs).x1[1]-y)*(ccv(cs).x2[1]-y)
+ (ccv(cs).x1[2]-z)*(ccv(cs).x2[2]-z) > -eps2) return 0;
}
cxx = cnx(ns); xs = cxx[0]-x; ys = cxx[1]-y; zs = cxx[2]-z;
cxx = cnx(n0); x0 = cxx[0]-x; y0 = cxx[1]-y; z0 = cxx[2]-z;
cxx = cnx(n1); x1 = cxx[0]-x; y1 = cxx[1]-y; z1 = cxx[2]-z;
cxx = cnx(n2); x2 = cxx[0]-x; y2 = cxx[1]-y; z2 = cxx[2]-z;
det0= x0*(y1*z2-z1*y2) + y0*(z1*x2-x1*z2) + z0*(x1*y2-y1*x2);
vol0= ibg2gVRel*det0 - ibg2gVAbs;
det = x0*(y1*zs-z1*ys) + y0*(z1*xs-x1*zs) + z0*(x1*ys-y1*xs);
if(det > vol0){
if(g3segment1line++ >= 1) return 0;
if(g3EDelete(cb,wzCellTypeSLine(tb,sb,0))) return 2;
return 0;
}
det = x1*(y2*zs-z2*ys) + y1*(z2*xs-x2*zs) + z1*(x2*ys-y2*xs);
if(det > vol0){
if(g3segment1line++ >= 1) return 0;
if(g3EDelete(cb,wzCellTypeSLine(tb,sb,1))) return 2;
return 0;
}
det = x2*(y0*zs-z0*ys) + y2*(z0*xs-x0*zs) + z2*(x0*ys-y0*xs);
if(det > vol0){
if(g3segment1line++ >= 1) return 0;
if(g3EDelete(cb,wzCellTypeSLine(tb,sb,2))) return 2;
return 0;
}
for(i=0;i<wzCellTypeSSize(ts,ss);i++){
if(n0==(nh=ccn(cs)[sh=wzCellTypeSPoint1(ts,ss,i)])) ss0=sh;
else if(n1==nh) ss1=sh;
else {wzAssert(n2==nh); ss2=sh;}
}
g0CGet(&cn,wzCellType3Tetrahedron); if(cn<0) return 0;
ccn(cn)[0] = nb;
ccn(cn)[1] = n1;
ccn(cn)[2] = n2;
ccn(cn)[3] = ns;
ccs(cn)[1] = cs;
ccs(cn)[2] = cb;
ccn(cb)[sb2] = ns;
ccn(cs)[ss1] = nb;
co=ccs(cb)[sb1]; ccs(cb)[sb1]=cs; ccs(cs)[ss]= co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cb) {ccs(co)[i]=cs;break;}
co=ccs(cb)[sb0]; ccs(cb)[sb0]=cn; ccs(cn)[3] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cb) {ccs(co)[i]=cn;break;}
co=ccs(cs)[ss2]; ccs(cs)[ss2]=cb; ccs(cb)[sb]= co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cs) {ccs(co)[i]=cb;break;}
co=ccs(cs)[ss0]; ccs(cs)[ss0]=cn; ccs(cn)[0] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cs) {ccs(co)[i]=cn;break;}
if(cRegion(cb)==ibg2NotFound)if(st[cb]==0){st[cb] = nosegment0; nosegment0 = cb;}
if(cRegion(cs)==ibg2NotFound)if(st[cs]==0){st[cs] = nosegment0; nosegment0 = cs;}
if(cRegion(cn)==ibg2NotFound){
if(st[cn]==0){st[cn] = nosegment0; nosegment0 = cn;}
}
g3cvol(cb); g3cvol(cs); g3cvol(cn);
return 1;
}
int wzDelaunayGenerator::g34Swap(int c1, int c2, int c3, int c4)
{int i,co,s12,s13,s1u,s1o,s21,s2u,s2o,s34,s3u,s3o,s43,s4u,s4o;
int n1,n2,n3,n4,nu,no;
wzFloat *cxx;
wzFloat x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4,xu,yu,zu,xo,yo,zo,det,vol0;
ibg2CellType t1=cct(c1);
for(i=0;i<4;i++){
if(ccs(c1)[i]==c2) {n1=ccn(c1)[i];s12=i;}
if(ccs(c2)[i]==c1) {n2=ccn(c2)[i];}
if(ccs(c1)[i]==c3) {n3=ccn(c1)[i];s13=i;}
if(ccs(c3)[i]==c1) {n4=ccn(c3)[i];}
}
if(ibg2Incorrect == ibg2iStatusOfLine(ibg2Geom,cnd(n1),cnd(n2))) return 0;
for(i=0;i<wzCellTypeSSize(t1,s12);i++){
if(wzCellTypeSSide(t1,s12,i)==s13){
nu = ccn(c1)[wzCellTypeSPoint1(t1,s12,i)];
no = ccn(c1)[wzCellTypeSPoint1(t1,s12,i+1)];
break;
}
}
cxx = cnx(n1); x = cxx[0]; y = cxx[1]; z = cxx[2];
cxx = cnx(n2); x2 = cxx[0]-x; y2 = cxx[1]-y; z2 = cxx[2]-z;
cxx = cnx(nu); xu = cxx[0]-x; yu = cxx[1]-y; zu = cxx[2]-z;
cxx = cnx(n4); x4 = cxx[0]-x; y4 = cxx[1]-y; z4 = cxx[2]-z;
cxx = cnx(no); xo = cxx[0]-x; yo = cxx[1]-y; zo = cxx[2]-z;
cxx = cnx(n3); x3 = cxx[0]-x; y3 = cxx[1]-y; z3 = cxx[2]-z;
vol0= -ibg2gVAbs;
det = x2*(yu*z4-zu*y4) + y2*(zu*x4-xu*z4) + z2*(xu*y4-yu*x4);
if(det > vol0) return 0;
det = x2*(y4*zo-z4*yo) + y2*(z4*xo-x4*zo) + z2*(x4*yo-y4*xo);
if(det > vol0) return 0;
det = x2*(yo*z3-zo*y3) + y2*(zo*x3-xo*z3) + z2*(xo*y3-yo*x3);
if(det > vol0) return 0;
det = x2*(y3*zu-z3*yu) + y2*(z3*xu-x3*zu) + z2*(x3*yu-y3*xu);
if(det > vol0) return 0;
for(i=0;i<4;i++){
if(ccn(c1)[i]==nu) s1u=i;
if(ccn(c1)[i]==no) s1o=i;
if(ccn(c2)[i]==n2) s21=i;
// if(ccn(c2)[i]==n3) s24=i;
if(ccn(c2)[i]==nu) s2u=i;
if(ccn(c2)[i]==no) s2o=i;
if(ccn(c3)[i]==n1) s34=i;
// if(ccn(c3)[i]==n4) s31=i;
if(ccn(c3)[i]==nu) s3u=i;
if(ccn(c3)[i]==no) s3o=i;
if(ccn(c4)[i]==n2) s43=i;
// if(ccn(c4)[i]==n4) s42=i;
if(ccn(c4)[i]==nu) s4u=i;
if(ccn(c4)[i]==no) s4o=i;
}
ccn(c1)[s1o] = n2;
ccn(c2)[s2u] = n1;
ccn(c3)[s3o] = n2;
ccn(c4)[s4u] = n1;
co=ccs(c1)[s1u]; ccs(c1)[s1u] = c2; ccs(c2)[s21] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==c1) {ccs(co)[i]=c2;break;}
co=ccs(c2)[s2o]; ccs(c2)[s2o] = c1; ccs(c1)[s12] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==c2) {ccs(co)[i]=c1;break;}
co=ccs(c3)[s3u]; ccs(c3)[s3u] = c4; ccs(c4)[s43] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==c3) {ccs(co)[i]=c4;break;}
co=ccs(c4)[s4o]; ccs(c4)[s4o] = c3; ccs(c3)[s34] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==c4) {ccs(co)[i]=c3;break;}
if(cRegion(c1)==ibg2NotFound)if(st[c1]==0){st[c1] = nosegment0; nosegment0 = c1;}
if(cRegion(c2)==ibg2NotFound)if(st[c2]==0){st[c2] = nosegment0; nosegment0 = c2;}
if(cRegion(c3)==ibg2NotFound)if(st[c3]==0){st[c3] = nosegment0; nosegment0 = c3;}
if(cRegion(c4)==ibg2NotFound)if(st[c4]==0){st[c4] = nosegment0; nosegment0 = c4;}
g3cvol(c1); g3cvol(c2); g3cvol(c3); g3cvol(c4);
return 1;
}
int wzDelaunayGenerator::g3EDelete(int cb, int eb)
{ibg2CellType tb=cct(cb),tl,tr;
int cl,cr,co,sbl,sbr,sbo,sbu,slb,slr,slo,slu,sl4,sr4,srb,sro,sru;
int nu,no,nl,nr,ns,i;
wzFloat *cxx;
wzFloat x,y,z,xl,yl,zl,xr,yr,zr,xo,yo,zo,xu,yu,zu,det1,det2,vol0;
int rc=0;
int c4=0;
beg:
cl=ccs(cb)[sbl=wzCellTypeESideL(tb,eb)]; if(ccdef(cl)) return 0; tl=cct(cl);
cr=ccs(cb)[sbr=wzCellTypeESideR(tb,eb)]; if(ccdef(cr)) return 0; tr=cct(cr);
nu=ccn(cb)[sbu=wzCellTypeEPoint1(tb,eb)];
no=ccn(cb)[sbo=wzCellTypeEPoint2(tb,eb)];
nl=ccn(cb)[sbl];
nr=ccn(cb)[sbr];
for(i=0;i<wzCellTypeSides(tl);i++){
if(ccs(cl)[i]==cb) {slb=i; ns=ccn(cl)[slb];}
else if(ccn(cl)[i]==nr) {sl4=i;
if(ccs(cl)[i]!=cr) {
/* four or more tetrahedrons aroung the given line:
for(srb=0;srb<4;srb++) if(ccs(cr)[srb]==cb) break;
if(g3SDelete(cb,sbl)){
cb=ccs(cr)[srb];
for(eb=0;eb<6;eb++){
if(ccn(cb)[wzCellTypeEPoint1(tb,eb)]==nu)
if(ccn(cb)[wzCellTypeEPoint2(tb,eb)]==no)
break;
if(ccn(cb)[wzCellTypeEPoint1(tb,eb)]==no)
if(ccn(cb)[wzCellTypeEPoint2(tb,eb)]==nu)
break;
}
wzAssert(eb<6);
rc=1;
goto beg;
}
for(slb=0;slb<4;slb++) if(ccs(cl)[slb]==cb) break;
if(g3SDelete(cb,sbr)){
cb=ccs(cl)[slb];
for(eb=0;eb<6;eb++){
if(ccn(cb)[wzCellTypeEPoint1(tb,eb)]==nu)
if(ccn(cb)[wzCellTypeEPoint2(tb,eb)]==no)
break;
if(ccn(cb)[wzCellTypeEPoint1(tb,eb)]==no)
if(ccn(cb)[wzCellTypeEPoint2(tb,eb)]==nu)
break;
}
wzAssert(eb<6);
rc=1;
goto beg;
}
/* special case: line with 4 Tetrahedrons
c4=ccs(cl)[sl4];
for(sr4=0;sr4<wzCellTypeSides(tr);sr4++){
if(ccn(cr)[sr4]==nl) {
if(ccs(cr)[sr4]!=c4){
/* five or more tetrahedrons around the line:
if(g3SDelete(cr,sr4)) return 2;
if(g3SDelete(cl,sl4)) return 2;
return rc;
}
break;
}
}
if(ccout(cl)){
if(ccdef(c4)) return rc;
if(ccout(cr)) return rc;
if(g34Swap(cb,cr,cl,c4)) return 1;
return rc;
}
if(ccout(cr)){
if(ccdef(c4)) return rc;
if(ccout(cl)) return rc;
if(g34Swap(cb,cl,cr,c4)) return 1;
return rc;
}
if(ccout(c4)) return rc;
if(g34Swap(cb,cl,cr,c4)) return 1;
if(g34Swap(cr,cb,c4,cl)) return 1;
return rc;
}
slr=i;
}
else if(ccn(cl)[i]==nu) {slu=i;}
else if(ccn(cl)[i]==no) {slo=i;}
else wzAssert(0);
}
if(ccout(cl)) return 0;
if(ccout(cr)) return 0;
cxx = cnx(ns); x = cxx[0]; y = cxx[1]; z = cxx[2];
cxx = cnx(nu); xu = cxx[0]-x; yu = cxx[1]-y; zu = cxx[2]-z;
cxx = cnx(nr); xr = cxx[0]-x; yr = cxx[1]-y; zr = cxx[2]-z;
cxx = cnx(no); xo = cxx[0]-x; yo = cxx[1]-y; zo = cxx[2]-z;
cxx = cnx(nl); xl = cxx[0]-x; yl = cxx[1]-y; zl = cxx[2]-z;
vol0= ibg2gVAbs;
det1 = xl*(yo*zr-zo*yr) + yl*(zo*xr-xo*zr) + zl*(xo*yr-yo*xr);
det2 = xr*(yu*zl-zu*yl) + yr*(zu*xl-xu*zl) + zr*(xu*yl-yu*xl);
if(det1 < ibg2gVRel*det2 + vol0) return rc;
if(det2 < ibg2gVRel*det1 + vol0) return rc;
for(i=0;i<wzCellTypeSides(tr);i++){
if(ccs(cr)[i]==cb) {srb=i; wzAssert(ns==ccn(cr)[srb]);}
else if(ccn(cr)[i]==nl) {wzAssert(ccs(cr)[i]==cl);}// srl=i;}
else if(ccn(cr)[i]==nu) {sru=i;}
else if(ccn(cr)[i]==no) {sro=i;}
else wzAssert(0);
}
ccn(cb)[sbo] = ns;
ccn(cl)[slu] = nl;
co=ccs(cb)[sbu]; ccs(cb)[sbu]=cl; ccs(cl)[slb] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cb) {ccs(co)[i]=cl;break;}
co=ccs(cl)[slo]; ccs(cl)[slo]=cb; ccs(cb)[sbl] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cl) {ccs(co)[i]=cb;break;}
co=ccs(cr)[sru]; ccs(cl)[slr] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cr) {ccs(co)[i]=cl;break;}
co=ccs(cr)[sro]; ccs(cb)[sbr] = co;
for(i=0;i<wzCellTypeSides(cct(co));i++) if(ccs(co)[i]==cr) {ccs(co)[i]=cb;break;}
if(cRegion(cb)==ibg2NotFound)if(st[cb]==0){st[cb] = nosegment0; nosegment0 = cb;}
if(cRegion(cl)==ibg2NotFound)if(st[cl]==0){st[cl] = nosegment0; nosegment0 = cl;}
g3cvol(cb); g3cvol(cl);
if(ccu(cr)<0) wzGridNFree(*this,-ccu(cr));
destroyCell(cr);
return 1;
}
int wzDelaunayGenerator::g3segment1(int c)
{ibg2CellType t=cct(c);
int rc,n,n0,e, *cn0=ccn(c);
if(t!=wzCellType3Tetrahedron) return 0;
if(cRegion(c)==ibg2Found) return 1;
g3segment1line = 1;
for(e=0;e<wzCellTypeLines(t);e++){
if(ibg2Incorrect != ibg2iStatusOfLine(ibg2Geom,cnd(cn0[wzCellTypeEPoint1(t,e)])
,cnd(cn0[wzCellTypeEPoint2(t,e)]))) continue;
if(g3EDelete(c,e)) return 1;
}
g3segment1line = 0;
n0 = -1;
for(n=0;n<wzCellTypePoints(t);n++){
if(ibg2pType(*cnd(cn0[n]))==ibg2SEdge){
if(n0==-1) n0=n;
else goto next;
}
if(ibg2pType(*cnd(cn0[n]))==ibg2SNode){
n0 = n; goto vert;
}
}
if(n0>=0){
rc = g3SDelete(c,n0);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -