ibgdshift.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 19:52:17.75	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDShift >>>  -	Intersection-Based Geometry Definition - Shift */
#define const
#include "ibgd.h"
#include "ibglib.h"
#include "ibgi.h"
#include "ibgd0.h"

typedef struct{
	ibGeometryRec	g;
	ibGeometry	gold;
 	ibgFloat	x[ibgDIM];
}ibgDShiftRec,*ibgDShift;

static int StatusOfEdge(ibGeometry g, ibgPoint *n1, ibgPoint *n2)
{ibgDShift sh = (ibgDShift) g;
 ibgFloat	x1[ibgDIM],x2[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){
	x1[d]		 = ibgpX(*n1)[d];
	x2[d]		 = ibgpX(*n2)[d];
 	ibgpX(*n1)[d]	-= sh->x[d];
 	ibgpX(*n2)[d]	-= sh->x[d];
 }
 rc = ibgiStatusOfEdge((sh->gold),n1,n2);
 for(d=0;d<ibgDIM;d++){ibgpX(*n1)[d]=x1[d];ibgpX(*n2)[d]=x2[d];}
 return rc;
}

static int RegionOfPoint(ibGeometry g, ibgPoint *nnew, const ibgPoint *nold)
{ibgDShift s = (ibgDShift ) g;
 ibgFloat xn[ibgDIM],xo[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){
	xn[d] = ibgpX(*nnew)[d];
	ibgpX(*nnew)[d] -= s->x[d];
 }
 if(nold){
	 for(d=0;d<ibgDIM;d++){
		xo[d] = ibgpX(*nold)[d];
		ibgpX(*nold)[d] -= s->x[d];
	 }
 }
 rc = ibgiRegionOfPoint((s->gold),nnew,nold);
 for(d=0;d<ibgDIM;d++){ibgpX(*nnew)[d] = xn[d];}
 if(nold){
	 for(d=0;d<ibgDIM;d++){ibgpX(*nold)[d] = xo[d];}
 }
 return rc;
}

static int FaceWithEdge(ibGeometry g, ibgPoint *nint,
	       	const ibgPoint *n1, const ibgPoint *n2)
{ibgDShift s = (ibgDShift ) g;
 ibgFloat xo1[ibgDIM],xo2[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){
	xo1[d] = ibgpX(*n1)[d];	ibgpX(*n1)[d] -= s->x[d];
	xo2[d] = ibgpX(*n2)[d];	ibgpX(*n2)[d] -= s->x[d];
 }
 rc = ibgiFaceWithEdge((s->gold),nint,n1,n2);
 for(d=0;d<ibgDIM;d++){
	ibgpX(*n1)[d] = xo1[d];
	ibgpX(*n2)[d] = xo2[d];
	ibgpX(*nint)[d] += s->x[d];}
 return rc;
}

static int LineWithTriangle(ibGeometry g, ibgPoint *nint, ibgPoint *nface,
     const ibgPoint *n1, const ibgPoint *n2, const ibgPoint *n3)
{ibgDShift s = (ibgDShift ) g;
 ibgFloat xob[ibgDIM],xo1[ibgDIM],xo2[ibgDIM],xo3[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){
	xob[d] = ibgpX(*nface)[d]; ibgpX(*nface)[d] -= s->x[d];
        xo1[d] = ibgpX(*n1)[d]; ibgpX(*n1)[d] -= s->x[d];
        xo2[d] = ibgpX(*n2)[d]; ibgpX(*n2)[d] -= s->x[d];
        xo3[d] = ibgpX(*n3)[d]; ibgpX(*n3)[d] -= s->x[d];
 }
 rc = ibgiLineWithTriangle((s->gold),nint,nface,n1,n2,n3);
 for(d=0;d<ibgDIM;d++){
	ibgpX(*nface)[d]= xob[d];
	ibgpX(*n1)[d]	= xo1[d];
	ibgpX(*n2)[d]	= xo2[d];
       	ibgpX(*n3)[d]	= xo3[d];
	ibgpX(*nint)[d]	+= s->x[d];
 }
 return rc;
}

static int NodeInTetrahedron(ibGeometry g, ibgPoint *nint, ibgPoint *nline,
       		     ibgPoint *n1, ibgPoint *n2, ibgPoint *n3, ibgPoint *n4)
{ibgDShift s = (ibgDShift ) g;
 ibgFloat xob[ibgDIM];
 ibgFloat xo1[ibgDIM],xo2[ibgDIM],xo3[ibgDIM],xo4[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){
	xob[d] = ibgpX(*nline)[d]; ibgpX(*nline)[d] -= s->x[d];
        xo1[d] = ibgpX(*n1)[d]; ibgpX(*n1)[d] -= s->x[d];
        xo2[d] = ibgpX(*n2)[d]; ibgpX(*n2)[d] -= s->x[d];
        xo3[d] = ibgpX(*n3)[d]; ibgpX(*n3)[d] -= s->x[d];
        xo4[d] = ibgpX(*n4)[d]; ibgpX(*n4)[d] -= s->x[d];
 }
 rc = ibgiNodeInTetrahedron((s->gold),nint,nline,n1,n2,n3,n4);
 for(d=0;d<ibgDIM;d++){
	ibgpX(*nline)[d]= xob[d];
	ibgpX(*n1)[d]	= xo1[d];
	ibgpX(*n2)[d]	= xo2[d];
	ibgpX(*n3)[d]	= xo3[d];
	ibgpX(*n4)[d]	= xo4[d];
	ibgpX(*nint)[d]	+= s->x[d];
 }
 return rc;
}

static void shiftFree(ibgPtObject sh)
{ibgdFree(((ibgDShift )sh)->gold);
 free(sh);
}

static ibGeometryClassRec	ibgDShiftClass;

ibGeometry ibgdShift(ibGeometry gold, ibgFloat x[ibgDIM])
{ibgDShift sh=(ibgDShift)malloc(sizeof(ibgDShiftRec));
 int d;
 for(d=0;d<ibgDIM;d++) sh->x[d]=x[d];
 ibgdInitialize((ibGeometry)sh,&ibgDShiftClass);
 sh->gold = gold;
 ibgDShiftClass.RegionOfPoint = RegionOfPoint;
 ibgDShiftClass.FaceWithEdge = FaceWithEdge;
 ibgDShiftClass.LineWithTriangle = LineWithTriangle;
 ibgDShiftClass.NodeInTetrahedron = NodeInTetrahedron;
 ibgDShiftClass.StatusOfEdge = StatusOfEdge;
 return (ibGeometry) sh;
}