ibgdgrid.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 19:57:32.20	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< ibgdGrid >>> - Intersection-Based Geometry Definition - with a Grid */
#define const
#include "ibgd.h"
#include "ibgg.h"
#include "ibgi.h"
#include "ibglib.h"
#include "ibgd0.h"
#include "ibgdefault.h"
#include "ibgoutput.h"

#define eps2 1.e-12
#define eps3 1.e-18
#define ibgpC(point,d)	ibgpX(point)[d]
#define ibgpIC(point)	ibgpI(point)[0]
#define ibgPointSIZE	sizeof(ibgPoint)

typedef struct{
	ibGeometryRec	g;
	ibGrid		*Grid;
	int		fanz;
	int		*flist;
}ibgDGridRec,*ibgDGrid;

/* Auf Adressen sollte nur "uber folgende Aufrufe zugegriffen werden: */
#define gnx(gg,i)   ((ibgFloat *) (ibgpX(ibgridPoint(*gg,i))))
#define gnf(gg,i)   ((ibgFloat *) (ibgpF(ibgridPoint(*gg,i))))
#define def(n)   ((n)>0)
#define cdef(gg,c)   ((c)> 0 && ibgridCellSegment(gg,c)> 0)
#define cund(gg,c)   ((c)<=0 || ibgridCellSegment(gg,c)<=0)
#define und(n)   ((n)<=0)

static int sig[2] = {1,-1};

static int RegionOfPoint(ibGeometry g, ibgPoint *nnew, const ibgPoint *nold)
{ibgDGrid grdf = (ibgDGrid ) g;
 ibGrid *gg  = grdf->Grid;
 int *gn0;
 int c0,c1,ch,*cs0,cv0,nl,i,i0,sl,f,fn,co;
 ibgCellType t,th;
 ibgFloat *x=ibgpX(*nnew);
 double vv[4],xx,yy,zz,x1,x2,x3,y1,y2,y3,z1,z2,z3,fv,vvs;
 ibgFloat* ff[IBGCNMAX];
 register ibgFloat *fj;
 xx = x[0]; yy = x[1]; zz = x[2];
#define cvmax 3
 cv0=0;
 if(nold==ibgNULL || ibgpIC(*nold)==0){
	for(c0=1;;c0++)	if(ibgridCellSegment(*gg,c0)>0) break;
 }else	c0=ibgpIC(*nold);
 switch(ibgcCODIM(t=ibgridCellType(*gg,c0))){
	case ibgSNode:
		co=ibgridCOutFirst(*gg,c0,t);
		c0=ibgridCOutside(*gg,co);
		t =ibgridCellType(*gg,c0);	/* kein break! */
	case ibgSLine:
	      	co=ibgridCOutFirst(*gg,c0,t);
		c0=ibgridCOutside(*gg,co);
		t =ibgridCellType(*gg,c0);	/* kein break! */
	case ibgSFace:
		co=ibgridCellLeft(*gg,c0,t);
		if(cdef(*gg,co))	c0=co;
		else			c0=ibgridCellRight(*gg,c0,t);
		t = ibgridCellType(*gg,c0);
 }
/* look at the current cell: */
gridRegionBegin:
 ibgassert(cdef(*gg,c0));
 nl=ibgcPoints(t=ibgridCellType(*gg,c0)); gn0=ibgridCPointList(*gg,c0);
 for(i=0;i<nl;i++) ff[i]=gnx(gg,gn0[i]);
 sl=ibgcSides(t);
 ibgassert(ibgcCODIM(t)==ibgSRegion);
 i0 = 0;
gridRegionContinue:
 switch(ibgcGDIM(t)){
  case 1:
   for(i=i0;i<sl;i++){
	if((vv[i]=sig[i]*(xx-ff[i][0]))<0) {goto gridRegionNext;}
   } break;
  case 2:
   for(i=i0;i<sl;i++){
	fj = ff[ibgcEPoint1(t,i)];
	x1 = fj[0] - xx; y1 = fj[1] - yy;
	fj = ff[ibgcEPoint2(t,i)];
	x2 = fj[0] - xx; y2 = fj[1] - yy;
	vv[i] = x1 * y2 - x2 * y1;
	if((vv[i]<-eps2)||((vv[i]<0)&&(cv0++<cvmax))) goto gridRegionNext;
   } break;
  case 3:
   for(i=i0;i<sl;i++){
	fj = ff[ibgcSPoint1(t,i,0)];
	x1 = fj[0] - xx; y1 = fj[1] - yy; z1 = fj[2] - zz;
	fj = ff[ibgcSPoint1(t,i,1)];
	x2 = fj[0] - xx; y2 = fj[1] - yy; z2 = fj[2] - zz;
	fj = ff[ibgcSPoint1(t,i,2)];
	x3 = fj[0] - xx; y3 = fj[1] - yy; z3 = fj[2] - zz;
	vv[i] = -(x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1));
	if((vv[i]<-eps3)||((vv[i]<0)&&(cv0++<cvmax))) goto gridRegionNext;
   } break;
 }
 goto gridRegionEnd;
