ibg2.cxx

Delaunay三角形的网格剖分程序

	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)
{wzCellType 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;
 wzCellType 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)
{wzCellType 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)
{wzCellType 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);
	if(rc==1) return 1;
    	if(rc==2) return 2;
 }
vert:
next:
 return 0;
}

int wzDelaunayGenerator::g3segment2(int c)
{wzCellType t=cct(c);
 if(t!=wzCellType3Tetrahedron) return 0;
 if(g3segment1(c)) return 1;
 if(g3CHalf(c)) return 1;
 return 0;
}

void wzDelaunayGenerator::MakeRegions(cogeometry g)
{
  int c,u;
 ibg2Geom = ibg2dCogeometry(g);
 ibg2tCopy(&(Grid.Topology),&(ibg2Geom->top));
 // unc=0; mec=0; mbc=0; mrc=0;
 nosegmentb = -1;
 nosegment0 = -1;
 nosegment1 = -1;
 nosegment2 = -1;
 wzCellTypeurrc   = 0;
 for(c=1;c<=Grid.cells.last();c++) st[c]=0;
 for(c=1;c<=Grid.cells.last();c++){
	if(ccund(c)) continue;
	ccu(c) = 0;
	if(cRegion(c)==ibg2NotFound){
	    st[c] = nosegmentb; nosegmentb = c; wzCellTypeurrc++;
	}
 }
 if(gdim()==2){
        nosegment0 = nosegmentb;
	do{
		while(nosegment0>0){
		      	c = nosegment0; nosegment0 = st[c]; st[c] = 0;
			if(cnosegment(ccu(c))){
			        g2segment(c);
			        if(ccu(c)<0){wzAssert(st[c]==0);
			        	st[c] = nosegment2; nosegment2 = c;
			        }
			}
		}
		if(nosegment1>0){