ibgdspline.c

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

/* last edit: Ilja Schmelzer -------------- 17-OCT-1994 19:41:59.20	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDSpline >>> - Intersection-Based Geometry Description

	Using a list of splines for boundary description

See "ibgdBound.c" for an example how to use the spline package.

The spline package was developed to allow the usage of existent geometry
descriptions with boundary grids or splines in IBGG. In general, we don't
recommend to use boundary splines for geometry description. Even a small
mismatch in the boundary definition or rounding errors can lead to global
errors in the  geometry description so that the grid generation is not stable.
We recommend to try small random shifts in the case of errors. */

/* The algorithm implemented here is not optimal --- for a (big) number n of
splines there will be O(n) behaviour for one request. Optimal are
quadtree/octree techniques with O(log n) behaviour. I will implement them if
they are necessary for someone. */

#include "ibg.h"
#include "ibgdefault.h"
#include "ibgi.h"
#include "ibglib.h"
#include "ibgd.h"
#include "ibgd0.h"
#include "ibgdspline.h"

#define big1 1.e6

static int Intersection(ibgdSplines* s, int *i,
		ibgFloat* xi,ibgFloat* x0,ibgFloat* x1)
{int r0 = ibgdSplineNone,rc,dim3,lx,ly,lz,size,j;
 ibgdSpline *sp;
 char* sc;
 ibgFloat xx0,xx1,yy0,yy1,zz0,zz1;
 if(s->dim==3)		dim3 = 1;
 else	 		dim3 = 0;
 if(x0[0]<x1[0]){	lx=1; xx0=x0[0]; xx1=x1[0];
 }else{			lx=0; xx0=x1[0]; xx1=x0[0];
 }
 if(x0[1]<x1[1]){	ly=1; yy0=x0[1]; yy1=x1[1];
 }else{			ly=0; yy0=x1[1]; yy1=x0[1];
 }
 if(dim3){
	 if(x0[2]<x1[2]){	lz=1; zz0=x0[2]; zz1=x1[2];
	 }else{			lz=0; zz0=x1[2]; zz1=x0[2];
	 }
 }
 size = s->size;
 sc   = s->spline;
 for(j=0;j<s->num;j++,sc+=size){
/*	#define sp ((ibgdSpline*) sc) */
	sp  = (ibgdSpline*) sc;
	if(xx0 > sp->xmax)     	continue;
	if(xx1 < sp->xmin)	continue;
	if(yy0 > sp->ymax)	continue;
	if(yy1 < sp->ymin)	continue;
	if(dim3){
		if(zz0 > sp->zmax)	continue;
		if(zz1 < sp->zmin)	continue;
	}
	if(rc = s->test((void*)sc,xi,x0,x1)){
		r0 = rc; *i = j;
		if(lx) xx1 = x1[0] = xi[0];	else xx0 = x1[0] = xi[0];
		if(ly) yy1 = x1[1] = xi[1];	else yy0 = x1[1] = xi[1];
		if(dim3){
			if(lz) zz1 = x1[2] = xi[2];
			else   zz0 = x1[2] = xi[2];
		}
	}
 }
 return r0;
}

static int RegionOfPoint(ibGeometry g, ibgPoint *nnew, const ibgPoint *nold)
{ibgdSplines* s = (ibgdSplines*) g;
 ibgSegment u0; ibgFloat x0[3],xi[3];
 int i,rc;
 if(nold && (ibgpType(*nold) == ibgSRegion)){
	u0 = ibgpSegment(*nold);
	x0[0] = ibgpX(*nold)[0];
	x0[1] = ibgpX(*nold)[1];
	x0[2] = ibgpX(*nold)[2];
 }else{
	u0 = s->outside;
	x0[0] = ibgpX(*nnew)[0];
	x0[1] = big1;
	x0[2] = ibgpX(*nnew)[2];
 }
 if(rc = Intersection(s,&i,xi,ibgpX(*nnew),x0)){
	if(rc==ibgdSplineLeft)
		u0 = ((ibgdSpline*)(s->spline+s->size*i))->left;
	else
		u0 = ((ibgdSpline*)(s->spline+s->size*i))->right;
 }
 ibgpSegment(*nnew)=u0;
 ibgpType(*nnew)=ibgSRegion;
 return ibgRegionFound;
}

static int FaceWithEdge(ibGeometry g, ibgPoint *nint,
	   		   const ibgPoint *n1, const ibgPoint *n2)
{ibgdSplines* s = (ibgdSplines*) g;
 ibgFloat x[3];
 int i,rc;
 x[0] = ibgpX(*n2)[0];
 x[1] = ibgpX(*n2)[1];
 x[2] = ibgpX(*n2)[2];
 if(Intersection(s,&i,ibgpX(*nint),ibgpX(*n1),x)){
	ibgpSegment(*nint) = ((ibgdSpline*)(s->spline+s->size*i))->face;
	ibgpType(*nint) = ibgSFace;
  	return  ibgFaceFound;
 }
 ibgpType(*nint) = ibgSNothing;
 return ibgNotFound;
}

static int Free(ibGeometry g)
{ibgdSplines* s = (ibgdSplines*) g;
 free(s->spline);
 return ibgSuccess;
}

static ibGeometryClassRec	ibgDSplineClass;

ibgdSplines* ibgdSplineListNew(int dim, int max, int size,
        	int (*test)(void* spline, ibgFloat* xi,
        	      	 ibgFloat* x1, ibgFloat* x2)
        )
{ibgdSplines *s = (ibgdSplines*) malloc(sizeof(ibgdSplines));
 s->dim  = dim;
 s->max  = max;
 s->num  = 0;
 s->size = size;
 s->test = test;
 s->spline	= (char*) malloc(size * max);
 s->outside	= 0;
 return s;
}


ibGeometry ibgdSplineList(ibgdSplines *spp)
{ibgDouble dx,dy,dz;
 int size,j,dim3;
 char *sc;
 if(spp->dim==3)	dim3 = 1;
 else	 		dim3 = 0;
 size = spp->size;
 sc   = spp->spline;
 for(j=0;j<spp->num;j++,sc+=size){
/*#define sp ((ibgdSpline*) sc) */
        ibgdSpline *sp = ((ibgdSpline*) sc);
	dx = 0.0001*(sp->xmax - sp->xmin) + 0.00001;
	sp->xmax += dx;
	sp->xmin -= dx;
	dy = 0.0001*(sp->ymax - sp->ymin) + 0.00001;
	sp->ymax += dy;
	sp->ymin -= dy;
	if(dim3){
		dz = 0.0001*(sp->zmax - sp->zmin) + 0.00001;
		sp->zmax += dz;
		sp->zmin -= dz;
	}
 }
 ibgdInitialize((ibGeometry)spp,&ibgDSplineClass);
 ibgDSplineClass.RegionOfPoint	= RegionOfPoint;
 ibgDSplineClass.FaceWithEdge= FaceWithEdge;
 ibgDSplineClass.Free		= Free;
 sc   = spp->spline;
 for(j=0;j<spp->num;j++,sc+=size){
        ibgdSpline *sp = ((ibgdSpline*) sc);
	if(sp->face != ibgDefaultFaceNr(sp->left,sp->right))
	    ibgtNewFace(&(spp->g.top),
			sp->face,ibgDefaultFaceNr(sp->left,sp->right));
 }
 return (ibGeometry)spp;
}