ibgdlinear.c

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

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

typedef struct{
	ibGeometryRec	g;
	ibGeometry	old;
 	ibgFloat	a[ibgDIM][ibgDIM],r[ibgDIM][ibgDIM];
}ibgDLinearRec,*ibgDLinear;

static int StatusOfEdge(ibGeometry g, ibgPoint *n1, ibgPoint *n2)
{ibgDLinear sl = (ibgDLinear ) 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]	= sl->r[d][0]*x1[0]
			+ sl->r[d][1]*x1[1]
			+ sl->r[d][2]*x1[2];
 	ibgpX(*n2)[d]	= sl->r[d][0]*x2[0]
			+ sl->r[d][1]*x2[1]
			+ sl->r[d][2]*x2[2];
 }
 rc = ibgiStatusOfEdge((sl->old),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, ibgPoint *nold)
{ibgDLinear sl = (ibgDLinear ) g;
 ibgFloat x0[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){x0[d]=ibgpX(*nnew)[d];}
 for(d=0;d<ibgDIM;d++){
 	ibgpX(*nnew)[d]	= sl->r[d][0]*x0[0]
			+ sl->r[d][1]*x0[1]
			+ sl->r[d][2]*x0[2];}
 rc = ibgiRegionOfPoint((sl->old),nnew,nold);
 for(d=0;d<ibgDIM;d++){ibgpX(*nnew)[d]=x0[d];}
 return rc;
}

static int FaceWithEdge(ibGeometry g, ibgPoint *nint,
					ibgPoint *n1, ibgPoint *n2)
{ibgDLinear sl = (ibgDLinear ) g;
 ibgFloat xo1[ibgDIM],xo2[ibgDIM],xon[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){xo1[d]=ibgpX(*n1)[d]; xo2[d]=ibgpX(*n2)[d];}
 for(d=0;d<ibgDIM;d++){
 	ibgpX(*n1)[d] = sl->r[d][0]*xo1[0]
		 + sl->r[d][1]*xo1[1]
		 + sl->r[d][2]*xo1[2];
	ibgpX(*n2)[d] = sl->r[d][0]*xo2[0]
		 + sl->r[d][1]*xo2[1]
		 + sl->r[d][2]*xo2[2];
 }
 rc = ibgiFaceWithEdge((sl->old),nint,n1,n2);
 for(d=0;d<ibgDIM;d++){
	ibgpX(*n1)[d]=xo1[d];
	ibgpX(*n2)[d]=xo2[d];
	xon[d]=ibgpX(*nint)[d];
 }
 for(d=0;d<ibgDIM;d++){
 	ibgpX(*nint)[d]  = sl->a[d][0]*xon[0]
			 + sl->a[d][1]*xon[1]
			 + sl->a[d][2]*xon[2];
 }
 return rc;
}

static int LineWithTriangle(ibGeometry g, ibgPoint *nint, ibgPoint *nface,
       	  	      ibgPoint *n1, ibgPoint *n2, ibgPoint *n3)
{ibgDLinear sl = (ibgDLinear ) g;
 ibgFloat xob[ibgDIM],xon[ibgDIM];
 ibgFloat xo1[ibgDIM],xo2[ibgDIM],xo3[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){xob[d]=ibgpX(*nface)[d];
 	xo1[d]=ibgpX(*n1)[d];
	xo2[d]=ibgpX(*n2)[d];
	xo3[d]=ibgpX(*n3)[d];
 }
 for(d=0;d<ibgDIM;d++){
	ibgpX(*nface)[d] = sl->r[d][0]*xob[0]
		 + sl->r[d][1]*xob[1]
		 + sl->r[d][2]*xob[2];
 	ibgpX(*n1)[d] = sl->r[d][0]*xo1[0]
		 + sl->r[d][1]*xo1[1]
		 + sl->r[d][2]*xo1[2];
	ibgpX(*n2)[d] = sl->r[d][0]*xo2[0]
		 + sl->r[d][1]*xo2[1]
		 + sl->r[d][2]*xo2[2];
	ibgpX(*n3)[d] = sl->r[d][0]*xo3[0]
		 + sl->r[d][1]*xo3[1]
		 + sl->r[d][2]*xo3[2];
 }
 rc = ibgiLineWithTriangle((sl->old),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];
	      xon[d]=ibgpX(*nint)[d];}
 for(d=0;d<ibgDIM;d++){
 	ibgpX(*nint)[d]   = sl->a[d][0]*xon[0]
		 	  + sl->a[d][1]*xon[1]
		 	  + sl->a[d][2]*xon[2];
 }
 return rc;
}

