ibgdsplitf.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 19:18:21.72	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDSplitFace >>> - Intersection-Based Geometry Definition
  				Splitting a Face
*/   
#define const
#include "ibglib.h"
#include "ibgd.h"
#include "ibgi.h"
#include "ibgd0.h"
#include "ibgdefault.h"

typedef struct _ibgSplitFace{
	ibGeometryRec	g;
	ibGeometry  	gold;
	ibgSegment		fold,fnew,line;
	ibgFloat	(*func)(ibgPtObject fThis, ibgPoint *x);
       	ibgPtObject	fThis;
	ibgFloat	epsilon;
}ibgSplitFaceRec,*ibgSplitFace;

static int RegionOfPoint(ibGeometry geom, ibgPoint *nnew, ibgPoint *nold)
{ibgSplitFace splt= (ibgSplitFace)geom;
 return ibgiRegionOfPoint(splt->gold,nnew,nold);
}

static int FaceWithEdge(ibGeometry g, ibgPoint *nint,
					ibgPoint *n1, ibgPoint *n2)
{ibgSplitFace splt = (ibgSplitFace ) g;
 int rc = ibgiFaceWithEdge(splt->gold,nint,n1,n2);
 if(ibgpSegment(*nint) == splt->fold){
 	if(splt->func(splt->fThis,nint)>0){
		ibgpSegment(*nint) = splt->fnew;
	}
 }
 return rc;
}

