ibg2.cxx

Delaunay三角形的网格剖分程序

/*
static const int ibg2ColBackground        = 0;
static const int ibg2ColError		 = 1;
static const int ibg2ColForeground        = 2;
static const int ibg2ColDefaultNode	 = 3;
static const int ibg2ColDefaultEdge	 = 4;
static const int ibg2ColDefaultFace	 = 5;
static const int ibg2ColDefaultBCondition = 5;	/* 6,7,8 for std. BC 
static const int ibg2ColDefaultRegion	 = 9;
static const int ibg2ColDefaultMaterial	 = 9;

static const int ibg2FileModNone  	= 0;
static const int ibg2FileModSimplex	= 1;

extern	wzGrid		*ibg2gCurrent;
extern  int	 (*(ibg2FileOutput[ibg2FileModMAX]))(wzGrid *g, FILE *file);

int ibg2ToSimplex(wzGrid *gg, FILE *file);

wzGrid *ibg2OtSimplex(FILE *file);

void  ibg2OutputInit();
void  ibg2Output(wzGrid*);
*/

wzGrid*		ibg2NoGrid=0;
wzName		ibg2NoName=0;
wzPoints       ibg2gExternaePoints;  

/*
int		ibg2gFace;
int		ibg2gEdge;
int		ibg2gNode;
wzFloat	ibg2gDmin[wzPointDim];
wzFloat	ibg2gBLWait,ibg2gBLHalf,ibg2gBLObtuse,ibg2gBctg;
wzFloat	ibg2gOref,ibg2gOinv,ibg2gOref2,ibg2gOmax;
wzFloat	ibg2gBV;
wzFloat	ibg2gPref;
int   		ibg2gPrefL;
*/
static wzPoint p0;


/*

void	ibgOldOutput(wzGrid* grid, char* file)
{
  wzNameFile = file;
  ibg2Output(grid);
  fflush(stdout);
}

*/
int ibg2OtSimplex(wzGrid* gg, FILE *file);

wzGrid *ibg2OtSimplex(FILE *file)
{
  wzGrid* gg=(wzGrid*) malloc(sizeof(wzGrid));
  ibg2OtSimplex(gg,file);
  return gg;
}

