ibg2.cxx

Delaunay三角形的网格剖分程序

#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "ibg2.hxx"

typedef struct _ibg2Topology	ibg2TopologyRec,*ibg2Topology;
const wzInteger IBG2TMAXFREGS=10;
struct _ibg2Topology{
	int		fregs;
	wzISegment*	freg[IBG2TMAXFREGS];
};
#define ibg2Segment wzISegment

#define ibg2pX(point)		((point).X)
#define ibg2pF(point)		((point).X+wzPointDim)
#define ibg2pI(point)		((point).I)
#define ibg2pSegment(point)	((point).S.index())
#define ibg2pType(point)	((point).S.codim())
#define ibg2pSetSegment(point,tt,uu)	((point).S = wzSegment((wzIndex)(uu),(wzDimension)(tt)))
#define ibg2SMAX 100
#define ibg2DefaultFaceNr(m1,m2) ((m1)>(m2)? (m1)*ibg2SMAX+(m2):(m2)*ibg2SMAX+(m1))
// old stuff ( to be removed in some future ): 


#define ibg2printf	printf
#define ibg2flush       fflush(stdout)

#define big1 1.e6
#define eps2 1.e-12
#define g0CGet(c,t)	(*(c) = createCell(t))
#define cnosegment(uu)	((uu)< 0)	/* undefined segment */
#define _nothing (-999999999)
#define _outside (0)


void	ibg2tInit(ibg2Topology t);
void	ibg2tCopy(ibg2Topology tnew, ibg2Topology told);
void	ibg2tFree(ibg2Topology t);
void ibg2tNewFace(ibg2Topology t, ibg2Segment snew, ibg2Segment old);
void ibg2tAppendFRegion(ibg2Topology t, ibg2Segment face, ibg2Segment reg);
int ibg2tStatusRV(ibg2Topology top, int u1, int u2);
int ibg2tStatusFV(ibg2Topology top, int u1, int u2);
int ibg2tStatusLV(ibg2Topology top, int u1, int u2);
int ibg2tStatusVV(ibg2Topology top, int u1, int u2);
int ibg2tStatusRF(ibg2Topology top, int u1, int u2);
int ibg2tStatusFF(ibg2Topology top, int u1, int u2);
int ibg2tStatusRL(ibg2Topology top, int u1, int u2);
int ibg2tStatusFL(ibg2Topology top, int u1, int u2);
int ibg2tStatusLL(ibg2Topology top, int u1, int u2);
int ibg2iStatusOfLine(ibg2Geometry geom, const wzPoint *n1, const wzPoint *n2);

void wzGridInit();


/*
extern wzFloat ibg2gDmin[wzPointDim];
extern wzFloat ibg2gDefaultRegionLength,ibg2gDefaultNodeLength;
extern wzFloat ibg2gDefaultFaceNormal,ibg2gDefaultFaceTangential;
extern wzFloat ibg2gDefaultEdgeNormal,ibg2gDefaultEdgeTangential;
void		ibg2gDefaultRefineRegion(ibg2PtObject This, wzPoint *point,
		   	      	wzFloat *length);
void		ibg2gDefaultRefinlFace(ibg2PtObject This, wzPoint *point,
	     	   	      	wzFloat *normal, wzFloat *tangential);
void		ibg2gDefaultRefineEdge(ibg2PtObject This, wzPoint *point,
		   	      	wzFloat *normal, wzFloat *tangential);
void		ibg2gDefaultRefinlNode(ibg2PtObject This, wzPoint *point,
		   	      	wzFloat *length);
wzBoolean	ibg2gDefaultRefineLine (ibg2PtObject This,
				wzPoint *n1, wzPoint *n2);
				*/
extern wzFloat ibg2gDelaunayBig, ibg2gDelaunayEps;
/*
extern wzFloat ibg2gBV;
extern wzFloat ibg2gBLWait,ibg2gBLHalf,ibg2gBLObtuse,ibg2gBctg;
extern wzFloat ibg2gOref,ibg2gOref2,ibg2gOinv,ibg2gOmax;
extern wzFloat ibg2gPref;
extern int	ibg2gPrefL;
*/
//wzGrid	*wzGridCopy(wzGrid *grid);
//wzGrid	*wzGridFree(wzGrid *gg);
void wzGridFunction(wzGrid *gg, int fnum,
	 	 	wzFloat	(*func)(ibg2PtObject fThis, wzPoint *x),
			ibg2PtObject	fThis);

