ibgoplane.cxx

Delaunay三角形的网格剖分程序

#include "ibgoctree.hxx"

ibgIndex ibgOctree::getPlane(ibgIndex e, ibgIndex s)
{
  ibgIndex bnu=lNodeDown[e];
  ibgIndex d  = lDirection(e);
  ibgIndex pi = lsPlane[d][s];
  ibgIndex pp = bnPlane[pi][bnu];
  if(!pp) return createSimplePlane(e,s);
  return pp;
}

ibgIndex ibgOctree::createSimplePlane(ibgIndex e, ibgIndex s)
{
  ibgIndex bnu=lNodeDown[e],nu=bnBase[bnu],bnur;
  ibgIndex bno=lNodeUp[e],  no=bnBase[bno],bnor;
  ibgIndex bnh,bnuh,bnoh,bnrh;
  ibgIndex nur,nor,nuro,nh=0,nuh=0,noh=0,nrh=0,nt;
  ibgIndex d = lDirection(e),rd,dd;
  ibgInteger cl,cr,cb;
  wzAssert(nUp[d][nu]==no);
  ibgIndex r = lsLine[d][s],p,p00,p10,p01,p11;
  if(is_undefined(nur=nNext[r][nu])){
    nh = nu;
    if(is_undefined(nu=nDown[d][nu])) return 0;
    if(is_undefined(nur=nNext[r][nu])) return 0;
    if(is_undefined(nor=nNext[r][no])) return 0;
  }else if(is_undefined(nor=nNext[r][no])){
    nh = no;
    if(is_undefined(no=nUp[d][no])) return 0;
    if(is_undefined(nor=nNext[r][no])) return 0;
  }
  if(is_undefined(nuro=nUp[d][nur])){
    nuh = nur;
    if(is_undefined(nur=nNext[r][nur])) return 0;
    if(is_undefined(nuro=nUp[d][nur])) return 0;
  }
  if(nuro==nor){
  }else if(nor==nUp[d][nuro]){
    nrh = nuro;
  }else if(nuro==nNext[r][nor]){
    noh = nor; nor = nuro;
  }else{
    noh = nor; nrh = nuro;
    if(is_undefined(nor=nUp[d][nuro])) return 0;
    if(nor != nNext[r][noh]) return 0;
  }
  if(isDown(r)){
    r = up(r);   
    nt=no;no=nor;nor=nt;
    nt=nu;nu=nur;nur=nt;
    nt=nh;nh=nrh;nrh=nt;
  }
  if(r<d){
    dd = r; r=d; d=dd;
    nt=nh;nh=nuh;nuh=nt;
    nt=no;no=nur;nur=nt;
    nt=noh;noh=nrh;nrh=nt;
  }
  wzAssert(d<r);
  dd = down(d);
  rd = down(r);
  if(!(bnu=nBoundary[nu]))	bnu=createSimpleBoundaryNode(nu);
  if(!(bno=nBoundary[no]))	bno=createSimpleBoundaryNode(no);
  if(!(bnur=nBoundary[nur]))	bnur=createSimpleBoundaryNode(nur);
  if(!(bnor=nBoundary[nor]))	bnor=createSimpleBoundaryNode(nor);
  p00=llPlane[d][r];  wzAssert(bnPlane[p00][bnu]==0);
  p10=llPlane[dd][r]; wzAssert(bnPlane[p10][bno]==0);
  p01=llPlane[d][rd]; wzAssert(bnPlane[p01][bnur]==0);
  p11=llPlane[dd][rd];wzAssert(bnPlane[p11][bnor]==0);
  if(PlaneSides){
    cl = nCell[llCell[pSide0[p00]][pSide1[p00]]][nu];
    cr = nCell[llCell[pSide1[p00]][pSide0[p00]]][nu];
    if(cl==cr){return 0;}
  }
  p = planes.create();
  bnPlane[p00][bnu]=p;	pNode[pReverse[p00]][p] = bnu;
  bnPlane[p10][bno]=p;  pNode[pReverse[p10]][p] = bno;
  bnPlane[p01][bnur]=p;	pNode[pReverse[p01]][p] = bnur;
  bnPlane[p11][bnor]=p;	pNode[pReverse[p11]][p] = bnor;
  pEdge[p] = 0;
  pType[p] = pOrtho[p00];
  if(PlaneSides){
    if(cl) if(cb=cBoundary[cl]) pCell[0][p] = cb; else pCell[0][p] = -cl; else pCell[0][p] = 0;
    if(cr) if(cb=cBoundary[cr]) pCell[1][p] = cb; else pCell[1][p] = -cr; else pCell[1][p] = 0;
  }
  if(nh){
    if(!(bnh=nBoundary[nh])) bnh=createSimpleBoundaryNode(nh);
    wzAssert(bnPlane[p00][bnh]==0);  bnPlane[p00][bnh]=p;
    wzAssert(bnPlane[p10][bnh]==0);  bnPlane[p10][bnh]=p;
  }
  if(nuh){
    if(!(bnuh=nBoundary[nuh])) bnuh=createSimpleBoundaryNode(nuh);
    wzAssert(bnPlane[p00][bnuh]==0);  bnPlane[p00][bnuh]=p;
    wzAssert(bnPlane[p01][bnuh]==0);  bnPlane[p01][bnuh]=p;
  }
  if(noh){
    if(!(bnoh=nBoundary[noh])) bnoh=createSimpleBoundaryNode(noh);
    wzAssert(bnPlane[p11][bnoh]==0);  bnPlane[p11][bnoh]=p;
    wzAssert(bnPlane[p10][bnoh]==0);  bnPlane[p10][bnoh]=p;
  }
  if(nrh){
    if(!(bnrh=nBoundary[nrh])) bnrh=createSimpleBoundaryNode(nrh);
    wzAssert(bnPlane[p11][bnrh]==0);  bnPlane[p11][bnrh]=p;
    wzAssert(bnPlane[p01][bnrh]==0);  bnPlane[p01][bnrh]=p;
  }
  return p;
}