#define buflen 200
int ibg2OtSimplex(wzGrid* gg, FILE *file)
{
  return 0;
  /*
int ccn,ccs,cc,cb,c,cl,cr,u,ul,ur,n,nn,rc,i,j,k;
 int points,rpoints,fpoints,cells,rcells,fcells;
 int t,cpoints,csides,num[9],gdim,xdim,fcellsguess;
 wzCellType ftyp,rtyp;
 float x,y,z;
 wzFloat xmin = big1,xmax = -big1;
 wzFloat ymin = big1,ymax = -big1;
 wzFloat zmin = big1,zmax = -big1;
 wzPoint *point;
 char buffer[buflen];
 fgets(buffer,buflen,file);
 if(!strcmp(buffer,"SimplexGrid 2.0")) return ibg2Error;
 fgets(buffer,buflen,file);
 while(buffer[0]=='#') fgets(buffer,buflen,file);
 if(!strcmp(buffer,"DIMENSION")) return ibg2Error;
 fgets(buffer,buflen,file);
 rc = sscanf(buffer,"%d",&gdim);	// the grid dimension 
 if(rc != 1) return ibg2Error;
 fgets(buffer,buflen,file);
 if(!strcmp(buffer,"NODES")) return ibg2Error;
 fgets(buffer,buflen,file);
 xdim = wzPointDim;
 rc = sscanf(buffer,"%d %d",&points,&xdim);	// the number of points 
 if(rc <= 1) return ibg2Error;

 gg->xdim()  	= wzPointDim;
 gg->gdim()      = gdim;
 gg->maxPoint 	= gg->lastPoint	= points;
 gg->Point	= (wzPoint*)malloc(sizeof(wzPoint)*(gg->maxPoint+10));
 gg->firstRegion= gg->firstCell	= gg->firstPoint	= 1;

 for(n=1;n<=points;n++){
	fgets(buffer,buflen,file);
	if(wzGridGDIM(*gg)==1){
	    	rc = sscanf(buffer,"%f",&x); y=0; z=0;
		if(rc != 1) goto readerror;
	}else if(wzGridGDIM(*gg)==2){
	    	rc = sscanf(buffer,"%f %f",&x,&y); z=0;
		if(rc != 2) goto readerror;
	}else if(wzGridGDIM(*gg)==3){
	    	rc = sscanf(buffer,"%f %f %f",&x,&y,&z);
		if(rc != 3) goto readerror;
	}
	point = &wzGridPoint(*gg,n);
	if(x<xmin)	xmin = x;
	if(y<ymin)	ymin = y;
	if(z<zmin)	zmin = z;
	if(x>xmax)	xmax = x;
	if(y>ymax)	ymax = y;
	if(z>zmax)	zmax = z;
	ibg2pX(*point)[0] = x;
	ibg2pX(*point)[1] = y;
	ibg2pX(*point)[2] = z;
 }
 wzGridXMin(*gg,0) = xmin;
 wzGridXMin(*gg,1) = ymin;
 wzGridXMin(*gg,2) = zmin;
 wzGridXMax(*gg,0) = xmax;
 wzGridXMax(*gg,1) = ymax;
 wzGridXMax(*gg,2) = zmax;


 fgets(buffer,buflen,file);
 if(!strcmp(buffer,"CELLS")) return ibg2Error;
 fgets(buffer,buflen,file);
 rc = sscanf(buffer,"%d",&rcells);	// the number of region cells 
 if(rc != 1) return ibg2Error;
 fcells = fcellsguess = rcells/2+30;
 gg->lastRegion	= rcells;
 gg->firstFace	= rcells+1;
 gg->lastFace	= gg->cells.last()	= gg->maxCell	= cells = rcells+fcells;
 gg->firstEdge	= gg->firstNode	= 1;
 gg->lastEdge	= gg->lastNode	= 0;
 if(wzGridGDIM(*gg)==1){
       	cpoints = 2*rcells + 1;
	csides = 2*rcells + 2*fcells + 1;
	rtyp = wzCellType1Line;	rpoints = 2;
	ftyp = wzCellType1Node;	fpoints = 1;
 }else if(wzGridGDIM(*gg)==2){
	cpoints = 3*rcells + 2*fcells + 1;
	csides = 3*rcells + 4*fcells + 1;
	rtyp = wzCellType2Triangle;	rpoints = 3;
	ftyp = wzCellType2Line;	fpoints = 2;
 }else if(wzGridGDIM(*gg)==3){
	cpoints = 4*rcells + 3*fcells + 1;
	csides = 4*rcells + 5*fcells + 1;
	rtyp = wzCellType3Tetrahedron;rpoints = 4;
     	ftyp = wzCellType3Triangle;	fpoints = 3;
 }		     
 gg->maxCPoint	= cpoints;
 gg->maxCSide	= csides;
 gg->CellType	= (wzCellType*)malloc((cells+1)*sizeof(wzCellType));
 gg->CellSegment = (ibg2Segment*)malloc((cells+1)*sizeof(ibg2Segment));
 gg->CPointFirst = (int*)malloc((cells+1)*sizeof(int));
 gg->CSideFirst = (int*)malloc((cells+1)*sizeof(int));
 gg->CPoint	= (int*)malloc((cpoints+1)*sizeof(int));
 gg->CSide	= (int*)malloc((csides+1)*sizeof(int));
 gg->lastCOut	= gg->maxCOut = 0;
 ccn=0; ccs=0;
 for(c=1;c<=rcells;c++){	// at first the rcells region cells: 
	fgets(buffer,buflen,file);
	rc = sscanf(buffer,"%d %d %d %d %d %d %d %d %d ",
		&num[0],&num[1],&num[2],&num[3],
		&num[4],&num[5],&num[6],&num[7],&num[8]);
	wzGridCellType(*gg,c) = rtyp;
	gg->CPointFirst[c] = ccn;
	gg->CSideFirst[c] = ccs;
	for(i=0;i<rpoints;i++){
		gg->CPoint[ccn++] = num[i];	// at first the points 
	}
	wzGridCellSegment(*gg,c) = num[i++];	// now the segment 
	for(j=0;j<rpoints;i++,j++){		// now the neighbour cells: 
		if(num[i]>0)		gg->CSide[ccs++] = num[i];
		else if(num[i]<0)	gg->CSide[ccs++] = rcells - num[i];
		else 			gg->CSide[ccs++] = 0;
	}
 }

 fgets(buffer,buflen,file);
 if(!strcmp(buffer,"FACES")) return ibg2Error;
 fgets(buffer,buflen,file);
 rc = sscanf(buffer,"%d",&fcells);	// the number of boundary face cells 
 if(rc != 1) return ibg2Error;
 gg->lastFace	= gg->cells.last()	= gg->maxCell	= cells = rcells+fcells;
 gg->firstEdge	= gg->firstNode	= 1;
 gg->lastEdge	= gg->lastNode	= 0;
 if(wzGridGDIM(*gg)==1){
       	cpoints = 2*rcells + 1;
	csides = 2*rcells + 2*fcells + 1;
 }else if(wzGridGDIM(*gg)==2){
	cpoints = 3*rcells + 2*fcells + 1;
	csides = 3*rcells + 4*fcells + 1;
 }else if(wzGridGDIM(*gg)==3){
	cpoints = 4*rcells + 3*fcells + 1;
	csides = 4*rcells + 5*fcells + 1;
 }		     
 gg->freeCPoint	= cpoints-1;
 gg->freeCSide	= csides-1;
 wzAssert(fcells<=fcellsguess);
  wenn das fehlschlaegt, reallocieren:
 gg->maxCPoint	= cpoints;
 gg->maxCSide	= csides;
 gg->CellType	= (wzCellType*)malloc((cells+1)*sizeof(wzCellType));
 gg->CellSegment = (ibg2Segment*)malloc((cells+1)*sizeof(ibg2Segment));
 gg->CPointFirst = (int*)malloc((cells+1)*sizeof(int));
 gg->CSideFirst = (int*)malloc((cells+1)*sizeof(int));
 gg->CPoint	= (int*)malloc((cpoints+1)*sizeof(int));
 gg->CSide	= (int*)malloc((csides+1)*sizeof(int));
  --------

 for(cb=1;cb<=fcells;c++,cb++){
	for(i=0;i<9;i++) num[i]=0;
	fgets(buffer,buflen,file);
	rc = sscanf(buffer,"%d %d %d %d %d %d %d %d %d ",
		&num[0],&num[1],&num[2],&num[3],
		&num[4],&num[5],&num[6],&num[7],&num[8]);
	gg->CellType[c] = ftyp;
	wzAssert(gg->CellType[c] == ftyp);
	gg->CPointFirst[c] = ccn;	// begin of point list 
	gg->CSideFirst[c] = ccs;	// begin of neighbour list 
	for(i=0;i<fpoints;i++){
	    	gg->CPoint[ccn++] = num[i];
	}
	u  = num[i++];		// segment  number 
	cl = num[i++];		// left  neighbour (cell or segment number) 
	cr = num[i++];		// right neighbour (cell or segment number) 
	for(j=0;j<fpoints;j++){
		gg->CSide[ccs++] = rcells+num[i++];
	}
	gg->CSide[ccs++] = cl;
	gg->CSide[ccs++] = cr;
	if(u==0){ // automatical generation of boundary face segment number 
	    if(cl>0) ul = wzGridCellSegment(*gg,cl);
	    else     ul = -cl;
	    if(cr>0) ur = wzGridCellSegment(*gg,cr);
	    else     ur = -cr;
	    u = ibg2DefaultFaceNr(ul,ur);
	}
	wzGridCellSegment(*gg,c) = u;
 }
 return ibg2Success;
readerror:
 wzGridFree(gg);
 return ibg2Error;
*/
}
void wzGridInit()
{wzCellTypeInit();}