#define wzGridGDIM(g)      ((g).gdim())
#define wzGridXDIM(g)      ((g).xdim())
#define wzGridFirstPoint(g)	(1)
#define wzGridLastPoint(g)	((g).Grid.points.last())
#define wzGridPoint(g,n)	((g).Grid.Point[n])
#define wzGridPointCell(g,n)	((g).Grid.PointCell[n])
#define wzGridFirstCell(g)   	(1)
#define wzGridLastCell(g)    	((g).Grid.cells.last())
#define wzGridCellType(g,c)  	((g).cct(c))
#define wzGridCellSegment(g,c)  ((g).ccu(c))
#define wzGridCelePoint(g,c,n)  ((g).ccn(c)[n])
#define wzGridCPointList(g,c)   ((g).ccn(c))
#define wzGridCellSide(g,c,s)	((g).ccs(c)[s])
#define wzGridCSideList(g,c)	((g).ccs(c))
#define wzGridNFree(g,n)	((g).Grid.points.destroy(n))
#define wzGridNGet(g,n)		((nh)=(g).Grid.points.create())
//#define wzGridHasPointCell(g)	((g).PointCell)
//#define wzGridHasSides(g)	((g).CSide)

/*............. list of "_outside" neighbours for cells with codim > 1:	*/
//#define wzGridHas_Outsides(g)	((g).COutside)
#define wzGridLastCOut(g)	((g).lastCOut)
#define wzGridCOutFirst(g,c,t)	((g).Data[(g).CellData[c]+wzCellTypeSExtF(t)])
#define wzGridCOutNext(g,o)	((g).COutNext[o])
#define wzGridCOutside(g,o)	((g).COutside[o])

/*........ begin and end for every geometrical type (codimension)   */
#define wzGridFirstRegionCell(g)  ((g).Grid.firstRegion)
#define wzGridLastRegionCell(g)   ((g).Grid.lastRegion)
#define wzGridFirstFaceCell(g)    ((g).Grid.firstFace)
#define wzGridLastFaceCell(g)     ((g).Grid.lastFace)
#define wzGridFirstEdgeCell(g)    ((g).Grid.firstEdge)
#define wzGridLastEdgeCell(g)     ((g).Grid.lastEdge)
#define wzGridFirstNodeCell(g)    ((g).Grid.firstNode)
#define wzGridLastNodeCell(g)     ((g).Grid.lastNode)

/*........ coordinates of the window for creating a geometry from the grid: */
#define wzGridXMin(g,i)        ((g).Grid.mincoord[i+1])
#define wzGridXMax(g,i)        ((g).Grid.maxcoord[i+1])

/*........ other data fields: */
#define wzGridHasCellData(g)	((g).Grid.CellData)
#define wzGridHasPointData(g)	((g).Grid.PointData)

#define wzGridPointData(g,n) ((char *)((g).Grid.PointData) + (g).NDataSize*(n))
#define wzGridCellData(g,c) ((char *)((g).Grid.CellData) + (g).CDataSize*(c))

//int	wzGridCNew(wzGrid *g, wzCellType t);
//int wzGridNRealloc(wzGrid*g);	/* to be used only in: */
//#define wzGridNNew(g)	  ((g)->lastPoint<(g)->maxPoint  ? ++(g)->lastPoint :\
//				     (wzGridNRealloc(g)  ? ++(g)->lastPoint : 0))
//int wzGridCORealloc();	/* to be used only in: */
//#define wzGridCONew(g)	  ((g)->lastCOut<(g)->maxCOut  ? ++(g)->lastCOut :\
	    			     (wzGridCORealloc(g) ? ++(g)->lastCOut : 0))

