ibgg.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 19:38:50.25	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
#define MODULENAME "ibgg"
#define PRINTLEVEL 2
#include "ibgg.h"
#include "ibglib.h"
#include "ibglimits.h"
#include "ibgc.h"
#include "ibgi.h"
#include "ibgd0.h"
#include "ibgdefault.h"
#include "ibgg0.h"
#include "ibgg0p.h"

#define cct(cc) ibgridCellType(ibgg,cc)
#define ccu(cc) ibgridCellSegment(ibgg,cc)
#define ccn(cc) ibgridCPointList(ibgg,cc)
#define ccs(cc)	ibgridCSideList(ibgg,cc)
#define ccdef(cc)	(ccu(cc)> 0)
#define ccout(cc)	(ccu(cc)==0)

#define ibggDelaunayNew 1
ibgPoints       ibggExternalPoints;
ibgPoints       pp;
int		ibggFace;
int		ibggLine;
int		ibggNode;
ibgFloat	ibggDmin[ibgDIM];
ibgFloat	ibggBLWait,ibggBLHalf,ibggBLObtuse,ibggBctg;
ibgFloat	ibggOref,ibggOinv,ibggOref2,ibggOmax;
ibgFloat	ibggBV;
ibgFloat	ibggPref;
int   		ibggPrefL;

ibGrid		ibgg;

ibGrid		ibgg0oldgrid;

int ibgridNRealloc(ibGrid *g)
{void *adr; int max;
 if(g->lastPoint<g->maxPoint)	return 1;
 ibgg = *g;
 max=3*ibgg.lastPoint/2 + 100;
 if((adr=realloc(ibgg.Point,max*sizeof(ibgPoint)))==ibgNULL){goto noverflow;}
 ibgg.Point = (ibgPoint *) adr;
 if(ibgg.PointCell!=ibgNULL){
 	if((adr=realloc(ibgg.PointCell,max*sizeof(int)))==ibgNULL)
	 	{goto noverflow;}
 	ibgg.PointCell = (int*)adr;
 }
 if(ibgridHasPointData(ibgg)){
	if((adr=realloc(ibgg.PointData,max*(ibgg.NDataSize)))==ibgNULL)
	 	{goto noverflow;}
	ibgg.PointData = adr;
 }
 ibgg.maxPoint=max-1;
 *g = ibgg;
 return 1;
noverflow:
 return 0;
}

int ibgridCNew(ibGrid *g, ibgCellType t)
{void *adr; int max,l;
 int c1;
 if(g->lastCell>=g->maxCell){
	max=3*g->lastCell/2 + 100;
	if((adr=realloc(g->CellType,max*sizeof(ibgCellType)))==ibgNULL)
		 	{g->lastCell--; ibgfatal; return 0;}
	g->CellType=(ibgCellType *)adr;
	if((adr=realloc(g->CellSegment,max*sizeof(int)))==ibgNULL)
		 	{g->lastCell--; ibgfatal; return 0;}
	g->CellSegment=(ibgSegment*)adr;
	if((adr=realloc(g->CPointFirst,max*sizeof(int)))==ibgNULL)
		 	{g->lastCell--; ibgfatal; return 0;}
	g->CPointFirst=(int*)adr;
	if((adr=realloc(g->CSideFirst,max*sizeof(int)))==ibgNULL)
		 	{g->lastCell--; ibgfatal; return 0;}
	g->CSideFirst=(int*)adr;
	if(ibgridHasCellData(*g)){
		if((adr=realloc(g->CellData,max*(g->CDataSize)))==ibgNULL)
			 	{g->lastCell--; ibgfatal; return 0;}
		g->CellData=adr;
	}
	g->maxCell = max-1;
 }
 c1 = ++(g->lastCell);
 (g->CPointFirst)[c1] = g->freeCPoint;
 g->freeCPoint += ibgcPoints(t);
 if(g->freeCPoint>=g->maxCPoint){
 	max=3*g->freeCPoint/2 + 100;
 	if((adr=realloc(g->CPoint,max*sizeof(int)))==ibgNULL)
 		{g->lastCell--; g->freeCPoint -= ibgcPoints(t);
 		 ibgfatal; return 0;}
  	g->CPoint=(int*)adr; g->maxCPoint=max-IBGCNMAX-1;
 }
 (g->CSideFirst)[c1] = g->freeCSide;
 switch(ibgcCODIM(t)){
 	case ibgSRegion:	g->freeCSide += (l=ibgcSides(t)); 	break;
 	case ibgSFace:	g->freeCSide += (l=ibgcSides(t) + 2);	break;
 	case ibgSLine:	g->freeCSide += (l=ibgcSides(t) + 1);	break;
 	case ibgSNode:	g->freeCSide += (l=ibgcSides(t) + 1);	break;
 }
 if(g->freeCSide>=g->maxCSide){
 	max=3*g->freeCSide/2 + 100;
 	if((adr=realloc(g->CSide,max*sizeof(int)))==ibgNULL)
 		{g->lastCell--; g->freeCPoint -= ibgcPoints(t);
 		 g->freeCSide -= l; ibgfatal; return 0;}
 	g->CSide=(int*)adr; g->maxCSide=max-IBGCSMAX-1;
 }
 g->CellType[c1] = t;
 g->CellSegment[c1] = 0;
 return c1;
}