/*	go to the next cell: */
gridRegionNext:
 c1=ibgridCSideList(*gg,c0)[i];
 if(cund(*gg,c1)){
	i0 = i+1; goto gridRegionContinue;
 }else if(ibgcCODIM(th=ibgridCellType(*gg,c1))==ibgSFace){
	if(c0==(ch=ibgridCellRight(*gg,c1,th)))	c1=ibgridCellLeft(*gg,c1,th);
	else					c1=ch;
	if(cund(*gg,c1)){
		i0 = i+1; goto gridRegionContinue;
 	}
 }
 c0=c1;
 goto gridRegionBegin;
/*	interpolate function values and return: */
gridRegionEnd:
#ifdef ibgFloatValues
 vvs=0;
 for(i=0;i<sl;i++) vvs 	 += vv[i];
 for(i=0;i<sl;i++) vv[i] /= vvs;
 gn0=ibgridCPointList(*gg,c0);
 for(f=0;f<grdf->fanz;f++){
	fn=grdf->flist[f];
	fv=0;
	for(i=0;i<sl;i++) fv += vv[i]*(ibgpF(ibgridPoint(*gg,gn0[i]))[fn]);
	ibgpF(*nnew)[fn]=fv;
 }
#endif
 ibgpIC(*nnew)	 = c0;
 ibgpSegment(*nnew) = ibgridCellSegment(*gg,c0);
 ibgpType(*nnew) = ibgSRegion;
 return ibgRegionFound;
}

