ibgg0.c

有限元学习研究用源代码(老外的),供科研人员参考

/* last edit: Ilja Schmelzer -------------- 24-MAR-1995 19:53:26.67	*/
/************************************************************************/
/*                                                                      */
/*  <<< I B G >>> - Intersection - Based Grid generation package 	*/
/*                                                                      */
/*  Version 1.1 by Ilja Schmelzer   [e-mail]       */
/*                                                                      */
/*  to be distributed under IBG license conditions (see "readme.ibg")	*/
/*                                                                      */
/************************************************************************/
#include "ibgi.h"
#include "ibgg0.h"
#include "ibgd0.h"
#include "ibgg.h"
#include "ibgdefault.h"
#include "ibglib.h"

#include "ibgg0p.h"

	ibGeometry	ibGeo;	/* current geometry */
	ibGrid0	ibgg0;	/* current internal grid */

#define ibgFound 0
#define ibgNotFound (-1)
#define ibgrefreg(n1,norm)		\
	(*ibggRef.RefineRegion)(ibggRef.pRefineRegion,n1,norm)
#define ibgrefface(n1,norm,tang)	\
	(*ibggRef.RefineFace)(ibggRef.pRefineFace,n1,norm,tang)
#define ibgrefline(n1,norm,tang)	\
	(*ibggRef.RefineLine)(ibggRef.pRefineLine,n1,norm,tang)
#define ibgrefvert(n1,norm)		\
	(*ibggRef.RefineNode)(ibggRef.pRefineNode,n1,norm)
#define ibgrefedge(n1,n2) (*ibggRef.RefineEdge)(ibggRef.pRefineEdge,n1,n2)

/* local creation procedures: */
static int    	g0BNew0();

/* local boundary shift procedures: */
static int	g0bshift(int nn, int firstcall);
static void	g0bjump(int n0);
static int	g0breg(int bb);

/* local refinement procedures: */
static void	g1nref(int nn);
static int	g1eref(int nu, unsigned ro);
static void	g2nref(int nn);
static int	g2eref(int nu, unsigned ro);
static int	g2elot(int nn, unsigned rr, unsigned ro);
static void	g3nref(int nn);
static int	g3eref(int nu, unsigned ro);
static int	g3elot(int nn, unsigned rr, unsigned ro);
static int	g3elo2(int nn, unsigned rr, unsigned ro);
static int	g3qref2(int qq, unsigned ro);
static int	g3qref1(int qq, unsigned ro);
static int	g3qref(int qq);

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

static ibgFloat ibgNorm,ibgtang;
/* <<< dimension-dependend arrays >>> */
/* dim == 3: */
#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}};

/* g3roq[r][o] is the octant of r, g3ror[r][o-1] and g3ror[r][o] */
/* g3rrq[r][o] is the positive octant of r,o;	*/
/* g3rrr[r][o] is the positive third direction to r,o;	*/

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

/* dim == 2: */
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}};

/* g2roq[r][o] is the quadrant of r and g2ror[r][o] */

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

int ibgg0getPoints(ibgPoints *p, ibGrid0 *g)
{
 ibgg0= *g; ibgg = *(g->gg);
 g0BNShift();
 g0QFree(g);
 p->spaceDimension = ibgg.spaceDimension;
 p->firstPoint     = 1;
 p->lastPoint      = ibgg.lastPoint;
 p->maxPoint       = ibgg.maxPoint;
 p->Point          = ibgg.Point;
 p->Near           = ibgg0.no;
}

