ibgdwindow.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 20:00:12.92	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDWindow >>> - Intersection-Based Geometry Definition
		1D-Continuation on the outside of a window
*/
#define const
#include "ibglib.h"
#include "ibgd.h"
#include "ibgi.h"
#include "ibgd0.h"
#include "ibgdefault.h"

typedef struct{
	ibGeometryRec	g;
   	ibGeometry	gold;
	ibgFloat	xmin[ibgDIM],xmax[ibgDIM];
}ibgWindowRec,*ibgWindow;

static int RegionOfPoint(ibGeometry geo, ibgPoint *nnew, ibgPoint *nold)
{ibgWindow wind= (ibgWindow )geo;
 ibgFloat xo[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){
	register x = xo[d] = ibgpX(*nnew)[d];
	if(x < wind->xmin[d]){
		ibgpX(*nnew)[d] = wind->xmin[d];
	}else if(x > wind->xmax[d]){
		ibgpX(*nnew)[d] = wind->xmax[d];
	}
 }
 rc = ibgiRegionOfPoint(wind->gold,nnew,nold);
 for(d=0;d<ibgDIM;d++){
	ibgpX(*nnew)[d] = xo[d];
 }
 return rc;
}

static int FaceWithEdge(ibGeometry geo, ibgPoint *nint,
					ibgPoint *n1, ibgPoint *n2)
{ibgWindow wind= (ibgWindow )geo;
 int rc,d,d1,d2;
 ibgFloat p1,p2,x1,x2,y1,y2,f1,f2,ff,xo1,xo2,v;
 ibgPoint nh;
 for(d=0;d<ibgDIM;d++){
	if(ibgpX(*n1)[d]<(v=wind->xmin[d])){
		if(ibgpX(*n2)[d]>v){
			goto fcut;
		}else{
			goto fproj;
		}
	}else if(ibgpX(*n2)[d]<v){
		if(ibgpX(*n1)[d]>v){
			goto fcut;
		}else{
			goto fproj;
		}
	}else if(ibgpX(*n1)[d]>(v=wind->xmax[d])){
		if(ibgpX(*n2)[d]<v){
			goto fcut;
		}else{
			goto fproj;
		}
	}else if(ibgpX(*n2)[d]>v){
		if(ibgpX(*n1)[d]<v){
			goto fcut;
		}else{
			goto fproj;
		}
	}
 }
 return ibgiFaceWithEdge(wind->gold,nint,n1,n2);
fproj:
 xo1 = ibgpX(*n1)[d];
 xo2 = ibgpX(*n2)[d];
 ibgpX(*n1)[d] = v;
 ibgpX(*n2)[d] = v;
 rc = ibgiFaceWithEdge(wind->gold,nint,n1,n2);
 ibgpX(*n1)[d] = xo1;
 ibgpX(*n2)[d] = xo2;
 d1 = (d+1)%ibgDIM;
 d2 = (d+2)%ibgDIM;
 x1 = ibgpX(*n1)[d1]-ibgpX(*nint)[d1];
 x2 = ibgpX(*n2)[d1]-ibgpX(*nint)[d1];
 y1 = ibgpX(*n1)[d2]-ibgpX(*nint)[d2];
 y2 = ibgpX(*n2)[d2]-ibgpX(*nint)[d2];
 f1 = x2*(x1-x2) + y2*(y1-y2);
 f2 = x1*(x1-x2) + y1*(y1-y2);
 ff = f1 - f2;
 if(ff==0){
	ff = f1 = 1; f2 = 0;
 }
 f1 /= ff; f2 /= ff;
 for(d=0;d<ibgDIM;d++)
 ibgpX(*nint)[d] = f1*xo1 + f2*xo2;
 return rc;
fcut:
 p2 = ibgpX(*n1)[d] - v;
 p1 = v - ibgpX(*n2)[d];
 p2 /= p2+p1; p1 /= p2+p1;
 for(d1=0;d1<ibgDIM;d1++){
	ibgpX(nh)[d1] = p1*ibgpX(*n1)[d1] + p2*ibgpX(*n2)[d1];
 }
 ibgiRegionOfPoint(wind->gold,&nh,n1);
 rc = FaceWithEdge(geo,nint,n1,&nh);
 if(rc==ibgNotFound){
	 rc = FaceWithEdge(geo,nint,&nh,n2);
 }
 return rc;
}

