ibggmake.c

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

/* last edit: Ilja Schmelzer -------------- 10-JUN-1994 15:05:53.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")	*/
/*                                                                      */
/************************************************************************/
#include "ibgi.h"
#include "ibgd0.h"
#include "ibggp.h"
#include "ibgg.h"
#include "ibgdefault.h"
#include "ibglib.h"

	ibGeometry	ibGeo;	/* current geometry */

#define ibgFound 0
#define ibgNotFound (-1)

static int g2swapm(int c, int i);
static void g2swapd(int c);
static int g3EDelete(int cb, int eb);

static ibgFloat ibggVAbs = eps3, ibggVRel = 0.05;
static int delaunay_on=0;
static int unc,mec,mbc,mrc;
static int g3half=1;
static int mst,*st,*st1,nosegmentb,nosegment0,nosegment1,nosegment2;
static int ibgcurrc;

/* defintion of segment of a correct element */
int ibgcRegion(int c)
{int n,nh,k,nn,m,u,*cn,i,hfirst=1,notff=0;
 ibgCellType t;
 ibgSegmentType tu;
 ibgPoint xx,*xa;
 if(ccu(c)<0){ ibgridNFree(ibgg,-ccu(c)); ccu(c) = 0;}
 cn=ccn(c); t=cct(c);
 for(n=0;n<ibgcPoints(t);n++){
	switch(cnt(nn=cn[n])){
    	case ibgSRegion:	if((m=cnu(nn)) <= 0){
					ccu(c) = ibgOutsideRegionNr;
					return ibgFound;
    				}
				xa = cnd(nn);
				goto mfound;
	case ibgSLine:		goto minside;
	case ibgSNode:	goto minside;
	}
 }
minside:
 for(k=0;k<ibgDIM;k++)  xx.xx[k] = 0;
 for(n=0;n<ibgcPoints(t);n++){nn=cn[n];
	if((cnu(nn)) <= 0){
	   	ccu(c) = ibgOutsideRegionNr;
		return ibgFound;
	}
 	for(k=0;k<ibgDIM;k++)  xx.xx[k] += cnx(nn)[k];
 }
 for(k=0;k<ibgDIM;k++)  xx.xx[k] /= ibgcPoints(t);
 ibgiRegionOfPoint(ibGeo,&xx,cnd(nn));
 m  = ibgpSegment(xx);
 xa = &xx;
mfound:
 ibgassert(m>0);
 for(n=0;n<ibgcPoints(t);n++){
 	u=cnu(nn=cn[n]);
	switch(cnt(nn)){
	case ibgSRegion:
		if(u <= 0) {ccu(c) = ibgOutsideRegionNr;
				       		return ibgFound;}
		if(u == m) continue;
	     	notff = nn;
	default:
	       	if(ibgiStatusOfEdge(ibGeo,xa,cnd(nn))!=ibgIncorrect) continue;
	     	notff = nn;
	}
 }
 if(notff==0){
	ccu(c)=m;
	return ibgFound;
 }
mnotfound:
 if(&xx == xa && hfirst){
	hfirst = 0;
  	/* may be a point near the dangerous corner will be better? */
	for(i=0;i<10;i++){
	  for(k=0;k<ibgDIM;k++)  xx.xx[k] = 0.5*(xx.xx[k]+ibgpX(*cnd(nn))[k]);
	  ibgiRegionOfPoint(ibGeo,&xx,cnd(nn));
	  if(ibgpSegment(xx)!=m){
		m = ibgpSegment(xx);
		goto mfound;
	  }
	}
 }
 ibgridNGet(ibgg,nh);
 *cnd(nh) = *xa;
 ccu(c) = -nh;
 return ibgNotFound;
}

/* definition of segment for an incorrect element */
void ibgcSetRegion(int c)
{int n = -ccu(c);
 ibgassert(n>0);
 ibgassert(n<=ibgg.lastPoint);
 ccu(c) = ibgpSegment(ibgridPoint(ibgg,n));
 ibgridNFree(ibgg,n);
}

