ibg2.cxx

Delaunay三角形的网格剖分程序

	    	cso=ccs(c2);cno=ccn(c2);
	 	for(s2=0;s2<wzCellTypeSides(t2);s2++) if(cso[s2]==co) break;
		wzAssert(s2<wzCellTypeSides(t2));
		for(p2=0;p2<wzCellTypeSSize(t2,s2);p2++)
			if(ccn(c2)[wzCellTypeSPoint1(t2,s2,p2)]==n0) break;
		co=c2;c2=cso[wzCellTypeSSide(t2,s2,p2)];
	    }
/* search of orientation in c2 
	    if(co==ccs(c2)[wzCellTypeSExtL(t2)]){
	    	for(s2=0;s2<wzCellTypeSides(t2);s2++)
			if(ccn(c2)[wzCellTypeSPoint1(t2,s2,0)]==n0) break;
		wzAssert(s2<wzCellTypeSides(t2));
		co=c2;c2=ccn(c2)[wzCellTypeSExtR(t2)];
	    }else if(co==ccn(c2)[wzCellTypeSExtR(t2)]){
	    	for(s2=0;s2<wzCellTypeSides(t2);s2++)
			if(ccn(c2)[wzCellTypeSPoint1(t2,s2,1)]==n0) break;
		wzAssert(s2<wzCellTypeSides(t2));
	      	co=c2;c2=ccn(c2)[wzCellTypeSExtL(t2)];
	    }else wzAssert(0);
	    if(first) {
		if(co==c0) if(ccu(c0)==ccu(c1)){/* typical - no new element 
			ccs(c0)[s0]=c1; ccs(c1)[s1]=c0; goto nexts0;
	     	}first=0;
/* creating a new ibg2SEdge-element: 
		if(gdim()==3) tl=wzCellType3Line; else tl=wzCellType2Node;
		cl  = createCell(tl);
		cln = ccn(cl); cct(cl) = tl;
		m0  = -1;
		for(i=0;i<wzCellTypePoints(tl);i++){
		      	cln[i] = ni = ccn(c0)[wzCellTypeSPoint1(t0,s0,i)];
			if(ibg2pType(wzGridPoint(*this,ni))==ibg2SEdge){
				if(m0== -1){
					m0 = ibg2pSegment(wzGridPoint(*this,ni));
				}else if(ibg2pSegment(wzGridPoint(*this,ni))!=m0){
					m0 = ibg2DefaultEdgeNr;
				}
			}
		}
		if(m0== -1) m0=ibg2DefaultEdgeNr;
		ccu(cl) = m0;
/* default space is only for one external neighbour wzCellTypeSExtF 
		//  		bs0 = wzGridCONew(&ibg2g); if(bs0==0) return;
		//  		bs1 = wzGridCONew(&ibg2g); if(bs1==0) return;
		ccs(cl)[wzCellTypeSExtF(tl)] = bs0;
		Grid.COutside[bs0] = c0; ccs(c0)[s0] = cl;
		Grid.COutside[bs1] = c1; ccs(c1)[s1] = cl;
		Grid.COutNext[bs0] = bs1;
	    }
	    if(co==c0) {Grid.COutNext[bs1] = bs0; break;}
	    //	    bso = bs1; bs1 = wzGridCONew(&ibg2g); if(bs1==0) return;
	    Grid.COutNext[bso] = bs1;
	    Grid.COutside[bs1] = co; ccs(co)[s2] = cl;
	}while(1);
nexts0:	;
    }
 }
 Grid.lastEdge=Grid.cells.last();
}

void wzDelaunayGenerator::MaklNodes(cogeometry g)
{int c0,c1,c2,cl,cl2,cs,n0;
 int *pn,bs0,bso;
 wzCellType t0,t1;
 int s0,s1,s2;
 if(gdim()<3){
	Grid.firstNode=1; Grid.lastNode=0;
	return;
 }
 ibg2Geom = ibg2dCogeometry(g);
 Grid.firstNode=Grid.cells.last()+1;
 pn=(int*)calloc(Grid.points.last()+1,sizeof(int));
 cl2=Grid.lastEdge;
 for(c0=Grid.firstEdge;c0<=cl2;c0++){t0=cct(c0);
     for(s0=0;s0<wzCellTypePoints(t0);s0++){
	n0=ccn(c0)[s0];
	if(pn[n0]==0){pn[n0]= -c0; continue;}
	if(pn[n0]< 0){
		c1 = -pn[n0];
		for(s1=0;s1<2;s1++) if(ccn(c1)[s1]==n0) break;
		wzAssert(s1<2);
		ccs(c1)[s1] = c0; ccs(c0)[s0] = c1;
		pn[n0] = c0; continue;
	}
	c1 = pn[n0]; t1=cct(c1);
	if(wzCellTypeCODIM(t1)==ibg2SEdge){
		cl = createCell(wzCellType3Node); if(cl==0) return;
		ccn(cl)[0]=n0;
		//		bs0 = wzGridCONew(&ibg2g); if(bs0==0) return;
		ccs(cl)[wzCellTypeSExtF(wzCellType3Node)] = bs0;
		Grid.COutside[bs0] = c0;
		Grid.COutNext[bs0] = bs0;
		for(s1=0;s1<2;s1++) if(ccn(c1)[s1]==n0) break;
		wzAssert(s1<2);
		c2=ccs(c1)[s1];
		for(s2=0;s2<2;s2++) if(ccn(c2)[s2]==n0) break;
		wzAssert(s2<2);
		ccs(c0)[s0]=cl; ccs(c1)[s1]=cl; ccs(c2)[s2]=cl;
		pn[n0] = cl;
	}else if(wzCellTypeCODIM(t1)==ibg2SNode){
	  //		bs0 = wzGridCONew(&ibg2g); if(bs0==0) return;
		ccs(c0)[s0]=c1;
		bso = ccs(cl)[wzCellTypeSExtF(wzCellType3Node)];
		Grid.COutside[bs0] = c0;
		Grid.COutNext[bs0] = Grid.COutNext[bso];
		Grid.COutNext[bso] = bs0;
	}else wzAssert(0);
     }
 }
 for(n0=1;n0<=Grid.points.last();n0++) if(pn[n0]<0) {
		cl = createCell(wzCellType3Node); if(cl==0) return;
		ccn(cl)[0]=n0;
		c0 = -pn[n0];
		if     (ccn(c0)[0]==n0) ccs(c0)[0]=cl;
		else if(ccn(c0)[1]==n0) ccs(c0)[1]=cl;
		else wzAssert(0);
		ccs(cl)[wzCellTypeSExtF(wzCellType3Node)]=c0;
		//		bs0 = wzGridCONew(&ibg2g); if(bs0==0) return;
		ccs(cl)[wzCellTypeSExtF(wzCellType3Node)] = bs0;
		Grid.COutside[bs0] = c0;
		Grid.COutNext[bs0] = bs0;
 }
 free(pn);
/* Test des Gitters auf Korrektheit 
 for(c0=1;c0<=Grid.cells.last();c0++) if(ccdef(c0)) {t0=cct(c0);
  	for(s0=0;s0<wzCellTypeSides(t0);s0++){
		cs=ccs(c0)[s0];
		wzAssert(cs>0);
		wzAssert(wzCellTypeCODIM(t0)<=wzCellTypeCODIM(cct(cs)));
		if(wzCellTypeCODIM(t0)==wzCellTypeCODIM(cct(cs))){
			wzAssert(ccu(c0)==ccu(cs));
		}else if(wzCellTypeCODIM(cct(cs))==ibg2SFace){
		}
	}
 }
 Grid.lastNode = Grid.cells.last();
}

*/
/*
void wzGridFunction(wzGrid *gg, int fnum,
	 	 	wzFloat	(*func)(ibg2PtObject fThis, wzPoint *x),
			ibg2PtObject	fThis)
{int nn;
 for(nn=wzGridFirstPoint(*gg);nn<=wzGridLastPoint(*gg);nn++){
	ibg2pF(wzGridPoint(*gg,nn))[fnum] = func(fThis,&wzGridPoint(*gg,nn));
 }
}
*/
//static wzFloat	ibg2gDminOld[3];