int ibgridCORealloc(ibGrid *gg)
{void *adr,*adrn,*adrs; int max;
 if(gg->lastCOut>=gg->maxCOut){
	max=3*gg->lastCOut/2 + 100;
	if(gg->maxCOut==0){
		if((adrn=malloc(max*sizeof(int)))==ibgNULL)
			       				{ibgfatal;return 0;}
		if((adrs=malloc(max*sizeof(int)))==ibgNULL)
			       				{ibgfatal;return 0;}
	}else{
		if((adrn=realloc(gg->COutNext,max*sizeof(int)))==ibgNULL)
			       				{ibgfatal;return 0;}
		if((adrs=realloc(gg->COutside,max*sizeof(int)))==ibgNULL)
			       		     		{ibgfatal;return 0;}
	}
 	gg->COutNext = (int *) adrn;
 	gg->COutside = (int *) adrs;
      	gg->maxCOut  = max-1;
 }
 return 1;
}

static int er1,er2;

static int gcellFace(int c, ibgSegment m0, ibgSegment m1)
{int bl,br,dl,dr,i,*cnn=ccn(c),b0,b1,first=1;
 ibgCellType t=cct(c);
 ibgPoint *n;
 b0 = ibgDefaultFaceNr(m0,m1);
 if(m1==ibgOutsideRegionNr) return b0;
 for(i=0;i<ibgcPoints(t);i++){
      	n = &(ibgridPoint(ibgg,cnn[i]));
	if(ibgpType(*n)==ibgSRegion){
			/* innerer Knoten auf dem Rand */
		er1++; /* kleine Fehlermeldung */
		return ibgError;
	}else if(ibgpType(*n)==ibgSFace){
		if((b1=ibgpSegment(*n))!=b0){
			if(first==0) {
			 	er2++;/* kleine Fehlermeldung */
				return ibgError;
			}
			if(b1>=ibgSMAX){
			 	er2++; /* kleine Fehlermeldung */
				return ibgError;
			}
			b0=b1;
			first=0; continue;
		}
	}
 }
 return b0;
}

static int gcellSetFace(int c, ibgSegment m0, ibgSegment m1)
{ibgCellType t=cct(c);
 return	ibgDefaultFaceNr(m0,m1);
}