wzName		ibg2OSNoName=0;
/*
wzGrid		ibg2g;
int ibg2g2find(wzFloat *xx, int c0);
*/

  /*
int wzGridNRealloc(wzGrid *g)
{
  wzAssert(0);
  return 0;
 if(g->lastPoint<g->maxPoint)	return 1;
 ibg2g = *g;
 max=3*lastPoint/2 + 100;
 if((adr=realloc(Point,max*sizeof(wzPoint)))==ibg2NULL){goto noverflow;}
 Point = (wzPoint *) adr;
 if(PointCell!=ibg2NULL){
 	if((adr=realloc(PointCell,max*sizeof(int)))==ibg2NULL)
	 	{goto noverflow;}
 	PointCell = (int*)adr;
 }
 if(wzGridHasPointData(*this)){
	if((adr=realloc(PointData,max*(NDataSize)))==ibg2NULL)
	 	{goto noverflow;}
	PointData = adr;
 }
 maxPoint=max-10;
 *g = ibg2g;
 return 1;
noverflow:
 return 0;
}
 */

static int er1,er2;
/*
int wzDelaunayGenerator::gceleFace(int c, ibg2Segment m0, ibg2Segment m1)
{
 int i,*cnn=ccn(c),b0,b1,first=1;
 wzCellType t=cct(c);
 wzPoint *n;
 b0 = ibg2DefaultFaceNr(m0,m1);
 if(m1==ibg2_OutsideRegionNr) return b0;
 for(i=0;i<wzCellTypePoints(t);i++){
      	n = &(wzGridPoint(*this,cnn[i]));
	if(ibg2pType(*n)==ibg2SRegion){
		er1++; 
		return ibg2Error;
	}else if(ibg2pType(*n)==ibg2SFace){
		if((b1=ibg2pSegment(*n))!=b0){
			if(first==0) {
			 	er2++;
				return ibg2Error;
			}
			if(b1>=ibg2SMAX){
			 	er2++;
				return ibg2Error;
			}
			b0=b1;
			first=0; continue;
		}
	}
 }
 return b0;
 return 0;
}
 */

