#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 asufficiently fine set of surface nodes for a nonconvex surface. Theproblem is that the Delaunay grid defines the convex hull of thisnon-convex surface. The idea is that, if we have an ~ homogeneous surface node set,surface triangles may be distinguished by the property that they havea 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 beenfilled. 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);}