ibgdbound.c

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

/* last edit: Ilja Schmelzer -------------- 20-OCT-1994 11:01:23.12	*/
/************************************************************************/
/*                                                                      */
/*  <<< 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")	*/
/*                                                                      */
/************************************************************************/
/* <<< IBGDBoundaryGrid >>> - Intersection-Based Geometry Description

		Using a simplicial boundary grid

We use here the spline package "ibgdSpline.h". So, this file can be used as an
example of a geometry definition using a set of splines to describe the
boundary.

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 very
stable. We recommend to use small random shifts if such problems occur. */

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

#define big1 1.e6

/* We use here the linear spline types ibgdSplineLine and ibgdSplineTriangle
defined in "ibgdSpline.h". For other spline types you have to define your own
spline types. The first component of the spline type must be equivalent to the
"visible" part of the spline type ibgdSpline (like a C++ derived class of the
basic class ibgdSpline). */

/* This is a function which initializes the spline data type. It has to contain
the initialization of the visible (basic class) part of the spline type. It can
be considered as the C++ constructor of the ibgdSplineLine class.
*/

void ibgdSplineLineNew(ibgdSplineLine *s,
	ibgFloat x1, ibgFloat y1,
	ibgFloat x2, ibgFloat y2,
	ibgSegment l, ibgSegment r, ibgSegment u)
{
 s->x1[0] = x1;
 s->x1[1] = y1;
 s->x2[0] = x2;
 s->x2[1] = y2;
 s->s.left = l;	        /* segment number of the left  region */
 s->s.right= r;	        /* segment number of the right region */
 s->s.face = u;	        /* segment number of the face itself */
 if(u==0) s->s.face = ibgDefaultFaceNr(r,l);
 else     s->s.face = u; 
/* Computation of a bounding rectangle (in 3D bounding cuboid): */
 if(x1<x2){
	s->s.xmin = x1; s->s.xmax = x2;
 }else{
	s->s.xmin = x2; s->s.xmax = x1;
 }
 if(y1<y2){
	s->s.ymin = y1; s->s.ymax = y2;
 }else{
	s->s.ymin = y2; s->s.ymax = y1;
 }
}

/* This function tests if the spline spl has an intersection with the segment
from v0 to v1. The coordinates of the first intersection from v0 have to be
returned in xi. The return values are ibgdSplineLeft / ibgdSplineRight
depending on the side of v0 or ibgdSplineNone if there is no intersection.
This function can be considered as a pure virtual function of the ibgdSpline
class which has to be implemented for every derived class like ibgdSplineLine,
with the first parameter as the this-pointer.
*/

int ibgdSplineLineTest(void* spl,
	   	ibgFloat* xi,ibgFloat* v0,ibgFloat* v1)
{ibgdSplineLine* s = (ibgdSplineLine*) spl;
 ibgDouble xx=v1[0]-v0[0],yy=v1[1]-v0[1];
 ibgDouble f1,f2,ff,gg,p,q,x10,x20,y10,y20;
 int rc;
 f1 = (x10 = s->x1[0]-v0[0]) * yy - (y10 = s->x1[1]-v0[1]) * xx;
 f2 = (x20 = s->x2[0]-v0[0]) * yy - (y20 = s->x2[1]-v0[1]) * xx;
 if(f1>0){
	if(f2>=0) return ibgdSplineNone;
	gg = x10*y20 - y10*x20;
	ff = f1 - f2;
	rc = ibgdSplineLeft;
 }else if(f1<0){
	if(f2<=0) return ibgdSplineNone;
	gg = y10*x20 - x10*y20;
	ff = f2 - f1;
	rc = ibgdSplineRight;
 }else	  	  return ibgdSplineNone;
 if(gg<=0)	  return ibgdSplineNone;
 if(gg>=ff)	  return ibgdSplineNone;
 p = (ff-gg)/ff;
 q =  gg    /ff;
 xi[0] = p * v0[0] + q * v1[0];
 xi[1] = p * v0[1] + q * v1[1];
 xi[2] = p * v0[2] + q * v1[2];
 return rc;
}


/* That's all we need to define a spline type. Now the implementation for the
3D case (spline type ibgdSplineTriangle): */