/* allocation routines: */
int g0CAlloc(int max)
{void *adr;int i;
 if((adr=malloc(max*sizeof(ibgCellType)))==ibgNULL) {goto overflow;}
 ibgg.CellType = (ibgCellType*) adr;
 if((adr=malloc(max*sizeof(ibgSegment)))==ibgNULL) {goto overflow;}
 ibgg.CellSegment = (ibgSegment*) adr;
 if((adr=malloc(max*sizeof(int)))==ibgNULL) {goto overflow;}
 ibgg.CPointFirst = (int*) adr;
 if((adr=malloc(max*sizeof(int)))==ibgNULL) {goto overflow;}
 ibgg.CSideFirst = (int*) adr;
 if((adr=malloc(2*max*sizeof(int)))==ibgNULL) {goto overflow;}
 ibgg.CPoint = (int*) adr;
 if((adr=malloc(2*max*sizeof(int)))==ibgNULL) {goto overflow;}
 ibgg.CSide = (int*) adr;
 if((adr=malloc(max*sizeof(tcv)))==ibgNULL) {goto overflow;}
 ibgg.CellData	= (tcv*) adr;
 ibgg.CDataSize	= sizeof(tcv);
 if((adr=malloc(ibgg.maxPoint*sizeof(int)))==ibgNULL) {goto overflow;}
 ibgg.PointCell	= (int*) adr;
 ibgg.maxCell	= max-1;
 ibgg.maxCSide	= 2*max-IBGCSMAX-1;
 ibgg.maxCPoint	= 2*max-IBGCNMAX-1;
 ibgg.maxCOut	= 0;
 ibgg.lastCOut	= 0;
 ibgg.freeCSide	= 0;
 ibgg.freeCPoint	= 0;
 ibgg.firstCell   = 1;
 ibgg.firstRegion = 1;
 ibgg.firstFace   = 1;
 ibgg.firstLine   = 1;
 ibgg.firstNode = 1;
 ibgg.lastCell   = 0;
 ibgg.lastRegion = 0;
 ibgg.lastFace   = 0;
 ibgg.lastLine   = 0;
 ibgg.lastNode = 0;
 for(i=0;i<IBGCTMAX;i++) ibgg.CFree[i] = ibgridNoCell;
 ibgg.NFree = ibgridNoPoint;
 return 1;
overflow:
 ibgfatal;
 return 0;
}

static int g0NAlloc(int last, int max)
{void *adr; int d;
 ibgg.lastPoint=last;
 if(last+10 >= max) max = last+100;
 max++;
 if((adr=malloc(max*sizeof(ibgPoint)))==ibgNULL) {goto overflow;}
 ibgg.Point = (ibgPoint *) adr;
 ibgg.CDataSize	= sizeof(tcv);
 ibgg0.nx = (char *) &(((ibgPoint *)adr)->xx[0]);
 ibgg0.nt = (char *) &(((ibgPoint *)adr)->tt);
 ibgg0.nu = (char *) &(((ibgPoint *)adr)->uu);
 ibgg0.nf = (char *) &(((ibgPoint *)adr)->ff[0]);
#ifdef GLEVEL
 ibgg0.nl = (char *) &(((ibgPoint *)adr)->ll[0]);
#endif
 if((adr=malloc(max*sizeof(int)))==ibgNULL)	{goto overflow;}
 ibgg0.nstack = (int*) adr;
 for(d=0;d<ibgg.gridDimension;d++){
 	if((adr=malloc(max*sizeof(int)))==ibgNULL)	{goto overflow;}
 	ibgg0.ndp[d] = (int*)adr;
 	ibgg0.ndr[d] = (int*)adr;
 	if((adr=malloc(max*sizeof(int)))==ibgNULL)	{goto overflow;}
 	ibgg0.ndn[d] = (int*)adr;
 	ibgg0.ndr[ibgg.gridDimension+d] = (int*)adr;
 }
 if((adr=malloc(max*sizeof(int)))==ibgNULL)		{goto overflow;}
 ibgg0.no = (int*)adr;
 ibgg.maxPoint=max-1;
 return 1;
overflow:
 ibgfatal;
 return 0;
}

/* reallocation of all parts of ibgg0: */

int g0NRealloc()
{void *adr; int max; int d;
 if(ibgg.lastPoint<ibgg.maxPoint)	return 1;
 if(ibgridNRealloc(&ibgg)==0)	return 0;
 ibgg0.nx = (char *) &(ibgg.Point->xx[0]);
 ibgg0.nt = (char *) &(ibgg.Point->tt);
 ibgg0.nu = (char *) &(ibgg.Point->uu);
 ibgg0.nf = (char *) &(ibgg.Point->ff[0]);
#ifdef GLEVEL
 ibgg0.nl = (char *) &(ibgg.Point->ll[0]);
#endif
 if((adr=realloc(ibgg0.nstack,(ibgg.maxPoint+1)*sizeof(int)))==ibgNULL)
        	       				{goto overflow;}
 ibgg0.nstack = (int*)adr;
 for(d=0;d<ibgg0.maxq;d++){
        if((adr=realloc(ibgg0.nq[d],(ibgg.maxPoint+1)*sizeof(int)))==ibgNULL)
        	       				{goto overflow;}
        ibgg0.nq[d] = (int*)adr;
 }
 for(d=0;d<ibgg0.maxd;d++){
        if((adr=realloc(ibgg0.ndp[d],(ibgg.maxPoint+1)*sizeof(int)))==ibgNULL)
        	       				{goto overflow;}
        ibgg0.ndp[d] = (int*)adr;
        ibgg0.ndr[d] = (int*)adr;
        if((adr=realloc(ibgg0.ndn[d],(ibgg.maxPoint+1)*sizeof(int)))==ibgNULL)
               	       				{goto overflow;}
        ibgg0.ndn[d] = (int*)adr;
        ibgg0.ndr[ibgg.gridDimension+d] = (int*)adr;
 }
 if(ibgg0.no!=ibgNULL){
        if((adr=realloc(ibgg0.no,(ibgg.maxPoint+1)*sizeof(int)))==ibgNULL)
        	       	       				{goto overflow;}
        ibgg0.no = (int*)adr;
 }
 if(ibgg0.nb!=ibgNULL){
        if((adr=realloc(ibgg0.nb,(ibgg.maxPoint+1)*sizeof(int)))==ibgNULL)
        	       	       				{goto overflow;}
        ibgg0.nb = (int*)adr;
 }
 return 1;
overflow:
 ibgfatal;
 return 0;
}