static int gmlev;

void g2segment(int c)
{int i,j,co,*cn=ccn(c);
 ibgCellType t=cct(c);
 ibgassert(t==ibgc2Triangle);
/* Korrektur des falschen Dreiecks: */
 for(i=0;i<3;i++){
	if(ibgIncorrect != ibgiStatusOfEdge(ibGeo,
			cnd(cn[ibgcSPoint1(t,i,0)]),
 			cnd(cn[ibgcSPoint1(t,i,1)]))) continue;
	co=ccs(c)[i]; if(ccout(co)) continue;
	if(gmlev==0){
		for(j=0;j<3;j++) if(ccs(co)[j]==c) break;
		if(ibgIncorrect != ibgiStatusOfEdge(ibGeo,cnd(cn[i]),cnd(ccn(co)[j])))
				   		if(g2swapm(c,i)) return;
	}else{
		gmlev--;
		if(g2swapm(c,i)) {gmlev++; return;}
		gmlev++;
	}
 }
}

static void g2swapd(int c)
{int *cn,*cno,*cs,*cso,*csb,nn,n0,n1,no,co,b,bo,i,i0,i1,j,j0,j1,k,m;
 ibgFloat x0,x1,xn,xo,y0,y1,yn,yo,sinn,sino,cosn,coso;
 ibgCellType t=cct(c);
 m=ccu(c); cs=ccs(c); cn=ccn(c);
 for(i=0;i<3;i++){
	if(m!=ccu((co=cs[i])))	continue;
 	cso=ccs(co); cno=ccn(co);
	if(cct(co)!=ibgc2Triangle)	continue;
	for(j=0;j<3;j++) if(cso[j]==c) goto cfound;
	ibgfatal; continue;
cfound:
	i0=ibgcSPoint1(t,i,0); i1=ibgcSPoint1(t,i,1);
	nn=cn[i]; n0=cn[i0]; n1=cn[i1]; no=cno[j];
	x0=cnx(n0)[0]; y0=cnx(n0)[1];
	x1=cnx(n1)[0]; y1=cnx(n1)[1];
	xn=cnx(nn)[0]; yn=cnx(nn)[1];
	xo=cnx(no)[0]; yo=cnx(no)[1];
	sinn  =	((x0-xn)*(y1-yn) - (x1-xn)*(y0-yn));	if(sinn<0) continue;
	sino  =	((x1-xo)*(y0-yo) - (x0-xo)*(y1-yo));	if(sino<0) continue;
	cosn =	((x0-xn)*(x1-xn) + (y0-yn)*(y1-yn));
	coso =	((x0-xo)*(x1-xo) + (y0-yo)*(y1-yo));
	if(cosn*sino+coso*sinn>-eps2) continue;
	j0=ibgcSPoint1(t,j,0); j1=ibgcSPoint1(t,j,1);
	b=cs[i0]; bo=cso[j0];
	cn[i1]=no; cno[j1]=nn;
	cs[i0]=co; cso[j0]=c;
	cs[i] =bo; cso[j] =b;
	if(b > 0){
		m=ibgcPoints(cct(b ));csb=ccs(b );
		for(k=0;k<m;k++) if(csb[k]==c )	{csb[k] =co; break;}
	}
	if(bo> 0){
		m=ibgcPoints(cct(bo));csb=ccs(bo);
		for(k=0;k<m;k++) if(csb[k]==co)	{csb[k] = c; break;}
  	}
	g2swapd(c);
	g2swapd(co);
	return;
 }
}

