wzmetric.cxx

Delaunay三角形的网格剖分程序

#include <math.h>
#include "wzmetric.hxx"

void wzDiagonalMetricSimple::getMetric(wzMetricData& data, const wzPoint& p) const
{
  int i;
  for(i=0;i<wzPointDim;i++) ((wzDiagonalMetricData*)(&data))->g_ii[i] = G_ii[i];
}

void wzDiagonalMetricSimple::getMetric(wzMetricData& data, const cogFlag1& f) const
{
  wzIndex i;
  wzFloat x,sq=0,l,oo[wzPointDim],o,t;
  for(i=0;i<wzPointDim;i++){
    x = oo[i] = f.p0[i] - f.po[i];
    sq += x*x;
  }
  l = sqrt(sq);
  if(l>0)
    for(i=0;i<wzPointDim;i++){
      o = oo[i]/l; t = 1-o;
      ((wzDiagonalMetricData*)(&data))->g_ii[i] = f_oo*o*o + f_tt*t*t;
    }
  else
    for(i=0;i<wzPointDim;i++){
      ((wzDiagonalMetricData*)(&data))->g_ii[i] = f_tt;
    }
}

wzFloat wzDiagonalMetricSimple::g_ii(const wzMetricData& data, 
		       const cogFlag1& f, wzIndex i) const
{
  return ((wzDiagonalMetricData*)(&data))->g_ii[i];
}

void wzMetricByFunction::getMetric(wzMetricData& data, const wzPoint& p) const
{
  int i;
  wzFloat f = (*Function)(p); 
  wzFloat g = 1.0/(f*f);
  for(i=0;i<wzPointDim;i++){
    if(g<G_ii[i])	((wzDiagonalMetricData*)(&data))->g_ii[i] = G_ii[i];
    else		((wzDiagonalMetricData*)(&data))->g_ii[i] = g;
  }
}

void wzMetricByFunction::getMetric(wzMetricData& data, const cogFlag1& f) const
{
  wzDiagonalMetricSimple::getMetric(data,f);
}

wzFloat wzMetric::distance(const wzMetricData& data,
			   const wzPoint& p, const wzPoint& q) const
{
  wzFloat d=0,x;
  int i;
  for(i=0;i<wzPointDim;i++){
    x=q[i]-p[i];
    d += x*x;
  }
  return d;
}

wzFloat wzMetric::distance(const wzMetricData& data,
			   const cogFlag1& f, const wzPoint& q) const
{
  wzFloat d=0,x; wzPoint& p=f.p1;
  int i;
  for(i=0;i<wzPointDim;i++){
    x=q[i]-p[i];
    d += x*x;
  }
  return d;
}

wzFloat wzDiagonalMetricSimple::distance(const wzMetricData& data,
					 const wzPoint& p, const wzPoint& q) const
{
  wzFloat d=0,x;
  int i;
  for(i=0;i<wzPointDim;i++){
    x=q[i]-p[i];
    d += ((wzDiagonalMetricData*)(&data))->g_ii[i]*x*x;
  }
  return d;
}

wzFloat wzDiagonalMetricSimple::distance(const wzMetricData& data,
					 const cogFlag1& f, const wzPoint& q) const
{
  wzFloat d=0,x; wzPoint& p=f.p1;
  int i;
  for(i=0;i<wzPointDim;i++){
    x=q[i]-p[i];
    d += ((wzDiagonalMetricData*)(&data))->g_ii[i]*x*x;
  }
  return d;
}

wzDiagonalMetricSimple::wzDiagonalMetricSimple(wzFloat dx, wzFloat dy, wzFloat dz)
{
  refinementGlobal(dx,dy,dz);
  faceRefinementGlobal(wzInfty,wzInfty);
}

void wzDiagonalMetricSimple::refinementGlobal(wzFloat dx, wzFloat dy, wzFloat dz)
{
  dx *= 1.0000001;
  dy *= 1.0000001;
  dz *= 1.0000001;
  if(dx==0) G_ii[0] = 0; else try{ G_ii[0] = 1/(dx*dx);}catch(...){G_ii[0] = 0;}
  if(dy==0) G_ii[1] = 0; else try{ G_ii[1] = 1/(dy*dy);}catch(...){G_ii[1] = 0;}
  if(dz==0) G_ii[2] = 0; else try{ G_ii[2] = 1/(dz*dz);}catch(...){G_ii[2] = 0;}
}

void wzDiagonalMetricSimple::faceRefinementGlobal(wzFloat o, wzFloat t)
{
  if(o<1.e-16) o=1.e-16;
  if(t<1.e-16) t=1.e-16;
  f_oo = 1/(o*o); 
  f_tt = 1/(t*t); 
}