wzdelaunay.cxx

Delaunay三角形的网格剖分程序

#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "ibg2.hxx"

void wzCellTypeInit(); 

ibgDelaunayMarker::ibgDelaunayMarker()
  :dsl((1+IBG2CSMAX)*sizeof(int))
  ,dsc(dsl.base)
  ,nsl(3*sizeof(int))
  ,nsc(nsl.base)
  ,nsco(nsc,1)
  ,nscs(nsc,2)
{
  for(int i=0;i<IBG2CSMAX;i++){
    dscl[i].reinitialize(dsc,1+i);
  }
}
void wzDelaunayGenerator::Delaunay(wzpoints plist, cogeometry cog)
{
  wzCellTypeInit();
  initializeDelaunayBox(plist->min,plist->max);
  include(plist);
  try{
    defineRegions(cog);
    createFaces(cog);
  }catch(...){;}
}

void wzDelaunayGenerator::Delaunay(wzpoints plist)
{
  wzCellTypeInit();
  initializeDelaunayBox(plist->min,plist->max);
  include(plist);
}

void wzDelaunayGenerator::include(wzpoints plist)
{
  switch(gdim()){
  case 0: Delaunay0(plist); break;
  case 1: Delaunay1(plist); break;
  case 2: Delaunay2(plist); break;
  case 3: Delaunay3(plist); break;
  }
}

#define ibg2pX(point)		((point).X)
#define ibg2pF(point)		((point).X+wzPointDim)
#define ibg2pI(point)		((point).I)
#define ibg2pSegment(point)	((point).S.index())
#define ibg2pType(point)	((point).S.codim())
#define ibg2pSetSegment(point,tt,uu)	((point).S = wzSegment((wzIndex)(uu),(wzDimension)(tt)))

const wzFloat	ibgDelaunayBig  = 10000.0;
const wzFloat	ibgDelaunayEps  = 2.e-6;
static wzFloat	ibgDelaunayEps2;
static wzFloat	ibgDelaunayEps3;

static int hcn3[16] = {1,3,2,5, 1,5,2,4, 1,4,2,6, 1,6,2,3};
static int hcs3[16] = {0,2,0,4, 0,3,0,1, 0,4,0,2, 0,1,0,3};
static int hcn2[ 6] = {1,3,2,  1,2,4};
static int hcs2[ 6] = {0,2,0,  0,0,1};
static int hcn1[ 2] = {1,2};
static int hcs1[ 2] = {0,0};
static int hcnum[4] = {0,1,2,4};
static int* hcn[ 4] = {0,hcn1,hcn2,hcn3};
static int* hcs[ 4] = {0,hcs1,hcs2,hcs3};
static wzCellType hctyp[4] = {wzCellType0Node,wzCellType1Line,wzCellType2Triangle,wzCellType3Tetrahedron};

void wzDelaunayGenerator::initializeDelaunayBox(wzFloat* min, wzFloat* max)
{
  int i,j,c,*n,*s,*hn,*hs;
  wzFloat m[wzPointDim],d0,dd;
  ibgSphere cvert;
  for(i=0;i<xdim();i++){
    Grid.mincoord[i+1] = min[i]; 
    Grid.maxcoord[i+1] = max[i];
  }
  if(gdim()==0) return;
  Grid.points.create(2*gdim());
  for(i=1;i<=hcnum[gdim()];i++){
    j = createCell(hctyp[gdim()]);
  }
  d0 = 0;
  for(i=0;i<xdim();i++){
    m[i] = (min[i]+max[i])/2; dd = max[i]-min[i]; if(dd>d0) d0=dd;
  }
  //  ibgDelaunayEps1 = d0*ibgDelaunayEps;
  ibgDelaunayEps2 = d0*ibgDelaunayEps*(max[0]-min[0])*ibgDelaunayEps;
  ibgDelaunayEps3 = 
    (max[0]-min[0])*ibgDelaunayEps*
    (max[1]-min[1])*ibgDelaunayEps*
    (max[2]-min[2])*ibgDelaunayEps;
  d0 *= ibgDelaunayBig;
  for(i=0;i<gdim();i++){
    for(j=0;j<xdim();j++){
      Grid.Point[2*i+1][j] = Grid.Point[2*i+2][j] = cvert.x1[j] = cvert.x2[j] = m[j];
    }
    Grid.Point[2*i+1][i] -= d0;
    Grid.Point[2*i+2][i] += d0;
    cnc(2*i+1) = cnc(2*i+2) = 1;
    cns(2*i+1,wzIsRegion,0);
    cns(2*i+2,wzIsRegion,0);
  }
  cvert.x1[0] -= d0;
  cvert.x2[0] += d0;
  hn = hcn[gdim()]; hs = hcs[gdim()];
  wzRangeLoop(Grid.cells,c){
    ccv(c)=cvert;
    n = ccn(c); s = ccs(c);
    for(i=0;i<wzCellTypepoints[hctyp[gdim()]];i++){*n++ = *hn++;}
    for(i=0;i<wzCellTypesides [hctyp[gdim()]];i++){*s++ = *hs++;}
  }
}