void ibgdSplineTriangleNew(ibgdSplineTriangle *s,
	ibgFloat x1, ibgFloat y1, ibgFloat z1,
	ibgFloat x2, ibgFloat y2, ibgFloat z2,
	ibgFloat x3, ibgFloat y3, ibgFloat z3,
	ibgSegment l, ibgSegment r, ibgSegment u)
{ibgFloat min,max;
 s->x1[0] = x1;
 s->x1[1] = y1;
 s->x1[2] = z1;
 s->x2[0] = x2;
 s->x2[1] = y2;
 s->x2[2] = z2;
 s->x3[0] = x3;
 s->x3[1] = y3;
 s->x3[2] = z3;
 s->s.left = l;
 s->s.right= r;
 s->s.face = u;
 min = x1; if(x2<min) min=x2; if(x3<min) min=x3; s->s.xmin = min;
 max = x1; if(x2>max) max=x2; if(x3>max) max=x3; s->s.xmax = max;
 min = y1; if(y2<min) min=y2; if(y3<min) min=y3; s->s.ymin = min;
 max = y1; if(y2>max) max=y2; if(y3>max) max=y3; s->s.ymax = max;
 min = z1; if(z2<min) min=z2; if(z3<min) min=z3; s->s.zmin = min;
 max = z1; if(z2>max) max=z2; if(z3>max) max=z3; s->s.zmax = max;
}


int ibgdSplineTriangleTest(void* This,
	   	ibgFloat* xi,ibgFloat* v0,ibgFloat* v1)
{ibgdSplineTriangle* s = (ibgdSplineTriangle*) This;
 ibgDouble xx=v1[0]-v0[0],yy=v1[1]-v0[1],zz=v1[2]-v0[2];
 ibgDouble x1=s->x1[0]-v0[0],y1=s->x1[1]-v0[1],z1=s->x1[2]-v0[2];
 ibgDouble x2=s->x2[0]-v0[0],y2=s->x2[1]-v0[1],z2=s->x2[2]-v0[2];
 ibgDouble x3=s->x3[0]-v0[0],y3=s->x3[1]-v0[1],z3=s->x3[2]-v0[2];
 ibgDouble f1,f2,f3,ff,gg,p,q;
 int rc;
 f1 = xx*(y1*z2-z1*y2) + yy*(z1*x2-x1*z2) + zz*(x1*y2-y1*x2);
 f2 = xx*(y2*z3-z2*y3) + yy*(z2*x3-x2*z3) + zz*(x2*y3-y2*x3);
 f3 = xx*(y3*z1-z3*y1) + yy*(z3*x1-x3*z1) + zz*(x3*y1-y3*x1);
 if(f1>0){
	if(f2<=0) return ibgdSplineNone;
	if(f3<=0) return ibgdSplineNone;
	gg = x1*(y2*z3-z2*y3) + y1*(z2*x3-x2*z3) + z1*(x2*y3-y2*x3);
	ff = f1 + f2 + f3;
        if(s->s.left <= 0) return ibgdSplineNone;
	rc = ibgdSplineLeft;
 }else if(f1<0){
	if(f2>=0) return ibgdSplineNone;
	if(f3>=0) return ibgdSplineNone;
	gg = - (x1*(y2*z3-z2*y3) + y1*(z2*x3-x2*z3) + z1*(x2*y3-y2*x3));
	ff = - (f1 + f2 + f3);
	if(s->s.right <= 0) return ibgdSplineNone;
	rc = ibgdSplineRight;
 }else	  	  return ibgdSplineNone;
 if(gg<=0)	  return ibgdSplineNone;
 if(gg>=ff)	  return ibgdSplineNone;
 p = (ff-gg)/ff;
 q =  gg    /ff;
 xi[0] = p * v0[0] + q * v1[0];
 xi[1] = p * v0[1] + q * v1[1];
 xi[2] = p * v0[2] + q * v1[2];
 return rc;
}

/* Here the constructor for the new geometry. It consist of the
allocation and initialization of the spline list and a call of ibgdSplineList
which creates the related geometry description */