static int FaceWithEdge(ibGeometry g, ibgPoint *nint,
	const ibgPoint *n1, const ibgPoint *n2)
{ibgDGrid grdf = (ibgDGrid )g;
 ibGrid *gg = grdf->Grid;
 int c0,co,c1,ch,c2,cb,cr,mm,nl,i,ii,io,in,sl,*gn0,*gs0,mark[IBGCSMAX],f,fn;
 ibgCellType t,th;
 double vv0,vv1,vv2,vv3,vvv,p0,p1,p2,vv[3],vvs,fv;
 double xb,xe,xbo,xeo,xx0,xn1,xn2,xx3,x,y,z;
 double yb,ye,ybo,yeo,yy0,yy1,yy2,yy3;
 double zb,ze,zbo,zeo,zz0,zz1,zz2,zz3;
 ibgFloat*	ff[IBGCNMAX];
 ibgFloat	*f0,*f1,*f2;
 register ibgFloat *fj;
 int cut,inv=0;
 c0 = ibgpIC(*n1);
 switch(ibgcCODIM(t=ibgridCellType(*gg,c0))){
	case ibgSNode:
		co=ibgridCOutFirst(*gg,c0,t);
		c0=ibgridCOutside(*gg,co);
		t =ibgridCellType(*gg,c0);	/* kein break! */
	case ibgSLine:
		co=ibgridCOutFirst(*gg,c0,t);
		c0=ibgridCOutside(*gg,co);
		t =ibgridCellType(*gg,c0);	/* kein break! */
	case ibgSFace:
		co=ibgridCellLeft(*gg,c0,t);
		if(cdef(*gg,co))	c0=co;
		else			c0=ibgridCellRight(*gg,c0,t);
		t = ibgridCellType(*gg,c0);
 }
 c2 = ibgpIC(*n2);
 xbo = xb = ibgpC(*n1,0);
 ybo = yb = ibgpC(*n1,1);
 zbo = zb = ibgpC(*n1,2);
 xeo = xe = ibgpC(*n2,0);
 yeo = ye = ibgpC(*n2,1);
 zeo = ze = ibgpC(*n2,2);
begx:
 xx0= xb-xe;	yy0= yb-ye;	zz0= zb-ze;
begc0:
 ibgassert(cdef(*gg,c0));
 ibgassert(ibgcCODIM(ibgridCellType(*gg,c0))==ibgSRegion);
 nl=ibgcPoints(t=ibgridCellType(*gg,c0)); gn0=ibgridCPointList(*gg,c0);
 for(i=0;i<nl;i++) ff[i]=gnx(gg,gn0[i]);
 sl=ibgcSides(t); gs0=ibgridCSideList(*gg,c0);
 for(i=0;i<sl;i++){
	if((ch=gs0[i])>0){
	 	switch(ibgcCODIM(th=ibgridCellType(*gg,ch))){
		case ibgSFace:
			if((c1=ibgridCellLeft(*gg,ch,th))<=0){
					{mark[i]=1;break;}
		  	}else if(c1==c0){
				if((c1=ibgridCellRight(*gg,ch,th))<=0)
					{mark[i]=1;break;}
			}
			if(ibgridCellSegment(*gg,c1)<=0)
					{mark[i]=1;break;}
		case ibgSRegion:
		default:
					{mark[i]=0;break;}
		}
	}else{	mark[i]=0;
	}
 }
 switch(ibgcGDIM(t)){
  case 1:
	if((vv0=(xe - ff[0][0]))<0) {ii=0; goto nextc0;}
	if((vv0=(ff[1][0] - xe))<0) {vvv = xb - ff[1][0]; ii=1; goto nextc0;}
	break;
  case 2:
   for(i=0;i<sl;i++){ii=i; io= -1;
begi2:
	fj  = ff[ibgcSPoint1(t,ii,0)];
 	xn1 = fj[0] - xe; yy1 = fj[1] - ye;
	fj  = ff[ibgcSPoint1(t,ii,1)];
 	xn2 = fj[0] - xe; yy2 = fj[1] - ye;
	vv0 = xn1*yy2 - xn2*yy1;
	if(vv0>=0) {mark[ii]=1; 
		if(io>-1) {ii=io; io= -1; goto begi2;}
		else continue;}
	vv1 = xx0*yy2 - xn2*yy0;
	if(vv1>0) {in=ibgcSSide(t,ii,1);
	       	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begi2;}
	}
	vv2 = xn1*yy0 - xx0*yy1;
	if(vv2>0) {in=ibgcSSide(t,ii,0);
	       	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begi2;}
	}
	goto nextc0;
   } break;
  case 3:
   for(i=0;i<sl;i++){ii=i; io= -1;
begi3:
	fj  = ff[ibgcSPoint1(t,ii,0)];
 	xn1 = fj[0] - xe; yy1 = fj[1] - ye; zz1 = fj[2] - ze;
	fj  = ff[ibgcSPoint1(t,ii,1)];
 	xn2 = fj[0] - xe; yy2 = fj[1] - ye; zz2 = fj[2] - ze;
	fj  = ff[ibgcSPoint1(t,ii,2)];
 	xx3 = fj[0] - xe; yy3 = fj[1] - ye; zz3 = fj[2] - ze;
	vv0	= xn1*(yy2*zz3-yy3*zz2)
		+ xn2*(yy3*zz1-yy1*zz3)
		+ xx3*(yy1*zz2-yy2*zz1);
	if(vv0<=0) {mark[ii]=1; 
		if(io>-1) {ii=io; io= -1; goto begi3;}
		else continue;}
	vv1	= xx0*(yy2*zz3-yy3*zz2)
		+ xn2*(yy3*zz0-yy0*zz3)
		+ xx3*(yy0*zz2-yy2*zz0);
	if(vv1<0) {in=ibgcSSide(t,ii,1);
		if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begi3;}
	}
	vv2	= xn1*(yy0*zz3-yy3*zz0)
		+ xx0*(yy3*zz1-yy1*zz3)
		+ xx3*(yy1*zz0-yy0*zz1);
	if(vv2<0) {in=ibgcSSide(t,ii,2);
		if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begi3;}
	}
	vv3	= xn1*(yy2*zz0-yy0*zz2)
		+ xn2*(yy0*zz1-yy1*zz0)
		+ xx0*(yy1*zz2-yy2*zz1);
	if(vv3<0) {in=ibgcSSide(t,ii,0);
		if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begi3;}
	}
	goto nextc0;
   } break;
 }
 if(ibgridCellSegment(*gg,c0)==ibgridCellSegment(*gg,c2)){
	ibgpType(*nint) = ibgSNothing;
	return ibgNotFound;
 }
 if(inv){
	ibgpC(*nint,0) = ibgpC(*n2,0);
	ibgpC(*nint,1) = ibgpC(*n2,1);
	ibgpC(*nint,2) = ibgpC(*n2,2);
	ibgpIC(*nint)  = c0;
	ibgpSegment(*nint)= ibgDefaultFaceNr(ibgridCellSegment(*gg,c0),
					  ibgridCellSegment(*gg,c2));
	ibgpType(*nint)= ibgSFace;
	return ibgFaceFound;
 }
 inv=1; 
 xbo = xb = ibgpC(*n2,0);
 ybo = yb = ibgpC(*n2,1);
 zbo = zb = ibgpC(*n2,2);
 xeo = xe = ibgpC(*n1,0);
 yeo = ye = ibgpC(*n1,1);
 zeo = ze = ibgpC(*n1,2);
 c0 = ibgpIC(*n2);
 c2 = ibgpIC(*n1);
 goto begx;