void ibgridMakeFaces(ibGrid *gg, ibGeometry g)
{int c0,co,c1,c2,cl,cl0,cl1,cl2,cs,cn,m,m0,m1,ni,n0,s,i,j,k,b0;
 int *cs0,*cs1,*cs2,*cn0,*cno,*cln,*css,*csn,*cnn,*cso,*pn,bs0,bs1,bso;
 ibgCellType t,ts,tn,tl,t0,t1,t2;
 int s0,s1,s2,p0,p1,p2;
 int bmindef=ibgDefaultFaceNr(1,0),bmax=0,*adr;
 ibgg = *gg;
 ibGeo = g;
 if(ibgg.gridDimension<1){
 	ibgg.firstFace=1; ibgg.lastFace=0;
	return;
 }
 er1 = er2 = 0;
 ibgg.firstFace= ibgg.lastCell+1;
/* creating ibgSFace-elements, their points, left and right neighbours: */
 cl0=ibgg.lastCell;
 for(c0=1;c0<=cl0;c0++) if(ccdef(c0)){
	if((m=ccu(c0))<=0) continue;
	t=cct(c0);
	for(s=0;s<ibgcSides(t);s++){
		cs=ccs(c0)[s];
		if(ccout(cs)) m1=ibgOutsideRegionNr; /* Au"senrand */
		else if(ibgcCODIM(cct(cs))!=ibgSRegion)	continue;
		else if(m<=(m1=ccu(cs)))		continue;
		switch(ibgcSSize(ts=cct(cs),s)){
			case 1:	tn=ibgc1Node;		break;
			case 2:	tn=ibgc2Edge;		break;
			case 3: tn=ibgc3Triangle;	break;
			case 4: tn=ibgc3Rectangle;	break;
		}
		ibgassert(ibgcPoints(tn)==ibgcSSize(t,s));
		cn = ibgridCNew(&ibgg,tn);
 		cct(cn)=tn;
		cnn=ccn(cn); csn=ccs(cn); css=ccs(cs); cn0=ccn(c0);
		for(i=0;i<ibgcPoints(tn);i++){
		      	cnn[i] = cn0[ibgcSPoint1(t,s,i)];
		}
		for(i=0;i<ibgcSides(tn);i++) csn[i]=0;
		csn[ibgcSExtR(tn)]=c0;
		csn[ibgcSExtL(tn)]=cs;
		ccs(c0)[s]=cn;
		for(i=0;i<ibgcSides(ts);i++) if(css[i]==c0) {css[i]=cn; break;}
		b0 = gcellFace(cn,m,m1);
		if(b0<0){ccu(cn)=b0=gcellSetFace(cn,m,m1);}
		else ccu(cn)=b0;
	}
 }
 ibgg.lastFace=ibgg.lastCell;
 *gg = ibgg;
}