static int g2swapm(int c, int i)
{int *cn,*cno,*cs,*cso,*csb,nn,n0,n1,no,co,b,bo,i0,i1,j,j0,j1,k,m;
 ibgFloat x,y,d;
 ibgCellType t=cct(c);
 cn=ccn(c); nn=cn[i];
 cs=ccs(c); co=cs[i]; cso=ccs(co); cno=ccn(co);
 if(ccout(co)) return 0;
 if(cct(co)!=ibgc2Triangle) ibgfatal;
 for(j=0;j<3;j++) if(cso[j]==c) goto cfound;
 ibgfatal; return 0;
cfound:
 i0=ibgcSPoint1(t,i,0); i1=ibgcSPoint1(t,i,1);
 n0=cn[i0]; n1=cn[i1]; no=cno[j];
 x=cnx(n0)[0]; y=cnx(n0)[1];
 d=(cnx(no)[0]-x)*(cnx(nn)[1]-y)-(cnx(nn)[0]-x)*(cnx(no)[1]-y);
 if(d<eps2) return 0;
 x=cnx(n1)[0]; y=cnx(n1)[1];
 d=(cnx(no)[0]-x)*(cnx(nn)[1]-y)-(cnx(nn)[0]-x)*(cnx(no)[1]-y);
 if(d>-eps2) return 0;
 j0=ibgcSPoint1(t,j,0); j1=ibgcSPoint1(t,j,1);
 b=cs[i0]; bo=cso[j0];
 cn[i1]=no; cno[j1]=nn;
 cs[i0]=co; cso[j0]=c;
 cs[i] =bo; cso[j] =b;
 if(b > 0){
	m=ibgcPoints(cct(b ));csb=ccs(b );
	for(k=0;k<m;k++) if(csb[k]==c )	{csb[k] =co; break;}
 }
 if(bo> 0){
	m=ibgcPoints(cct(bo));csb=ccs(bo);
	for(k=0;k<m;k++) if(csb[k]==co)	{csb[k] = c; break;}
 }
 if(ibgcRegion(c)==ibgNotFound)	g2segment(c);
 else				g2swapd(c);
 if(ibgcRegion(co)==ibgNotFound)g2segment(co);
 else				g2swapd(co);
 return co;
}
/* Berechnen und Eintragen der Daten der Umkugel: */
static float g3sdelmin=0.1, g3sdelmax=0.9;
int g3cvol(int c)
{int *cn0;
 ibgFloat *cx0,*cx1,*cx2,*cx3;
 ibgFloat x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3,l1,l2,l3,a1,a2,a3,mx,my,mz,det;
 cn0=ccn(c);
 cx0=cnx(*cn0++); cx1=cnx(*cn0++); cx2=cnx(*cn0++); cx3=cnx(*cn0++);
 x  = cx0[0];   y  = cx0[1];   z  = cx0[2];
 x1 = cx1[0]-x; y1 = cx1[1]-y; z1 = cx1[2]-z;
 x2 = cx2[0]-x; y2 = cx2[1]-y; z2 = cx2[2]-z;
 x3 = cx3[0]-x; y3 = cx3[1]-y; z3 = cx3[2]-z;
 l1=x1*x1 + y1*y1 + z1*z1;
 l2=x2*x2 + y2*y2 + z2*z2;
 l3=x3*x3 + y3*y3 + z3*z3;
 a1=y2*z3 - y3*z2;
 a2=y3*z1 - y1*z3;
 a3=y1*z2 - y2*z1;
 det=x1*a1 + x2*a2 + x3*a3;
 if(det>-eps3) return 0; /* Entartung */
 mx=(l1*a1+l2*a2+l3*a3)/det;
 my=(l1*(z2*x3-z3*x2)+l2*(z3*x1-z1*x3)+l3*(z1*x2-z2*x1))
 		/det;
 mz=(l1*(x2*y3-x3*y2)+l2*(x3*y1-x1*y3)+l3*(x1*y2-x2*y1))
 		/det;
 ccv(c).x1[0]=x;
 ccv(c).x1[1]=y;
 ccv(c).x1[2]=z;
 ccv(c).x2[0]=x+mx;
 ccv(c).x2[1]=y+my;
 ccv(c).x2[2]=z+mz;
 return 1;
}

