ibgdsplitr.c

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

/* last edit: Ilja Schmelzer -------------- 24-MAR-1995 17:51:53.58	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDSplitRegion >>> - Intersection-Based Geometry Definition
			      	Splitting a Region
*/
#define const
#include "ibglib.h"
#include "ibgd.h"
#include "ibgi.h"
#include "ibglib.h"
#include "ibgd0.h"
#include "ibgt.h"
#include "ibgdefault.h"

typedef struct _ibgSplitRegion{
	ibGeometryRec	g;
	ibGeometry	gold;
	ibgTopology	told;
	ibgSegment		rold,rnew,face;
	ibgFloat	(*func)(ibgPtObject fThis, ibgPoint *x);
	ibgPtObject	fThis;
	ibgFloat	epsilon;
}ibgSplitRegionRec,*ibgSplitRegion;

static int RegionOfPoint(ibGeometry geom, ibgPoint *nnew, const ibgPoint *nold)
{ibgSplitRegion splt= (ibgSplitRegion)geom;
 ibgiRegionOfPoint(splt->gold,nnew,nold);
 if( ibgpSegment(*nnew) == splt->rold)
 	if(splt->func(splt->fThis,nnew) > 0){
		ibgpSegment(*nnew) = splt->rnew;
	}
 return ibgRegionFound;
}

static ibgSegment splitChangeFaceNr(ibgSplitRegion splt,ibgSegment u)
{register int u0,u1;
 if((u0=u/ibgSMAX) == 0)	return u;
 u1 = u%ibgSMAX;
 if(splt->rold == u0){
 	if(splt->rnew > u1)	return ibgSMAX*splt->rnew+u1;
 	if(splt->rnew ==u1)	ibgfatal;
 	else			return ibgSMAX*u1+splt->rnew;
 }
 if(splt->rold == u1){
 	if(splt->rnew > u0)	return ibgSMAX*splt->rnew+u0;
 	if(splt->rnew ==u0)	ibgfatal;
 	else			return ibgSMAX*u0+splt->rnew;
 }
 ibgfatal;
 return u;
}

static ibgSegment splitCreateLineNr(ibgSplitRegion splt,ibgSegment u)
{int u0,u1;
 if((u0=u/ibgSMAX) == 0){
	ibgassert(ibgIncorrect !=ibgtStatusRF(&splt->g.top,splt->rold,u));
	if(ibgCorrect != ibgtStatusRF(&splt->g.top,splt->rnew,u)){
	      	return ibgDefaultLineNr2(splt->rold,splt->rnew);
	}else{
	      	return u;
	}
 }
 u1 = u%ibgSMAX;
 ibgassert((splt->rold == u0)||(splt->rold == u1));
 if(splt->rnew > u0) return ibgSMAX*(ibgSMAX*splt->rnew+u0)+u1;
 if(splt->rnew ==u0) return u;
 if(splt->rnew > u1) return ibgSMAX*(ibgSMAX*u0+splt->rnew)+u1;
 if(splt->rnew ==u1) return u;
 else                return ibgSMAX*(ibgSMAX*u0+u1)+splt->rnew;
}

static ibgSegment splitChangeLineNr(ibgSplitRegion splt,ibgSegment u)
{int u0,u1,u2,u3;
 if((u0=u/ibgSMAX) == 0)	return u;
 u1 = u0/ibgSMAX;
 if(u1 == 0){
	u1 = u %ibgSMAX;
	if(splt->rold == u0){
		if(splt->rnew > u1)	return ibgSMAX*splt->rnew+u1;
		if(splt->rnew ==u1)	ibgfatal;
		else			return ibgSMAX*u1+splt->rnew;
	}
	if(splt->rold == u1){
		if(splt->rnew > u0)	return ibgSMAX*splt->rnew+u0;
		if(splt->rnew ==u0)	ibgfatal;
		else			return ibgSMAX*u0+splt->rnew;
	}
 }
 u2 = u0%ibgSMAX;
 u3 = u %ibgSMAX;
 if(splt->rold == u1){
	if(splt->rnew > u2) return ibgSMAX*(ibgSMAX*splt->rnew+u2)+u3;
	if(splt->rnew ==u2) return (ibgSMAX*u2)+u3;
	if(splt->rnew > u3) return ibgSMAX*(ibgSMAX*u2+splt->rnew)+u3;
	if(splt->rnew ==u3) return (ibgSMAX*u2)+u3;
	else		   return ibgSMAX*(ibgSMAX*u2+u3)+splt->rnew;
 }
 if(splt->rold == u2){
	if(splt->rnew > u1) return ibgSMAX*(ibgSMAX*splt->rnew+u1)+u3;
	if(splt->rnew ==u1) return (ibgSMAX*u1)+u3;
	if(splt->rnew > u3) return ibgSMAX*(ibgSMAX*u1+splt->rnew)+u3;
	if(splt->rnew ==u3) return (ibgSMAX*u1)+u3;
	else		   return ibgSMAX*(ibgSMAX*u1+u3)+splt->rnew;
 }
 if(splt->rold == u3){
 	if(splt->rnew > u1) return ibgSMAX*(ibgSMAX*splt->rnew+u1)+u2;
	if(splt->rnew ==u1) return (ibgSMAX*u1)+u2;
	if(splt->rnew > u2) return ibgSMAX*(ibgSMAX*u1+splt->rnew)+u2;
	if(splt->rnew ==u2) return (ibgSMAX*u1)+u2;
	else		   return ibgSMAX*(ibgSMAX*u1+u2)+splt->rnew;
 }
 ibgfatal;
 return u;
}

