ibgregion.cxx

Delaunay三角形的网格剖分程序

#include "wzgrid.hxx"
#include "cogwarnings.hxx"

// if >0, the region number for errorneous region cells;
wzIndex ibgErrorneousRegion=19;
//const  wzIndex defineRegionsWithInnerPoint = 1;
const  wzIndex defineRegionsWithNodes = 2;
static wzIndex defineRegionsMode = defineRegionsWithNodes;
//static defineRegionsMode = defineRegionsWithInnerPoint;

void wzDelaunayGenerator::defineRegions(cogeometry g)
{
  int c;
  for(c=1;c<=Grid.cells.last();c++){
    if(ccund(c)) continue;
    defineRegionOfCell(c,g);
  }
  Grid.firstRegion = 1;
  Grid.lastRegion  = Grid.cells.last();
}

void wzDelaunayGenerator::defineRegionOfCell(int c, cogeometry g)
{
  int n,k,nn,m=0,m1;
  wzPoint xx;
  wzCellType t = cct(c);
  int npoints = wzCellTypePoints(t);
  int *cn=ccn(c);
  for(n=0;n<npoints;n++){
    switch(cnt(nn=cn[n])){
    case wzIsRegion:
      if((m1=cnu(nn)) <= 0){ccu(c) = 0; return;}
      if(m){
	if(m1 != m) m = -1;
      }else{
	m = m1;
      }
    }
  }
  if(m>0 && (defineRegionsMode==defineRegionsWithNodes)){ccu(c)=m; return;}
  if(m<0){
    cogInfoBadRegionCell();
    if(ibgErrorneousRegion){
      ccu(c) = ibgErrorneousRegion; return;
    }
  }
  for(n=0;n<npoints;n++){
    nn=cn[n]; for(k=0;k<wzPointDim;k++) xx.X[k] += cnx(nn)[k];
  }
  for(k=0;k<wzPointDim;k++)  xx.X[k] /= npoints;
  g->setStartpoint(*cnd(nn)); g->Point(xx);
  ccu(c)=xx.segment().index();
}