static int g3segment1edge;
int g3CHalf(int c)
{
 ibgCellType t=cct(c);
 int s1,s2,sr,sl,so,sh,c1,co,cu,ce,cn,cr,n1,n2,nh,e,i,rc;
 int *cn0=ccn(c);
 ibgFloat *cxx;
 ibgFloat x,y,z,xh,yh,zh,x2,y2,z2,d;
 if(t!=ibgc3Tetrahedron) return 0;
 for(e=0;e<ibgcEdges(t);e++){
	if(ibgIncorrect != ibgiStatusOfEdge(ibGeo,cnd(cn0[ibgcEPoint1(t,e)])
			   	       ,cnd(cn0[ibgcEPoint2(t,e)]))) continue;
	n1=cn0[s1=ibgcEPoint1(t,e)];
	n2=cn0[s2=ibgcEPoint2(t,e)];
	if(cnt(n1)!=ibgSRegion) continue;
	if(cnt(n2)!=ibgSRegion) continue;
	sr=ibgcESideR(t,e);
	sl=ibgcESideL(t,e);
	goto newpoint;
 }
/*
 for(e=0;e<ibgcEdges(t);e++){
	if(ibgIncorrect != ibgiStatusOfEdge(ibGeo,cnd(cn0[ibgcEPoint1(t,e)])
				       ,cnd(cn0[ibgcEPoint2(t,e)]))) continue;
	n1=cn0[s1=ibgcEPoint1(t,e)];
	n2=cn0[s2=ibgcEPoint2(t,e)];
	if(cnt(n1)!=ibgSRegion) {
 		if(cnt(n2)!=ibgSRegion) continue;
		nh=n1; n1=n2; n2=nh;
		sh=s1; s1=s2; s2=sh;
		sr=ibgcESideL(t,e);
		sl=ibgcESideR(t,e);
	} else {sr=ibgcESideR(t,e);
	 	sl=ibgcESideL(t,e);
	}
	goto newpoint;
 }
 for(e=0;e<ibgcEdges(t);e++){
	if(ibgIncorrect != ibgiStatusOfEdge(ibGeo,cnd(cn0[ibgcEPoint1(t,e)])
				       ,cnd(cn0[ibgcEPoint2(t,e)]))) continue;
	n1=cn0[s1=ibgcEPoint1(t,e)];
	n2=cn0[s2=ibgcEPoint2(t,e)];
	ibgfatal;
	sr=ibgcESideR(t,e);
	sl=ibgcESideL(t,e);
	goto newpoint;
 }
 ibgfatal;
*/
 return 0;
newpoint:
 ibgridNGet(ibgg,nh);
 rc = ibgiFaceWithEdge(ibGeo,cnd(nh),cnd(n1),cnd(n2));
 if(rc != ibgFaceFound) {
	ibgfatal;
	return 0;
 }
/* test if nh is near n1 or n2: */
 cxx = cnx(n1); x  = cxx[0];	y  = cxx[1];	z  = cxx[2];
 cxx = cnx(nh); xh = cxx[0]-x;	yh = cxx[1]-y;	zh = cxx[2]-z;
 cxx = cnx(n2); x2 = cxx[0]-x;	y2 = cxx[1]-y;	z2 = cxx[2]-z;
 if((d=(xh*x2+yh*y2+zh*z2)/(x2*x2+y2*y2+z2*z2))<g3sdelmin){
	cxx = cnx(n1);
	x  = cxx[0];	y  = cxx[1];	z  = cxx[2];
       	*cnd(n1) = *cnd(nh);
	cxx[0] = x;	cxx[1] = y;	cxx[2] = z;
	ibgridNFree(ibgg,nh);
	if(ibgcRegion(c)==ibgNotFound) return 0;
	return 1;
 }else if(d>g3sdelmax){
	cxx = cnx(n2);
	x  = cxx[0];	y  = cxx[1];	z  = cxx[2];
	*cnd(n2) = *cnd(nh);
	cxx[0] = x;	cxx[1] = y;	cxx[2] = z;
	ibgridNFree(ibgg,nh);
	if(ibgcRegion(c)==ibgNotFound) return 0;
	return 1;
 }
/* ibgc3Tetrahedron-specials used: side == point not in this side */
 c1 = c; co = 0;
 cr = ccs(c1)[sr];
 do{
	g0CGet(&cn,ibgc3Tetrahedron);
	if(cn<=0) {ibgfatal; return 0;}
	ccn(cn)[s1]=nh;
	ccn(cn)[s2]=n2;
	ccn(cn)[sr]=ccn(c1)[sr];
	ccn(cn)[sl]=ccn(c1)[sl];
	ccn(c1)[s2]=nh;
	ccs(cn)[s1]=ce=ccs(c1)[s1];
	for(i=0;i<ibgcSides(cct(ce));i++)
		if(ccs(ce)[i]==c1) {ccs(ce)[i]=cn; break;}
	ccs(cn)[s2]=c1;
	if(co>0){
		ccs(cn)[sl]=co;
		ccs(co)[so]=cn;
	}
 	ccs(c1)[s1]=cn;
	if(ibgcRegion(c1)==ibgNotFound){
		if(st[c1]==0){st[c1] = nosegment0; nosegment0 = c1;}
	}
	if(ibgcRegion(cn)==ibgNotFound){
		if(cn>=mst){mst=(3*mst)/2;
			st1=(int*)realloc(st,mst*sizeof(int));
	 		if(st1!=ibgNULL) st=st1; else ibgfatal;
		}
		if(st[cn]==0){ st[cn] = nosegment0; nosegment0 = cn;}
	}
	g3cvol(c1);	g3cvol(cn);
/* going to the next element: */
	co = cn; so = sr; cu = c1; c1 = ccs(c1)[sr];
	s1=s2=sr=sl= -1;
	for(i=0;i<4;i++){
		if     (ccn(c1)[i]==n1)	s1=i;
		else if(ccn(c1)[i]==n2) s2=i;
		else if(ccs(c1)[i]==cu) sl=i;
		else  			sr=i;
	}
	ibgassert(s1>=0);	if(c1!=c) ibgassert(s2>=0); ibgassert(sr>=0);	ibgassert(sl>=0);
 }while(c1!=c);
 cn = ccs(c1)[s1];
 ccs(cn)[sl]=co;
 ccs(co)[so]=cn;
 return 1;
}