static int LastUsedCell;

int wzDelaunayGenerator::g1find(wzFloat *xx, int c0)
{
  int *cn0;
  if(ccu(c0)<0){
    c0 = 1;
    while (ccu(c0)<0) c0++;
  }
beg:
  wzAssert(ccu(c0)>=0);
  cn0=ccn(c0);
  if(xx[0]<cnx(cn0[0])[0]) {c0=ccs(c0)[0]; goto beg;} 
  if(xx[0]>cnx(cn0[1])[0]) {c0=ccs(c0)[1]; goto beg;} 
  return LastUsedCell = c0;
}

int wzDelaunayGenerator::g2find(wzFloat *xx, int c0)
{
  int nl,*cn0;
  wzCellType t;
  int i;
  wzFloat vv,x1,x2,y1,y2;
  wzFloat* ff[IBG2CSMAX];
  register wzFloat *fj;
  if(ccund(c0)) c0 = LastUsedCell;
  if(ccund(c0)){
    wzRangeLoop(Grid.cells,c0){
      if(ccund(c0)) continue;
      break;
    }
  }
beg:
  wzAssert(!ccund(c0));
  wzAssert(ccu(c0)>=0);
  nl=wzCellTypePoints(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[wzCellTypeSPoint1(t,i,0)];
    x1 = fj[0]-xx[0]; y1=fj[1]-xx[1];
    fj = ff[wzCellTypeSPoint1(t,i,1)];
    x2 = fj[0]-xx[0]; y2=fj[1]-xx[1];
    vv = x1 * y2 - x2 * y1;
    if(vv<-ibgDelaunayEps2) {c0=ccs(c0)[i]; goto beg;}
  }
  return LastUsedCell = c0;
}

int wzDelaunayGenerator::g3find(wzFloat *xx, int c0)
{
  register int nl,*cn0;
  int i,sl;
  wzCellType t;
  long double vv,x1,x2,x3,y1,y2,y3,z1,z2,z3,xx0=xx[0],xx1=xx[1],xx2=xx[2];
  wzFloat* ff[IBG2CNMAX];
  register wzFloat *fj;
  if(ccund(c0)) c0 = LastUsedCell;
  if(ccund(c0)){
    wzRangeLoop(Grid.cells,c0){
      if(ccund(c0)) continue;
      break;
    }
  }
beg:
  wzAssert(!ccund(c0));
  wzAssert(ccu(c0)>=0);
  nl=wzCellTypePoints(t=cct(c0)); cn0=ccn(c0);
  for(i=0;i<nl;i++) ff[i]=cnx(cn0[i]);
  sl=wzCellTypeSides(t);
  for(i=0;i<sl;i++){
    fj = ff[wzCellTypeSPoint1(t,i,0)];
    x1 = ((long double)fj[0])-xx0;y1=((long double)fj[1])-xx1;z1=((long double)fj[2])-xx2;
    fj = ff[wzCellTypeSPoint1(t,i,1)];
    x2 = ((long double)fj[0])-xx0;y2=((long double)fj[1])-xx1;z2=((long double)fj[2])-xx2;
    fj = ff[wzCellTypeSPoint1(t,i,2)];
    x3 = ((long double)fj[0])-xx0;y3=((long double)fj[1])-xx1;z3=((long double)fj[2])-xx2;
    vv = x1 * (y2*z3-y3*z2) + x2 * (y3*z1-y1*z3) + x3 * (y1*z2-y2*z1);
    if(vv>ibgDelaunayEps3) {c0=ccs(c0)[i]; goto beg;}
  }
  return LastUsedCell = c0;
}

