cog1d.cxx

Delaunay三角形的网格剖分程序

#include "cog1d.hxx"

void Cogeometry1D::setLast(cogIndex l) const
{
 ((Cogeometry1D *) this)->last = l;
}

Cogeometry1D::Cogeometry1D(cogFloat xmin, cogFloat xmax, 
		     cogIndex n,
		     cogRegion first,
		     cogRegion outside)
:last(0)
{cogFace f0; cogIndex i,r0 = first;
 num=0; r(0)=outside;
 if(n==0){after(xmin,first,f0); return;}
 cogFloat delta = (xmax-xmin)/n;
 if(delta > 0){
	for(i=0;i<n;i++){after(xmin+i*delta,cogRegion(r0+i),f0);}
	after(xmax,outside,f0);
 }else if(delta<0){
	for(i=0;i<n;i++){before(xmin+i*delta,cogRegion(r0+i),f0);}
	before(xmax,outside,f0);
 }else{
	before(xmax,first,f0);
 }
}

Cogeometry1D::Cogeometry1D(cogRegion r0,
	cogFloat x1, cogRegion r1,
        cogFloat x2, cogRegion r2,
	cogFloat x3, cogRegion r3,
	cogFloat x4, cogRegion r4,
	cogFloat x5, cogRegion r5,
	cogFloat x6, cogRegion r6,
	cogFloat x7, cogRegion r7,
	cogFloat x8, cogRegion r8,
	cogFloat x9, cogRegion r9)
:last(0)
{cogFace f0;
 num=0; r(0)=r0;
 if(x1==cogNoFloat) {return;} after(x1,r1,f0);
 if(x2==cogNoFloat) {return;} after(x2,r2,f0);
 if(x3==cogNoFloat) {return;} after(x3,r3,f0);
 if(x4==cogNoFloat) {return;} after(x4,r4,f0);
 if(x5==cogNoFloat) {return;} after(x5,r5,f0);
 if(x6==cogNoFloat) {return;} after(x6,r6,f0);
 if(x7==cogNoFloat) {return;} after(x7,r7,f0);
 if(x8==cogNoFloat) {return;} after(x8,r8,f0);
 if(x9==cogNoFloat) {return;} after(x9,r9,f0);
}

Cogeometry1D::Cogeometry1D(cogRegion r0,
        cogFloat x1, cogRegion r1, cogFace f1,
	cogFloat x2, cogRegion r2, cogFace f2,
	cogFloat x3, cogRegion r3, cogFace f3,
	cogFloat x4, cogRegion r4, cogFace f4,
	cogFloat x5, cogRegion r5, cogFace f5,
	cogFloat x6, cogRegion r6, cogFace f6,
	cogFloat x7, cogRegion r7, cogFace f7,
	cogFloat x8, cogRegion r8, cogFace f8,
	cogFloat x9, cogRegion r9, cogFace f9)
:last(0)
{num=0; r(0)=r0;
 if(x1==cogNoFloat) {return;} after(x1,r1,f1);
 if(x2==cogNoFloat) {return;} after(x2,r2,f2);
 if(x3==cogNoFloat) {return;} after(x3,r3,f3);
 if(x4==cogNoFloat) {return;} after(x4,r4,f4);
 if(x5==cogNoFloat) {return;} after(x5,r5,f5);
 if(x6==cogNoFloat) {return;} after(x6,r6,f6);
 if(x7==cogNoFloat) {return;} after(x7,r7,f7);
 if(x8==cogNoFloat) {return;} after(x8,r8,f8);
 if(x9==cogNoFloat) {return;} after(x9,r9,f9);
}

Cogeometry1D::Cogeometry1D(cogIndex n, cogFloat *xx, cogRegion* rr)
:num(0),last(0)
{int i; cogFace def;
 r(0)=rr[num];
 for(i=0;i<n;i++){
	before(xx[i],rr[i],def);
 }
}

