wzgridcont.cxx

Delaunay三角形的网格剖分程序

#include "wzgrid.hxx"

wzgrid wzGridContinuation3D(wzgrid grid2D, wzFloat dz, wzIndex nz)
{
  wzIndex nu,no,o0,o1,o2,or,i,j,nj,cj,cd,os;
  wzCellType ot;
  wzSegment s;
  wzInteger *ocn,*ccn;
  wzFloat z;
  wzPoint *point,*pbase;
  wzgrid gg = new wzGrid(3,3);
  wzIndex n = grid2D->points.length();
  wzIndex c = grid2D->cells.length();
  gg->points.create((nz+1)*n);
  gg->cells.create( 3*nz*c);
  gg->data.create(6*3*nz*c);
  for(j=0;j<=nz;j++){
    nj = n*j; z = dz*j;
    wzRangeLoop(grid2D->points,i){
      point = &gg->Point[++nj];
      pbase = &grid2D->Point[i];
      point->x() = pbase->x();
      point->y() = pbase->y();
      point->z() = z;
      point->segment() = pbase->segment();
    }
  }
  cj = 0; cd = 1;
  wzRangeLoop(grid2D->cells,i){
    ot  = grid2D->CellType[i];
    if(ot != wzCellType2Triangle) continue;
    os  = grid2D->CellSegment[i];
    ocn = &grid2D->Data[grid2D->CellData[i]];
    if(ocn[0]<ocn[1]) 
      if     (ocn[1]<ocn[2]){o0=0;o1=1;o2=2;or=0;}
      else if(ocn[0]<ocn[2]){o0=0;o1=2;o2=1;or=1;}
      else                  {o0=2;o1=0;o2=1;or=0;}
    else
      if     (ocn[1]>ocn[2]){o0=2;o1=1;o2=0;or=1;}
      else if(ocn[0]>ocn[2]){o0=1;o1=2;o2=0;or=0;}
      else                  {o0=1;o1=0;o2=2;or=1;}
    wzAssert(ocn[o0]<ocn[o1]);
    wzAssert(ocn[o1]<ocn[o2]);
    wzAssert(ocn[o0]>0);
    wzAssert(ocn[o2]<=n);
    wzAssert(or<2);
    nu = 0;
    for(j=0;j<nz;j++){
      no = nu+n;
      // first tetradehron:
      cj++; cd += 8;
      gg->CellType[cj] = wzCellType3Tetrahedron;
      gg->CellData[cj] = cd;
      gg->CellSegment[cj] = os;
      ccn = &gg->Data[cd];
      ccn[0] = nu + ocn[o0];
      ccn[1] = nu + ocn[o1];
      ccn[2] = nu + ocn[o2];
      ccn[3] = no + ocn[o2];
      // second tetradehron:
      cj++; cd += 8;
      gg->CellType[cj] = wzCellType3Tetrahedron;
      gg->CellData[cj] = cd;
      gg->CellSegment[cj] = os;
      ccn = &gg->Data[cd];
      ccn[0] = nu + ocn[o0];
      ccn[1] = nu + ocn[o1];
      ccn[2] = no + ocn[o1];
      ccn[3] = no + ocn[o2];
      // third tetradehron:
      cj++; cd += 8;
      gg->CellType[cj] = wzCellType3Tetrahedron;
      gg->CellData[cj] = cd;
      gg->CellSegment[cj] = os;
      ccn = &gg->Data[cd];
      ccn[0] = nu + ocn[o0];
      ccn[1] = no + ocn[o0];
      ccn[2] = no + ocn[o1];
      ccn[3] = no + ocn[o2];
      nu=no;
    }
  }
  return gg;
}