static int NodeInTetrahedron(ibGeometry g, ibgPoint *nint, ibgPoint *nline,
       		      ibgPoint *n1, ibgPoint *n2, ibgPoint *n3, ibgPoint *n4)
{ibgDLinear sl = (ibgDLinear ) g;
 ibgFloat xob[ibgDIM],xon[ibgDIM];
 ibgFloat xo1[ibgDIM],xo2[ibgDIM],xo3[ibgDIM],xo4[ibgDIM];
 int d,rc;
 for(d=0;d<ibgDIM;d++){xob[d]=ibgpX(*nline)[d];
 	xo1[d]=ibgpX(*n1)[d];
	xo2[d]=ibgpX(*n2)[d];
	xo3[d]=ibgpX(*n3)[d];
	xo4[d]=ibgpX(*n4)[d];
 }
 for(d=0;d<ibgDIM;d++){
	ibgpX(*nline)[d] = sl->r[d][0]*xob[0]
		 + sl->r[d][1]*xob[1]
		 + sl->r[d][2]*xob[2];
 	ibgpX(*n1)[d] = sl->r[d][0]*xo1[0]
		 + sl->r[d][1]*xo1[1]
		 + sl->r[d][2]*xo1[2];
	ibgpX(*n2)[d] = sl->r[d][0]*xo2[0]
		 + sl->r[d][1]*xo2[1]
		 + sl->r[d][2]*xo2[2];
	ibgpX(*n3)[d] = sl->r[d][0]*xo3[0]
		 + sl->r[d][1]*xo3[1]
		 + sl->r[d][2]*xo3[2];
	ibgpX(*n4)[d] = sl->r[d][0]*xo4[0]
		 + sl->r[d][1]*xo4[1]
		 + sl->r[d][2]*xo4[2];
 }
 rc = ibgiNodeInTetrahedron((sl->old),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];
	      xon[d]=ibgpX(*nint)[d];}
 for(d=0;d<ibgDIM;d++){
 	ibgpX(*nint)[d]   = sl->a[d][0]*xon[0]
		 	  + sl->a[d][1]*xon[1]
		 	  + sl->a[d][2]*xon[2];
 }
 return rc;
}

static void Free(ibgPtObject sl)
{
 ibgdFree(((ibgDLinear )sl)->old);
}

static ibGeometryClassRec	ibgDLinearClass;

ibGeometry ibgdLinearTrans(ibGeometry	gold,
				ibgFloat a[ibgDIM][ibgDIM])
{
 ibgDLinear sl=(ibgDLinear)malloc(sizeof(ibgDLinearRec));
 double det;
 int d,e;
 ibgdInitialize((ibGeometry)sl,&ibgDLinearClass);
 det = a[0][0]*(((double)a[1][1])*a[2][2]-((double)a[1][2])*a[2][1])
     + a[0][1]*(((double)a[1][2])*a[2][0]-((double)a[1][0])*a[2][2])
     + a[0][2]*(((double)a[1][0])*a[2][1]-((double)a[1][1])*a[2][0]);
 for(d=0;d<ibgDIM;d++) for(e=0;e<ibgDIM;e++) sl->a[d][e]=a[d][e];
 sl->r[0][0]=(((double)a[1][1])*a[2][2]-((double)a[1][2])*a[2][1])/det;
 sl->r[0][1]=(((double)a[1][2])*a[2][0]-((double)a[1][0])*a[2][2])/det;
 sl->r[0][2]=(((double)a[1][0])*a[2][1]-((double)a[1][1])*a[2][0])/det;
 sl->r[1][0]=(((double)a[2][1])*a[0][2]-((double)a[2][2])*a[0][1])/det;
 sl->r[1][1]=(((double)a[2][2])*a[0][0]-((double)a[2][0])*a[0][2])/det;
 sl->r[1][2]=(((double)a[2][0])*a[0][1]-((double)a[2][1])*a[0][0])/det;
 sl->r[2][0]=(((double)a[0][1])*a[1][2]-((double)a[0][2])*a[1][1])/det;
 sl->r[2][1]=(((double)a[0][2])*a[1][0]-((double)a[0][0])*a[1][2])/det;
 sl->r[2][2]=(((double)a[0][0])*a[1][1]-((double)a[0][1])*a[1][0])/det;
 sl->old = gold;
 ibgDLinearClass.RegionOfPoint = RegionOfPoint;
 ibgDLinearClass.FaceWithEdge = FaceWithEdge;
 ibgDLinearClass.LineWithTriangle = LineWithTriangle;
 ibgDLinearClass.NodeInTetrahedron = NodeInTetrahedron;
 ibgDLinearClass.StatusOfEdge = StatusOfEdge;
 return (ibGeometry) sl;
}