nextc0:
 cb=ibgridCSideList(*gg,c0)[ii]; t =ibgridCellType(*gg,cb);
 if(ibgcCODIM(t)==ibgSRegion) {
	if(cb==c2) {return ibgNotFound;}
	c0=cb; goto begc0;
 }
 if(inv){
	if(c0==(cr=ibgridCellLeft(*gg,cb,t))) cr=ibgridCellRight(*gg,cb,t);
	else ibgassert(c0==ibgridCellRight(*gg,cb,t));
	if(ibgridCellSegment(*gg,cr)!=ibgridCellSegment(*gg,c2)) {c0=cr; goto begc0;}
 }
 c0=cb;
 switch(ibgcGDIM(t)){
  case 1:
	if(ii)	vvv = xb - ff[0][0];
	else	vvv = ff[1][0] - xb;
	if(vvv<0) vvv=0;
	break;
  case 2:
	vvv	= (xn1-xx0)*(yy2-yy0) - (xn2-xx0)*(yy1-yy0);
	if(vvv<0) vvv=0;
	break;
  case 3:
	yy1 -= yy0; yy2 -= yy0; yy3 -= yy0;
	zz1 -= zz0; zz2 -= zz0; zz3 -= zz0;
	vvv	= (xn1-xx0)*(yy2*zz3-yy3*zz2)
		+ (xn2-xx0)*(yy3*zz1-yy1*zz3)
  	  	+ (xx3-xx0)*(yy1*zz2-yy2*zz1);
	if(vvv>0) vvv=0;
	break;
 }
 p1=vv0/(vv0-vvv); p2=vvv/(vvv-vv0);
 ibgpC(*nint,0) = p1*xbo + p2*xeo;
 ibgpC(*nint,1) = p1*ybo + p2*yeo;
 ibgpC(*nint,2) = p1*zbo + p2*zeo;
 ibgpIC(*nint)	= c0;