#define ibg2SaveRegionMode  2

	ibg2Geometry	ibg2Geom;	/* current geometry */

#define ibg2Found 0
#define ibg2NotFound (-1)

static wzFloat ibg2gVAbs = 1.e-18, ibg2gVRel = 0.05;
static int delaunay_on=0;
static int nosegmentb,nosegment0,nosegment1,nosegment2;
static int wzCellTypeurrc;

/* defintion of segment of a correct element 
int wzDelaunayGenerator::cRegion(int c)
{int n,nh,k,nn,m,u,*cn,i,hfirst=1,notff=0;
 wzCellType t;
 wzPoint xx,*xa;
 if(ccu(c)<0){ wzGridNFree(*this,-ccu(c)); ccu(c) = 0;}
 cn=ccn(c); t=cct(c);
 for(n=0;n<wzCellTypePoints(t);n++){
	switch(cnt(nn=cn[n])){
    	case ibg2SRegion:	if((m=cnu(nn)) <= 0){
					ccu(c) = ibg2_OutsideRegionNr;
					return ibg2Found;
    				}
				xa = cnd(nn);
				goto mfound;
	case ibg2SEdge:
                 		goto minside;
	case ibg2SNode:	        goto minside;
	}
 }
minside:
 for(k=0;k<wzPointDim;k++)  xx.X[k] = 0;
 for(n=0;n<wzCellTypePoints(t);n++){nn=cn[n];
	if((cnu(nn)) <= 0){
	   	ccu(c) = ibg2_OutsideRegionNr;
		return ibg2Found;
	}
 	for(k=0;k<wzPointDim;k++)  xx.X[k] += cnx(nn)[k];
 }
 for(k=0;k<wzPointDim;k++)  xx.X[k] /= wzCellTypePoints(t);
 ibg2iRegionOfPoint(ibg2Geom,&xx,cnd(nn));
 m  = ibg2pSegment(xx);
 xa = &xx;
mfound:
 wzAssert(m>0);
 for(n=0;n<wzCellTypePoints(t);n++){
 	u=cnu(nn=cn[n]);
	switch(cnt(nn)){
	case ibg2SRegion:
		if(u <= 0) {ccu(c) = ibg2_OutsideRegionNr;
				       		return ibg2Found;}
		if(u == m) continue;
	     	notff = nn;
	default:
	       	if(ibg2iStatusOfLine(ibg2Geom,xa,cnd(nn))!=ibg2Incorrect) continue;
	     	notff = nn;
	}
 }
 if(notff==0){
	ccu(c)=m;
	return ibg2Found;
 }
 if(&xx == xa && hfirst){
	hfirst = 0;
  	/* may be a point near the dangerous corner will be better? 
	for(i=0;i<10;i++){
	  for(k=0;k<wzPointDim;k++)  xx.X[k] = 0.5*(xx.X[k]+ibg2pX(*cnd(nn))[k]);
	  ibg2iRegionOfPoint(ibg2Geom,&xx,cnd(nn));
	  if(ibg2pSegment(xx)!=m){
		m = ibg2pSegment(xx);
		goto mfound;
	  }
	}
 }
 wzGridNGet(*this,nh);
 *cnd(nh) = *xa;
 ccu(c) = -nh;
 return ibg2NotFound;
}
*/
/* definition of segment for an incorrect element 
void wzDelaunayGenerator::cSetRegion(int c)
{int n = -ccu(c);
 wzAssert(n>0);
 wzAssert(n<=Grid.points.last());
 ccu(c) = ibg2pSegment(wzGridPoint(*this,n));
 wzGridNFree(*this,n);
}
*/
static int gmlev;
/*
void wzDelaunayGenerator::g2segment(int c)
{int i,j,co,*cn=ccn(c);
 wzCellType t=cct(c);
 wzAssert(t==wzCellType2Triangle);
/* Korrektur des falschen Dreiecks: 
 for(i=0;i<3;i++){
	if(ibg2Incorrect != ibg2iStatusOfLine(ibg2Geom,
			cnd(cn[wzCellTypeSPoint1(t,i,0)]),
 			cnd(cn[wzCellTypeSPoint1(t,i,1)]))) continue;
	co=ccs(c)[i]; if(ccout(co)) continue;
	if(gmlev==0){
		for(j=0;j<3;j++) if(ccs(co)[j]==c) break;
		if(ibg2Incorrect != ibg2iStatusOfLine(ibg2Geom,cnd(cn[i]),cnd(ccn(co)[j])))
				   		if(g2swapm(c,i)) return;
	}else{
		gmlev--;
		if(g2swapm(c,i)) {gmlev++; return;}
		gmlev++;
	}
 }
}
void wzDelaunayGenerator::g2swapd(int c)
{int *cn,*cno,*cs,*cso,*csb,nn,n0,n1,no,co,b,bo,i,i0,i1,j,j0,j1,k,m;
 wzFloat x0,x1,xn,xo,y0,y1,yn,yo,sinn,sino,cosn,coso;
 wzCellType t=cct(c);
 m=ccu(c); cs=ccs(c); cn=ccn(c);
 for(i=0;i<3;i++){
	if(m!=ccu((co=cs[i])))	continue;
 	cso=ccs(co); cno=ccn(co);
	if(cct(co)!=wzCellType2Triangle)	continue;
	for(j=0;j<3;j++) if(cso[j]==c) goto cfound;
	wzAssert(0); continue;
cfound:
	i0=wzCellTypeSPoint1(t,i,0); i1=wzCellTypeSPoint1(t,i,1);
	nn=cn[i]; n0=cn[i0]; n1=cn[i1]; no=cno[j];
	x0=cnx(n0)[0]; y0=cnx(n0)[1];
	x1=cnx(n1)[0]; y1=cnx(n1)[1];
	xn=cnx(nn)[0]; yn=cnx(nn)[1];
	xo=cnx(no)[0]; yo=cnx(no)[1];
	sinn  =	((x0-xn)*(y1-yn) - (x1-xn)*(y0-yn));	if(sinn<0) continue;
	sino  =	((x1-xo)*(y0-yo) - (x0-xo)*(y1-yo));	if(sino<0) continue;
	cosn =	((x0-xn)*(x1-xn) + (y0-yn)*(y1-yn));
	coso =	((x0-xo)*(x1-xo) + (y0-yo)*(y1-yo));
	if(cosn*sino+coso*sinn>-eps2) continue;
	j0=wzCellTypeSPoint1(t,j,0); j1=wzCellTypeSPoint1(t,j,1);
	b=cs[i0]; bo=cso[j0];
	cn[i1]=no; cno[j1]=nn;
	cs[i0]=co; cso[j0]=c;
	cs[i] =bo; cso[j] =b;
	if(b > 0){
		m=wzCellTypePoints(cct(b ));csb=ccs(b );
		for(k=0;k<m;k++) if(csb[k]==c )	{csb[k] =co; break;}
	}
	if(bo> 0){
		m=wzCellTypePoints(cct(bo));csb=ccs(bo);
		for(k=0;k<m;k++) if(csb[k]==co)	{csb[k] = c; break;}
  	}
	g2swapd(c);
	g2swapd(co);
	return;
 }
}
int wzDelaunayGenerator::g2swapm(int c, int i)
{int *cn,*cno,*cs,*cso,*csb,nn,n0,n1,no,co,b,bo,i0,i1,j,j0,j1,k,m;
 wzFloat x,y,d;
 wzCellType t=cct(c);
 cn=ccn(c); nn=cn[i];
 cs=ccs(c); co=cs[i]; cso=ccs(co); cno=ccn(co);
 if(ccout(co)) return 0;
 if(cct(co)!=wzCellType2Triangle) wzAssert(0);
 for(j=0;j<3;j++) if(cso[j]==c) goto cfound;
 wzAssert(0); return 0;
cfound:
 i0=wzCellTypeSPoint1(t,i,0); i1=wzCellTypeSPoint1(t,i,1);
 n0=cn[i0]; n1=cn[i1]; no=cno[j];
 x=cnx(n0)[0]; y=cnx(n0)[1];
 d=(cnx(no)[0]-x)*(cnx(nn)[1]-y)-(cnx(nn)[0]-x)*(cnx(no)[1]-y);
 if(d<eps2) return 0;
 x=cnx(n1)[0]; y=cnx(n1)[1];
 d=(cnx(no)[0]-x)*(cnx(nn)[1]-y)-(cnx(nn)[0]-x)*(cnx(no)[1]-y);
 if(d>-eps2) return 0;
 j0=wzCellTypeSPoint1(t,j,0); j1=wzCellTypeSPoint1(t,j,1);
 b=cs[i0]; bo=cso[j0];
 cn[i1]=no; cno[j1]=nn;
 cs[i0]=co; cso[j0]=c;
 cs[i] =bo; cso[j] =b;
 if(b > 0){
	m=wzCellTypePoints(cct(b ));csb=ccs(b );
	for(k=0;k<m;k++) if(csb[k]==c )	{csb[k] =co; break;}
 }
 if(bo> 0){
	m=wzCellTypePoints(cct(bo));csb=ccs(bo);
	for(k=0;k<m;k++) if(csb[k]==co)	{csb[k] = c; break;}
 }
 if(cRegion(c)==ibg2NotFound)	g2segment(c);
 else				g2swapd(c);
 if(cRegion(co)==ibg2NotFound)  g2segment(co);
 else				g2swapd(co);
 return co;
}
/* Berechnen und Eintragen der Daten der Umkugel:
static float g3sdelmin=0.1, g3sdelmax=0.9;

int wzDelaunayGenerator::g3cvol(int c)
{int *cn0;
 wzFloat *cx0,*cx1,*cx2,*cx3;
 wzFloat x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,l1,l2,l3,a1,a2,a3,mx,my,mz,det;
 cn0=ccn(c);
 cx0=cnx(*cn0++); cx1=cnx(*cn0++); cx2=cnx(*cn0++); cx3=cnx(*cn0++);
 x  = cx0[0];   y  = cx0[1];   z  = cx0[2];
 x1 = cx1[0]-x; y1 = cx1[1]-y; z1 = cx1[2]-z;
 x2 = cx2[0]-x; y2 = cx2[1]-y; z2 = cx2[2]-z;
 x3 = cx3[0]-x; y3 = cx3[1]-y; z3 = cx3[2]-z;
 l1=x1*x1 + y1*y1 + z1*z1;
 l2=x2*x2 + y2*y2 + z2*z2;
 l3=x3*x3 + y3*y3 + z3*z3;
 a1=y2*z3 - y3*z2;
 a2=y3*z1 - y1*z3;
 a3=y1*z2 - y2*z1;
 det=x1*a1 + x2*a2 + x3*a3;
 if(det>-ibgDelaunayEps3) return 0; /* Entartung 
 mx=(l1*a1+l2*a2+l3*a3)/det;
 my=(l1*(z2*x3-z3*x2)+l2*(z3*x1-z1*x3)+l3*(z1*x2-z2*x1))
 		/det;
 mz=(l1*(x2*y3-x3*y2)+l2*(x3*y1-x1*y3)+l3*(x1*y2-x2*y1))
 		/det;
 ccv(c).x1[0]=x;
 ccv(c).x1[1]=y;
 ccv(c).x1[2]=z;
 ccv(c).x2[0]=x+mx;
 ccv(c).x2[1]=y+my;
 ccv(c).x2[2]=z+mz;
 return 1;
}

static int g3segment1line;

int wzDelaunayGenerator::g3CHalf(int c)
{
 wzCellType 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;