ibgomark.cxx

Delaunay三角形的网格剖分程序

#include "cogenoctree.hxx"
#include "ibgoctree.hxx"
#include "cogwarnings.hxx"

const wzIndex CogenOctree::UnderdogRemoval = 0x01;
const wzIndex CogenOctree::LongSideIntersectionRemoval = 0x02;
const wzIndex CogenOctree::IrregularRemoval = 0x04;
const wzIndex CogenOctree::FullRemoval = 0xff;

const wzIndex ibgoUnderFace = 0x2;
const wzIndex ibgoOnNode    = 0x2;
const wzIndex ibgoIrregular = 0x4;
const wzIndex ibgoSupports  = 0x8;
const wzIndex ibgoUnhide    = ~ibgoHidden;
const wzIndex ibgoDontHide  = ibgoSupports | ibgoOnNode;

void CogenOctree::removeLongSideIntersections(wzBoolean on)
{
  if(on){
    RegionNodeRemoval |= LongSideIntersectionRemoval;
    FaceNodeRemoval   |= LongSideIntersectionRemoval;
  }else{
    RegionNodeRemoval &= ~LongSideIntersectionRemoval;
    FaceNodeRemoval   &= ~LongSideIntersectionRemoval;
  }
}

void CogenOctree::removeIrregularNodes(wzBoolean on)
{
  if(on){
    RegionNodeRemoval |= IrregularRemoval;
  }else{
    RegionNodeRemoval &= ~IrregularRemoval;
  }
}

void CogenOctree::removeUnderdogs(wzFloat fac)
{
  if(fac){
    RegionNodeRemoval |= UnderdogRemoval;
    FaceNodeRemoval   |= UnderdogRemoval;
    EdgeNodeRemoval   |= UnderdogRemoval;
    VertexNodeRemoval |= UnderdogRemoval;
  }else{
    RegionNodeRemoval &= ~UnderdogRemoval;
    FaceNodeRemoval   &= ~UnderdogRemoval;
    EdgeNodeRemoval   &= ~UnderdogRemoval;
    VertexNodeRemoval &= ~UnderdogRemoval;
  }
}


wzFloat ibgOctree::LongSideIntersectionFactor = 0.45;
// wzFloat ibgOctree::LongSideIntersectionFactor = 0.25;


void CogenOctree::hideForbiddenFaceDirections(ibgoctree gg) const
{
  if(!suppressFaceIntersection[3]) return;
  ibgOctree& g = *gg;
  wzIndex ff,ll,d;
  wzRangeLoop(g.faces,ff){
    ll = g.fLine[ff];	 	
    d  = g.lDirection(ll);
    if(suppressFaceIntersection[d]){
      g.fType[ ff] = ibgoHidden;
    }
  }
}


void CogenOctree::hideLongSideIntersections(ibgoctree gg) const
{
  wzBoolean omitFace;
  if(suppressFaceIntersection[3]) return;
  if(!(RegionNodeRemoval & LongSideIntersectionRemoval)) return;
  if(FaceNodeRemoval & LongSideIntersectionRemoval) omitFace = 1; 
  ibgOctree& g = *gg;
  wzInteger nu,no,nru,nro,noo,nuu;
  wzIndex d,u,s,r,dr;
  wzIndex ff,fs,ll,ls,bu,bo; //bru,bro,bou,boo;
  wzFloat dx,dy,dx2,dy2;
  wzBoolean isIrregular,isShort,Long1,Long2,nosupp[]={1,1,1,1};
  wzRangeLoop(g.faces,ff){
    if(g.fType[ff]&ibgoDontHide) continue;
    ll = g.fLine[ff];	 	
    d  = g.lDirection(ll);	u  = g.down(d);
    bo = g.lNodeUp[ll];		no = g.bnBase[bo];
    bu = g.lNodeDown[ll];	nu = g.bnBase[bu];
    dx = g.LongSideIntersectionFactor*(g.nPoint[no][d] - g.nPoint[nu][d]);
    dx2 = g.chart->g_ii(g.fMetric[ff],g.fPoint[ff],d)*dx*dx;
    isShort=1;isIrregular=0;
    for(s=0;s<g.LineSides;s++){
      r = g.lsLine[d][s]; dr = g.Up(r);
      if(g.is_defined(nru=g.nNext[r][nu]))
	if(g.is_defined(nro=g.nNext[r][no])){
	  dy = g.nPoint[nru][dr] - g.nPoint[nu][dr];
	  dy2 = g.chart->g_ii(g.fMetric[ff],g.fPoint[ff],dr)*dy*dy;
	  if(dy2<dx2){
	    if(isShort){
	      isShort=0;Long1=dr;Long2=3;nosupp[0]=nosupp[1]=nosupp[2]=1;
	    }else{
	      if(Long1 != dr) Long2=dr;
	    }
	    if((ls=g.bnNext[r][bu])&&(fs=g.lFaceUp[ls])){
	      (g.fType[fs] |= ibgoSupports) &= ibgoUnhide; nosupp[dr]=0;
	    }
	    if((ls=g.bnNext[r][bo])&&(fs=g.lFaceUp[ls])){
	      (g.fType[fs] |= ibgoSupports) &= ibgoUnhide; nosupp[dr]=0;
	    }
	  }
	}else isIrregular=1;
      else isIrregular=1;
    }
    if(isIrregular) continue;
    if(isShort) continue;
    if(nosupp[Long1] && nosupp[Long2]) continue;
    g.bnType[bu] = ibgoHidden;
    g.bnType[bo] = ibgoHidden;
    if(omitFace) g.fType[ ff] = ibgoHidden;
  }
}

