ibgi.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 18:59:17.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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGI >>> - Intersection-Based Geometry Interface */

#include "ibgi.h"
#include "ibglib.h"
#include "ibgd0.h"

ibgPoint ibgNoPoint;

int	ibgiRegionOfPoint(ibGeometry g, ibgPoint *nnew, ibgPoint *nold)
{register int rc = (g)->Class->RegionOfPoint(g,nnew,nold);
 if(rc != ibgNotFound){
	ibgassert(ibgpType(*nnew) == ibgSRegion);
 }else{
	ibgassert(ibgpType(*nnew) == ibgSNothing);
 }
 return rc;
}

int	ibgiFaceWithEdge  (ibGeometry g, ibgPoint *nint, ibgPoint *n1, ibgPoint *n2)
{register int rc = (g)->Class->FaceWithEdge(g,nint,n1,n2);
 if(rc != ibgNotFound){
	ibgassert(ibgpType(*nint) == ibgSFace);
 }else{
	ibgassert(ibgpType(*nint) == ibgSNothing);
 }
 return rc;
}

int	ibgiLineWithTriangle  (ibGeometry g, ibgPoint *nint, ibgPoint *nface,
	     	     	ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{register int rc = (g)->Class->LineWithTriangle(g,nint,nface,n1,n2,n3);
 if(rc == ibgLineFound){
	ibgassert(ibgpType(*nint) == ibgSLine);
 }else if(rc != ibgNotFound){
	ibgassert(ibgpType(*nint) == ibgSFace);
 }else{
	ibgassert(ibgpType(*nint) == ibgSNothing);
 }
 return rc;
}

int	ibgiNodeInTetrahedron(ibGeometry g, ibgPoint *nint, ibgPoint *nline,
			ibgPoint *n1, ibgPoint *n2, ibgPoint *n3, ibgPoint *n4)
{register int rc = (g)->Class->NodeInTetrahedron(g,nint,nline,n1,n2,n3,n4);
 if(rc == ibgNodeFound){
	ibgassert(ibgpType(*nint) == ibgSNode);
 }else if(rc != ibgNotFound){
	ibgassert(ibgpType(*nint) == ibgSLine);
 }else{
	ibgassert(ibgpType(*nint) == ibgSNothing);
 }
 return rc;
}

int	ibgiLineWithRectangle(ibGeometry g, ibgPoint *nint, ibgPoint *nface,
	     	     	ibgPoint *n1, ibgPoint *n2, ibgPoint *n3, ibgPoint *n4)
{int i,s,rc;
 ibgPoint nold;
 s = -1;
 for(i=0;i<g->LoopBreakRectangle;i++){
	switch(s){
	case -1: rc = ibgiLineWithTriangle(g,nint,nface,n1,n2,n3);
		switch(rc){
			case ibgiOnSide12: s = ibgiOnSide12; break;
			case ibgiOnSide23: s = ibgiOnSide23; break;
			case ibgiOnSide31: s = -2; break;
		}break;
	case -2: rc = ibgiLineWithTriangle(g,nint,&nold,n1,n3,n4);
		switch(rc){
			case ibgiOnSide12: s = -3; break;
			case ibgiOnSide23: s = ibgiOnSide34; break;
			case ibgiOnSide31: s = ibgiOnSide41; break;
		}break;
	case -3: rc = ibgiLineWithTriangle(g,nint,&nold,n3,n1,n2);
		switch(rc){
			case ibgiOnSide12: s = -2; break;
			case ibgiOnSide23: s = ibgiOnSide12; break;
			case ibgiOnSide31: s = ibgiOnSide23; break;
		}break;
	}
	if(rc==ibgLineFound) return rc;
	if(s > 0)	 return s;
	nold = *nint;
 }
 ibgpType(*nint) = ibgSNothing;
 ibgpSegment(*nint) = 0;
 return ibgNotFound;
}

int	ibgiNodeInHexahedron(ibGeometry g, ibgPoint *nint, ibgPoint *nline,
			ibgPoint *n1, ibgPoint *n2, ibgPoint *n3, ibgPoint *n4,
	    		ibgPoint *n5, ibgPoint *n6, ibgPoint *n7, ibgPoint *n8)
{int i,s,rc;
 ibgPoint nold;
 register dd;
 register ibgFloat *c0 = ibgpX(*nline);
 register ibgFloat *c1 = ibgpX(*n1), *c2 = ibgpX(*n3), *c3 = ibgpX(*n7);
 register ibgFloat x1=c1[0]-c0[0],x2=c2[0]-c0[0],x3=c3[0]-c0[0];
 register ibgFloat y1=c1[1]-c0[1],y2=c2[1]-c0[1],y3=c3[1]-c0[1];
 register ibgFloat z1=c1[2]-c0[2],z2=c2[2]-c0[2],z3=c3[2]-c0[2];
 dd = x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1);
 if(dd>0) s = -1;
 else	  s = -2;
 for(i=0;i<g->LoopBreakHexahedron;i++){
	switch(s){
	case -1: rc = ibgiNodeInTetrahedron(g,nint,nline,n1,n2,n3,n7);
		switch(rc){
			case ibgiOnSide123: s = ibgiOnSide1234; break;
			case ibgiOnSide124: s = -26; break;
			case ibgiOnSide134: s = -34; break;
			case ibgiOnSide234: s = ibgiOnSide2367; break;
		}break;
	case -2: rc = ibgiNodeInTetrahedron(g,nint,nline,n1,n3,n4,n7);
		switch(rc){
			case ibgiOnSide123: s = ibgiOnSide1234; break;
			case ibgiOnSide124: s = -32; break;
			case ibgiOnSide134: s = -48; break;
			case ibgiOnSide234: s = ibgiOnSide3478; break;
		}break;
	case -23: rc = ibgiNodeInTetrahedron(g,nint,&nold,n7,n2,n1,n3);
		switch(rc){
			case ibgiOnSide123: s = -26; break;
			case ibgiOnSide124: s = ibgiOnSide2367; break;
			case ibgiOnSide134: s = -34; break;
			case ibgiOnSide234: s = ibgiOnSide1234; break;
		}break;
	case -32: rc = ibgiNodeInTetrahedron(g,nint,&nold,n1,n3,n7,n2);
		switch(rc){
			case ibgiOnSide123: s = -34; break;
			case ibgiOnSide124: s = ibgiOnSide1234; break;
			case ibgiOnSide134: s = -26; break;
			case ibgiOnSide234: s = ibgiOnSide2367; break;
		}break;
	case -34: rc = ibgiNodeInTetrahedron(g,nint,&nold,n7,n3,n1,n4);
		switch(rc){
			case ibgiOnSide123: s = -32; break;
			case ibgiOnSide124: s = ibgiOnSide3478; break;
			case ibgiOnSide134: s = -48; break;
			case ibgiOnSide234: s = ibgiOnSide1234; break;
		}break;
	case -43: rc = ibgiNodeInTetrahedron(g,nint,&nold,n1,n4,n7,n3);
		switch(rc){
			case ibgiOnSide123: s = -48; break;
			case ibgiOnSide124: s = ibgiOnSide1234; break;
			case ibgiOnSide134: s = -32; break;
			case ibgiOnSide234: s = ibgiOnSide3478; break;
		}break;
	case -48: rc = ibgiNodeInTetrahedron(g,nint,&nold,n7,n4,n1,n8);
		switch(rc){
			case ibgiOnSide123: s = -43; break;
			case ibgiOnSide124: s = ibgiOnSide3478; break;
			case ibgiOnSide134: s = -85; break;
			case ibgiOnSide234: s = ibgiOnSide1458; break;
		}break;
	case -84: rc = ibgiNodeInTetrahedron(g,nint,&nold,n1,n8,n7,n4);
		switch(rc){
			case ibgiOnSide123: s = -85; break;
			case ibgiOnSide124: s = ibgiOnSide1458; break;
			case ibgiOnSide134: s = -43; break;
			case ibgiOnSide234: s = ibgiOnSide3478; break;
		}break;
	case -85: rc = ibgiNodeInTetrahedron(g,nint,&nold,n7,n8,n1,n5);
		switch(rc){
			case ibgiOnSide123: s = -84; break;
			case ibgiOnSide124: s = ibgiOnSide5678; break;
			case ibgiOnSide134: s = -56; break;
			case ibgiOnSide234: s = ibgiOnSide1458; break;
		}break;
	case -58: rc = ibgiNodeInTetrahedron(g,nint,&nold,n1,n5,n7,n8);
		switch(rc){
			case ibgiOnSide123: s = -56; break;
			case ibgiOnSide124: s = ibgiOnSide1458; break;
			case ibgiOnSide134: s = -84; break;
			case ibgiOnSide234: s = ibgiOnSide5678; break;
		}break;
	case -56: rc = ibgiNodeInTetrahedron(g,nint,&nold,n7,n5,n1,n6);
		switch(rc){
			case ibgiOnSide123: s = -58; break;
			case ibgiOnSide124: s = ibgiOnSide5678; break;
			case ibgiOnSide134: s = -62; break;
			case ibgiOnSide234: s = ibgiOnSide1256; break;
		}break;
	case -65: rc = ibgiNodeInTetrahedron(g,nint,&nold,n1,n6,n7,n5);
		switch(rc){
			case ibgiOnSide123: s = -62; break;
			case ibgiOnSide124: s = ibgiOnSide1256; break;
			case ibgiOnSide134: s = -58; break;
			case ibgiOnSide234: s = ibgiOnSide5678; break;
		}break;
	case -62: rc = ibgiNodeInTetrahedron(g,nint,&nold,n7,n6,n1,n2);
		switch(rc){
			case ibgiOnSide123: s = -65; break;
			case ibgiOnSide124: s = ibgiOnSide2367; break;
			case ibgiOnSide134: s = -23; break;
			case ibgiOnSide234: s = ibgiOnSide1256; break;
		}break;
	case -26: rc = ibgiNodeInTetrahedron(g,nint,&nold,n1,n2,n7,n6);
		switch(rc){
			case ibgiOnSide123: s = -23; break;
			case ibgiOnSide124: s = ibgiOnSide1256; break;
			case ibgiOnSide134: s = -65; break;
			case ibgiOnSide234: s = ibgiOnSide2367; break;
		}break;
	}
	if(rc==ibgNodeFound) return rc;
	if(s > 0)	 return s;
 }
 ibgpType(*nint) = ibgSNothing;
 ibgpSegment(*nint) = 0;
 return ibgNotFound;
}