ibg2.cxx

Delaunay三角形的网格剖分程序

			continue;
		}
		break;
	}while(1);
 }else if(gdim()==3){
	while(nosegmentb>0){
                c = nosegmentb; nosegmentb = st[c]; st[c] = 0;
	        if(ccund(c)) continue;
	        while(cnosegment(ccu(c))){
	            u = g3segment1(c);
	            if(st[c]!=0) break;
	            if(u==0){
	        	st[c] = nosegment0; nosegment0 = c;
	                break;
	            }
	        }
	}
	do{
		while(nosegment0>0){
		    if(wzCellTypeurrc<0) break;
		    c = nosegment0; nosegment0 = st[c]; st[c] = 0;
		    if(ccund(c)) continue;
		    while(cnosegment(ccu(c))){
			wzCellTypeurrc--;
			u = g3segment1(c);
			if(st[c]!=0) break;
			if(u==0){
			    st[c] = nosegment1; nosegment1 = c;
			    break;
			}
		    }
		}
		if(nosegment1>0){
	  	        c = nosegment1; nosegment1 = st[c]; st[c] = 0;
	   	        if(ccund(c)) continue;
		        while(cnosegment(ccu(c))){
		        	u = g3segment2(c);
		        	if(st[c]!=0) break;
		             	if(u==0){
		        		st[c] = nosegment2; nosegment2 = c;
		        		break;
		        	}
		        }
		        continue;
		}
		break;
	}while(1);
 }
 while(nosegment2>0){
      	c = nosegment2; nosegment2 = st[c]; st[c] = 0;
	if(ccund(c)) continue;
   	if(cnosegment(ccu(c))){
		cSetRegion(c);
	}
 }
 Grid.firstRegion = 1;
 Grid.lastRegion  = Grid.cells.last();
}
*/

#define ibg2Found 0
#define ibg2NotFound (-1)
#define ibg2refreg(n1,norm)		\
	(*ibg2gRef.RefineRegion)(ibg2gRef.pRefineRegion,n1,norm)
#define ibg2refface(n1,norm,tang)	\
	(*ibg2gRef.RefinlFace)(ibg2gRef.pRefinlFace,n1,norm,tang)
#define ibg2refedge(n1,norm,tang)	\
	(*ibg2gRef.RefineEdge)(ibg2gRef.pRefineEdge,n1,norm,tang)
#define ibg2refvert(n1,norm)		\
	(*ibg2gRef.RefinlNode)(ibg2gRef.pRefinlNode,n1,norm)
#define ibg2refline(n1,n2) (*ibg2gRef.RefineLine)(ibg2gRef.pRefineLine,n1,n2)

int g3EDelete(int cb, int eb);

/* dimension-dependend data */
#define MAXR1 2
#define MAXO1 0

#define MAXR2 4
#define MAXO2 2
#define MAXQ2 0

#define MAXR3 6
#define MAXO3 4
#define MAXQ3 8

#define MAXP3 5	/* P: the number after O */

