📄 ibggdel.c
字号:
/* last edit: Ilja Schmelzer -------------- 24-MAR-1995 19:53:26.67 */
/************************************************************************/
/* */
/* <<< I B G >>> - Intersection - Based Grid generation package */
/* */
/* Version 1.1 by Ilja Schmelzer schmelzer@iaas-berlin.d400.de */
/* */
/* to be distributed under IBG license conditions (see "readme.ibg") */
/* */
/************************************************************************/
#include "ibgi.h"
#include "ibgg0.h"
#include "ibgd0.h"
#include "ibgg.h"
#include "ibgdefault.h"
#include "ibglib.h"
#include "ibgg0p.h"
ibgFloat ibggDelaunayBig = 1.e3;
ibgFloat ibggDelaunayEps = 1.e-6;
int ibggDelaunay0(ibGrid *g, ibgPoints *list);
int ibggDelaunay1(ibGrid *g, ibgPoints *list);
int ibggDelaunay2(ibGrid *g, ibgPoints *list);
int ibggDelaunay3(ibGrid *g, ibgPoints *list);
int ibggDelaunay1(ibGrid *g, ibgPoints *list){return 0;}
/* int ibggDelaunay3(ibGrid *g, ibgPoints *list){return 0;} */
int ibgAppendPoint(ibgPoints *list, ibgPoint *p)
{int l,l0,i;
if(l=(list->lastPoint + 1) > list->maxPoint){
void* adr;
l = 2*l+50;
adr = realloc(list->Point,(l+2)*sizeof(ibgPoint));
if(adr) list->Point = (ibgPoint*) adr;
else return ibgError;
adr = realloc(list->Near,(l+2)*sizeof(int));
if(adr) list->Near = (int*) adr;
else return ibgError;
list->maxPoint = l+1;
}
l = ++list->lastPoint;
list->Point[l] = *p;
list->Near[l] = 0;
return ibgSuccess;
}
int ibgAppendPoints(ibgPoints *list, ibgPoints *rest)
{int l,l0,i;
if(l=(list->lastPoint + rest->lastPoint - rest->firstPoint + 1)
> list->maxPoint){void* adr;
adr = realloc(list->Point,(l+2)*sizeof(ibgPoint));
if(adr) list->Point = (ibgPoint*) adr;
else return ibgError;
adr = realloc(list->Near,(l+2)*sizeof(int));
if(adr) list->Near = (int*) adr;
else return ibgError;
list->maxPoint = l+1;
}
l = list->lastPoint;
l0 = l - rest->firstPoint + 1;
for(i=rest->firstPoint;i<=rest->lastPoint; i++){l++;
list->Point[l] = rest->Point[i];
if(rest->Near) list->Near[l] = rest->Near[i] + l0;
else list->Near[l] = 0;
}
list->lastPoint = l;
return ibgSuccess;
}
int ibggDelaunay(ibGrid *g, int dim, ibgPoints *list)
{
switch(dim){
case 0: return ibggDelaunay0(g,list);
case 1: return ibggDelaunay1(g,list);
case 2: return ibggDelaunay2(g,list);
case 3: return ibggDelaunay3(g,list);
}
}
/* Delaunay-Triangulierung des Gitters */
static ibgFloat big;
static int ibgg2find(ibgFloat *xx, int c0)
{register int nl,*cn0;
ibgCellType t;
int i,sl;
ibgFloat vv,x1,x2,y1,y2;
ibgFloat* ff[IBGCSMAX];
register ibgFloat *fj;
if(ccu(c0)<0){
c0 = 1;
while (ccu(c0)<0) c0++;
}
beg:
ibgassert(ccu(c0)>=0);
nl=ibgcPoints(t=cct(c0)); cn0=ccn(c0);
for(i=0;i<nl;i++) ff[i]=cnx(cn0[i]);
for(i=0;i<nl;i++){
fj = ff[ibgcSPoint1(t,i,0)];
x1 = fj[0]-xx[0]; y1=fj[1]-xx[1];
fj = ff[ibgcSPoint1(t,i,1)];
x2 = fj[0]-xx[0]; y2=fj[1]-xx[1];
vv = x1 * y2 - x2 * y1;
if(vv<-eps2) {c0=ccs(c0)[i]; goto beg;}
}
return c0;
}
/* Hilfsdaten f"ur Startsituation */
static int hcn2[6] = {1,3,2, 1,2,4};
static int hcs2[6] = {0,2,0, 0,0,1};
static int determ=1,dslast=0,nslast=0;
int ibggDelaunay2(ibGrid *g, ibgPoints *list)
{int n,b,s,s0,sh,ms,i,i0,ih,j,d,si,sl,del;
int n0,nn,no,nt,nh,nl,c0,c1,cc,ct,ch,co,cl;
int *cs0,*cs1,*csi,*csj,*csh,*cn0,*cni,*cnj,*cnh,*ss,*se;
ibgFloat *cx,*cx1,*cx2,*cx3;
ibgFloat x,y,z,xm,ym,zm,dx,dy,dz,dd,x1,x2,x3,y1,y2,y3,z1,z2,z3,vv;
ibgFloat l1,l2,l3,s1,s2,s3,h1,h2,h3,n1,n2,n3,det,a1,a2,a3,mx,my,mz;
ibgCellType t,t0,t1,th,tt;
ibgFloat xmin[3],xmax[3],xbig,ybig,eps;
ibgSegmentType tg; ibgPoint nd; int bb;
tcv v,cvert;
#define DSMAX 500
int dsc[DSMAX],dscl[DSMAX][IBGCSMAX],dsl,dsi,ds0,dsh;
#define NSMAX 600
int nsc[NSMAX],nscs[NSMAX],nsco[NSMAX],nsl,nsi;
ibgg.spaceDimension = list->spaceDimension;
ibgg.firstPoint = list->firstPoint;
ibgg.lastPoint = list->lastPoint;
ibgg.maxPoint = list->maxPoint;
ibgg.Point = list->Point;
/* Speicherbereitstellung: */
nl = ibgg.lastPoint;
ibgg.lastPoint += 4; if(ibgridNRealloc(&ibgg)==0) return ibgError;
g0CAlloc(2*ibgg.lastPoint);
/* Berechnen des Rahmens */
xmin[0] = xmax[0] = cnx(1)[0];
xmin[1] = xmax[1] = cnx(1)[1];
for(i=2;i<=nl;i++){
if(cnx(i)[0] < xmin[0]) xmin[0] = cnx(i)[0];
if(cnx(i)[0] > xmax[0]) xmax[0] = cnx(i)[0];
if(cnx(i)[1] < xmin[1]) xmin[1] = cnx(i)[1];
if(cnx(i)[1] > xmax[1]) xmax[1] = cnx(i)[1];
}
/* Aufbau der Startsituation:
4 unendlich ferne Knoten vom Mittelpunkt aus am Ende dazu.
Sie bilden 2 Dreiecke um die x-Achse herum.
*/
xm=0.5*(xmin[0]+xmax[0]); dx=(xmax[0]-xmin[0]);
ym=0.5*(xmin[1]+xmax[1]); dy=(xmax[1]-xmin[1]);
dd=dx*dx+dy*dy;
xbig = dx * ibggDelaunayBig;
ybig = dy * ibggDelaunayBig;
if(xbig > ybig) big = xbig;
else big = ybig;
eps = (dx + dy) * ibggDelaunayEps;
for(i=nl+1;i<=ibgg.lastPoint;i++){
cnx(i)[0]=xm; cnx(i)[1]=ym;
cnt(i)=ibgSRegion; cnu(i)=outside;
}
cnx(nl+1)[0] -= big; cnx(nl+2)[0] += big;
cnx(nl+3)[1] -= big; cnx(nl+4)[1] += big;
cvert.x1[0]=xm+xbig; cvert.x1[1]=ym;
cvert.x2[0]=xm-xbig; cvert.x2[1]=ym;
for(i=1;i<3;i++){
cct(i)=ibgc2Triangle; ccu(i)=0;
ccv(i)=cvert;
ibgg.CPointFirst[i] = 3*i-2;
ibgg.CSideFirst[i] = 3*i-2;
}
ibgg.freeCPoint=ibgg.freeCSide=7; ibgg.lastCell=2;
for(i=1;i<7;i++){
ibgg.CPoint[i] = nl+hcn2[i-1];
ibgg.CSide[i] = hcs2[i-1];
}
if((dslast>=DSMAX)||(dslast<=1)) dslast=DSMAX-1;
/* Zyklus "uber die Knoten: */
cnc(0)=1; nl=ibgg.lastPoint-4;
for(n0=1;n0<=nl;n0++){
/* Suche nach einer Zelle, die n0 enth"alt: */
n0zyk: if(list->Near){ /* Der nahegelegene Knoten */
if((no=list->Near[n0])==0) no = n0-1;
}else {no=n0-1;}
c0=ibgg2find(cnx(n0),cnc(no));
x=cnx(n0)[0]; y=cnx(n0)[1];
/* Anlegen folgender Listen:
dsc[dsl]: Zellen, die durch den Knoten n0 zerst"ort werden;
nsc[nsl]: neue Zellen;
*/
nsl=0; dsl=1; dsi=1; ccu(c0)= -dsl; dsc[dsl++]=c0;
while(dsi<dsl){
c0 = dsc[dsi];
cs0 = ccs(c0); ms=ibgcSides((t0=cct(c0)));
cn0 = ccn(c0);
for(s=0;s<ms;s++){
cc= *(cs0++);
if(cc==0) {del=0; goto deltest;}
if(ccu(cc)<0) {dscl[dsi][s]= -1; continue;}
v = ccv(cc); t = cct(cc);
if((vv= (v.x1[0]-x)*(v.x2[0]-x)
+ (v.x1[1]-y)*(v.x2[1]-y)) < -eps2) del=1;
else del=0;
if(dsl>=dslast) del=0;
deltest: if(del){
/* n0 in Umkreis -> Zelle cc wird zerst"ort */
ccu(cc) = -dsl; dsc[dsl++] = cc;
dscl[dsi][s] = -1;}
else{
/* Anstelle von c0 entsteht neue Nachbarzelle c1 von cc */
g0CGet(&c1,ibgc2Triangle);
if(c1<0) {ibgfatal; return ibgError;}
cni=ccn(c1); csi=ccs(c1);
cn0 = ccn(c0);
cni[0]=n0; cnc(n0)=c1; cni++; cnj=cni;
csi[0]=cc; csi++;
for(i=0;i<2;i++){
cni[i] = nt = cn0[ibgcSPoint1(t0,s,i)];
cnc(nt) = c1; csi[i] = -1;
}
/* Berechnen und Eintragen der Daten des Umkreises: */
cx1=cnx(*cnj++); cx2=cnx(*cnj++);
x1=cx1[0]-x; y1=cx1[1]-y;
x2=cx2[0]-x; y2=cx2[1]-y;
l1=x1*x1 + y1*y1;
l2=x2*x2 + y2*y2;
det=x1*y2 - x2*y1;
if((det<eps2)){ /* Entartung */
ibgassert(determ);
g0CDel(c1);
/* cc has to be deleted, but it's neighbours don't know this if they are
deleted before (v-test is ok for cc). test deleted neighbours for this effect*/
csh=ccs(cc);
for(sh=0;sh<ibgcSides(t);sh++){
if(!(dsh= -ccu(ch=csh[sh]))) continue;
if(ch==c0) continue;
if(dsh>dsi) continue;
for(ih=0;ih<ibgcSides(cct(ch));ih++){
if(ccs(ch)[ih]==cc){
cl=(dscl[ih])[dsh];
ibgassert(cl>0);
g0CDel(cl);
(dscl[ih])[dsh] = -1;
}
}
}
del=1; goto deltest;
}
mx=(l1*y2 - l2*y1)/det;
my=(l2*x1 - l1*x2)/det;
ccv(c1).x1[0]=x;
ccv(c1).x1[1]=y;
ccv(c1).x2[0]=x+mx;
ccv(c1).x2[1]=y+my;
cvert = ccv(c1);
/* Eintragen von c1 in die nsc- und dsc-Listen: */
nsc[nsl]=c1; nscs[nsl]=s; nsco[nsl]=c0;
nsl++; if(nsl>=NSMAX) ibgfatal;
dscl[dsi][s]=c1;
}
}
dsi++;
}
/* Herstellung der Verbindungen der neuen Zellen untereinander: */
for(nsi=0;nsi<nsl;nsi++){ /* Zyklus "uber neue Zellen */
nsizyk: c1 = nsc[nsi]; cs1= ccs(c1); /* neue Zelle */
c0 =nsco[nsi]; s0 = nscs[nsi]; /* alte Zelle */
ds0= -ccu(c0);
t0 =cct(c0); t1=cct(c1); cn0=ccn(c0); cs0=ccs(c0);
/* Eintragen des Zeigers von cc auf c1: */
if((cc=cs1[0])>0){
se=(ss=ccs(cc))+ibgcSides(t);
for(;ss<se;ss++) if(ss[0]==c0)
{ss[0]=c1; break;}
}
ccu(c1)=0;
for(i0=0;i0<2;i0++){
if(cs1[ibgcSSide(t1,0,i0)]>0) continue;
nh=cn0[ibgcSPoint1(t0,s0,i0)];
ct=dscl[ds0][ibgcSSide(t0,s0,i0)];
if(ct<=0){
/* "Ubergang zur Nachbarzelle */
ch=c0;th=t0;sh=s0;ih=i0;csh=cs0;
do{ co =ch;
ch =csh[ibgcSSide(th,sh,ih)];
dsh= -ccu(ch); ibgassert(dsh>0);
th =cct(ch);
cnh=ccn(ch);
csh=ccs(ch);
for(sh=0;;sh++) {
if(csh[sh]==co) break;
ibgassert(sh<4);}
for(ih=0;;ih++) {
if(cnh[ibgcSPoint1(th,sh,ih)]==nh)
break;
ibgassert(ih<2);}
ct=dscl[dsh][ibgcSSide(th,sh,ih)];
}while(ct<=0);
}
cs1[ibgcSSide(t1,0,i0)]=ct;
cnj = ccn(ct); csj = ccs(ct); tt=cct(ct);
for(j=0;j<4;j++) if(cnj[ibgcSPoint1(tt,0,j)]==nh){
csj[ibgcSSide(tt,0,j)]= c1; break;
}
}
}
/* Eintragen der zerst"orten Zellen in die free-Listen */
for(dsi=1;dsi<dsl;dsi++){
cc=dsc[dsi]; g0CDel(cc);
}
/* Versuche, Nachbarzellen zusammenzukleben:
for(nsi=1;nsi<nsl;nsi++){
cc=nsc[nsi];
switch(t=cct(cc)){
case ibgc2Triangle:
for(i=0;i<4;i++){
cn=
} */
}
/* Wegschmei"sen der alten Felder: */
ibgassert(ibgg.gridDimension==2);
*g = ibgg;
return ibgSuccess;
}
/*
static int ibgg1New(ibGrid0 *g0)
{ibGrid *gg;
ibgPoint nd; ibgSegmentType t; int nn,bb,c,d,na,ne;
ibgFloat xm,dx,big;
int n,c0,*cn0,*cs0,nc;
ibgg.spaceDimension = list->spaceDimension;
ibgg.firstPoint = list->firstPoint;
ibgg.lastPoint = list->lastPoint;
ibgg.maxPoint = list->maxPoint;
ibgg.Point = list->Point;
g0CAlloc(ibgg.lastPoint+10);
ibgg.lastCell=ibgg.lastPoint+1;
c0=0; cn0=ibgg.CPoint; cs0=ibgg.CSide; nc=0; na=1;
for(n=1;n<=ibgg.lastPoint;n++){
c0++;
*(cn0++) = n;
*(cs0++) = cndr(n,1);
if((*(cn0++) = *(cs0++) = cndr(n,0))==outside) ne=n;
ibgg.CPointFirst[c0] = ibgg.CSideFirst[c0] = nc;
nc += 2;
cct(c0) = ibgc1Edge;
ccu(c0) = 0;
}
ibgassert(cndn(na,0)==outside);
ibgassert(cndp(ne,0)==outside);
ibgg.lastPoint += 2; if(g0NRealloc()==0) return ibgError;
xm=0.5*(xmin(0)+xmax(0)); dx=(xmax(0)-xmin(0));
big=0.01*dx*dx/ibGeo->Delta;
cnx(ibgg.lastPoint-1)[0] = xm-big;
cnx(ibgg.lastPoint )[0] = xm+big;
cnu(ibgg.lastPoint-1) = nothing;
cnu(ibgg.lastPoint ) = nothing;
cnt(ibgg.lastPoint-1) = ibgSRegion;
cnt(ibgg.lastPoint ) = ibgSRegion;
c0++;
ibgg.CPointFirst[c0] = ibgg.CSideFirst[c0] = nc;
nc += 2;
cct(c0) = ibgc1Edge;
ccn(c0)[1] = na;
ccs(c0)[1] = na; ccs(na)[0] = c0;
ccn(c0)[0] = ibgg.lastPoint-1; ccs(c0)[0] = 0;
ccn(ne)[1] = ibgg.lastPoint; ccs(ne)[1] = 0;
ibgassert(nc<ibgg.maxCPoint);
ibgassert(nc<ibgg.maxCSide);
ibgassert(c0<ibgg.maxCell);
ibgg.freeCPoint = ibgg.freeCSide = nc;
ibgg.lastCell = c0;
ibgassert(c0==ibgg.lastPoint-1);
*g = ibgg;
return ibgSuccess;
}
*/
int ibggDelaunay0(ibGrid *g, ibgPoints *list)
{int n;
int *n0 = list->Near;
ibgg.spaceDimension = list->spaceDimension;
ibgg.firstPoint = list->firstPoint;
ibgg.lastPoint = list->lastPoint;
ibgg.maxPoint = list->maxPoint;
ibgg.Point = list->Point;
g0CAlloc(ibgg.lastPoint);
ibgg.lastCell = ibgg.lastPoint;
for(n=ibgg.firstPoint;n<=ibgg.lastPoint;n++){
ibgg.CPointFirst[n] = n;
ibgg.CSideFirst[n] = 0;
ibgg.CPoint[n] = n;
cct(n) = ibgc0Node;
ccu(n) = 0;
}
*g = ibgg;
return ibgSuccess;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -