cognavigate.cxx

Delaunay三角形的网格剖分程序

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

// in this file, the functions for navigation among the flags of the
// same intersection point are implemented:
/*
void	Cogeometry::otherSide(cogFlag1& f, const cogFlag1& o) const;
void	Cogeometry::otherSide(cogFlag2& f, const cogFlag2& o) const;
void	Cogeometry::otherSide(cogFlag3& f, const cogFlag3& o) const;
void	Cogeometry::nextFace (cogFlag2& f, const cogFlag2& o) const;
void	Cogeometry::otherLine(cogFlag3& f, const cogFlag3& o) const;
*/
void	Cogeometry::otherSide(cogFlag1& f, const cogFlag1& o) const
{
  f.p1 = o.p1;
  f.p0 = o.po;
  f.po = o.p0;
}
void	Cogeometry::otherSide(cogFlag2& f, const cogFlag2& o) const
{
  f.p2 = o.p2;
  f.p1 = o.p1;
  f.p0 = o.po;
  f.po = o.p0;
}
void	Cogeometry::otherSide(cogFlag3& f, const cogFlag3& o) const
{
  f.p3 = o.p3;
  f.p2 = o.p2;
  f.p1 = o.p1;
  f.p0 = o.po;
  f.po = o.p0;
}

void	Cogeometry::nextFace (cogFlag3& f, const cogFlag3& o) const
{
  cogFlag2 o2(o.p2,o.p1,o.p0,o.po);
  nextFace(f,o2);
  f.p3 = o.p3; 
}

void   Cogeometry::nextFace(cogFlag2& f, const cogFlag2& o) const
{
  cogIndex i,rc;
  cogFloat x[3],y[3],z[3],u[3],l;
  f.p2 = o.p2;
  for(i=0;i<3;i++){
    x[i] = o.p2[i]-o.p1[i];
    y[i] = o.po[i]-o.p1[i];
  }
  z[0] = x[1]*y[2]-x[2]*y[1];
  z[1] = x[2]*y[0]-x[0]*y[2];
  z[2] = x[0]*y[1]-x[1]*y[0];
  l = 1.1*sqrt(z[0]*z[0]+z[1]*z[1]+z[2]*z[2]);
  u[0] = (z[1]*x[2]-z[2]*x[1])/l;
  u[1] = (z[2]*x[0]-z[0]*x[2])/l;
  u[2] = (z[0]*x[1]-z[1]*x[0])/l;
  wzAssert(u[0]*y[0]+u[1]*y[1]+u[2]*y[2]>0);
  wzPoint po(o.p2[0]+x[0],o.p2[1]+x[1],o.p2[2]+x[2]);
  wzPoint p0(o.p1[0]+u[0],o.p1[1]+u[1],o.p1[2]+u[2]);
  wzPoint p1(  po[0]+u[0],  po[1]+u[1],  po[2]+u[2]);
  wzPoint p2(  po[0]-u[0],  po[1]-u[1],  po[2]-u[2]);
  wzPoint p3(o.p1[0]-u[0],o.p1[1]-u[1],o.p1[2]-u[2]);
  Point(p0); rc = Line(f,cogLine(o.po,p0));	if(rc==cogRCFaceFound) return;
  Point(p1); rc = Line(f,cogLine(p0,p1));	if(rc==cogRCFaceFound) return;
  Point(p2); rc = Line(f,cogLine(p1,p2));	if(rc==cogRCFaceFound) return;
  Point(p3); rc = Line(f,cogLine(p2,p3));	if(rc==cogRCFaceFound) return;
  rc = Line(f,cogLine(p3,o.p0));		if(rc==cogRCFaceFound) return;
  wzAssert(0);
}

void	Cogeometry::otherLine(cogFlag3& f, const cogFlag3& o) const
{
  // not yet implemented
  wzAssert(0);
}