static int LineWithTriangle(ibGeometry g, ibgPoint *nint, ibgPoint *nface,
		    	ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{ibgSplitFace splt = (ibgSplitFace ) g;
 ibgPoint n12,n23,n31,nold,nfold;
 ibGeometry go=splt->gold;
 ibgSegment fu;
 int mf,m1,m2,m3,m4,b1,i,d,s,s0,sn,rc;
 ibgFloat fi,f0,f1,f2,f3,p12,p23,p31,q12,q23,q31,dd;
 f0 = (splt->func)(splt->fThis,nface);
 nfold = *nface;
 if((fu=ibgpSegment(*nface))==splt->fnew){
/* potential error if an old part of fnew exists!!! */
	if(f0>0) ibgpSegment(nfold) = splt->fold;
 }else if(fu==splt->fold){;
 }else{
	return ibgiLineWithTriangle(go,nint,&nfold,n1,n2,n3);
 }
 s0 = f0>0 ? 1 : -1; f0 *= s0;
 rc = ibgiLineWithTriangle(go,nint,&nfold,n1,n2,n3);
 fi = (splt->func)(splt->fThis,nint);
 if(rc != ibgLineFound){
	if(ibgpSegment(*nint)==splt->fold){
		if(fi>0){
			ibgpSegment(*nint) = splt->fnew;
		}
	}
 }
 fi *= s0;
 if(fi >= splt->epsilon)	return rc;
 if(rc == ibgLineFound)	return rc;
 if(fi >= -splt->epsilon){
	ibgpSegment(*nint) = splt->line;
	ibgpType(*nint) = ibgSLine;
	return ibgLineFound;
 }
 f1 = (splt->func)(splt->fThis,n1);
 f2 = (splt->func)(splt->fThis,n2);
 f3 = (splt->func)(splt->fThis,n3);
 if(f1*f2>=0){
	p12 = q12 = 0.5;
 }else{	p12 = f2/(f2-f1);
	q12 = f1/(f1-f2);
 }
 if(f2*f3>=0){
	p23 = q23 = 0.5;
 }else{	p23 = f3/(f3-f2);
	q23 = f2/(f2-f3);
 }
 if(f3*f1>=0){
	p31 = q31 = 0.5;
 }else{	p31 = f1/(f1-f3);
	q31 = f3/(f3-f1);
 }
 for(d=0;d<ibgDIM;d++){
	ibgpX(n12)[d]=(p12*ibgpX(*n1)[d]+q12*ibgpX(*n2)[d]);
	ibgpX(n23)[d]=(p23*ibgpX(*n2)[d]+q23*ibgpX(*n3)[d]);
	ibgpX(n31)[d]=(p31*ibgpX(*n3)[d]+q31*ibgpX(*n1)[d]);
 }
 ibgiRegionOfPoint(go,&n12,n1);
 ibgiRegionOfPoint(go,&n23,n2);
 ibgiRegionOfPoint(go,&n31,n3);
 dd=0;
 for(d=0;d<ibgDIM;d++) {
	dd += (ibgpX(*nface)[d] - ibgpX(n12)[d])
	    * (ibgpX(   n12)[d] - ibgpX(*n1)[d]);
 }
 if(dd<0) s= -1; else s= -2;
 for(i=0;i<ibgdLSteps;i++){
	switch(s){
	case -1:
		sn=LineWithTriangle(g,nint,nface,n1,&n12,&n31);
		switch(sn){
		case ibgiOnSide12:	s = ibgiOnSide12; break;
		case ibgiOnSide23:	s = -10; break;
		case ibgiOnSide31:	s = ibgiOnSide31; break;
		}break;
	case -2:
		sn=LineWithTriangle(g,nint,nface,&n12,n2,&n23);
		switch(sn){
		case ibgiOnSide12:	s = ibgiOnSide12; break;
		case ibgiOnSide23:	s = ibgiOnSide23; break;
		case ibgiOnSide31:	s = -20; break;
		}break;
	case -23:
		sn=LineWithTriangle(g,nint,&nold,&n12,&n31,n1);
		switch(sn){
		case ibgiOnSide12:	s = -10; break;
		case ibgiOnSide23:	s = ibgiOnSide31; break;
		case ibgiOnSide31:	s = ibgiOnSide12; break;
		}break;
	case -31:
		sn=LineWithTriangle(g,nint,&nold,&n23,&n12,n2);
		switch(sn){
		case ibgiOnSide12:	s = -20; break;
		case ibgiOnSide23:	s = ibgiOnSide12; break;
		case ibgiOnSide31:	s = ibgiOnSide23; break;
		}break;
	case -12:
		sn=LineWithTriangle(g,nint,&nold,&n31,&n23,n3);
		switch(sn){
		case ibgiOnSide12:	s = -30; break;
		case ibgiOnSide23:	s = ibgiOnSide23; break;
		case ibgiOnSide31:	s = ibgiOnSide31; break;
		}break;
	case -10:
		sn=LineWithTriangle(g,nint,&nold,&n31,&n12,&n23);
		switch(sn){
		case ibgiOnSide12:	s = -23; break;
		case ibgiOnSide23:	s = -31; break;
		case ibgiOnSide31:	s = -12; break;
		}break;
	case -20:
		sn=LineWithTriangle(g,nint,&nold,&n12,&n23,&n31);
		switch(sn){
		case ibgiOnSide12:	s = -31; break;
		case ibgiOnSide23:	s = -12; break;
		case ibgiOnSide31:	s = -23; break;
		}break;
	case -30:
		sn=LineWithTriangle(g,nint,&nold,&n23,&n31,&n12);
		switch(sn){
		case ibgiOnSide12:	s = -12; break;
		case ibgiOnSide23:	s = -23; break;
		case ibgiOnSide31:	s = -31; 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 g)
{ibgSplitFace splt = (ibgSplitFace) g;
 ibgdFree(splt->gold);
 return ibgSuccess;
}

static ibGeometryClassRec	ibgSplitFaceClass;

ibGeometry ibgdSplitFace(ibGeometry 	gold,
			ibgSegment		fnew,
			ibgSegment		fold,
			ibgFloat	(*func)(ibgPtObject fThis, ibgPoint *n),
			ibgPtObject	fThis,
	    		ibgFloat	Dfunc)
{ibgSplitFace splt=(ibgSplitFace) malloc(sizeof(ibgSplitFaceRec));
 ibgdInitialize((ibGeometry)splt,&ibgSplitFaceClass);
 ibgtCopy(&(splt->g.top),&(gold->top));
 ibgtNewFace(&(splt->g.top),fnew,fold);
 splt->gold=gold;
 splt->fold=fold; splt->fnew=fnew;
 if(splt->fold >= ibgSMAX){
	splt->line=ibgDefaultLineNr2(splt->fold/ibgSMAX,splt->fold%ibgSMAX);
 }else{
	splt->line=ibgDefaultLineNr;
 }
 splt->func=func; splt->fThis=fThis;
 splt->epsilon = 0.5*splt->g.Delta*Dfunc;
 ibgSplitFaceClass.RegionOfPoint = RegionOfPoint;
 ibgSplitFaceClass.FaceWithEdge = FaceWithEdge;
 ibgSplitFaceClass.LineWithTriangle = LineWithTriangle;
 ibgSplitFaceClass.Free = Free;
 return (ibGeometry)splt;
}