static int SplitRegionZero(ibgPoint *nh, ibgPoint *n1, ibgPoint *n2,
			   	ibgSplitRegion splt)
{double ff,f1,f2,delta; ibgPoint c1,c2; int d,ic,dc;
 f1=(*splt->func)(splt->fThis,n1);
 f2=(*splt->func)(splt->fThis,n2);
 if(f1*f2>0)			{*nh= *n1;	return ibgError;}
 if(f1==0)			{*nh= *n1;	return ibgFaceFound;}
 if(f2==0)			{*nh= *n2;	return ibgFaceFound;}
 c1= *n1; c2= *n2; ic=0; dc=1;
 do{
     if(dc){dc=0; delta=0;
	for(d=0;d<ibgDIM;d++){
		register double dn	= (ibgpX(c1)[d] - ibgpX(c2)[d]);
		delta 			+= dn>0 ? dn : -dn;
		ibgpX(*nh)[d]		= (f2*ibgpX(c1)[d]
					-  f1*ibgpX(c2)[d]) / (f2-f1);
	}
     }else{ dc=1; delta=0;
	for(d=0;d<ibgDIM;d++){
		register double dn 	= (ibgpX(c1)[d] - ibgpX(c2)[d]);
		delta 			+= dn>0 ? dn : -dn;
		ibgpX(*nh)[d] 		= 0.5*(ibgpX(c1)[d] + ibgpX(c2)[d]);
	}
     }
     ibgiRegionOfPoint(splt->gold,nh,&c1);
     if(delta < splt->g.Delta)		return ibgFaceFound;
     if(ibgpSegment(*nh) != splt->rold)	return ibgNotFound; /* another region */
     ff=(splt->func)(splt->fThis,nh);
     if (ff > splt->epsilon)			ff += splt->epsilon;
     else if (ff < -splt->epsilon)		ff -= splt->epsilon;
     else				return ibgFaceFound;
     if(ff*f1>0){f1=ff;for(d=0;d<ibgDIM;d++) ibgpX(c1)[d]=ibgpX(*nh)[d];}
     else       {f2=ff;for(d=0;d<ibgDIM;d++) ibgpX(c2)[d]=ibgpX(*nh)[d];}
 }while(1);				return ibgFaceFound;
}