//	void CogenOctree::hideIrregular(ibgoctree gg, cogeometry g) const;
// Hides nodes in boundary environment which are not regular

void CogenOctree::hideIrregular(ibgoctree gg) const
{
  ibgOctree& g = *gg;
  wzInteger nu,no,nru,nro,noo,nou;
  wzIndex d,dr,u,s,t,o,r;
  wzIndex ff,ll,bu,bo,bru,bro,bou,boo;
  wzIndex rHide = RegionNodeRemoval & IrregularRemoval;
  // even if we should not remove them, they should be marked;
  wzRangeLoop(g.faces,ff){
    ll = g.fLine[ff];	 	
    d  = g.lDirection(ll);	u  = g.down(d);
    bo = g.lNodeUp[ll];		no = g.bnBase[bo];
    bu = g.lNodeDown[ll];	nu = g.bnBase[bu];
    for(s=0;s<g.LineSides;s++){
      r = g.lsLine[d][s]; dr = g.Up(r);
      if(g.is_nothing(nru=g.nNext[r][nu])){
	if(rHide) g.bnType[bu] = ibgoHidden;
	g.bnType[bu] |= ibgoIrregular;
      }else if(g.is_defined(nru)){
	bru = g.nBoundary[nru]; 
	if(!bru) bru = g.createSimpleBoundaryNode(nru);
	if(g.is_nothing(g.nNext[d][nru])){
	  if(rHide) g.bnType[bru] = ibgoHidden;
	  g.bnType[bru] |= ibgoIrregular;
	}
	o = g.llOrtho[d][r];
	for(t=0;t<g.PlaneSides;t++){
	  if(g.is_nothing(nou=g.nNext[o][nru])){
	    if(rHide) g.bnType[bru] = ibgoHidden;
	    g.bnType[bru] |= ibgoIrregular;
	  }else if(g.is_defined(nou)){
	    if(g.is_nothing(g.nNext[d][nou])){
	      bou = g.nBoundary[nou];
	      if(!bou) bou = g.createSimpleBoundaryNode(nou);
	      if(rHide) g.bnType[bou] = ibgoHidden;
	      g.bnType[bou] |= ibgoIrregular;
	    }
	  }
	  o = g.Other(o);
	}
      }
      if(g.is_nothing(nro=g.nNext[r][no])){
	if(rHide) g.bnType[bo] = ibgoHidden;
	g.bnType[bo] |= ibgoIrregular;
      }else if(g.is_defined(nro)){
	bro = g.nBoundary[nro];
	if(!bro) bro = g.createSimpleBoundaryNode(nro);
	if(g.is_nothing(g.nNext[u][nro])){
	  if(rHide) g.bnType[bro] = ibgoHidden;
	  g.bnType[bro] |= ibgoIrregular;
	}
	o = g.llOrtho[d][r];
	for(t=0;t<g.PlaneSides;t++){
	  if(g.is_nothing(noo=g.nNext[o][nro])){
	    if(rHide) g.bnType[bro] = ibgoHidden;
	    g.bnType[bro] |= ibgoIrregular;
	  }else if(g.is_defined(noo)){
	    boo = g.nBoundary[noo];
	    if(g.is_nothing(g.nNext[u][noo])){
	      if(!boo) boo = g.createSimpleBoundaryNode(noo);
	      if(rHide) g.bnType[boo] = ibgoHidden;
	      g.bnType[boo] |= ibgoIrregular;
	    }
	  }
	  o = g.Other(o);
	}
      }
    }
  }  
}