void ibgridMakeLines(ibGrid *gg, ibGeometry g)
{int c0,co,c1,c2,cl,cl1,cl2,cs,cn,m,m0,ni,n0,s,i,j,k,b0;
 int *cs0,*cs1,*cs2,*cn0,*cno,*cln,*css,*csn,*cnn,*cso,*pn,bs0,bs1,bso;
 ibgCellType t,ts,tn,tl,t0,t1,t2;
 int s0,s1,s2,p0,p1,p2,first,bl,br,dl,dr;
 int bmindef=ibgDefaultFaceNr(1,0),bmax=0,*adr;
 if(ibgg.gridDimension<2){
	ibgg.firstFace=1; ibgg.lastFace=0;
	return;
 }
 ibgg  = *gg;
 ibGeo = g;
 ibgg.firstLine=ibgg.lastCell+1;
 for(c0=ibgg.firstFace;c0<=ibgg.lastFace;c0++){t0=cct(c0);
   for(s0=0;s0<ibgcSides(t0);s0++){
	if((cs0=ccs(c0))[s0]!=0) continue;
	n0=(cn0=ccn(c0))[ibgcSPoint1(t0,s0,0)];
/* search of the first (after c0) ibgSFace-element c1: */
	co=c0; c1=cs0[ibgcSExtR(t0)];
	while(ibgcCODIM(t1=cct(c1))==ibgSRegion){
		cso=ccs(c1);cno=ccn(c1);
	 	for(s1=0;s1<ibgcSides(t1);s1++) if(cso[s1]==co) break;
		ibgassert(s1<ibgcSides(t1));
		for(p1=0;p1<ibgcSSize(t1,s1);p1++)
			if(cno[ibgcSPoint1(t1,s1,p1)]==n0) break;
	    	co=c1;c1=cso[ibgcSSide(t1,s1,p1)];
	}
/* search of orientation in c1 */
	if(co==(cs1=ccs(c1))[ibgcSExtL(t1)]){
		for(s1=0;s1<ibgcSides(t1);s1++)
			if(ccn(c1)[ibgcSPoint1(t1,s1,0)]==n0) break;
		ibgassert(s1<ibgcSides(t1));
		co=c1;c2=cs1[ibgcSExtR(t1)];
	}else if(co==cs1[ibgcSExtR(t1)]){
		for(s1=0;s1<ibgcSides(t1);s1++)
			if(ccn(c1)[ibgcSPoint1(t1,s1,1)]==n0) break;
		ibgassert(s1<ibgcSides(t1));
	      	co=c1;c2=cs1[ibgcSExtL(t1)];
	}else ibgfatal;
	ibgassert(c0!=c1);
	first=1;
	do{
/* search of the next ibgSFace-element c2: */
    	    while(ibgcCODIM(t2=cct(c2))==ibgSRegion){
	    	cso=ccs(c2);cno=ccn(c2);
	 	for(s2=0;s2<ibgcSides(t2);s2++) if(cso[s2]==co) break;
		ibgassert(s2<ibgcSides(t2));
		for(p2=0;p2<ibgcSSize(t2,s2);p2++)
			if(ccn(c2)[ibgcSPoint1(t2,s2,p2)]==n0) break;
		co=c2;c2=cso[ibgcSSide(t2,s2,p2)];
	    }
/* search of orientation in c2 */
	    if(co==ccs(c2)[ibgcSExtL(t2)]){
	    	for(s2=0;s2<ibgcSides(t2);s2++)
			if(ccn(c2)[ibgcSPoint1(t2,s2,0)]==n0) break;
		ibgassert(s2<ibgcSides(t2));
		co=c2;c2=ccs(c2)[ibgcSExtR(t2)];
	    }else if(co==ccs(c2)[ibgcSExtR(t2)]){
	    	for(s2=0;s2<ibgcSides(t2);s2++)
			if(ccn(c2)[ibgcSPoint1(t2,s2,1)]==n0) break;
		ibgassert(s2<ibgcSides(t2));
	      	co=c2;c2=ccs(c2)[ibgcSExtL(t2)];
	    }else ibgfatal;
	    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 ibgSLine-element: */
		if(ibgg.gridDimension==3) tl=ibgc3Edge; else tl=ibgc2Node;
		cl  = ibgridCNew(&ibgg,tl);
		cln = ccn(cl); cct(cl) = tl;
		m0  = -1;
		for(i=0;i<ibgcPoints(tl);i++){
		      	cln[i] = ni = ccn(c0)[ibgcSPoint1(t0,s0,i)];
			if(ibgpType(ibgridPoint(ibgg,ni))==ibgSLine){
				if(m0== -1){
					m0 = ibgpSegment(ibgridPoint(ibgg,ni));
				}else if(ibgpSegment(ibgridPoint(ibgg,ni))!=m0){
					m0 = ibgDefaultLineNr;
				}
			}
		}
		if(m0== -1) m0=ibgDefaultLineNr;
		ccu(cl) = m0;
/* default space is only for one external neighbour ibgcSExtF */
  		bs0 = ibgridCONew(&ibgg); if(bs0==0) return;
  		bs1 = ibgridCONew(&ibgg); if(bs1==0) return;
		ccs(cl)[ibgcSExtF(tl)] = bs0;
		ibgg.COutside[bs0] = c0; ccs(c0)[s0] = cl;
		ibgg.COutside[bs1] = c1; ccs(c1)[s1] = cl;
		ibgg.COutNext[bs0] = bs1;
	    }
	    if(co==c0) {ibgg.COutNext[bs1] = bs0; break;}
	    bso = bs1; bs1 = ibgridCONew(&ibgg); if(bs1==0) return;
	    ibgg.COutNext[bso] = bs1;
	    ibgg.COutside[bs1] = co; ccs(co)[s2] = cl;
	}while(1);
nexts0:	;
    }
 }
 ibgg.lastLine=ibgg.lastCell;
 *gg = ibgg;
}
void ibgridMakeNodes(ibGrid *gg, ibGeometry g)
{int c0,co,c1,c2,cl,cl0,cl2,cs,cn,m,m0,ni,n0,s,i,j,k,b0;
 int *cs0,*cs1,*cs2,*cn0,*cno,*cln,*css,*csn,*cnn,*cso,*pn,bs0,bs1,bso;
 ibgCellType t,ts,tn,tl,t0,t1,t2;
 int s0,s1,s2,p0,p1,p2,first,bl,br,dl,dr;
 int bmindef=ibgDefaultFaceNr(1,0),bmax=0,*adr;
 if(ibgg.gridDimension<3){
	ibgg.firstNode=1; ibgg.lastNode=0;
	return;