ibginside.cxx

Delaunay三角形的网格剖分程序

#include "ibggenerator.hxx"
#include "ibgoctree.hxx"
#include "ibg2.hxx"
#include <stdlib.h>
#include <time.h>
double tmrCPUTime(void);   /* seconds */
double tmrRealTime(void);  /* seconds */
char * tmrDate(void);

extern int ibgNumberOfFaceCells;
void ibgDistinguishOutsideInside(wzgrid g, wzRegion out, wzRegion ins, wzRegion undef);

/*

The following algorithm assumes a Delaunay grid created by a
sufficiently fine set of surface nodes for a nonconvex surface. The
problem is that the Delaunay grid defines the convex hull of this
non-convex surface. 

The idea is that, if we have an ~ homogeneous surface node set,
surface triangles may be distinguished by the property that they have
a small diameter.  Starting from outside, we fill the region inside 
through boundary triangles which are too large.

The algorithm depends on the correct choice for epsilon. Thus,
different epsilons have to be tried until the whole domain has been
filled. The last try before seems the best.

*/

static int OnlyEpsilonMode=0;

void ibgDistinguishOutsideInside(wzgrid g, wzFloat eps0, 
				 wzRegion out, wzRegion ins, wzRegion undef)
{
  wzArray<int> stack(g->cells);
  int i,s,m,m0,m1,c,c0,cs,*ccn,*ccs,*csn,stb,stl,st1=0,st0=0,l,n1,n2;
  wzRegion r0,r1;
  wzFloat edge[10],maxedge,dd,maxsedge;
  wzRangeLoop(g->cells,c){
    if(g->ccund(c)) continue;
    g->ccu(c) = 0;
  }
  wzRangeLoop(g->cells,c0){
    if(g->ccund(c)) continue;
    if(g->ccu(c0) != 0) continue;
    stl = 1; stb = 0; stack[stl] = c0;
    g->ccu(c0) = c0;
    while(stb < stl){
      stb++;
      c = stack[stb];
      ccs = g->ccs(c); 
      ccn = g->ccn(c);
      for(l=0;l<wzCellTypesides[wzCellType3Tetrahedron];l++){
	edge[l] = 0;
      }
      maxedge = 0;
      for(l=0;l<wzCellTypelines[wzCellType3Tetrahedron];l++){
	n1 = ccn[wzCellTypeepoint1[wzCellType3Tetrahedron][l]];
	n2 = ccn[wzCellTypeepoint2[wzCellType3Tetrahedron][l]];
	dd = cogLine(g->Point[n1],g->Point[n2]).diameter();
	if(dd>edge[wzCellTypeesidel[wzCellType3Tetrahedron][l]]){
	  edge[wzCellTypeesidel[wzCellType3Tetrahedron][l]] = dd;
	}
	if(dd>edge[wzCellTypeesider[wzCellType3Tetrahedron][l]]){
	  edge[wzCellTypeesider[wzCellType3Tetrahedron][l]] = dd;
	}
	if(dd>maxedge) maxedge = dd;
      }
      for(l=0;l<wzCellTypesides[wzCellType3Tetrahedron];l++){
	cs = ccs[l];
	if(cs<=0) continue;
	if(g->ccu(cs)) continue;
	if(edge[l]>eps0) goto addside;
	if(OnlyEpsilonMode) continue;
	csn = g->ccn(cs); maxsedge=0;
	for(s=0;s<wzCellTypelines[wzCellType3Tetrahedron];s++){
	  n1 = csn[wzCellTypeepoint1[wzCellType3Tetrahedron][s]];
	  n2 = csn[wzCellTypeepoint2[wzCellType3Tetrahedron][s]];
	  dd = cogLine(g->Point[n1],g->Point[n2]).diameter();
	  if(dd>maxsedge) maxsedge=dd;
	}
	if(maxsedge<eps0){
	  if(maxsedge>edge[l] && edge[l] >= maxedge) goto addside;
	  if(maxedge>edge[l] && edge[l] > 0.9999*maxsedge) goto addside;
	}
	continue;
      addside:
	stl++; stack[stl] = cs; g->ccu(cs) = c0;
      }
    }
    if(stl>st1){
      if(stl>st0){
	st1 = st0; m1 = m0; st0 = stl; m0 = c0;
      }else{
	st1 = stl; m1 = c0;
      }
    }
  }
  m = 0;
  wzRangeLoop(g->cells,c){
    if(g->ccund(c)) continue;
    ccn = g->ccn(c);
    for(i=0;i<4;i++) if(ccn[i]<7){ m = g->ccu(c); goto found;}
  }
found:
  if(m==m0){
    r0 = out; r1 = ins;
  }else if (m==m1){
    r1 = out; r0 = ins;
  }else{
    wzAssert(0);
  }
  wzRangeLoop(g->cells,c){
    if(g->ccund(c)) continue;
    m = g->ccu(c);
    if(m==m0){
      g->ccu(c) = r0;
    }else if(m==m1){
      g->ccu(c) = r1;
    }else{
      g->ccu(c) = undef;
    }
  }
  g->lastRegion = g->cells.last();
}

#include <math.h>

void ibgDistinguishOutsideInside(wzgrid g, wzRegion out, wzRegion ins, wzRegion undef)
{
  //  wzFloat epsilon[20];
  wzFloat eps0 = g->cgmax(0)-g->cgmin(0);
  wzFloat epsf = 0.5;
  wzIndex i,c,m,mu,mi,mo,mi0=0;
  OnlyEpsilonMode = 1;
  wzFloat epsbest = eps0;
  for(i=0;i<25;i++){
    ibgDistinguishOutsideInside(g,eps0,out,ins,undef);
    mo=mi=mu=0;
    if(i>10)  OnlyEpsilonMode = 0;
    wzRangeLoop(g->cells,c){
      if(g->ccund(c)) continue;
      m = g->ccu(c);
      if(m==ins) mi++;
      if(m==undef) mu++;
      if(m==out) mo++;
    }
    if(mi>=mi0){
      epsbest = eps0;
      eps0 *= epsf; mi0 = mi; 
    }else{
      eps0 /= epsf; epsf = sqrt(epsf); 
      eps0 *= epsf;
    }
  }
  ibgDistinguishOutsideInside(g,epsbest,out,ins,out);
}