cogview.cxx

Delaunay三角形的网格剖分程序

#include <math.h>
#include <stdio.h>
extern "C"{
#include <GL/gl.h>
#include "gltools/glwin.h"
}

#include "wzoutput.hxx"
#include "wzview.hxx"

static void cogViewCallback(wzWindow win, void *data);
static void DrawTriangle(double *p1, double *p2, double *p3);
static cogeometry geom;
static wzFloat *sh  = wzCurrentData::shift;
static wzFloat *fac = wzCurrentData::factor;

void cogeometryView(wzWindow win, wzIndex i, 
		    cogeometry g, char* name, char* title)
{
  wzwDefineObject(win,i,
		  name[0],name,title,"cogview.html#Keys",
		  (wzwCallback)cogViewCallback,*(void**)&g);
}

static void DrawTriangle(double *p1, double *p2, double *p3)
{
  const int on12=0,non12=3,non23=4,non13=5,n123=6,nno=7;
  int i,face;
  wzColor col;
  cogPoint x1(sh[0]+p1[0]*fac[0],sh[1]+p1[1]*fac[1],sh[2]+p1[2]*fac[2]);
  cogPoint x2(sh[0]+p2[0]*fac[0],sh[1]+p2[1]*fac[1],sh[2]+p2[2]*fac[2]);
  cogPoint x3(sh[0]+p3[0]*fac[0],sh[1]+p3[1]*fac[1],sh[2]+p3[2]*fac[2]);
  cogPoint c1,c2,c3,c0,co,mm;
  cogIndex rc;
  cogFlag1 f1(c1,c0,co),f2(c2,c0,co),f3(c3,c0,co);
  geom->Point(x1); geom->Point(x2); geom->Point(x3);
  face = nno;
  if(x1.segment() != x2.segment()){
    face= on12;
    rc=geom->Line(f3,cogLine(x1,x2));
    if(rc!=cogRCFaceFound) wzOutput::Default<<"not found\n";
    for(i=0;i<3;i++){c3[i] -= sh[i];  c3[i] /= fac[i];}
  }
  if(x1.segment() != x3.segment()){
    if(face==on12) {
      if (x2.segment() == x3.segment()) face = non23;
      else                              face = n123;
    }else 	face=non12;
    rc = geom->Line(f2,cogLine(x1,x3));
    if(rc!=cogRCFaceFound) wzOutput::Default<<"not found\n";
    for(i=0;i<3;i++){c2[i] -= sh[i];  c2[i] /= fac[i];}
  }
  if(x2.segment() != x3.segment()){
    if(face==on12) face=non13;
    rc=geom->Line(f1,cogLine(x3,x2));
    if(rc!=cogRCFaceFound) wzOutput::Default<<"not found\n";
    for(i=0;i<3;i++){c1[i] -= sh[i];  c1[i] /= fac[i];}
  }
  switch (face){
  case non13:
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p3[0],p3[1],p3[2]);
    glVertex3f(p1[0],p1[1],p1[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(p3[0],p3[1],p3[2]);
    col=x2.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(p2[0],p2[1],p2[2]);
    break;
  case non23:
    col=x2.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p2[0],p2[1],p2[2]);
    glVertex3f(p3[0],p3[1],p3[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    col=x2.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(p2[0],p2[1],p2[2]);
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(p1[0],p1[1],p1[2]);
    break;
  case non12:
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p1[0],p1[1],p1[2]);
    glVertex3f(p2[0],p2[1],p2[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(p1[0],p1[1],p1[2]);
    col=x3.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(p3[0],p3[1],p3[2]);
    break;
  case n123:
    mm[0] = (c1[0]+c2[0]+c3[0])/3;
    mm[1] = (c1[1]+c2[1]+c3[1])/3;
    mm[2] = (c1[2]+c2[2]+c3[2])/3;
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p1[0],p1[1],p1[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(mm[0],mm[1],mm[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    col=x2.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p2[0],p2[1],p2[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(mm[0],mm[1],mm[2]);
    glVertex3f(c3[0],c3[1],c3[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    col=x3.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p3[0],p3[1],p3[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(mm[0],mm[1],mm[2]);
    glVertex3f(c1[0],c1[1],c1[2]);
    glVertex3f(c2[0],c2[1],c2[2]);
    break;
  default:
    col=x1.segment().region().color();
    glColor3f(col.r()/255.0,col.g()/255.0,col.b()/255.0);
    glVertex3f(p1[0],p1[1],p1[2]);
    glVertex3f(p2[0],p2[1],p2[2]);
    glVertex3f(p3[0],p3[1],p3[2]);
    break;
  }
}

static void cogViewCallback(wzWindow win, void *data)
{
  double a[3],aa,x1[3],x2[3],x3[3],x4[3],xm[3],xp[3],b,c,d,f[3];
  int n,dy,dz,dir,hold,triang;
  geom = *(cogeometry *) &data;
  glBegin(GL_TRIANGLES);
  wzwGetVolume(win,&xm[0],&xp[0],&xm[1],&xp[1],&xm[2],&xp[2]);
  wzwGetPlane(win,&a[0],&a[1],&a[2],&d);
  wzwGetDrawLevel(win,&hold,&triang);
  n = (int) (sqrt((double)triang)/2); if(n<=2) n=2;
  if      (a[0]>0.5)  {dir=0;dy=1;dz=2;aa=a[0];}
  else if (a[1]>0.5)  {dir=1;dy=0;dz=2;aa=a[1];}
  else                {dir=2;dy=0;dz=1;aa=a[2];}
  f[dy]= (xp[dy]-xm[dy])/n;	// h[dy]=f[dy]/2;
  f[dz]= (xp[dz]-xm[dz])/n;	// h[dz]=f[dz]/2;
  b=a[dy]; c=a[dz]; aa *= xp[dir]-xm[dir];
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++){
      x1[dy]=xm[dy]+f[dy]*i;     x1[dz]=xm[dz]+f[dz]*j;
      x2[dy]=xm[dy]+f[dy]*(i+1); x2[dz]=xm[dz]+f[dz]*j;
      x3[dy]=xm[dy]+f[dy]*i;     x3[dz]=xm[dz]+f[dz]*(j+1);
      x4[dy]=xm[dy]+f[dy]*(i+1); x4[dz]=xm[dz]+f[dz]*(j+1);
      x1[dir]=-(d+b*x1[dy]+c*x1[dz]);
      x2[dir]=-(d+b*x2[dy]+c*x2[dz]);
      x3[dir]=-(d+b*x3[dy]+c*x3[dz]);
      x4[dir]=-(d+b*x4[dy]+c*x4[dz]);
      DrawTriangle(x1,x2,x3);
      DrawTriangle(x3,x2,x4);
    }
  glEnd();
  geom = 0;
}