cogenoctree.cxx

Delaunay三角形的网格剖分程序

#include "ibgoctree.hxx"
#include "ibg2.hxx"
#include <stdlib.h>
#include "cogenoctree.hxx"
#include "ibggenerator.hxx"
#include "cogwarnings.hxx"

wzFloat  CogenOctree::epsilon = 0;
wzFloat  CogenOctree::relativeMinimalDistance[3];

wzBoolean  CogenOctree::SimpleBinaryIntersectionForLine = 1;

const wzIndex cogFullBoxType	= 0;
const wzIndex cogWedgeType	= 1;

cogenOctree	CogenOctree::cogenoctree() const {return (CogenOctree *) this;}

CogenOctree::CogenOctree()
  :blist(sizeof(CogenOctreeBox))
  ,flist(sizeof(CogenOctreeFace))
  ,rlist(sizeof(CogenOctreeRefineItem))
  ,box(blist.base)
  ,face(flist.base)
  ,refine(rlist.base)
  ,boundaryLevel(4)
  ,RegionNodeRemoval(0)
  ,FaceNodeRemoval(0)
  ,EdgeNodeRemoval(0)
  ,VertexNodeRemoval(0)
{
  setBox();
  relativeMinimalDistance[0] = 1.e-3;
  relativeMinimalDistance[1] = 1.e-3;
  relativeMinimalDistance[2] = 1.e-3;
  // defining an epsilon close enough to the rounding error:
  epsilon = 1.0;
  for(int i=0;i<200;i++){
    if(((wzFloat)(1.0+epsilon))<=(wzFloat)1.0) break;
    epsilon /= 2;
  }
  // changing this epsilon seems necessary for nonlinear coordinates
  // with possibly much higher inaccuracy;
  epsilon *= 10;
  // enforceRegionNodeRemoval();
  //  enforceFaceNodeRemoval();
  //  omitRegionNodeRemoval();
  //  omitFaceNodeRemoval();
  removeUnderdogs(wzTrue);
  removeIrregularNodes(wzTrue);
  removeLongSideIntersections(wzTrue);
}

void CogenOctree::setBorder(
			    wzFloat xmin, wzFloat xmax,
			    wzFloat ymin, wzFloat ymax,
			    wzFloat zmin, wzFloat zmax)
{
  min[0]=xmin;  max[0]=xmax;
  min[1]=ymin;  max[1]=ymax;
  min[2]=zmin;  max[2]=zmax;
  for(int i=0;i<wzPointDim;i++){
    if(min[i] != -wzInfty) addPlane(i,min[i]);
    if(max[i] !=  wzInfty) addPlane(i,max[i]);
  }
}

void CogenOctree::addPlane(wzIndex dir, wzFloat x)
{
  if(x>=min[dir] && x<=max[dir]) wzTensorGrid::add(dir,x);
}

void CogenOctree::setFaceInBox(wzFace f,
			      wzFloat xmin,wzFloat xmax,
			      wzFloat ymin,wzFloat ymax,
			      wzFloat zmin,wzFloat zmax)
{
  wzFloat Min[3]={xmin,ymin,zmin};
  wzFloat Max[3]={xmax,ymax,zmax};
  setFaceInBox(f,Min,Max);
}

void CogenOctree::setFaceInWedge(wzFace f, wzWedge dir,
				 wzFloat xmin,wzFloat xmax,
				 wzFloat ymin,wzFloat ymax,
				 wzFloat zmin,wzFloat zmax)
{
  wzFloat Min[3]={xmin,ymin,zmin};
  wzFloat Max[3]={xmax,ymax,zmax};
  setFaceInWedge(f,dir,Min,Max);
}