static int g0BAlloc(int max)
{void *adr;
 if((adr=malloc(max*sizeof(ibgPoint)))==ibgNULL)	{goto overflow;}
 ibgg0.bPoint = (ibgPoint *) adr;
 ibgg0.bx = (char *) &(((ibgPoint *)adr)->xx[0]);
 ibgg0.bt = (char *) &(((ibgPoint *)adr)->tt);
 ibgg0.bu = (char *) &(((ibgPoint *)adr)->uu);
 ibgg0.bf = (char *) &(((ibgPoint *)adr)->ff[0]);
 if((adr=malloc(max*sizeof(int)))==ibgNULL)	{goto overflow;}
 ibgg0.bn = (int *) adr;
 if((adr=malloc(max*sizeof(int)))==ibgNULL)	{goto overflow;}
 ibgg0.bstack = (int *) adr;
 if((adr=calloc((ibgg.maxPoint+1),sizeof(int)))==ibgNULL){goto overflow;}
 ibgg0.nb = (int *) adr;
 ibgg0.bfree	= -1;
 ibgg0.lastBound	= 0;
 ibgg0.maxBound	= max-1;
 return 1;
overflow:
 ibgfatal;
 return 0;
}

static int g0BRealloc()
{void *adr; int max;
 if(ibgg0.lastBound<ibgg0.maxBound) return 1;
 max=3*ibgg0.lastBound/2;
 if((adr=realloc(ibgg0.bPoint,max*sizeof(ibgPoint)))==ibgNULL)	{goto overflow;}
 ibgg0.bPoint = (ibgPoint *) adr;
 ibgg0.bx = (char *) &(((ibgPoint *)adr)->xx[0]);
 ibgg0.bt = (char *) &(((ibgPoint *)adr)->tt);
 ibgg0.bu = (char *) &(((ibgPoint *)adr)->uu);
 ibgg0.bf = (char *) &(((ibgPoint *)adr)->ff[0]);
 if((adr=realloc(ibgg0.bn,max*sizeof(int)))==ibgNULL)	{goto overflow;}
 ibgg0.bn = (int *) adr;
 if((adr=realloc(ibgg0.bstack,max*sizeof(int)))==ibgNULL)	{goto overflow;}
 ibgg0.bstack = (int *) adr;
 ibgg0.maxBound=max-1;
 return 1;
overflow:
 ibgfatal;
 return 0;
}

int g0QFree()
{int q,d;
 if(ibgg0.maxq==0) return 0;
 ibgg0.maxq = 0;
 for(q=0;q<ibgg0.maxq;q++){
        if(ibgg0.qn[q] != ibgNULL)	free(ibgg0.qn[q]);
        if(ibgg0.nq[q] != ibgNULL)	free(ibgg0.nq[q]);
 }
 for(d=0;d<ibgg0.maxd;d++){
        if(ibgg0.qdp[d] != ibgNULL)	free(ibgg0.qdp[d]);
        if(ibgg0.qdn[d] != ibgNULL)	free(ibgg0.qdn[d]);
        if(ibgg0.ndp[d] != ibgNULL)	free(ibgg0.ndp[d]);
        if(ibgg0.ndn[d] != ibgNULL)	free(ibgg0.ndn[d]);
 }
 if(ibgg0.qn0 != ibgNULL)		free(ibgg0.qn0);
 ibgg0.maxd = 0;
 ibgg0.maxq = 0;
 return 1;
}

