cog.cxx

Delaunay三角形的网格剖分程序

#include <math.h>
#include "cog.hxx"
wzRegion 	Cogeometry::InitialDefaultRegion = wzRegion(1);
wzFace		Cogeometry::InitialDefaultFace	 = wzFace(1);
Cogeometry::Cogeometry(cogFloat delta)
  :DefaultRegion(InitialDefaultRegion)
  ,DefaultFace  (InitialDefaultFace)
{
  setDelta(delta);
  specialRefinement = 0;
  cogLinkPrintCommands();
  Startpoint = 0;
  floatvalues = 0;
}

CogeometryDecorator::CogeometryDecorator(cogeometry content)
  :Cogeometry(content->Delta)
  ,original(content)
{
  floatvalues = content->floatValues();
}

void Cogeometry::setDelta(cogFloat delta)
{
  Delta=delta;
  ShiftDelta = 0.1*Delta;
  Delta2D = 32*Delta;
  Delta3D = 128*Delta;
  singleIntersectionLimit2D = 10*Delta2D;
}



cogFloat cogTetrahedron::volume() const
{
  cogFloat x[3],y[3],z[3]; int i;
  for(i=0;i<3;i++){
    x[i] = c1[i]-c0[i];
    y[i] = c2[i]-c0[i];
    z[i] = c3[i]-c0[i];
  }
  return 
    x[0]*(y[1]*z[2]-y[2]*z[1])
    + x[1]*(y[2]*z[0]-y[0]*z[2])
    + x[2]*(y[0]*z[1]-y[1]*z[0]);
}

cogFloat cogTetrahedron::diameter() const
{
  cogFloat d=0;
  cogIndex i;
  for(i=0;i<3;i++){
    if(c3[i]>c0[i]){
      if(c2[i]>c3[i]){
	if(c1[i]>c2[i])		{d += c1[i]-c0[i];}
	else if(c1[i]>c0[i]) 	{d += c2[i]-c0[i];}
	else 			{d += c2[i]-c1[i];}
      }else if(c2[i]<c0[i]){
	if(c1[i]>c3[i])		{d += c1[i]-c2[i];}
	else if(c1[i]>c2[i]) 	{d += c3[i]-c2[i];}
	else 			{d += c3[i]-c1[i];}
      }else{
	if(c1[i]>c3[i])		{d += c1[i]-c0[i];}
	else if(c1[i]>c0[i]) 	{d += c3[i]-c0[i];}
	else 			{d += c3[i]-c1[i];}
      }
    }else{
      if(c2[i]>c0[i]){
	if(c1[i]>c2[i])		{d += c1[i]-c3[i];}
	else if(c1[i]>c3[i]) 	{d += c2[i]-c3[i];}
	else 			{d += c2[i]-c1[i];}
      }else if(c2[i]<c3[i]){
	if(c1[i]>c0[i])		{d += c1[i]-c2[i];}
	else if(c1[i]>c2[i]) 	{d += c0[i]-c2[i];}
	else 			{d += c0[i]-c1[i];}
      }else{
	if(c1[i]>c0[i])		{d += c1[i]-c3[i];}
	else if(c1[i]>c3[i]) 	{d += c0[i]-c3[i];}
	else 			{d += c0[i]-c1[i];}
      }
    }
  }
  return d;
}

void cogLine::barycentre(cogPoint& p) const
{
  switch(c0.spaceDimension()){
  case 3:
    p[2] = (c0[2]==c1[2])? c0[2]: (c0[2]+c1[2])/2;
  case 2:
    p[1] = (c0[1]==c1[1])? c0[1]: (c0[1]+c1[1])/2;
  case 1:
    p[0] = (c0[0]==c1[0])? c0[0]: (c0[0]+c1[0])/2;
  }
}

void cogTriangle::barycentre(cogPoint& p) const
{
 p[0] = (c0[0]+c1[0]+c2[0])/3;
 p[1] = (c0[1]+c1[1]+c2[1])/3;
 if(c0.spaceDimension()>2) p[2] = (c0[2]+c1[2]+c2[2])/3;
}

void cogTetrahedron::barycentre(cogPoint& p) const
{
 p[0] = (c0[0]+c1[0]+c2[0]+c3[0])/4;
 p[1] = (c0[1]+c1[1]+c2[1]+c3[1])/4;
 p[2] = (c0[2]+c1[2]+c2[2]+c3[2])/4;
}