ibGeometry ibgdBoundaryGrid(ibGrid *gg,ibgSegment outside)
{ibgdSplines *s;
 ibgdSplineLine *s2;
 ibgdSplineTriangle *s3;
 int c,t,u, reg;
 int fm,num,cl,cr,*cn;
 ibgPoint *n1,*n2,*n3;

/* a very small shift of the geometry to avoid exact match between grid points
and splines. Such a match usually leads to incorrect geometries because of
rounding errors. */
 ibgFloat xs = 0.1732e-6;
 ibgFloat ys = 0.1132e-6;
 ibgFloat zs = 0.1347e-6;
/* different modi of the ibgg grid format - with or without region cells */
 if(ibgridLastRegionCell(*gg)>0) reg = 1; else reg = 0;
/* an upper bound for the number of splines (for memory allocation): */
 fm = ibgridLastFaceCell(*gg) - ibgridFirstFaceCell(*gg) + 2;
/* initialization of a cycle over all face cells with counter num: */
 num=0;
 switch(ibgridGDIM(*gg)){
 case 1: ibgfatal;
	break;
 case 2:
/* creation of the spline list, space allocation for the splines: */
	s  = ibgdSplineListNew(2,fm,sizeof(ibgdSplineLine),ibgdSplineLineTest);
/* this is the array of splines we have to initialize now: */
	s2 = (ibgdSplineLine*) s->spline;
	foribgFaceCells(*gg,c,t,u){
		if(t != ibgc2Edge){ibgfatal; continue;}
/* left and right neighbour cell in a region: */
	       	cl = ibgridCellLeft(*gg,c,t);
		cr = ibgridCellRight(*gg,c,t);
		if(reg){
/* the region number of these cells. For a pure boundary grid (reg==0)
cl and cr are already these region numbers. */
		        if(cl>0) cl = ibgridCellSegment(*gg,cl);
		        if(cr>0) cr = ibgridCellSegment(*gg,cr);
		}else{
		        if(cl<0) cl = -cl;
			if(cr<0) cr = -cr;
		}
/* the points of the segment: */
		cn = ibgridCPointList(*gg,c);
		n1 = &ibgridPoint(*gg,cn[0]);
		n2 = &ibgridPoint(*gg,cn[1]);
/* now we initialize the spline with number num: */
		ibgdSplineLineNew(&s2[num],
		        ibgpX(*n1)[0]+xs, ibgpX(*n1)[1]+ys,
		        ibgpX(*n2)[0]+xs, ibgpX(*n2)[1]+ys,
		        cl,cr,u);
		num++;
	}endibgFaceCells(*gg,c,t,u);
	break;
 case 3:
/* analogical for the 3D case (spline type ibgdSplineTriangle): */
	s  = ibgdSplineListNew(3,fm,sizeof(ibgdSplineTriangle),
		        			ibgdSplineTriangleTest);
	s3 = (ibgdSplineTriangle*) s->spline;
	foribgFaceCells(*gg,c,t,u){
		if(t != ibgc3Triangle){ibgfatal; continue;}
	       	cl = ibgridCellLeft(*gg,c,t);
		cr = ibgridCellRight(*gg,c,t);
		if(reg){
		        if(cl>0) cl = ibgridCellSegment(*gg,cl);
		        if(cr>0) cr = ibgridCellSegment(*gg,cr);
		}else{
		        if(cl<0) cl = -cl;
			if(cr<0) cr = -cr;
		}
		cn = ibgridCPointList(*gg,c);
		n1 = &ibgridPoint(*gg,cn[0]);
		n2 = &ibgridPoint(*gg,cn[1]);
		n3 = &ibgridPoint(*gg,cn[2]);
		ibgdSplineTriangleNew(&s3[num],
		        ibgpX(*n1)[0]+xs, ibgpX(*n1)[1]+ys, ibgpX(*n1)[2]+zs,
		        ibgpX(*n2)[0]+xs, ibgpX(*n2)[1]+ys, ibgpX(*n2)[2]+zs,
		        ibgpX(*n3)[0]+xs, ibgpX(*n3)[1]+ys, ibgpX(*n3)[2]+zs,
		        cl,cr,u);
		num++;
	}endibgFaceCells(*gg,c,t,u);
	break;
 }
/* the number of splines in the list: */
 s->num = num;
 s->outside = outside;

/* Now we create the geometry description defined by the spline list.
Remark that this routine assumes that the splines in the list are already
initialized. For example, a search tree for fast search may be created in this
function. */

 return ibgdSplineList(s);
}