cogenrefine.cxx

Delaunay三角形的网格剖分程序

#include "cogenrefine.hxx"

static const wzFloat MinDelta = 1.e-16;
static const wzIndex nRefPoint = 0;
static const wzIndex nRefBox = 1;
//static const wzIndex nRefPoint2 = 2;
//static const wzIndex nRefPoint3 = 3;
static const wzIndex nRefRegion = 9;
static const wzIndex fRefFace = 10;
static const wzIndex fRefBoundary = 11;
static const wzIndex fRefBetween = 12;
  
cogenRefine	CogenRefine::cogenrefine() const {return (CogenRefine*) this;}

CogenRefinementFunctionAnisotropic::CogenRefinementFunctionAnisotropic
(cogenerator gen, wzIndex f, wzIndex d)
  :CogenDecorator(gen)
  ,oldRefinement(gen->refinement())
  ,off(f),dim(d)
{
  wzPoint p; if (p.floatValues() < f+d) p.setFloatValues(f+d);
}

CogenRefinementFunctionIsotropic::CogenRefinementFunctionIsotropic
(cogenerator gen, wzIndex f)
  :CogenDecorator(gen)
  ,oldRefinement(gen->refinement())
  ,off(f)
{
  wzPoint p; if (p.floatValues() < f+1) p.setFloatValues(f+1);
}

CogenRefine::CogenRefine(cogenerator gen)
  :CogenDecorator(gen)
  ,nList(sizeof(CogenRefineItem))
  ,nItem(nList.base)
  ,fList(sizeof(CogenRefineItem))
  ,fItem(fList.base)
  ,oldRefinement(gen->refinement())
  ,PointMode(nRefPoint)
{
}

void CogenRefine::setPointModeRectangular()
{PointMode = nRefBox;}
void CogenRefine::setPointModeSpherical()  
{PointMode = nRefPoint;}

// isotropic variants:

wzIndex CogenRefine::refinementIsotropicInRegion(wzFloat d, wzRegion r)
{return refinementInRegion(d,d,d,r);}

wzIndex CogenRefine::refinementIsotropicNearBoundary(wzFloat d, wzRegion r)
{return refinementNearBoundary(d,d,d,r);}

wzIndex CogenRefine::refinementIsotropicBetween(wzFloat d, wzRegion r1, wzRegion r2)
{return refinementBetween(d,d,d,r1,r2);}

wzIndex CogenRefine::refinementIsotropicNearFace(wzFloat d, wzFace f)
{return refinementNearFace(d,d,d,f);}

wzIndex CogenRefine::refinementIsotropicNearPoint(wzFloat d, wzFloat n,
						    wzFloat xx, wzFloat yy, wzFloat zz)
{return refinementNearPoint(d,d,d,n,n,n,xx,yy,zz);}

wzIndex CogenRefine::refinementInRegion(wzFloat dx, wzFloat dy, wzFloat dz, wzRegion s)
{
  wzIndex i = nList.create();
  CogenRefineItem& r = nItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = nRefRegion;
  r.s    = s;
  return i;
}

wzIndex CogenRefine::refinementNearFace(wzFloat dx, wzFloat dy, wzFloat dz, wzFace s)
{
  wzIndex i = fList.create();
  CogenRefineItem& r = fItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = fRefFace;
  r.s    = s;
  return i;
}

wzIndex CogenRefine::refinementNearBoundary(wzFloat dx, wzFloat dy, wzFloat dz, wzRegion s)
{
  wzIndex i = fList.create();
  CogenRefineItem& r = fItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = fRefBoundary;
  r.s    = s;
  return i;
}

wzIndex CogenRefine::refinementBetween(wzFloat dx, wzFloat dy, wzFloat dz, 
				       wzRegion s1, wzRegion s2)
{
  wzIndex i = fList.create();
  CogenRefineItem& r = fItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = fRefBetween;
  r.s    = s1;
  r.so   = s2;
  return i;
}

wzIndex CogenRefine::refinementNearPoint(wzFloat dx, wzFloat dy, wzFloat dz,
					 wzFloat nx, wzFloat ny, wzFloat nz,
					 wzFloat xx, wzFloat yy, wzFloat zz)
{
  wzIndex i = nList.create();
  CogenRefineItem& r = nItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  if(nx<dx) nx=dx;
  if(ny<dy) ny=dy;
  if(nz<dz) nz=dz;
  r.near[0] = nx;
  r.near[1] = ny;
  r.near[2] = nz;
  r.x[0] = xx;
  r.x[1] = yy;
  r.x[2] = zz;
  r.typ  = PointMode;
  return i;
}