/*

enum q3richt{nnn,pnn,npn,nnp,npp,pnp,ppn,ppp,ooo};
enum n3richt{x3p,y3p,z3p,x3n,y3n,z3n,o3o};
static enum q3richt g3rrq[MAXR3][MAXR3]=
{{ooo,ppp,pnp,ooo,pnn,ppn},{ppn,ooo,ppp,npp,ooo,npn},{ppp,npp,ooo,nnp,pnp,ooo},
 {ooo,npn,npp,ooo,nnp,nnn},{pnp,ooo,nnp,nnn,ooo,pnn},{pnn,ppn,ooo,npn,nnn,ooo}};

static enum q3richt g3roq[MAXR3][MAXP3]=
{{ppp,ppn,pnn,pnp,ppp},{ppp,npp,npn,ppn,ppp},{ppp,pnp,nnp,npp,ppp},
 {npp,npn,nnn,nnp,npp},{pnp,nnp,nnn,pnn,pnp},{ppn,pnn,nnn,npn,ppn}};

static enum n3richt g3ror[MAXR3][MAXP3]=
{{y3p,z3n,y3n,z3p,y3p},{z3p,x3n,z3n,x3p,z3p},{x3p,y3n,x3n,y3p,x3p},
 {y3p,z3n,y3n,z3p,y3p},{z3p,x3n,z3n,x3p,z3p},{x3p,y3n,x3n,y3p,x3p}};

static enum n3richt g3rrr[MAXR3][MAXR3]=
{{o3o,z3p,y3n,o3o,z3n,y3p},{z3n,o3o,x3p,z3p,o3o,x3n},{y3p,x3n,o3o,y3n,x3p,o3o},
 {o3o,z3n,y3p,o3o,z3p,y3n},{z3p,o3o,x3n,z3n,o3o,x3p},{y3n,x3p,o3o,y3p,x3n,o3o}};


enum q3richt g3qother[MAXQ3]= {ppp,npp,pnp,ppn,pnn,npn,nnp,nnn};

enum q2richt{nn2,pn2,np2,pp2};
enum n2richt{x2p,y2p,x2n,y2n};
static enum q2richt g2roq[MAXR2][MAXO2]=
{{pp2,pn2},{np2,pp2},{np2,nn2},{nn2,pn2}};

static enum n2richt g2ror[MAXR2][MAXO2]=
{{y2p,y2n},{x2p,x2n},{y2p,y2n},{x2p,x2n}};


#if MAXQ2 > 0
     	static enum q2richt g2qother[MAXQ2]= {pp2,np2,pn2,nn2};
#endif

*/


typedef struct ibg2ExtData0 ibg2ExtData;

/* initialization, reading the file with the external parameters */
ibg2ExtData *ibg2ExtInit(char* filename);

/* getting the coarse grid information from the external parameters */
wzInteger  ibg2ExtReg(ibg2ExtData *r, wzInteger d);
wzInteger  ibg2ExtNum(ibg2ExtData *r, wzInteger d);
wzFloat* ibg2ExtX(ibg2ExtData *r, wzInteger d);

/* the refinement criteria using the external parameters */

void ibg2ExtRegion(ibg2PtObject This, wzPoint *x,  wzFloat *l);
int  ibg2ExtLine  (ibg2PtObject This, wzPoint *x1, wzPoint *x2);

	/* data  of the circumscribing sphere of the cell */


	/* cell data */

/* allocation routines: */
  /*
int g0CAlloc(wzGrid* g, int max)
{
 void *adr;int i;
 if((adr=malloc(max*sizeof(wzCellType)))==ibg2NULL) {goto overflow;}
 g->CellType = (wzCellType*) adr;
 if((adr=malloc(max*sizeof(ibg2Segment)))==ibg2NULL) {goto overflow;}
 g->CellSegment = (ibg2Segment*) adr;
 if((adr=malloc(max*sizeof(int)))==ibg2NULL) {goto overflow;}
 g->CPointFirst = (int*) adr;
 if((adr=malloc(max*sizeof(int)))==ibg2NULL) {goto overflow;}
 g->CSideFirst = (int*) adr;
 if((adr=malloc(2*max*sizeof(int)))==ibg2NULL) {goto overflow;}
 g->CPoint = (int*) adr;
 if((adr=malloc(2*max*sizeof(int)))==ibg2NULL) {goto overflow;}
 g->CSide = (int*) adr;
 if((adr=malloc(max*sizeof(ibgSphere)))==ibg2NULL) {goto overflow;}
 g->CellData	= (ibgSphere*) adr;
 g->CDataSize	= sizeof(ibgSphere);
 if((adr=malloc(g->maxPoint*sizeof(int)))==ibg2NULL) {goto overflow;}
 g->PointCell	= (int*) adr;
 g->maxCell	= max-1;
 g->maxCSide	= 2*max-IBG2CSMAX-1;
 g->maxCPoint	= 2*max-IBG2CNMAX-1;
 g->maxCOut	= 0;
 g->lastCOut	= 0;
 g->freeCSide	= 0;
 g->freeCPoint	= 0;
 g->firstCell   = 1;
 g->lastCell   = 0;
 g->NFree = wzGridNoPoint;
 return 1;
}
 */