Cogeometry1D::Cogeometry1D(cogIndex n, cogFloat *xx,
        			cogRegion* rr, cogFace *ff)
:num(0),last(0)
{int i;
 r(0)=rr[num];
 for(i=0;i<n;i++){
        before(xx[i],rr[i],ff[i]);
 }
}

void Cogeometry1D::before(cogFloat xx, cogRegion b, cogFace face)
{int i;
 r(num+1)=r[num];
 if(num==0 || xx>x[num-1]){x(num)=xx;r[num]=b;f(num)=face;num++;return;}
 for(i=num-1;i>0;i--){
        r[i+1] = r[i];
        if(xx>x[i-1]){x[i]=xx; r[i]=b; f[i]=face; return;}
        x[i]=x[i-1]; f[i]=f[i-1];
 }
 x[0]=xx;r[1]=r[0];r[0]=b;f[0]=face;
}

void Cogeometry1D::after(cogFloat xx, cogRegion a, cogFace face)
{int i;
 r(num+1)=r[num];
 if(num==0 || xx>x[num-1]){x(num)=xx;r[num+1]=a;f(num)=face;num++;return;}
 for(i=num-1;i>0;i--){
        r[i+1] = r[i];
        if(xx>x[i-1]){x[i]=xx; r[i+1]=a; f[i]=face; return;}
        x[i]=x[i-1]; f[i]=f[i-1];
 }
 x[0]=xx;r[1]=a;f[0]=face;
}

cogIndex   Cogeometry1D::Point(cogPoint& s) const
{int i=last;
 cogFloat xx=s[0];
 while(xx>=x[i]){
        i++;
        if(i==num){
        	setLast(num-1); s.segment()=r[num]; return cogRCRegionFound;
        }
 } // now xx<x[i];
 while(i>0){
        i--;
        if(xx>=x[i]){
                setLast(i); s.segment()=r[i+1]; return cogRCRegionFound;
        }
 }
 setLast(0); s.segment()=r[0]; return cogRCRegionFound;
}

cogIndex   Cogeometry1D::Line (cogFlag1& f1, const cogLine& s) const
{cogFloat x0=s.c0[0],x1=s.c1[0];
 cogIndex li; cogFloat xi;
 cogPoint p0(x0); Point(p0); // set last so that x[last]<=x0<x[last+1] or x0<x[0];
 if(x0<x1){
   if(x0<x[0]){
     li = 0;
   }else{
     li = last+1;
   }
   if(li==num || x1<x[li]-Delta){ return cogRCFaceNotFound;}
   if(x1<x[li]+Delta){
     if(s.c0.segment()==s.c1.segment()) return cogRCFaceNotFound;
   }
   f1.p1.affine(((xi= x[li])-x0)/(x1-x0),s.c1,s.c0); f1.po = f1.p0 = f1.p1;
   f1.p1[0] = xi;
   f1.p0[0] = xi-ShiftDelta;
   f1.po[0] = xi+ShiftDelta;
   f1.p1.segment() = f[li];
   f1.p0.segment() = r[li];
   f1.po.segment() = r[li+1];
 }else if(x1<x0){
   li = last;
   if(li==0 || x1>x[li]+Delta){ return cogRCFaceNotFound;}
   if(x1>x[li]-Delta){
     if(s.c0.segment()==s.c1.segment()) return cogRCFaceNotFound;
   }
   f1.p1.affine(((xi= x[li])-x0)/(x1-x0),s.c1,s.c0); f1.po = f1.p0 = f1.p1;
   f1.p1[0] = xi;
   f1.p0[0] = xi+ShiftDelta;
   f1.po[0] = xi-ShiftDelta;
   f1.p1.segment() = f[li];
   f1.p0.segment() = r[li+1];
   f1.po.segment() = r[li];
 }else{
   return cogRCFaceNotFound;
 }
 return cogRCFaceFound;
}