static int g0QAlloc(int last, int max)
{void *adr; int q,d;
 ibgg0.lastQuad = 0;
 ibgg0.maxQuad  = 0;
 if(ibgg0.maxq==0) return 1;
 if(last >= max) max = last+1;
 max++;
 for(q=0;q<ibgg0.maxq;q++){
        if((adr=malloc(max*sizeof(int)))==ibgNULL)
        	{g0QFree(); goto overflow; }
        ibgg0.qn[q] = (int *) adr;
        if((adr=malloc(max*sizeof(int)))==ibgNULL)
        	{g0QFree(); goto overflow; }
        ibgg0.nq[q] = (int *) adr;
 }
/*
 for(d=0;d<ibgg.gridDimension;d++){
        if((adr=malloc(max*sizeof(int)))==ibgNULL)
        	{g0QFree(); goto overflow; }
        ibgg0.qdp[d] = adr;
        ibgg0.qdr[d] = adr;
        if((adr=malloc(max*sizeof(int)))==ibgNULL)
        	{g0QFree(); goto overflow; }
        ibgg0.qdn[d] = adr;
        ibgg0.qdr[ibgg.gridDimension+d] = adr;
 }
 if((adr=malloc(max*sizeof(int)))==ibgNULL)
        	{g0QFree(); goto overflow; }
 ibgg0.qn0 = adr;
*/
 ibgg0.lastQuad = last;
 ibgg0.maxQuad  = max-1;
 return 1;
overflow:
 ibgfatal;
 return 0;
}

static int g0QRealloc()
{void *adr; int q,d,max;
 if(ibgg0.lastQuad<ibgg0.maxQuad) return 1;
 max = 3*ibgg0.lastQuad/2;
 for(q=0;q<ibgg0.maxq;q++){
        if((adr=realloc(ibgg0.qn[q],max*sizeof(int)))==ibgNULL){goto overflow;}
        ibgg0.qn[q] = (int *) adr;
 }
/*
 for(d=0;d<ibgg.gridDimension;d++){
        if((adr=realloc(ibgg0.qdp[d],max*sizeof(int)))==ibgNULL){goto overflow;}
        ibgg0.qdp[d] = adr;
        ibgg0.qdr[d] = adr;
        if((adr=realloc(ibgg0.qdn[d],max*sizeof(int)))==ibgNULL){goto overflow;}
        ibgg0.qdn[d] = adr;
        ibgg0.qdr[ibgg.gridDimension+d] = adr;
 }
 if((adr=realloc(ibgg0.qn0,max*sizeof(int)))==ibgNULL)	{goto overflow;}
 ibgg0.qn0 = adr;
*/
 ibgg0.maxQuad=max-1;
 return 1;
overflow:
 ibgfatal;
 return 0;
}

static int g0BNew0()
{int b;
 if((b=ibgg0.bfree)>0){
        ibgassert(cbn(b)== -100000);
        ibgg0.bfree = ibgg0.bstack[b];
        ibgassert((ibgg0.bfree>0)?(cbn(ibgg0.bfree)== -100000):1);
        return b;
 }
 if(ibgg0.lastBound < ibgg0.maxBound){
        cbn(b = ++ibgg0.lastBound) = 0;
        return b;
 }
 if(g0BRealloc()){
        cbn(b = ++ibgg0.lastBound) = 0;
        return b;
 }
 return 0;
}

int g0NNew0()
{int n;
 if((n=ibgg0.nfree)>0){
        ibgg0.nfree = ibgg0.nstack[n]; return n;
 }
 if(ibgg.lastPoint < ibgg.maxPoint){
        return ++ibgg.lastPoint;
 }
 if(g0NRealloc()){
        return ++ibgg.lastPoint;
 }
 return 0;
}

/* shift of point coordinates to boundary coordinates, free of boundary data */
int g0BNShift()
{int nn,bb;
 ibgPoint nd;
 ibgSegmentType t;
/* Speicherung der boundary-Koordinaten auf den Ursprungsknoten: */
 for(nn=1;nn<=ibgg.lastPoint;nn++) if(cnb(nn)!=0){
        bb = cnb(nn);
        *cnd(nn) = *cbd(bb);
 }
 free(ibgg0.nb); free(ibgg0.bn); free(ibgg0.bPoint);