void nwmax();
wzFloat ibg2gDelaunayBig = 1.e3;
wzFloat ibg2gDelaunayEps = 1.e-6;


int   ibg2AppendPoint(wzpoints list, wzPoint *p)
{
  int l;
  l  = list->create();
  list->Point[l] = *p;
  list->Near[l]  = 0;
  return ibg2Success;
}

int   ibg2AppendPoints(wzpoints list, wzpoints rest)
{
  int l,i0,i;
  i0 = list->create(rest->length());
  l  = list->last();
  wzRangeLoop(*rest,i){
    list->Point[i0] = rest->Point[i];
    if(rest->Near[i0])	list->Near[i0] = rest->Near[i] + l;
    else           	list->Near[i0] = 0;
    i0++;
  }
  return ibg2Success;
}

#define ibg2Found 0
#define ibg2NotFound (-1)
#define ibg2refreg(n1,norm)		\
	(*ibg2gRef.RefineRegion)(ibg2gRef.pRefineRegion,n1,norm)
#define ibg2refface(n1,norm,tang)	\
	(*ibg2gRef.RefinlFace)(ibg2gRef.pRefinlFace,n1,norm,tang)
#define ibg2refedge(n1,norm,tang)	\
	(*ibg2gRef.RefineEdge)(ibg2gRef.pRefineEdge,n1,norm,tang)
#define ibg2refvert(n1,norm)		\
	(*ibg2gRef.RefinlNode)(ibg2gRef.pRefinlNode,n1,norm)
#define ibg2refline(n1,n2) (*ibg2gRef.RefineLine)(ibg2gRef.pRefineLine,n1,n2)

int g3EDelete(int cb, int eb);

/* Hilfsdaten f"ur Startsituation */
//static int hcn3[16] =	{1,3,2,5, 1,5,2,4, 1,4,2,6, 1,6,2,3};
//static int hcs3[16] =	{0,2,0,4, 0,3,0,1, 0,4,0,2, 0,1,0,3};
int ibg2g3find(wzFloat *xx, int c0);


#define eps2 1.e-12