int   wzDelaunayGenerator::Delaunay0(wzpoints list)
{
  if(int l = list->last()){
    int c = createCell(wzCellType0Node);
    int n = Grid.points.append();
    ccn(c)[0]=n; cnc(n) = c;
    wzPoint& p = Grid.Point[ccn(c)[0]];
    p = list->Point[list->first()];
  }
  return wzSuccess;
}

int   wzDelaunayGenerator::Delaunay1(wzpoints list)
{
  int n,n0,c0,c1,cs,no,nb,nshift;
  n0 = nb = Grid.points.append(list->length());
  nshift = nb - list->first(); n0--;
  wzRangeLoop(*list,n){
    n0++;
    Grid.Point[n0]=list->Point[n];
    // Suche nach einer Zelle, die n0 enth"alt: 
    no = list->Near[n] + nshift;
    if(no<nb) no = n0-1;
    c0=g1find(cnx(n0),cnc(no));
    c1 = createCell(wzCellType1Line);
    ccs(c1)[1] = cs = ccs(c0)[1];
    ccs(c0)[1] = ccs(cs)[0] = c1;
    ccs(c1)[0] = c0;  
    ccn(c1)[1] = ccn(c0)[1];
    ccn(c1)[0] = ccn(c0)[1] = n0;
    cnc(n0)    = c1;
    ccu(c1)    = ccu(c0)    = 0;
  }
  return wzSuccess;
}

int   wzDelaunayGenerator::Delaunay2(wzpoints list)
{
  int n,s,s0,sh,ms,i,i0,ih,j,del;
  int n0,no,nb,nt,nshift,nh,c0,c1,cc,ct,ch,co,cl;
  int *cs0,*cs1,*csi,*csj,*csh,*cn0,*cni,*cnj,*cnh,*ss,*se;
  wzFloat *cx1,*cx2;
  wzFloat x,y,x1,x2,y1,y2;
  wzFloat l1,l2,det,mx,my;
  wzCellType t,t0,t1,th,tt;
  wzPoint nd;
  ibgSphere v;
  int dsi,ds0,dsh,nsi;
  LastUsedCell=0;
  n0 = nb = Grid.points.append(list->length());
  nshift = nb - list->first(); n0--;
  wzRangeLoop(*list,n){
    n0++;
    Grid.Point[n0]=list->Point[n];
    // Suche nach einer Zelle, die n0 enth"alt: 
    no = list->Near[n] + nshift;
    if(no<nb) no = n0-1;
    c0=g2find(cnx(n0),cnc(no));
    x=cnx(n0)[0]; y=cnx(n0)[1];
    nsl=0; dsl=0; dsi=1; dsc[++dsl]=c0; ccu(c0)= -dsl;
    while(dsi<=dsl){
      c0  = dsc[dsi];
      cs0 = ccs(c0); ms=wzCellTypeSides((t0=cct(c0)));
      cn0 = ccn(c0);
      for(s=0;s<ms;s++){
	cc= cs0[s];
	if(cc<=0) {del=0; goto deltest;}
	if(ccu(cc)<0) {(dscl[s])[dsi]= -1; continue;}
	v = ccv(cc); t = cct(cc);
	if((    (v.x1[0]-x)*(v.x2[0]-x)
		+   (v.x1[1]-y)*(v.x2[1]-y)) < ibgDelaunayEps2)	del=1;
	else							del=0;
      deltest:
	if(del){
	  dsc[++dsl] = cc; ccu(cc) = -dsl;