ibgg.c

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

 }
 ibgg  = *gg;
 ibGeo = g;
 ibgg.firstNode=ibgg.lastCell+1;
 pn=(int*)calloc(ibgg.lastPoint+1,sizeof(int));
 cl2=ibgg.lastLine;
 for(c0=ibgg.firstLine;c0<=cl2;c0++){t0=cct(c0);
     for(s0=0;s0<ibgcPoints(t0);s0++){
	n0=ccn(c0)[s0];
	if(pn[n0]==0){pn[n0]= -c0; continue;}
	if(pn[n0]< 0){
		c1 = -pn[n0];
		for(s1=0;s1<2;s1++) if(ccn(c1)[s1]==n0) break;
		ibgassert(s1<2);
		ccs(c1)[s1] = c0; ccs(c0)[s0] = c1;
		pn[n0] = c0; continue;
	}
	c1 = pn[n0]; t1=cct(c1);
	if(ibgcCODIM(t1)==ibgSLine){
		cl = ibgridCNew(&ibgg,ibgc3Node); if(cl==0) return;
		ccn(cl)[0]=n0;
		bs0 = ibgridCONew(&ibgg); if(bs0==0) return;
		ccs(cl)[ibgcSExtF(ibgc3Node)] = bs0;
		ibgg.COutside[bs0] = c0;
		ibgg.COutNext[bs0] = bs0;
		for(s1=0;s1<2;s1++) if(ccn(c1)[s1]==n0) break;
		ibgassert(s1<2);
		c2=ccs(c1)[s1];
		for(s2=0;s2<2;s2++) if(ccn(c2)[s2]==n0) break;
		ibgassert(s2<2);
		ccs(c0)[s0]=cl; ccs(c1)[s1]=cl; ccs(c2)[s2]=cl;
		pn[n0] = cl;
	}else if(ibgcCODIM(t1)==ibgSNode){
		bs0 = ibgridCONew(&ibgg); if(bs0==0) return;
		ccs(c0)[s0]=c1;
		bso = ccs(cl)[ibgcSExtF(ibgc3Node)];
		ibgg.COutside[bs0] = c0;
		ibgg.COutNext[bs0] = ibgg.COutNext[bso];
		ibgg.COutNext[bso] = bs0;
	}else ibgfatal;
     }
 }
 for(n0=1;n0<=ibgg.lastPoint;n0++) if(pn[n0]<0) {
		cl = ibgridCNew(&ibgg,ibgc3Node); if(cl==0) return;
		ccn(cl)[0]=n0;
		c0 = -pn[n0];
		if     (ccn(c0)[0]==n0) ccs(c0)[0]=cl;
		else if(ccn(c0)[1]==n0) ccs(c0)[1]=cl;
		else ibgfatal;
		ccs(cl)[ibgcSExtF(ibgc3Node)]=c0;
		bs0 = ibgridCONew(&ibgg); if(bs0==0) return;
		ccs(cl)[ibgcSExtF(ibgc3Node)] = bs0;
		ibgg.COutside[bs0] = c0;
		ibgg.COutNext[bs0] = bs0;
 }
 free(pn);
/* Test des Gitters auf Korrektheit */
 for(c0=1;c0<=ibgg.lastCell;c0++) if(ccdef(c0)) {t0=cct(c0);
  	for(s0=0;s0<ibgcSides(t0);s0++){
		cs=ccs(c0)[s0];
		ibgassert(cs>0);
		ibgassert(ibgcCODIM(t0)<=ibgcCODIM(cct(cs)));
		if(ibgcCODIM(t0)==ibgcCODIM(cct(cs))){
			ibgassert(ccu(c0)==ccu(cs));
		}else if(ibgcCODIM(cct(cs))==ibgSFace){
		}
	}
 }
 ibgg.lastNode = ibgg.lastCell;
 *gg = ibgg;
}


void ibGridFunction(ibGrid *gg, int fnum,
	 	 	ibgFloat	(*func)(ibgPtObject fThis, ibgPoint *x),
			ibgPtObject	fThis)
{int nn;
 for(nn=ibgridFirstPoint(*gg);nn<=ibgridLastPoint(*gg);nn++){
	ibgpF(ibgridPoint(*gg,nn))[fnum] = func(fThis,&ibgridPoint(*gg,nn));
 }
}