int g3SDelete(int cb, int sb)
{ibgCellType tb=cct(cb),ts;
 int ub,us,un,cs,cn,co,ss,ss0,ss1,ss2,sb0,sb1,sb2,sh,nb,n0,n1,n2,ns,nh,i;
 ibgFloat *cxx;
 ibgFloat x,y,z,x0,y0,z0,x1,y1,z1,x2,y2,z2,xs,ys,zs,det,det0,vol0;
 cs=ccs(cb)[sb]; ts=cct(cs);
 if(ccout(cs)) return 0;
 for(ss=0;ss<ibgcSides(ts);ss++) if(ccs(cs)[ss]==cb) break;
 nb  = ccn(cb)[sb];
 ns  = ccn(cs)[ss];
 n0  = ccn(cb)[sb0=ibgcSPoint1(tb,sb,0)];
 n1  = ccn(cb)[sb1=ibgcSPoint1(tb,sb,1)];
 n2  = ccn(cb)[sb2=ibgcSPoint1(tb,sb,2)];
 if(ibgIncorrect == ibgiStatusOfEdge(ibGeo,cnd(nb),cnd(ns))) return 0;
 cxx = cnx(nb); x  = cxx[0];	 y  = cxx[1];	  z  = cxx[2];
 if(delaunay_on){			/* Umkugeltest:  */
	 if((ccv(cs).x1[0]-x)*(ccv(cs).x2[0]-x)
	  + (ccv(cs).x1[1]-y)*(ccv(cs).x2[1]-y)
	  + (ccv(cs).x1[2]-z)*(ccv(cs).x2[2]-z) > -eps2) return 0;