ibgosplit.cxx

Delaunay三角形的网格剖分程序

#include "ibgoctree.hxx"

/* basic splitting functions (pure splitting, correct data structures
 * before and after splitting, no refinement strategies)

ibgiCell ibgOctree3D::splitCell(ibgiCell qq, unsigned ro);
ibgIndex ibgOctree::splitBoundaryLine(ibgIndex nu, ibgIndex dir, ibgIndex nh);

 */

ibgIndex ibgOctree::splitBoundaryLine(ibgiNode nu, ibgIndex dir, ibgiNode nh)
{
  ibgIndex i,lo,lu,bnh,bnu;
  wzInteger ff,fo;
  nBoundary[nh] = 0;
  bnu = nBoundary[nu]; 
  if(bnu==0)	{return 0;}
  lo = bnUp[dir][bnu];	
  if(lo==0)	{return 0;}
  bnh = createSimpleBoundaryNode(nh);
  lu  = bnUp[dir][bnu] = bnDown[dir][bnh] = lines.create();
  bnUp[dir][bnh]= lo;
  lNodeDown[lo] = bnh;
  lNodeDown[lu] = bnu;
  lNodeUp[lu]   = bnh;
  lType[lu] = lType[lo];
  for(i=0;i<LineSides;i++) lPlane[i][lu]  = lPlane[i][lo];
// does not split planes!!
  lFaceUp[lu] = lFaceDown[lu] = 0;
  ff = lFaceDown[lo];
  if(ff && fPoint[ff][dir]<nPoint[nh][dir]){
    lFaceDown[lu] = ff; fLine[ff] = lu; fo = ff; 
    while((ff = -fUp[fo])>0){
      if(fPoint[ff][dir]>=nPoint[nh][dir]){
	fDown[ff] = fUp[fo] = 0;
	break;
      }
      fLine[ff] = lu; fo = ff;
    }
    lFaceUp[lu] = fo; 
    if(ff<0)  lFaceDown[lo] = 0;
    else      lFaceDown[lo] = ff;
    if(ff==0) lFaceUp[lo] = 0;
    //    wzAssert(fPointUp[fo].segment() == nPoint[nh].segment());
  }else{
    //    wzAssert(nPoint[nh].segment() == nPoint[nu].segment());
  }
  return lu;
}

ibgiCell ibgOctree3D::splitCell(ibgiCell qq, unsigned ro)
{
  int nu,no,nh,n0,n1,n2,n3,nr,nl,nur,nhr,nor,nurr;//,nt;
  int ou,oo,oh;
  int qn,o,w,eo,eu,fo,fu,go,gu;
  int rr,ru=Other(ro);
  
  if(is_undefined(qq)) return 0;
  nu = cNode[(lsCell[ru][3])][qq];  no = cNode[(int)lsCell[ro][3]][qq];
  for(o=0;o<LineSides;o++){
    ou = nu; nu = cNode[lsCell[ru][o]][qq];
    oo = no; no = cNode[lsCell[ro][o]][qq];
    if(no==nNext[ro][nu])	return 0;
    wzAssert(no==nNext[ro][nNext[ro][nu]]);
    if((nor=nNext[lsLine[ro][o]][oo])==no) continue;
    wzAssert(no==nNext[lsLine[ro][o]][nor]);
    if((nur=nNext[lsLine[ru][o]][ou])==nu) continue;
    wzAssert(nu==nNext[lsLine[ro][o]][nur]);
    if(nor==nNext[ro][nur])	return 0;
    wzAssert(is_undefined(nt=nNext[ro][nur])||nor==nNext[ro][nt]);
  }
  
  qn = cells.create();
  cBoundary[qn]=0;
  eu = lsCell[ru][3]; eo = lsCell[ro][3];
  nu = cNode[eu][qq]; nh = nNext[ro][nu]; wzAssert(is_defined(nh));
  no = cNode[eo][qq]; wzAssert(nNext[ro][nh] == no);
  fu = cReverse[eo]; fo = cReverse[eu];
  for(o=0;o<LineSides;o++){
    ou = nu; gu = fu; oh = nh; go = fo;
    eu = lsCell[ru][o];  eo = lsCell[ro][o];
    fu = cReverse[eo]; fo = cReverse[eu];
    nu = cNode[eu][qq];  nh = nNext[ro][nu]; wzAssert(is_defined(nh));
    no = cNode[eo][qq];  wzAssert(nNext[ro][nh] == no);
    cNode[eu][qn] = nu;  cNode[eo][qn] = nh; cNode[eu][qq] = nh;
    nCell[fo][nu] = qn;  nCell[fu][nh] = qn; wzAssert(nCell[fo][nh]==qq);
    rr = lsLine[ru][o];
    if((nhr=nNext[rr][oh])!=nh) {
      wzAssert(is_defined(nhr));
      wzAssert(nNext[rr][nhr]==nh);
      nCell[fu][nhr]=qn; nCell[gu][nhr]=qn;
    }
    if((nur=nNext[rr][ou])!=nu) {
      nCell[fo][nur]=qn; nCell[go][nur]=qn;
      if(o>1)continue;
      if(is_undefined(nurr=nNext[lsLine[ru][o+1]][nur])) continue;
      for(w=0;w<LineSides;w++){
	if(nCell[lsCell[ro][w]][nurr]==qq)
	  nCell[lsCell[ro][w]][nurr]= qn;
      }
    }
  }
  n0 = cNode[lsCell[ru][0]][qq];
  n1 = cNode[lsCell[ru][1]][qq];
  n2 = cNode[lsCell[ru][2]][qq];
  n3 = cNode[lsCell[ru][3]][qq];
  if(n1==(nr=nNext[lsLine[ru][1]][n0])){
    if(n2==(nr=nNext[lsLine[ru][2]][n1])) return qn;
    if(n0==(nl=nNext[lsLine[ru][0]][n3])) return qn;
    nNext[lsLine[ru][3]][nr] = nl;
    nNext[lsLine[ru][1]][nl] = nr;
  }else if(n3==(nl=nNext[lsLine[ru][3]][n2])){
    if(n2==(nr=nNext[lsLine[ru][2]][n1])) return qn;
    if(n0==(nl=nNext[lsLine[ru][0]][n3])) return qn;
    nNext[lsLine[ru][3]][nr] = nl;
    nNext[lsLine[ru][1]][nl] = nr;
  }else  {
    nNext[lsLine[ru][2]][nr] = nl;
    nNext[lsLine[ru][0]][nl] = nr;
  }
  return qn;
}