ibGrid *ibGridCopy(ibGrid *grid)
{int i;
 ibGrid *gg=(ibGrid*)malloc(sizeof(ibGrid));
 *gg = ibgg = *grid;
 gg->Point	= (ibgPoint*)malloc(IBGPOINTSIZE*(gg->maxPoint+1));
 gg->CellType	= (ibgCellType*)malloc(sizeof(ibgCellType)*(gg->maxCell+1));
 gg->CellSegment	= (ibgSegment*)malloc(sizeof(int)*(gg->maxCell+1));
 gg->CPointFirst	= (int*)malloc(sizeof(int)*(gg->maxCell+1));
 gg->CSideFirst	= (int*)malloc(sizeof(int)*(gg->maxCell+1));
 gg->CPoint	= (int*)malloc(sizeof(int)*(gg->maxCPoint+1));
 gg->CSide	= (int*)malloc(sizeof(int)*(gg->maxCSide+1));
 gg->COutNext	= (int*)malloc(sizeof(int)*(gg->maxCOut+1));
 gg->COutside	= (int*)malloc(sizeof(int)*(gg->maxCOut+1));
 for(i=gg->firstPoint;i<=gg->lastPoint;i++){
	gg->Point[i] = ibgg.Point[i];
 }
 for(i=gg->firstCell;i<=gg->lastCell;i++){
	gg->CellType[i]		= ibgg.CellType[i];
	gg->CellSegment[i]	= ibgg.CellSegment[i];
	gg->CPointFirst[i]	= ibgg.CPointFirst[i];
	gg->CSideFirst[i]	= ibgg.CSideFirst[i];
 }
 for(i=0;i<gg->freeCPoint;i++){
		gg->CPoint[i]	= ibgg.CPoint[i];
 }
 for(i=0;i<gg->freeCSide;i++){
		gg->CSide[i]	= ibgg.CSide[i];
 }
 for(i=0;i<=gg->lastCOut;i++){
		gg->COutside[i]	= ibgg.COutside[i];
		gg->COutNext[i]	= ibgg.COutNext[i];
 }
 return gg;
}

ibGrid *ibGridFree(ibGrid *gg)
{
 if((gg->Point)		!= ibgNULL) free(gg->Point);
 if((gg->CellType)	!= ibgNULL) free(gg->CellType);
 if((gg->CellSegment)	!= ibgNULL) free(gg->CellSegment);
 if((gg->CPointFirst)	!= ibgNULL) free(gg->CPointFirst);
 if((gg->CSideFirst)	!= ibgNULL) free(gg->CSideFirst);
 if((gg->CPoint)		!= ibgNULL) free(gg->CPoint);
 if((gg->CSide)		!= ibgNULL) free(gg->CSide);
 if((gg->COutNext)	!= ibgNULL) free(gg->COutNext);
 if((gg->COutside)	!= ibgNULL) free(gg->COutside);
 if((gg->CellData)	!= ibgNULL) free(gg->CellData);
 if((gg->PointData)	!= ibgNULL) free(gg->PointData);
 if((gg->PointCell)	!= ibgNULL) free(gg->PointCell);
 free(gg);
 return ibgNULL;
}

ibggRefPar ibggRef;

static ibgFloat	ibggDminOld[3];

ibGrid	*ibGridGenerate(ibGeometry g,
	       	unsigned maxPoint,
      		unsigned nx, ibgFloat *x, unsigned fx,
		unsigned ny, ibgFloat *y, unsigned fy,
	    	unsigned nz, ibgFloat *z, unsigned fz,
	     	void (*RefineRegion)(ibgPtObject This, ibgPoint *point,
			   	      	ibgFloat *length),
		ibgPtObject pRefineRegion,
		void (*RefineFace)(ibgPtObject This, ibgPoint *point,
				      	ibgFloat *normal, ibgFloat *tangential),
	     	ibgPtObject pRefineFace,
	     	void (*RefineLine)(ibgPtObject This, ibgPoint *point,
			   	      	ibgFloat *normal, ibgFloat *tangential),
	     	ibgPtObject pRefineLine,
	     	void (*RefineNode)(ibgPtObject This, ibgPoint *point,
			   	      	ibgFloat *length),
	     	ibgPtObject pRefineNode,
		ibgBoolean (*RefineEdge) (ibgPtObject This,
					ibgPoint *n1, ibgPoint *n2),
		ibgPtObject pRefineEdge
	)
{ibGrid *gg; int i,rc;
 ibGrid0 *g0;
 clock_t time0=clock();
 clock_t timeo,timen;
 long int timed;
 ibgg0SetRefPar(&ibggRef,
		RefineRegion,pRefineRegion,
		RefineFace,pRefineFace,
		RefineLine,pRefineLine,
		RefineNode,pRefineNode,
		RefineEdge,pRefineEdge);
 if(nx<1)	nx=1;
 if(ny<1)	ny=1;
 if(nz<1)	nz=1;
/* if here is access violation the input is incorrect: */
 for(i=1;i<nx;i++) if(x[i]<=x[i-1]+g->Delta) goto incorrinput;
 for(i=1;i<ny;i++) if(y[i]<=y[i-1]+g->Delta) goto incorrinput;
 for(i=1;i<nz;i++) if(z[i]<=z[i-1]+g->Delta) goto incorrinput;
 if(nx==1){
	if(ny>1)	goto incorrinput;
	if(nz>1)	goto incorrinput;
 }else{
	if(g->Delta/(x[nx-1]-x[0]) < 10*IBG_FLT_EPSILON) goto deltatoosmall;
 }
 if(ny==1){
	if(nz>1)	goto incorrinput;
 }else{
	if(g->Delta/(y[ny-1]-y[0]) < 10*IBG_FLT_EPSILON) goto deltatoosmall;
 }
 if(nz==1){;
 }else{
	if(g->Delta/(z[nz-1]-z[0]) < 10*IBG_FLT_EPSILON) goto deltatoosmall;
 }
 ibgmsgreset;
 ibgmessage(ibgMIStart);
 ibgmsgreset;
 g0 = ibgg0CoarseGrid(maxPoint,nx,x,fx,ny,y,fy,nz,z,fz);
 if(g0==ibgNULL)	goto nogrid;
/*
 for(i=0;i<ibgDIM;i++){
	ibggDminOld[i] = ibggDmin[i];
	if(ibggDmin[i]>ibggDefaultLineNormal)ibggDmin[i]=ibggDefaultLineNormal;
 }
*/
 timed = (1000*(long int)((timeo=clock()) - time0))/CLOCKS_PER_SEC;
 ibgmsgreset;
 ibgmessage(ibgMIRefine);
 rc = ibgg0RegionOfPoint(g0,g);
 if(rc==ibgError){
	ibgg0Free(g0);
	goto nogrid;
 }

 rc = ibgg0Refine(g0,g,&ibggRef);
 if(rc==ibgError){
	goto nogrid;
 }
 timed = (1000*(long int)((timen=clock()) - timeo))/CLOCKS_PER_SEC;
 if(ibgMStatus == ibgMIRefine) ibgmsg(" %d ms;",timed);
 timeo=timen;
 ibgmsgreset;
 ibgmessage(ibgMIBShift);
 rc = ibgg0Boundary(g0,g,&ibggRef,maxPoint/2);
 if(rc==ibgError){
	goto nogrid;
 }
 timed = (1000*(long int)((timen=clock()) - timeo))/CLOCKS_PER_SEC;
 if(ibgMStatus == ibgMIBShift)	ibgmsg(" %d ms;",timed);
 timeo = timen;
 ibgmsgreset;
 ibgmessage(ibgMIDelaun);
 if (ibggDelaunayNew) {
     rc = ibgg0getPoints(&pp,g0);
     rc = ibgAppendPoints(&pp,&ibggExternalPoints);
     gg = (ibGrid*) malloc(sizeof(ibGrid));
     rc = ibggDelaunay(gg,g0->gg->gridDimension,&pp);
 }else{
     rc = ibgg0Delaunay(g0);
     if(rc==ibgError){
         goto nogrid;
     }
     timed = (1000*(long int)((timen=clock()) - timeo))/CLOCKS_PER_SEC;
     if(ibgMStatus == ibgMIDelaun)	ibgmsg(" %d ms;",timed);
     timeo = timen;
     ibgmsgreset;
     ibgmessage(ibgMIBGridG);
     ibgflush;
/*  gg = ibgg0Correct(g0); */
     gg = g0->gg;
     if(gg==ibgNULL){
         goto nogrid;
     }
 }
 ibgridMakeRegions(gg,g);
 ibgridMakeFaces(gg,g);
 ibgridMakeLines(gg,g);
 ibgridMakeNodes(gg,g);
/*
 for(i=0;i<ibgDIM;i++)	ibggDmin[i] = ibggDminOld[i];
*/
 timed = (1000*(long int)((timen=clock()) - timeo))/CLOCKS_PER_SEC;
 if(ibgMStatus == ibgMIBGridG)	ibgmsg(" %d ms;",timed);
 ibgmsgreset;
 ibgmessage(ibgMISuccess);
 timed = (1000*(long int)((timen=clock()) - time0))/CLOCKS_PER_SEC;
 ibgmsg(" with %d points, %d cells, in %d ms",
       	ibgg.lastPoint,ibgg.lastCell,timed);
 ibgg0oldgrid = *gg;
 ibgmsgreset;
 return gg;
incorrinput:
 ibgmessage(ibgMEIncCoarse);
 goto nogrid;
deltatoosmall:
 ibgmessage(ibgMESmallDelta);
 goto nogrid;
nogrid:
 ibgmsgreset;
 ibgmessage(ibgMINoGrid);
 return ibgNULL;
}

/* default refinement parameters and routines: */
ibgFloat ibggDefaultRegionLength,ibggDefaultNodeLength;
ibgFloat ibggDefaultFaceNormal,ibggDefaultFaceTangential;
ibgFloat ibggDefaultLineNormal,ibggDefaultLineTangential;

void ibggDefaultRefineRegion(ibgPtObject This, ibgPoint *point,
			   	      	ibgFloat *length)
{*length = ibggDefaultRegionLength;}

ibgBoolean ibggDefaultRefineEdge (ibgPtObject This, ibgPoint *n1, ibgPoint *n2)
{return ibgFalse;}

void ibggDefaultRefineFace(ibgPtObject This, ibgPoint *nface,
			   	      	ibgFloat *normal, ibgFloat *tangential)
{*normal	= ibggDefaultFaceNormal;
 *tangential  	= ibggDefaultFaceTangential;
}

void ibggDefaultRefineLine(ibgPtObject This, ibgPoint *nline,
			   	      	ibgFloat *normal, ibgFloat *tangential)
{*normal	= ibggDefaultLineNormal;
 *tangential  	= ibggDefaultLineTangential;
}

void ibggDefaultRefineNode(ibgPtObject This, ibgPoint *point,
			   	      	ibgFloat *length)
{*length = ibggDefaultNodeLength;}

void ibGridInit()
{ibgcInit();
/* initialization of parameters from the documentation: */
 ibggDmin[0]= 0.001;
 ibggDmin[1]= 0.001;
 ibggDmin[2]= 0.001;
 ibggDefaultRegionLength	= ibgInfty;
 ibggDefaultFaceNormal		= ibgInfty;
 ibggDefaultFaceTangential	= ibgInfty;
 ibggDefaultLineNormal		= 0.01;
 ibggDefaultLineTangential	= ibgInfty;
 ibggDefaultNodeLength	= ibgInfty;
/* initialization of other parameters: */
 ibggExternalPoints.Point = malloc(100*sizeof(ibgPoint));
 ibggExternalPoints.Near  = malloc(100*sizeof(int));
 ibggExternalPoints.maxPoint = 99;
 ibggExternalPoints.lastPoint = 0;
 ibggExternalPoints.firstPoint = 1;
 ibggFace=1;
 ibggLine=1;
 ibggNode=1;
 ibggBLWait=1.01;
 ibggBLHalf=2.01;
 ibggBLObtuse=1.01;
 ibggBctg=0.1;
 ibggOref=0.49;
 ibggOmax=0.0;
 ibggBV=1.01;	/* internal: 1 + rounding error */
 ibggPref=2.01;
 ibggPrefL=1;
/* defining dependend parameters: */
 ibggOinv=1.0/ibggOref;
 ibggOref2=ibggOref*ibggOref;
}