static int LineWithTriangle(ibGeometry geo, ibgPoint *nint, ibgPoint *nface,
		   	     	ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{ibgWindow wind = (ibgWindow ) geo;
 ibgPoint nh,nh1,nh2,nh3,nold;
 ibgFloat x1,x2,x3,y1,y2,y3,f1,f2,f3,ff,xo1,xo2,xo3,xof,v,p1,p2;
 int	s12,s23,s31,s,sn,d,d1,d2,rc,i;
 for(d=0;d<ibgDIM;d++){
	if(ibgpX(*n1)[d]<(v=wind->xmin[d])){
		if(ibgpX(*n2)[d]>v){
			nh1 = *n1; nh2 = *n2; nh3 = *n3;
			if(ibgpX(*nface)[d]<v) s = -1; else s = -2;
			s12=ibgiOnSide12;s23=ibgiOnSide23;s31=ibgiOnSide31;
			goto lcut;
		}else if(ibgpX(*n3)[d]>v){
			nh1 = *n3; nh2 = *n1; nh3 = *n2; s = -3;
			s12=ibgiOnSide31;s23=ibgiOnSide12;s31=ibgiOnSide23;
			goto lcut;
		}else{
			goto lproj;
		}
	}else if(ibgpX(*n2)[d]<v){
		if(ibgpX(*n1)[d]>v){
			nh1 = *n1; nh2 = *n2; nh3 = *n3;
			if(ibgpX(*nface)[d]>v) s = -1; else s = -2;
			s12=ibgiOnSide12;s23=ibgiOnSide23;s31=ibgiOnSide31;
			goto lcut;
		}else if(ibgpX(*n3)[d]>v){
			nh1 = *n2; nh2 = *n3; nh3 = *n1; s = -4;
			s12=ibgiOnSide23;s23=ibgiOnSide31;s31=ibgiOnSide12;
			goto lcut;
		}else{
			goto lproj;
		}
	}else if(ibgpX(*n3)[d]<v){
		if(ibgpX(*n1)[d]>v){
			nh1 = *n3; nh2 = *n1; nh3 = *n2; s = -3;
			s12=ibgiOnSide31;s23=ibgiOnSide12;s31=ibgiOnSide23;
			goto lcut;
		}else if(ibgpX(*n2)[d]>v){
			nh1 = *n2; nh2 = *n3; nh3 = *n1; s = -4;
			s12=ibgiOnSide23;s23=ibgiOnSide31;s31=ibgiOnSide12;
			goto lcut;
		}else{
			goto lproj;
		}
	}else if(ibgpX(*n1)[d]>(v=wind->xmax[d])){
		if(ibgpX(*n2)[d]<v){
			nh1 = *n1; nh2 = *n2; nh3 = *n3;
			if(ibgpX(*nface)[d]>v) s = -1; else s = -2;
			s12=ibgiOnSide12;s23=ibgiOnSide23;s31=ibgiOnSide31;
			goto lcut;
		}else if(ibgpX(*n3)[d]<v){
			nh1 = *n3; nh2 = *n1; nh3 = *n2; s = -3;
			s12=ibgiOnSide31;s23=ibgiOnSide12;s31=ibgiOnSide23;
			goto lcut;
		}else{
			goto lproj;
		}
	}else if(ibgpX(*n2)[d]>v){
		if(ibgpX(*n1)[d]<v){
			nh1 = *n1; nh2 = *n2; nh3 = *n3;
			if(ibgpX(*nface)[d]<v) s = -1; else s = -2;
			s12=ibgiOnSide12;s23=ibgiOnSide23;s31=ibgiOnSide31;
			goto lcut;
		}else if(ibgpX(*n3)[d]<v){
			nh1 = *n2; nh2 = *n3; nh3 = *n1; s = -4;
			s12=ibgiOnSide23;s23=ibgiOnSide31;s31=ibgiOnSide12;
			goto lcut;
		}else{
			goto lproj;
		}
	}else if(ibgpX(*n3)[d]>v){
		if(ibgpX(*n1)[d]<v){
			nh1 = *n3; nh2 = *n1; nh3 = *n2; s = -3;
			s12=ibgiOnSide31;s23=ibgiOnSide12;s31=ibgiOnSide23;
			goto lcut;
		}else if(ibgpX(*n2)[d]<v){
			nh1 = *n2; nh2 = *n3; nh3 = *n1; s = -4;
			s12=ibgiOnSide23;s23=ibgiOnSide31;s31=ibgiOnSide12;
			goto lcut;
		}else{
			goto lproj;
		}
	}
 }
 return ibgiLineWithTriangle(wind->gold,nint,nface,n1,n2,n3);
lproj:
 xo1 = ibgpX(*n1)[d];
 xo2 = ibgpX(*n2)[d];
 xo3 = ibgpX(*n3)[d];
 xof = ibgpX(*nface)[d];
 ibgpX(*n1)[d] = v;
 ibgpX(*n2)[d] = v;
 ibgpX(*n3)[d] = v;
 ibgpX(*nface)[d] = v;
 rc = ibgiLineWithTriangle(wind->gold,nint,nface,n1,n2,n3);
 d1 = (d+1)%ibgDIM;
 d2 = (d+2)%ibgDIM;
 x1 = ibgpX(*n1)[d1]-ibgpX(*nint)[d1];
 x2 = ibgpX(*n2)[d1]-ibgpX(*nint)[d1];
 x3 = ibgpX(*n3)[d1]-ibgpX(*nint)[d1];
 y1 = ibgpX(*n1)[d2]-ibgpX(*nint)[d2];
 y2 = ibgpX(*n2)[d2]-ibgpX(*nint)[d2];
 y3 = ibgpX(*n3)[d2]-ibgpX(*nint)[d2];
 f1 = x2*y3 - x3*y2;
 f2 = x3*y1 - x1*y3;
 f3 = x1*y2 - x2*y1;
 if((ff=f1+f2+f3)==0){
 	f1 = x2*x3 + y2*y3;
 	f2 = x3*x1 + y3*y1;
 	f3 = x1*x2 + y1*y2;
 	if(f1<f2){
		if(f2<f3){
			f2=0; ff = (f3-f1); f1 *= -1;
		}else if(f1<f3){
		 	f3=0; ff = (f1-f2); f2 *= -1;
		}else{
		 	f1=0; ff = (f2-f3); f3 *= -1;
		}
	}else{
		if(f2>f3){
			f2=0; ff = (f3-f1); f1 *= -1;
		}else if(f1>f3){
		 	f3=0; ff = (f1-f2); f2 *= -1;
		}else{
		 	f1=0; ff = (f2-f3); f3 *= -1;
		}
 	}
	if(ff==0){
		f1=ff=1; f2=f3=0;
	}
 }
 f1 /= ff;	f2 /= ff;	f3 /= ff;
 ibgpX(*n1)[d] = xo1;
 ibgpX(*n2)[d] = xo2;
 ibgpX(*n3)[d] = xo3;
 ibgpX(*nface)[d] = xof;
 ibgpX(*nint)[d]  = f1*xo1 + f2*xo2 + f3*xo3;
 return rc;
lcut:
 p2 = ibgpX(nh1)[d] - v;
 p1 = v - ibgpX(nh2)[d];
 p2 /= p2+p1; p1 /= p2+p1;
 for(d1=0;d1<ibgDIM;d1++){
	ibgpX(nh)[d1] = p1*ibgpX(nh1)[d1] + p2*ibgpX(nh2)[d1];
 }
 for(i=0;i<ibgdLSteps;i++){
	switch(s){
	case -1:
		sn=LineWithTriangle(geo,nint,nface,&nh1,&nh,&nh3);
		switch(sn){
		case ibgiOnSide12:	s = s12; break;
		case ibgiOnSide23:	s = -6; break;
		case ibgiOnSide31:	s = s31; break;
		}break;
	case -2:
		sn=LineWithTriangle(geo,nint,nface,&nh,&nh2,&nh3);
		switch(sn){
		case ibgiOnSide12:	s = s12; break;
		case ibgiOnSide23:	s = s23; break;
		case ibgiOnSide31:	s = -5; break;
		}break;
	case -3:
		sn=LineWithTriangle(geo,nint,nface,&nh2,&nh3,&nh);
		switch(sn){
		case ibgiOnSide12:	s = s23; break;
		case ibgiOnSide23:	s = -5; break;
		case ibgiOnSide31:	s = s12; break;
		}break;
	case -4:
		sn=LineWithTriangle(geo,nint,nface,&nh3,&nh1,&nh);
		switch(sn){
		case ibgiOnSide12:	s = s31; break;
		case ibgiOnSide23:	s = -6; break;
		case ibgiOnSide31:	s = s12; break;
		}break;
	case -5:
		sn=LineWithTriangle(geo,nint,&nold,&nh,&nh3,&nh1);
		switch(sn){
		case ibgiOnSide12:	s = -6; break;
		case ibgiOnSide23:	s = s31; break;
		case ibgiOnSide31:	s = s12; break;
		}break;
	case -6:
		sn=LineWithTriangle(geo,nint,&nold,&nh3,&nh,&nh2);
		switch(sn){
		case ibgiOnSide12:	s = -5; break;
		case ibgiOnSide23:	s = s12; break;
		case ibgiOnSide31:	s = s23; break;
		}break;
	}
	if(sn== ibgLineFound){
		return ibgLineFound;
  	}else if(s>0){
	      	return s;
	}
	nold = *nint;
 }
 ibgpSegment(*nint) = ibgDefaultLineNr;
 ibgpType(*nint) = ibgSNothing;
 return ibgNotFound;
}

static int Free(ibGeometry geo)
{ibgWindow wind = (ibgWindow) geo;
 ibgdFree(wind->gold);
 return ibgSuccess;
}

static ibGeometryClassRec	ibgWindowClass;

ibGeometry ibgdWindow(ibGeometry gold,
		ibgFloat	xmin[ibgDIM],
		ibgFloat	xmax[ibgDIM])
{ibgWindow wind=(ibgWindow)malloc(sizeof(ibgWindowRec));
 int d;
 ibgdInitialize((ibGeometry)wind,&ibgWindowClass);
 wind->gold=gold;
 for(d=0;d<ibgDIM;d++){
	wind->xmin[d] = xmin[d] + wind->g.Delta;
	wind->xmax[d] = xmax[d] - wind->g.Delta;
 }
 ibgWindowClass.RegionOfPoint = RegionOfPoint;
 ibgWindowClass.FaceWithEdge = FaceWithEdge;
 ibgWindowClass.LineWithTriangle = LineWithTriangle;
 ibgWindowClass.Free = Free;
 return (ibGeometry) wind;
}