#ifdef ibgFloatValues
 gn0=ibgridCPointList(*gg,c0); sl=ibgcPoints(t);
 switch(ibgcGDIM(t)){
  case 1: vv[0]=1.0;
	break;
  case 2: x = p1*xb + p2*xe; y = p1*yb + p2*ye;
	f0=gnx(gg,gn0[0]); f1=gnx(gg,gn0[1]);
	xx0=f1[0]-f0[0]; yy0=f1[1]-f0[1];
	vv[0]=(f1[0]-x)*xx0 + (f1[1]-y)*yy0;
	vv[1]=(x-f0[0])*xx0 + (y-f0[1])*yy0;
	vvs  = vv[0]+vv[1];
	vv[0] /= vvs; vv[1] /= vvs;
	break;
  case 3:  x = p1*xb + p2*xe; y = p1*yb + p2*ye; z = p1*zb + p2*ze;
	f0=gnx(gg,gn0[0]); f1=gnx(gg,gn0[1]); f2=gnx(gg,gn0[2]);
	xx0=(f1[1]-f0[1])*(f2[2]-f0[2]) - (f1[2]-f0[2])*(f2[1]-f0[1]);
	yy0=(f1[2]-f0[2])*(f2[0]-f0[0]) - (f1[0]-f0[0])*(f2[2]-f0[2]);
	zz0=(f1[0]-f0[0])*(f2[1]-f0[1]) - (f1[1]-f0[1])*(f2[0]-f0[0]);
	vv[0] = (f1[0]-x)*(yy0*(f2[2]-z) - zz0*(f2[1]-y))
	      + (f1[1]-y)*(zz0*(f2[0]-x) - xx0*(f2[2]-z))	
	      + (f1[2]-z)*(xx0*(f2[1]-y) - yy0*(f2[0]-x));	
	vv[1] = (f2[0]-x)*(yy0*(f0[2]-z) - zz0*(f0[1]-y))
	      + (f2[1]-y)*(zz0*(f0[0]-x) - xx0*(f0[2]-z))	
	      + (f2[2]-z)*(xx0*(f0[1]-y) - yy0*(f0[0]-x));	
	vv[2] = (f0[0]-x)*(yy0*(f1[2]-z) - zz0*(f1[1]-y))
	      + (f0[1]-y)*(zz0*(f1[0]-x) - xx0*(f1[2]-z))	
	      + (f0[2]-z)*(xx0*(f1[1]-y) - yy0*(f1[0]-x));	
	vvs   = vv[0]+vv[1]+vv[2];
	vv[0] /= vvs; vv[1] /= vvs; vv[2] /= vvs;
	break;
 }
 for(f=0;f<grdf->fanz;f++){
	fn=grdf->flist[f];
	fv=0;
	for(i=0;i<sl;i++) fv += vv[i]*(ibgpF(ibgridPoint(*gg,gn0[i]))[fn]);
	ibgpF(*nint)[fn]=fv;
 }
#endif
 ibgpSegment(*nint) = ibgridCellSegment(*gg,c0);
 ibgpType(*nint) = ibgSFace;
 return ibgFaceFound;
}

static int LineWithTriangle(ibGeometry g, ibgPoint *nint, ibgPoint *nface,
			ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{ibgDGrid grdf = (ibgDGrid )g;
 ibGrid *gg=grdf->Grid;
 int c0;
 goto gridLineError;
/*
 z23 = (x3-xp)*(y2-yp)-(y3-yp)*(x2-xp);
 if(z23<0){    goto} 	 Verlassen des Vierecks */
 /* drinne -> zum naechsten Nachbarn */

gridLineEnd:
 ibgpSegment(*nint) = ibgridCellSegment(*gg,c0);
 ibgpType(*nint) = ibgSLine;
 return ibgLineFound;
/* external return --- the standard routine: */
gridLineError:
 return ibgDefaultLineWithTriangle(g,nint,nface,n1,n2,n3);
}

static int Free(ibGeometry g)
{
 ibGridFree(((ibgDGrid )g)->Grid);
 return ibgSuccess;
}

static int stdlist[ibgFloatValues];
static ibGeometryClassRec	ibgDGridClass;

ibGeometry ibgdGrid(ibGrid *gg, int fanz, int *flist)
{ibgDGrid grdf; int i;
 if(gg->lastRegion==0) return ibgdBoundaryGrid(gg,ibgRSpace);
 grdf = (ibgDGrid) malloc(sizeof(ibgDGridRec));
 grdf->Grid = gg;
 grdf->fanz = fanz;
 if(flist==ibgNULL){
	ibgassert(fanz<=ibgFloatValues);
  	for(i=0;i<fanz;i++) stdlist[i]=i;
	grdf->flist = stdlist;
 }else{
	grdf->flist = flist;
	for(i=0;i<fanz;i++) ibgassert(flist[i]<ibgFloatValues);
 }
 ibgdInitialize((ibGeometry)grdf,&ibgDGridClass);
 ibgtCopy(&(grdf->g.top),&(gg->Topology));
 ibgDGridClass.RegionOfPoint = RegionOfPoint;
 ibgDGridClass.FaceWithEdge = FaceWithEdge;
 ibgDGridClass.LineWithTriangle = LineWithTriangle;
 ibgDGridClass.Free = Free;
 return (ibGeometry) grdf;
}