/*
int	ibg2ToYDL(wzGrid *gg, FILE* file);
int     ibg2ToOFF(wzGrid *gg, FILE* file);
void	ibg2ToGMS_Init();
int	ibg2ToGMS_X11(wzGrid *g);
int	ibg2ToGMS_PS (wzGrid *g, FILE *file);
int	ibg2ToGMS3D_PS (wzGrid *g, FILE *file);
int	ibg2ToGMS3D_X11 (wzGrid *g);

static const int ibg2MIsolator		= 2;
static const int ibg2MConductor		= 3;
static const int ibg2MSemiconductor	= 4;
static const int ibg2USiO2		= 2;
static const int ibg2UAl			= 3;
static const int ibg2USilicon		= 4;
static const int ibg2UNi2O3    		= 5;
static const int ibg2UAg			= 6;

int		ibg2Material[ibg2SMAX];
int		ibg2BCondition[ibg2SMAX];
int		ibg2ColRegion[ibg2SMAX];
int		ibg2CoeFace[ibg2SMAX];
int		ibg2ColEdge[ibg2SMAX];
int		ibg2ColNode[ibg2SMAX];
int		ibg2ColMaterial[ibg2SMAX];
int		ibg2ColBCondition[ibg2SMAX];
wzBoolean	ibg2OSSpecial[ibg2SMAX];
char*		wzNameRegion[ibg2SMAX];
char*		wzNamlFace[ibg2SMAX];
char*		wzNameEdge[ibg2SMAX];
char*		wzNamlNode[ibg2SMAX];
char*		wzNameMaterial[ibg2SMAX];
char*		wzNameBCondition[ibg2SMAX];

ibg2Segment		ibg2OSList[ibg2OSMaxLength];
wzFloat	ibg2WX[wzPointDim];
wzFloat	ibg2WAlpha,ibg2WBeta;
wzFloat	ibg2WXMin[wzPointDim];
wzFloat	ibg2WXMax[wzPointDim];
ibg2Color	ibg2ColorOf[ibg2ColorMAX];
char*		ibg2ColorName[ibg2ColorMAX];
wzBoolean	ibg2OutputType[IBG2STMAX];
char*		wzSName[IBG2STMAX] = {"IsRegion","IsFace","IsEdge","IsNode",
					""};
int 		(*ibg2BCDefault)(ibg2Segment u1, ibg2Segment u2) = ibg2BCDefaultDirichlet;
wzBoolean 	ibg2MOn;
wzBoolean	ibg2BCOn;
ibg2SegmentType	ibg2OSType;
wzBoolean	ibg2OS_Outside;
int		ibg2OSMode;
int		ibg2OSListLength;
char*		ibg2OSName;
int		ibg2IdPoint;
int		ibg2IdCell;
wzBoolean	ibg2Window;
int		ibg2ColorFree;
int		ibg2ColorCurrent;
int		ibg2ColDoubleBoundary;
int		ibg2ColSingleBoundary;
int		ibg2Col_OutsideBoundary;
ibg2Segment		ibg2OSSegment;
wzBoolean	ibg2ResetGrid;
char*		wzNameApplication;
char*		wzNameFile;
//int  		ibg2FileMode;
wzGrid		*ibg2gCurrent;
int  	 	ibg2FileModFree;
char*	 	ibg2FileModName[ibg2FileModMAX];
char*	 	ibg2FillType[ibg2FileModMAX];
int	 	(*(ibg2FileOutput[ibg2FileModMAX]))(wzGrid *g, FILE *file);
//static 	int	ibg2readcolors();

void	ibg2OutputInit()
{ibg2Segment u; int color,fmod;
 ibg2BCDefault	= ibg2BCDefaultDirichlet;
 ibg2MOn		= ibg2False;
 ibg2BCOn	= ibg2False;
 ibg2OSType	= ibg2SRegion;
 ibg2OS_Outside	= ibg2False;
 ibg2OSMode	= ibg2OSModSpecial;
 ibg2OSListLength	= 0;
 ibg2OSName	= ibg2OSNoName;
 ibg2IdPoint	= ibg2IdModOff;
 ibg2IdCell	= ibg2IdModOff;
 ibg2Window	= ibg2False;
 ibg2ColorFree	= 0;
 ibg2ColorCurrent	= -1;
 ibg2ColDoubleBoundary	= ibg2ColModMax;
 ibg2ColSingleBoundary	= ibg2ColModSelf;
 ibg2Col_OutsideBoundary	= ibg2ColModSelf;
 ibg2OSSegment		= ibg2RSpace;
 ibg2ResetGrid 		= ibg2True;
 if(wzNameApplication == ibg2NULL) wzNameApplication	= "ibg-test";
 wzNameFile		= "ibg-test";
 ibg2FileMode		= ibg2FileModNone;
 ibg2FileModFree		= 0;
 ibg2gCurrent		= ibg2NoGrid;

 for(u=0;u<ibg2SMAX;u++){
	ibg2Material[u]		= u;
	ibg2BCondition[u]	= u;
	wzNameRegion[u]	= ibg2OSNoName;
	wzNamlFace[u]		= ibg2OSNoName;
	wzNameMaterial[u]	= ibg2OSNoName;
	wzNameBCondition[u]	= ibg2OSNoName;
	ibg2ColRegion[u]	= ibg2ColDefaultRegion;
	ibg2CoeFace[u]		= ibg2ColDefaultFace;
	ibg2ColEdge[u]		= ibg2ColDefaultEdge;
	ibg2ColNode[u]		= ibg2ColDefaultNode;
	ibg2ColMaterial[u]	= ibg2ColDefaultRegion;
	ibg2ColBCondition[u]	= ibg2ColDefaultFace;
	ibg2OSSpecial[u]		= ibg2True;
 }
 for(u=0;u<IBG2STMAX;u++)	ibg2OutputType[u] = ibg2False;
 wzNameRegion[0]	= "Unknown";
 wzNamlFace[0]		= "Unknown";
 wzNameMaterial[0]	= "Unknown";
 wzNameBCondition[0]	= "Unknown";
 wzNameMaterial[ibg2MAir]		= "Air";
 wzNameRegion[ibg2RSpace]		= "Space";
 wzNameBCondition[ibg2BCDirichlet] 	= "Dirichlet";
 wzNameBCondition[ibg2BCNeumann]	= "Neumann";
 wzNameBCondition[ibg2BCMixed]		= "Mixed";
 // ibg2readcolors();
 color = ibg2DeclareColor("ColBackground",	"#fdf5e6"); /* OldLace	
 wzAssert(color == ibg2ColBackground);
 color = ibg2DeclareColor("ColError",		"#ff0000"); /* red	
 wzAssert(color == ibg2ColError);
 color = ibg2DeclareColor("ColForeground",	"#000000"); /* black	
 wzAssert(color == ibg2ColForeground);
 color = ibg2DeclareColor("ColDefaultNode",	"#b22222"); /* firebrick  
 wzAssert(color == ibg2ColDefaultNode);
 color = ibg2DeclareColor("ColDefaultEdge",	"#ff6347"); /* tomato	
 wzAssert(color == ibg2ColDefaultEdge);
 color = ibg2DeclareColor("ColDefaultFace",	"#e9967a"); /* DarkSalmon 
 wzAssert(color == ibg2ColDefaultFace);
 wzAssert(color == ibg2ColDefaultBCondition);
 ibg2ColBCondition[ibg2BCDirichlet]    	=
 	 ibg2DeclareColor("ColDirichlet", 	"#32cd32"); /* LimeGreen 
 ibg2ColBCondition[ibg2BCNeumann]		=
	 ibg2DeclareColor("ColNeumann",		"#4169e1"); /* RoyalBlue 
 ibg2ColBCondition[ibg2BCMixed]		=
	 ibg2DeclareColor("ColMixed",		"#20b2aa"); /* LightSeaGreen 
 color = ibg2DeclareColor("ColDefaultRegion",	"#a0522d"); /* sienna	
 wzAssert(color == ibg2ColDefaultRegion);
 wzAssert(color == ibg2ColDefaultMaterial);
 ibg2ColMaterial[ibg2MAir] 	       	=  ibg2ColRegion[ibg2RSpace]	=
	 ibg2DeclareColor("ColAir",  		"#c0c0c0");	/* gray	
 fmod =	ibg2DeclareFileMode("None",	"",	0);
 wzAssert(fmod == ibg2FileModNone);
 ibg2FileMode =	ibg2DeclareFileMode("Simplex",	"sg",	ibg2ToSimplex);
 wzAssert(ibg2FileMode == ibg2FileModSimplex);
}

void	ibg2Output(wzGrid *g)
{char c,filename[500];
 int i,rc,notype=1;
 FILE *file;
 int d;
 if(g==ibg2NULL) return;
 ibg2gCurrent = g;
 if(ibg2ResetGrid){
	ibg2ResetGrid = ibg2False;
	for(d=0;d<wzPointDim;d++){
		ibg2WXMin[d] = wzGridXMin(*g,d);
		ibg2WXMax[d] = wzGridXMax(*g,d);
		ibg2WX[d]    = wzGridXMax(*g,d) - wzGridXMin(*g,d);
		if(ibg2WX[d]<= 0.0) ibg2WX[d] = 1.0;
	}
 	ibg2Window = ibg2True;
	if(wzGridGDIM(*g) == 3){
	      	ibg2OutputType[ibg2SRegion]	= ibg2False;
		ibg2OutputType[ibg2SFace]		= ibg2True;
	}else{
	      	ibg2OutputType[ibg2SRegion]	= ibg2True;
		ibg2OutputType[ibg2SFace]		= ibg2False;
	}
 }
 ibg2Control();
 if((ibg2FileMode <= ibg2FileModNone)||
    (ibg2FileMode >= ibg2FileModFree)) {
     ibg2message(