void CogenRefine::getMetric(wzMetricData& data, const cogFlag1& f) const
{
  wzIndex it,i;
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  oldRefinement->getMetric(data,f);
  //  wzDiagonalMetricSimple::getMetric(data,f);
  wzRangeLoop(fList,it){
    const CogenRefineItem& r = fItem[it];
    switch(r.typ){
    case fRefFace:
      if(r.s != f.p1.segment()) continue; 
      break;
    case fRefBetween:
      if(r.s == f.p0.segment()){
	if(r.so != f.po.segment()) continue; 
      }else if(r.so == f.p0.segment()){
	if(r.s != f.po.segment()) continue; 
      }else{
	continue;
      }
      break;
    case fRefBoundary:
      if(r.s != f.p0.segment()){
	if(r.s != f.po.segment()) continue; 
      }
      break;
    default:
      wzAssert(0); continue;
    }
    for(i=0;i<wzPointDim;i++){
      if(d.g_ii[i]<r.dist[i]) d.g_ii[i] = r.dist[i];
    }
  }
}

void CogenRefine::getMetric(wzMetricData& data, const wzPoint& p) const
{
  wzIndex it,i;
  wzFloat x,f,g;
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  oldRefinement->getMetric(data,p);
  for(i=0;i<wzPointDim;i++) 
    if(d.g_ii[i] < G_ii[i]) d.g_ii[i] = G_ii[i];
  wzRangeLoop(nList,it){
    const CogenRefineItem& r = nItem[it];
    switch(r.typ){
    case nRefRegion:
      if(r.s == p.segment()){
	for(i=0;i<wzPointDim;i++){
	  if(d.g_ii[i]<r.dist[i]) d.g_ii[i] = r.dist[i];
	}
      }
      break;
    case nRefBox:
      f = 1;
      for(i=0;i<wzPointDim;i++){
	x = (r.x[i] - p[i])/r.near[i];
	x *= x; x *= x; x *= x;
	f *= 1/(1+x*x);
      }
      for(i=0;i<wzPointDim;i++){
	g = f*r.dist[i];
	if(d.g_ii[i]<g) d.g_ii[i] = g;
      }
      break;
    case nRefPoint:
      f = 1;
      for(i=0;i<wzPointDim;i++){
	x = (r.x[i] - p[i])/r.near[i];
	f += x*x;
      }
      for(i=0;i<wzPointDim;i++){
	g = r.dist[i]/f;
	if(d.g_ii[i]<g) d.g_ii[i] = g;
      }
      break;
    }
  }
}

void CogenRefinementFunctionAnisotropic::getMetric
(wzMetricData& data, const wzPoint& p) const
{
  oldRefinement->getMetric(data,p);
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  for(wzIndex i=0;i<dim;i++){
    wzFloat dx = p.f(off+i),g;
    if(dx<MinDelta) dx=MinDelta; 
    g = 1.0/(dx*dx);
    if(d.g_ii[i]<g) d.g_ii[i] = g;
  }
}

void CogenRefinementFunctionAnisotropic::getMetric
(wzMetricData& data, const cogFlag1& f) const
{
  oldRefinement->getMetric(data,f);
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  for(wzIndex i=0;i<dim;i++){
    wzFloat dx = f.p1.f(off+i),g;
    if(dx<MinDelta) dx=MinDelta; 
    g = 1.0/(dx*dx);
    if(d.g_ii[i]<g) d.g_ii[i] = g;
  }
}

void CogenRefinementFunctionIsotropic::getMetric
(wzMetricData& data, const wzPoint& p) const
{
  oldRefinement->getMetric(data,p);
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  wzFloat dx = p.f(off),g;
  if(dx<MinDelta) dx=MinDelta; 
  g = 1.0/(dx*dx);
  for(wzIndex i=0;i<wzPointDim;i++){
    if(d.g_ii[i]<g) d.g_ii[i] = g;
  }
}

void CogenRefinementFunctionIsotropic::getMetric
(wzMetricData& data, const cogFlag1& f) const
{
  oldRefinement->getMetric(data,f);
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  wzFloat dx = f.p1.f(off),g;
  if(dx<MinDelta) dx=MinDelta; 
  g = 1.0/(dx*dx);
  for(wzIndex i=0;i<wzPointDim;i++){
    if(d.g_ii[i]<g) d.g_ii[i] = g;
  }
}