/*
int wzDelaunayGenerator::gcellSetFace(int c, ibg2Segment m0, ibg2Segment m1)
{wzCellType t=cct(c);
 return	ibg2DefaultFaceNr(m0,m1);
}

void wzDelaunayGenerator::MaklFaces(cogeometry g)
{int c0,cl0,cs,cn,m,m1,s,i,b0;
 int *cn0,*css,*csn,*cnn;
 wzCellType t,ts,tn;
 ibg2Geom = ibg2dCogeometry(g);
 if(gdim()<1){
 	Grid.firstFace=1; Grid.lastFace=0;
	return;
 }
 er1 = er2 = 0;
 Grid.firstFace= Grid.cells.last()+1;
 cl0=Grid.cells.last();
 for(c0=1;c0<=cl0;c0++) if(ccdef(c0)){
	if((m=ccu(c0))<=0) continue;
	t=cct(c0);
	for(s=0;s<wzCellTypeSides(t);s++){
		cs=ccs(c0)[s];
		if(ccout(cs)) m1=ibg2_OutsideRegionNr; /* Au"senrand 
		else if(wzCellTypeCODIM(cct(cs))!=ibg2SRegion)	continue;
		else if(m<=(m1=ccu(cs)))		continue;
		switch(wzCellTypeSSize(ts=cct(cs),s)){
			case 1:	tn=wzCellType1Node;		break;
			case 2:	tn=wzCellType2Line;		break;
			case 3: tn=wzCellType3Triangle;	break;
			case 4: tn=wzCellType3Rectangle;	break;
		}
		wzAssert(wzCellTypePoints(tn)==wzCellTypeSSize(t,s));
		cn = createCell(tn);
 		cct(cn)=tn;
		cnn=ccn(cn); csn=ccs(cn); css=ccs(cs); cn0=ccn(c0);
		for(i=0;i<wzCellTypePoints(tn);i++){
		      	cnn[i] = cn0[wzCellTypeSPoint1(t,s,i)];
		}
		for(i=0;i<wzCellTypeSides(tn);i++) csn[i]=0;
		cnn[wzCellTypeSExtR(tn)]=c0;
		cnn[wzCellTypeSExtL(tn)]=cs;
		ccs(c0)[s]=cn;
		for(i=0;i<wzCellTypeSides(ts);i++) if(css[i]==c0) {css[i]=cn; break;}
		b0 = gceleFace(cn,m,m1);
		if(b0<0){ccu(cn)=b0=gcellSetFace(cn,m,m1);}
		else ccu(cn)=b0;
	}
 }
 Grid.lastFace = Grid.cells.last();
}
void wzDelaunayGenerator::MakeEdges(cogeometry g)
{int c0,co,c1,c2,cl,m0,ni,n0,i;
 int *cn0,*cno,*cn1,*cln,*cso,bs0,bs1,bso;
 wzCellType tl,t0,t1,t2;
 int s0,s1,s2,p1,p2,first;
 if(gdim()<2){
	Grid.firstFace=1; Grid.lastFace=0;
	return;
 }
 ibg2Geom = ibg2dCogeometry(g);
 Grid.firstEdge=Grid.cells.last()+1;
 for(c0=Grid.firstFace;c0<=Grid.lastFace;c0++){t0=cct(c0);
   for(s0=0;s0<wzCellTypeSides(t0);s0++){
	if((ccs(c0))[s0]!=0) continue;
	n0=(cn0=ccn(c0))[wzCellTypeSPoint1(t0,s0,0)];
/* search of the first (after c0) ibg2SFace-element c1: 
	co=c0; c1=cn0[wzCellTypeSExtR(t0)];
	while(wzCellTypeCODIM(t1=cct(c1))==ibg2SRegion){
		cso=ccs(c1);cno=ccn(c1);
	 	for(s1=0;s1<wzCellTypeSides(t1);s1++) if(cso[s1]==co) break;
		wzAssert(s1<wzCellTypeSides(t1));
		for(p1=0;p1<wzCellTypeSSize(t1,s1);p1++)
			if(cno[wzCellTypeSPoint1(t1,s1,p1)]==n0) break;
	    	co=c1;c1=cso[wzCellTypeSSide(t1,s1,p1)];
	}
/* search of orientation in c1 
	if(co==(cn1=ccn(c1))[wzCellTypeSExtL(t1)]){
		for(s1=0;s1<wzCellTypeSides(t1);s1++)
			if(cn1[wzCellTypeSPoint1(t1,s1,0)]==n0) break;
		wzAssert(s1<wzCellTypeSides(t1));
		co=c1;c2=cn1[wzCellTypeSExtR(t1)];
	}else if(co==cn1[wzCellTypeSExtR(t1)]){
		for(s1=0;s1<wzCellTypeSides(t1);s1++)
			if(ccn(c1)[wzCellTypeSPoint1(t1,s1,1)]==n0) break;
		wzAssert(s1<wzCellTypeSides(t1));
	      	co=c1;c2=cn1[wzCellTypeSExtL(t1)];
	}else wzAssert(0);
	wzAssert(c0!=c1);
	first=1;
	do{
/* search of the next ibg2SFace-element c2: 
    	    while(wzCellTypeCODIM(t2=cct(c2))==ibg2SRegion){