static int FaceWithEdge(ibGeometry This, ibgPoint *nint,
			const ibgPoint *n1, const ibgPoint *n2)
{ibgSplitRegion splt = ((ibgSplitRegion)This);
 ibgPoint nh,nb; int i,rc;
 if(ibgpSegment(*n1)==splt->rold){
	if(ibgpSegment(*n2)==splt->rnew){
		nb = *n2;
		ibgiRegionOfPoint(splt->gold,&nb,n1);
	       	if(splt->rold==ibgpSegment(nb))
			{ibgpSegment(nb)=splt->rnew;
			 rc=SplitRegionZero(&nh,n1,&nb,splt);}
		else	{rc=ibgNotFound; nh= *n2;}
	}else		{rc=ibgNotFound; nh= *n2;}
	while(rc==ibgNotFound){
		rc=ibgiFaceWithEdge(splt->gold,nint,n1,&nh);
		if(rc==ibgNotFound) return ibgNotFound;
		if((splt->func)(splt->fThis,nint)<=0)	return rc;
		rc=SplitRegionZero(&nh,n1,nint,splt);
	}
	*nint = nh;
	ibgpType(*nint) = ibgSFace;
	ibgpSegment(*nint) = splt->face;
	return ibgFaceFound;
 }else if(ibgpSegment(*n1)==splt->rnew){
	nh = *n1; rc = ibgiRegionOfPoint(splt->gold,&nh,n1);
 	if((rc==ibgError) || (splt->rnew != ibgpSegment(nh))){
		nb = nh;
		if(ibgpSegment(*n2) == splt->rold){
			rc = SplitRegionZero(&nh,&nb,n2,splt);
		}else{	nh = *n2; rc = ibgiRegionOfPoint(splt->gold,&nh,n2);
		 	if((rc==ibgError)) ibgpSegment(nh)=splt->rold;
			rc = ibgNotFound;
		}
	}else{
		rc = ibgiFaceWithEdge(splt->gold,nint,n1,n2);
		if(rc==ibgNotFound) return ibgNotFound;
		ibgassert(ibgIncorrect!=ibgtStatusRF(splt->told,splt->rnew,
							      ibgpSegment(*nint)));
		if(ibgIncorrect==ibgtStatusRF(splt->told,splt->rold,ibgpSegment(*nint)))
						return ibgFaceFound;
		if((splt->func)(splt->fThis,nint)<=0)	return ibgFaceFound;
		nb = *nint;
  		if(ibgpSegment(*n2)==splt->rold){
		 	rc = SplitRegionZero(&nh,&nb,n2,splt);
		}else{	rc = ibgNotFound; nh = *n2;
		}
	}
	while(rc==ibgNotFound){
		rc = ibgiFaceWithEdge(splt->gold,nint,&nb,&nh);
		if(rc==ibgNotFound) return ibgNotFound;
		if((splt->func)(splt->fThis,nint)>0){
		    if(ibgCorrect==ibgtStatusRF(splt->told,splt->rold
						,ibgpSegment(*nint))){
		    	ibgpSegment(*nint)=splitChangeFaceNr(splt,ibgpSegment(*nint));
		    }
		    return  ibgFaceFound;
		}
		rc=SplitRegionZero(&nh,&nb,nint,splt);
	}
	*nint = nh;
	ibgpType(*nint) = ibgSFace;
	ibgpSegment(*nint) = splt->face;
	return ibgFaceFound;
 }
 else if(ibgpSegment(*n2)==splt->rnew){
	nh = *n2; rc = ibgiRegionOfPoint(splt->gold,&nh,n2);
 	if(rc == ibgError) ibgpSegment(nh) = splt->rold;
      	rc = ibgiFaceWithEdge(splt->gold,nint,n1,&nh);
 }
 else	rc = ibgiFaceWithEdge(splt->gold,nint,n1,n2);
 if(rc==ibgNotFound) return ibgNotFound;

 if(ibgCorrect==ibgtStatusRF(splt->told,splt->rold,ibgpSegment(*nint))){
      	if((splt->func)(splt->fThis,nint)>0){
  		ibgpSegment(*nint) = splitChangeFaceNr(splt,ibgpSegment(*nint));
	}
 }
 return  ibgFaceFound;
}


#define ibgSplitFound	4567

static int SplitRegionLine0(ibgSplitRegion splt,
		      	ibgPoint *nint, ibgPoint *nface,
			ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{ibgPoint n12,n23,n31,nold;
 ibGeometry g=splt->gold;
 int i,d,s,sn,rc;
 ibgFloat fi,f0,f1,f2,f3,p12,p23,p31,q12,q23,q31,dd,s0;
 f0 = (splt->func)(splt->fThis,nface);
 s0 = f0>0 ? 1 : -1; f0 *= s0;
 rc = ibgiLineWithTriangle(splt->gold,nint,nface,n1,n2,n3);
 fi = s0*(splt->func)(splt->fThis,nint);
 if(fi >= splt->epsilon)	return rc;
 if(rc == ibgLineFound)	return rc;
 if(fi >= -splt->epsilon)return ibgSplitFound;
 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(p12<0.1)	{p12=0.1; q12=0.9;}
 if(q12<0.1)	{q12=0.1; p12=0.9;}
 if(f2*f3>=0){
	p23 = q23 = 0.5;