const int hideNothing = 0;
const int hideRefined = 1;
const int hideFew = 2;
//const int hideNearlyAll = 3;
const int hideAll = 4;

//static int	ibgBasicMode = hideRefined;
static int	ibgBasicMode = hideFew;
static float	ibgOnLongSideButClose = 0.1;
static float	ibgFaceVeryClose = 0.01;
static float	ibgFaceRelativeClose = 0.1;
static float	ibgFaceRelativeFar   = 0.9;
static float	ibgOtherFaceTooClose = 0.1;
static float	ibgEdgeTooClose = 0.3;
static float	ibgEdgeTooClosePar = 0.5;

void CogenOctreeOldmark(ibgoctree gg, cogeometry cog)
{
  ibgInteger bn,n,i,rc,ri,rb,rx,fl,ff,fr,d,dr,r,ls,ll,e,s,mark,dim,border;
  ibgInteger nu,no,nt,nr,bu,bo,ou,pe;
  const ibgIndex Default=0,FarFaceType=1,FaceType=2,CloseFaceType=3,OnFace=4,RefinedNode=5; 
  ibgFloat dx,dy,df,de,dp,de1,dp1,rt,rf,xf,yf,yr,xi,di,fi,db,fb,xu,xo;
  cogInfoOmittedFaces.report_only();
  cogInfoOmittedNodes.report_only();
  wzRangeLoop(gg->bnodes,bn){
    n = gg->bnBase[bn];
    rc = ri = rb = rx = Default; dim=border=0;
    for(i=0;i<gg->LinesUp;i++){
      if(gg->is_outside(no=gg->nUp  [i][n])){
	border++;
	if((fl=gg->bnDown[i][bn]) && (ff=gg->lFaceUp[fl])){
	  db = (xi=gg->nPoint[n][i])-gg->fPoint[ff][i];
	  wzAssert(db>-0.0001);
	  if(db<cog->Delta){
	    gg->fPoint[ff][i] = gg->nPoint[n][i];
	    gg->fPointUp[ff][i] += db;
	    gg->fPointDown[ff][i] += db;
	    rx = OnFace; continue;
	  }
	  fb = db/(xi-gg->nPoint[gg->nDown[i][n]][i]);
	  if(fb<ibgFaceRelativeClose)	{rb = CloseFaceType;}
	  else if(rb<FaceType){
	    if(fb<ibgFaceRelativeFar)	{rb = FaceType;}
	    else 			{rb = FarFaceType;}
	  }else if(rb<FarFaceType)	{rb = FarFaceType;}
	}
      }else if(gg->is_outside(nu=gg->nDown[i][n])){
	border++;
	if((fl=gg->bnUp[i][bn]) && (ff=gg->lFaceDown[fl])){
	  db = gg->fPoint[ff][i]-(xi=gg->nPoint[n][i]);
	  wzAssert(db>-0.0001);
	  if(db<cog->Delta){
	    gg->fPoint[ff][i] = gg->nPoint[n][i]; 
	    gg->fPointUp[ff][i] -= db;
	    gg->fPointDown[ff][i] -= db;
	    rx = OnFace; continue;
	  }
	  fb = db/(gg->nPoint[gg->nUp[i][n]][i]-xi);
	  if(fb<ibgFaceRelativeClose)	{rb = CloseFaceType;}
	  else if(rb<FaceType){
	    if(fb<ibgFaceRelativeFar)	{rb = FaceType;}
	    else 			{rb = FarFaceType;}
	  }else if(rb<FarFaceType)	{rb = FarFaceType;}
	}
      }else if(gg->is_undefined(no)){
	if(rx==Default) rx = RefinedNode;
      }else if(gg->is_undefined(nu)){
	if(rx==Default) rx = RefinedNode;
      }else{
	dim++;
	if((fl=gg->bnDown[i][bn]) && (ff=gg->lFaceUp[fl])){
	  di = (xi=gg->nPoint[n][i])-gg->fPoint[ff][i];
	  wzAssert(di>-0.0001);
	  if(di<cog->Delta){
	    gg->fPoint[ff][i] = gg->nPoint[n][i];
	    gg->fPointUp[ff][i] += di;
	    gg->fPointDown[ff][i] += di;
	    rx = OnFace; continue;
	  }
	  fi = di/(xi-gg->nPoint[gg->nDown[i][n]][i]);
	  if(fi<ibgFaceRelativeClose)	{ri = CloseFaceType;}