void CogenOctree::computeWedgeData(wzFloat *a, wzFloat *d,
				   wzWedge dir,
				   wzFloat *xmin, wzFloat *xmax)
{
  wzFloat x[3];
  switch(dir){
  case wzWedgeRightUp:
    x[0]=xmin[0]; x[1]=xmax[1]; x[2]=xmin[2];
    a[0]=xmax[1]-xmin[1]; a[1]=xmax[0]-xmin[0]; a[2]=0;
    break;
  case wzWedgeRightDown:
    x[0]=xmin[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=xmax[1]-xmin[1]; a[1]=xmin[0]-xmax[0]; a[2]=0;
    break;
  case wzWedgeLeftUp:
    x[0]=xmin[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=xmin[1]-xmax[1]; a[1]=xmax[0]-xmin[0]; a[2]=0;
    break;
  case wzWedgeLeftDown:
    x[0]=xmin[0]; x[1]=xmax[1]; x[2]=xmin[2];
    a[0]=xmin[1]-xmax[1]; a[1]=xmin[0]-xmax[0]; a[2]=0;
    break;
  case wzWedgeRightNear:
    x[0]=xmax[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=xmax[2]-xmin[2]; a[1]=0; a[2]=xmax[0]-xmin[0];
    break;
  case wzWedgeRightFar:
    x[0]=xmin[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=xmax[2]-xmin[2]; a[1]=0; a[2]=xmin[0]-xmax[0];
    break;
  case wzWedgeLeftNear:
    x[0]=xmin[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=xmin[2]-xmax[2]; a[1]=0; a[2]=xmax[0]-xmin[0];
    break;
  case wzWedgeLeftFar:
    x[0]=xmax[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=xmin[2]-xmax[2]; a[1]=0; a[2]=xmin[0]-xmax[0];
    break;
  case wzWedgeUpNear:
    x[0]=xmin[0]; x[1]=xmax[1]; x[2]=xmin[2];
    a[0]=0; a[1]=xmax[2]-xmin[2]; a[2]=xmax[1]-xmin[1];
    break;
  case wzWedgeUpFar:
    x[0]=xmin[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=0; a[1]=xmax[2]-xmin[2]; a[2]=xmin[1]-xmax[1];
    break;
  case wzWedgeDownNear:
    x[0]=xmin[0]; x[1]=xmin[1]; x[2]=xmin[2];
    a[0]=0; a[1]=xmin[2]-xmax[2]; a[2]=xmax[1]-xmin[1];
    break;
  case wzWedgeDownFar:
    x[0]=xmin[0]; x[1]=xmax[1]; x[2]=xmin[2];
    a[0]=0; a[1]=xmin[2]-xmax[2]; a[2]=xmin[1]-xmax[1];
    break;
  }
  *d = 0;
  for(int i=0;i<3;i++){
    *d 	-= x[i]*a[i];
  }
}

void CogenOctree::setFaceInWedge(wzFace f, wzWedge dir,
				 wzFloat *xmin, wzFloat *xmax)
{
  int i;
  for(i=0;i<3;i++){
    if(xmax[i]<min[i]){cogInfoEmptyFace();}
    if(xmin[i]>max[i]){cogInfoEmptyFace();}
  }
  wzIndex l = flist.create();
  face[l].face = f;
  face[l].type = cogWedgeType;
  face[l].part = dir;
  computeWedgeData(face[l].a,&(face[l].d),dir,xmin,xmax);
  face[l].d    += epsilon;
  for(i=0;i<3;i++){
    if(xmin[i]<min[i])	face[l].min[i] = min[i];
    else		face[l].min[i] = xmin[i];
    if(xmax[i]>max[i])	face[l].max[i] = max[i];
    else		face[l].max[i] = xmax[i];
  }
}

void CogenOctree::setFaceInBox(wzFace f, wzFloat *xmin, wzFloat *xmax)
{
  int i;
  for(i=0;i<3;i++){
    if(xmax[i]<min[i]){cogInfoEmptyFace();}
    if(xmin[i]>max[i]){cogInfoEmptyFace();}
  }
  wzIndex l = flist.create();
  face[l].type = cogFullBoxType;
  face[l].face = f;
  for(i=0;i<3;i++){
    if(xmin[i]<min[i])	face[l].min[i] = min[i];
    else		face[l].min[i] = xmin[i];
    if(xmax[i]>max[i])	face[l].max[i] = max[i];
    else		face[l].max[i] = xmax[i];
    face[l].a[i] = 0;
  }
  face[l].d = 0;
}

wzIndex CogenOctree::addBox(wzRegion r,
			    wzFloat xmin,wzFloat xmax,
			    wzFloat ymin,wzFloat ymax,
			    wzFloat zmin,wzFloat zmax)
{
  wzFloat Min[3]={xmin,ymin,zmin};
  wzFloat Max[3]={xmax,ymax,zmax};
  return addBox(r,Min,Max);
}

wzIndex CogenOctree::addWedge(wzRegion r, wzWedge dir,
			    wzFloat xmin,wzFloat xmax,
			    wzFloat ymin,wzFloat ymax,
			    wzFloat zmin,wzFloat zmax)
{
  wzFloat Min[3]={xmin,ymin,zmin};
  wzFloat Max[3]={xmax,ymax,zmax};
  return addWedge(r,dir,Min,Max);
}

wzIndex CogenOctree::addBox(wzRegion r,wzFloat *xmin, wzFloat *xmax)
{
  int i;
  for(i=0;i<3;i++){
    if(xmax[i]<min[i]){cogInfoEmptyRegion();}
    if(xmin[i]>max[i]){cogInfoEmptyRegion();}
  }
  wzIndex l = blist.create();
  box[l].type   = cogFullBoxType;
  box[l].region = r;
  for(i=0;i<3;i++){
    if(xmin[i]<min[i])	box[l].min[i] = min[i];
    else		box[l].min[i] = xmin[i];
    if(xmax[i]>max[i])	box[l].max[i] = max[i];
    else		box[l].max[i] = xmax[i];
    box[l].a[i] = 0;
    box[l].mino[i] = box[l].min[i];
    box[l].maxo[i] = box[l].max[i];
  }
  box[l].d = box[l].dorig = 0;
  return l;
}

wzIndex CogenOctree::addWedge(wzRegion r, wzWedge dir, wzFloat *xmin, wzFloat *xmax)
{
  int i;
  for(i=0;i<3;i++){
    if(xmax[i]<min[i]){cogInfoEmptyRegion();}
    if(xmin[i]>max[i]){cogInfoEmptyRegion();}
  }
  wzIndex l = blist.create();
  box[l].type   = cogWedgeType;
  box[l].part   = dir;
  box[l].region = r;
  computeWedgeData(box[l].a,&(box[l].d),dir,xmin,xmax);
  box[l].dorig  = box[l].d;
  for(i=0;i<3;i++){
    if(xmin[i]<min[i])	box[l].min[i] = min[i];
    else		box[l].min[i] = xmin[i];
    if(xmax[i]>max[i])	box[l].max[i] = max[i];
    else		box[l].max[i] = xmax[i];
    box[l].mino[i] = box[l].min[i];
    box[l].maxo[i] = box[l].max[i];
  }
  return l;
}

static const wzFloat dfac = 1+1.e-7;
static const wzFloat dadd = 1/(wzInfty*wzInfty);

wzIndex CogenOctree::refinementInBox(wzFloat dx, wzFloat dy, wzFloat dz, wzIndex b)
{
  wzIndex r = rlist.create(); wzFloat dd;
  refine[r].box = b;
  dd = dfac*dx*dx + dadd;  refine[r].dist[0] = 1/dd;
  dd = dfac*dy*dy + dadd;  refine[r].dist[1] = 1/dd;
  dd = dfac*dz*dz + dadd;  refine[r].dist[2] = 1/dd;
  return r;
}

void CogenOctree::setFaceOfBox(wzFace f, wzIndex b)
{
  CogenOctreeBox& bb = *(box+b);
  wzAssert(bb.type == cogFullBoxType);
  setFaceInBox(f,bb.min,bb.max);
}

void CogenOctree::setFaceOfWedge(wzFace f, wzIndex w)
{
  CogenOctreeBox& ww = *(box+w);
  wzAssert(ww.type == cogWedgeType);
  setFaceInWedge(f,(wzWedge)ww.part,ww.min,ww.max);
}

void CogenOctree::setFaceOfWedgeSide(wzFace f, wzIndex w)
{setFaceOfWedge(f,w);}  // incorrect simplification;

void CogenOctree::setFaceOfWedgeSide(wzFace f, wzIndex w, wzIndex dir)
{setFaceOfBoxSide(f,w,dir);}

void CogenOctree::setFaceOfBoxSide(wzFace f, wzIndex b, wzIndex dir)
{
  CogenOctreeBox& B = *(box+b);
  switch(dir){
  case 0:setFaceInBox(f,B.max[0],B.max[0],B.min[1],B.max[1],B.min[2],B.max[2]);break;
  case 1:setFaceInBox(f,B.min[0],B.max[0],B.max[1],B.max[1],B.min[2],B.max[2]);break;
  case 2:setFaceInBox(f,B.min[0],B.max[0],B.min[1],B.max[1],B.max[2],B.max[2]);break;
  case 3:setFaceInBox(f,B.min[0],B.min[0],B.min[1],B.max[1],B.min[2],B.max[2]);break;
  case 4:setFaceInBox(f,B.min[0],B.max[0],B.min[1],B.min[1],B.min[2],B.max[2]);break;
  case 5:setFaceInBox(f,B.min[0],B.max[0],B.min[1],B.max[1],B.min[2],B.min[2]);break;
  }
}

// Epsilong should be something absolute, like 1.e-7.  Renormalization
// in agreement with the size of the grid will be done below.

void CogenOctree::setEpsilon(wzFloat epsilonNew)
{
  epsilon = epsilonNew;
  if(!checked){
    return;
  }
  int i,l;
  wzFloat delta;
  if(max[0]>min[0]) delta = max[0]-min[0];
  if(max[1]>min[1] && max[1]<min[1]+delta) delta = max[1]-min[1];
  if(max[2]>min[2] && max[2]<min[2]+delta) delta = max[2]-min[2];
  wzRangeLoop(blist,l){
    CogenOctreeBox& rr = box[l];
    for(i=0;i<3;i++){
      if(rr.mino[i]>=0)