#define forwzPoints(g,in,t,u,point) \
	for(in=wzGridFirstPoint(g); in<=wzGridLastPoint(g); in++){\
		point = &wzGridPoint(g,in);\
		u = ibg2pSegment(*point);\
		t = ibg2pType(*point);\
		if((t == ibg2SRegion) && (u <= 0)) continue;
#define endwzPoints(g,in,t,u,point)	}
/*............ cycle over cells of a grid g:
	parameter: ic - cell number, u - segment of the cell, t - type of the cell
*/
#define foribg2Cells(g,ic,t,u)\
		for(ic=wzGridFirstCell(g);ic<=wzGridLastCell(g);ic++)\
		if(t=wzGridCellType(g,ic),(u=wzGridCellSegment(g,ic))>0)
#define endibg2Cells(g,ic,t,u)
#define foribg2RegionCells(g,ic,t,u)	\
	for(ic=wzGridFirstRegionCell(g);ic<=wzGridLastRegionCell(g);ic++)\
	if((u=wzGridCellSegment(g,ic))>0)\
	if(wzCellTypeCODIM(t=wzGridCellType(g,ic))==ibg2SRegion)
#define endibg2RegionCells(g,ic,t,u)
#define foribg2FaceCells(g,ic,t,u)	\
	for(ic=wzGridFirstFaceCell(g);ic<=wzGridLastFaceCell(g);ic++)\
	if((u=wzGridCellSegment(g,ic))>0)\
	if(wzCellTypeCODIM(t=wzGridCellType(g,ic))==ibg2SFace)
#define endibg2FaceCells(g,ic,t,u)
#define foribg2EdgeCells(g,ic,t,u)	\
	for(ic=wzGridFirstEdgeCell(g);ic<=wzGridLastEdgeCell(g);ic++)\
		if(wzCellTypeCODIM(t=wzGridCellType(g,ic))==ibg2SEdge)\
			if((u=wzGridCellSegment(g,ic))>0)
#define endibg2EdgeCells(g,ic,t,u)
#define foribg2NodeCells(g,ic,t,u)	\
	for(ic=wzGridFirstNodeCell(g);ic<=wzGridLastNodeCell(g);ic++)\
		if(wzCellTypeCODIM(t=wzGridCellType(g,ic))==ibg2SNode)\
			if((u=wzGridCellSegment(g,ic))>0)
#define endibg2NodeCells(g,ic,t,u)

			  /*
extern wzGrid  ibg2g;
static wzInteger* ccn(wzInteger cc) {return wzGridCPointList(*this,cc);}
static wzFloat* cnx(wzInteger nn){return ibg2pX(wzGridPoint(*this,nn));}
static ibg2Segment ccu(wzInteger cc) {return wzGridCellSegment(*this,cc);}
static wzCellType cct(wzInteger cc) {return wzGridCellType(*this,cc);}
static wzInteger* ccs(wzInteger cc) {return wzGridCSideList(*this,cc);}
*/

//int wzGridCNew(wzGrid *g, wzCellType t);

#define	wzGridCFree(g,c)	(g).destroyCell(c)
#define wzGridCGet(g,c,t)	((*c) = (g).createCell(t))
/* ......... macros connected with the CFree-field:
#define	wzGridNoCell	123456789
#define	wzGridCFree(g,c)	\
	(wzGridCellSegment(g,c) = -(g).CFree[wzGridCellType(g,c)],\
	(g).CFree[wzGridCellType(g,c)] = c,\
	wzGridCellType(g,c) = wzCellType0Nothing)
#define wzGridCGet(g,c,t)		\
	((*(c)=(g).CFree[t])!=wzGridNoCell )?(	\
		(g).CFree[t] = - wzGridCellSegment(g,*(c)),	\
		wzGridCellType(g,*(c)) = (t),	\
		wzGridCellSegment(g,*(c)) =  0 	\
	):(	*(c) = wzGridCNew(&(g),t)  	\
	)
 */

/* ......... macros connected with the NFree-field: */
#define	wzGridNoPoint	12345678

extern	int	ibg2MsgLevel;
extern	int	ibg2MsgOn;
extern	int	ibg2MStatus;


void	ibg2Message(int err, char* file, int edge);

#define ibg2message(x)	ibg2Message(x,__FILE__,__LINE__)
#define ibg2msg		if(ibg2MsgOn) ibg2printf
#define ibg2msgreset	ibg2message(ibg2MNone)
#define ibg2warning	ibg2message( ibg2MWInternal)

/* list of messages for external users (they have to be documented) */
#define ibg2MNone	  0

#define ibg2MIInit	  1
#define ibg2MIQuit	  2
#define ibg2MIStart	  3
#define ibg2MISuccess	  4
#define ibg2MIRefine	  5
#define ibg2MIBShift	  6
#define ibg2MIDelaun	  7
#define ibg2MIbg2ridG	  8
#define ibg2MINoGrid	  9 
#define ibg2MIWriteFile	 10 

#define ibg2MInformation  19

#define ibg2MWGridCoarse	 20
#define ibg2MWBoundMod	 21
#define ibg2MWBoundInc	 22
#define ibg2MWGridInc     23
#define ibg2MWSmallGeom	 24
#define ibg2MWUncGeom	 25
#define ibg2MWIncFFace	 26
#define ibg2MWInternal	 27

#define ibg2MWarning      49

#define ibg2MESmallDelta	 50
#define ibg2MEMaxFRegs	 51
#define ibg2MEOverflow	 52
#define ibg2MEIncCoarse	 53
#define ibg2MELoadError   54
#define ibg2MEIncGeom	 55
#define ibg2MEDemo	 56
#define ibg2MEIncFileMode 57
#define ibg2MEWriteFile   58
#define ibg2MEOpenFile    59
#define ibg2MEWrongFormat 60
#define ibg2MEInternal	 80

#define ibg2MError        99 

/*
extern	wzGrid		ibg2g;
*/

#define g0ibig 123456789

#undef cnl
#define cndef(nn)	((nn)> 0)     	/* defined point */
#define cnund(nn)	((nn)<=0)	/* undefined point */
#define cn_outside(nn)	((nn)==0)	/* point _outside */
#define cn_nothing(nn)	((nn)< 0)	/* undefined point */

#define rother(r)	((r) < gdim() ? (r) + gdim() : (r) - gdim())
#define rposit(r)	((r) < gdim() ? (r)         : (r) - gdim())
#define rnegat(r)	((r) < gdim() ? (r) + gdim() : (r)        )
#define rispos(r)	((r) < gdim())
#define dnegat(d)	((d) + gdim())
/*
#define g0NNew		g0NNew0()
#define cct(cc) wzGridCellType(*this,cc)
#define ccu(cc) wzGridCellSegment(*this,cc)
#define ccn(cc) wzGridCPointList(*this,cc)
#define ccs(cc)	wzGridCSideList(*this,cc)
*/
#define ccdef(cc)	(ccu(cc)> 0)
#define ccout(cc)	(ccu(cc)==0)

ibg2Geometry ibg2dGrid(wzGrid *gg, int fanz, int *flist);
void ibg2ControlInit();
void ibg2Control();

/*
void	ibg2OutputInit();
void	ibg2Output(wzGrid *g);
wzBoolean	ibg2OSelected(ibg2Segment u);
wzBoolean	ibg2OSelectedRegion(ibg2Segment u);
wzBoolean	ibg2OSelectedFace(ibg2Segment u, ibg2Segment left, ibg2Segment right);
wzBoolean	ibg2OSelectedEdge(ibg2Segment u);
wzBoolean	ibg2OSelectedNode(ibg2Segment u);
wzBoolean	ibg2IdentPoint(int n, char* ident);
wzBoolean	ibg2IdentCell(int c, char* ident);
wzBoolean	ibg2InWindow(wzPoint *n);
#define	ibg2IsInWindow(n)	(ibg2Window ? ibg2InWindow(n) : ibg2True )
wzBoolean	ibg2WindowCoord(wzPoint *n, wzFloat *x);
int	ibg2ColorRegion(ibg2Segment u);
int	ibg2ColorFace(ibg2Segment u, ibg2Segment left, ibg2Segment right);
int	ibg2ColorEdge(ibg2Segment u);
int	ibg2ColorNode(ibg2Segment u);
wzBoolean	ibg2OSetColor(int id, char* colorname);
int	ibg2DeclareColor(char* name, char* value);
int	ibg2DeclareFileMode(char* mode, char* filetype,
	   	   	int (*output)(wzGrid *g, FILE *file)
	);
int wzGridSave          (wzGrid *gg, FILE *file);
int wzGridWrite         (wzGrid *gg, FILE *file);
int ibg2ToSimplex	(wzGrid *gg, FILE *file);
wzGrid *wzGridLoad      (FILE *file);
wzGrid *wzGridRead      (FILE *file);
wzGrid *ibg2OtSimplex	(FILE *file);
extern	int		ibg2Material[ibg2SMAX];
extern	wzBoolean	ibg2MOn;
extern	int		ibg2BCondition[ibg2SMAX];
extern	wzBoolean	ibg2BCOn;
static const int	ibg2RSpace       = 1;
static const int	ibg2MAir	        = 1;
static const int	ibg2BCDirichlet	= 1;
static const int	ibg2BCNeumann	= 2;
static const int	ibg2BCMixed	= 3;
extern	int	 (*ibg2BCDefault)(ibg2Segment u1, ibg2Segment u2);
int	ibg2BCDefaultDirichlet(ibg2Segment u1, ibg2Segment u2);
int	ibg2BCDefaultNeumann(ibg2Segment u1, ibg2Segment u2);
int	ibg2BCDefaultMixed(ibg2Segment u1, ibg2Segment u2);
int	ibg2BCDefaultMax(ibg2Segment u1, ibg2Segment u2);
int	ibg2BCDefaultMin(ibg2Segment u1, ibg2Segment u2);
extern	char*		wzSName[IBG2STMAX];
extern	char*		wzNameRegion[ibg2SMAX];
extern	char*		wzNamlFace[ibg2SMAX];
extern	char*		wzNameEdge[ibg2SMAX];
extern	char*		wzNamlNode[ibg2SMAX];
extern	char*		wzNameMaterial[ibg2SMAX];
extern	char*		wzNameBCondition[ibg2SMAX];
extern	char*		wzNameApplication;
extern	ibg2SegmentType	ibg2OSType;
extern	wzBoolean	ibg2OutputType[IBG2STMAX];
extern	wzBoolean	ibg2OS_Outside;
extern	int		ibg2OSMode;
*/
typedef enum{
	ibg2OSModAll,	/* all segments */
	ibg2OSModOne,	/* only segment ibg2OSSegment */
	ibg2OSModSpecial,	/* segments with ibg2OSSpecial[u] */
	ibg2OSModMat,	/* segments by type (non-zero material/boundarytype) */
	ibg2OSModList,	/* the segments from the list. */
	ibg2OSModName,	/* the segments named ibg2OSName */
	ibg2OSModNamed,	/* the segments with names */
	ibg2OSMod_Nothing
}ibg2OSModType;

extern	ibg2Segment		ibg2OSSegment;
extern	wzBoolean	ibg2OSSpecial[ibg2SMAX];
extern	int		ibg2OSListLength;
#define	ibg2OSMaxLength	100
extern	ibg2Segment		ibg2OSList[ibg2OSMaxLength];
extern	wzName		ibg2OSName;
extern  wzName		ibg2OSNoName;

extern	wzBoolean	ibg2ResetGrid;
extern	int	ibg2IdPoint;	/* which type of identifier for points */
extern	int	ibg2IdCell;	/* which type of identifier for cells */

typedef enum{
	ibg2IdModOff,
	ibg2IdModOn,
	ibg2IdModTyp,
	ibg2IdModSegment,
	ibg2IdModMaterial,
	ibg2IdModNType,
	ibg2IdModNSegment,
	ibg2IdModNMaterial,
	ibg2IdMod_Nothing
}ibg2IdModType;
extern	wzBoolean	ibg2Window;		/* allows to avoid tests */
extern	wzFloat	ibg2WAlpha,ibg2WBeta;
extern	wzFloat	ibg2WX[wzPointDim];
extern	wzFloat	ibg2WXMin[wzPointDim];
extern	wzFloat	ibg2WXMax[wzPointDim];
extern	int	ibg2ColRegion[ibg2SMAX];
extern	int	ibg2CoeFace[ibg2SMAX];
extern	int	ibg2ColEdge[ibg2SMAX];
extern	int	ibg2ColNode[ibg2SMAX];
extern	int	ibg2ColMaterial[ibg2SMAX];
extern	int	ibg2ColBCondition[ibg2SMAX];

typedef	struct{
	unsigned	r;
	unsigned	g;
	unsigned	b;
}ibg2Color;

static const int ibg2ColOn	= (0xffff);

#define	ibg2ColorMAX	0xff
extern	ibg2Color	ibg2ColorOf[ibg2ColorMAX];
extern	char*		ibg2ColorName[ibg2ColorMAX];
extern	int		ibg2ColorFree;
extern	int		ibg2ColorCurrent;
extern	int	ibg2ColDoubleBoundary;	/* above regions are selected */
extern	int	ibg2ColSingleBoundary;	/* one region is selected */
extern	int	ibg2Col_OutsideBoundary;	/* one region is _outside */
typedef enum{
ibg2ColModZero,	/* zero color value */
  ibg2ColModMax,	/* the maximum of left and right */
  ibg2ColModMin,	/* the minimum of left and right */
  ibg2ColModSelf,	/* the selected region */
  ibg2ColModOther,	/* the region which is not selected */
  ibg2ColModBC,	/* the boundary condition */
  ibg2ColMod_Nothing
}ibg2ColModType;