cogengrid.cxx

Delaunay三角形的网格剖分程序

#include "cogengrid.hxx"

#include "wzwarning.hxx"
static wzCounter cogInfoPointCalls(wzCounter::report_severity,
	      "CogenGrid::Point called","CogenGrid::Point calls: <>\n");
static wzCounter cogInfoPointCallCells(wzCounter::report_severity,
	      "CogenGrid::Point cell tested","CogenGrid::Point cells tested: <>\n");
static wzCounter cogInfoLineCalls(wzCounter::report_severity,
	      "CogenGrid::line called","CogenGrid::line calls: <>\n");
static wzCounter cogInfoFaceNotFound(wzCounter::report_severity,
	      "CogenGrid::face not found","CogenGrid::faces not found: <>\n");
static wzCounter cogInfoReverse(wzCounter::report_severity,
	      "CogenGrid::reverse used","CogenGrid::reverse used: <>\n");
static wzCounter cogInfoReverse2(wzCounter::report_severity,
	      "CogenGrid::reverse 2 used","CogenGrid::reverse 2 used: <>\n");

const wzIndex wzCellOfPoint=0;
static int sig[2] = {1,-1};

CogenGrid::CogenGrid(wzgrid g):grid(g),lastCellUsed(0)
{
  setDelta(1.e-4);
  floatvalues=grid->floatValues();
  addBox(grid->cgmin(0),grid->cgmax(0),
	 grid->cgmin(1),grid->cgmax(1),
	 grid->cgmin(2),grid->cgmax(2));
}

cogenGrid CogenGrid::cogengrid()
{return this;}

void CogenGrid::setDelta(cogFloat delta)
{
  Cogeometry::setDelta(delta);
  epsilon1 = delta;
  epsilon2 = epsilon1*epsilon1/2.0;
  epsilon3 = epsilon1*epsilon2/3.0;
}

wzIndex CogenGrid::Point(wzPoint& p0) const
{
  int *cn0;
  int c0,c1,ch,cv0,nl,i,i0,sl,co;
  wzCellType t,th;
  wzFloat *x= p0.X;
  double vv[4],xx,yy,zz,x1,x2,x3,y1,y2,y3,z1,z2,z3;
  wzFloat* ff[IBG2CNMAX];
  register wzFloat *fj;
  cogInfoPointCalls++;
  xx = x[0]; yy = x[1]; zz = x[2];
#define cvmax 3
  cv0=0;
  if(Startpoint == 0){
    if(lastCellUsed == 0){
      wzRangeLoop(grid->cells,c0){
	if(grid->ccdef(c0)) break;
      }
    }else{
      c0 = lastCellUsed;
    }
  }else{
    c0=Startpoint->I[wzCellOfPoint]; setStartpoint(0);
  }
  switch(wzCellTypeCODIM(t=grid->cct(c0))){
    /*
      case wzIsNode:
      co=wzGridCOutFirst(*gg,c0,t);
      c0=wzGridCOutside(*gg,co);
      t =cct(c0);       
      case wzIsEdge:
      co=wzGridCOutFirst(*gg,c0,t);
      c0=wzGridCOutside(*gg,co);
      t =cct(c0);
      */
  case wzIsFace:
    co=grid->ccleft(c0,t);
    if(grid->ccdef(co))	c0=co;
    else		c0=grid->ccright(c0,t);
    t = grid->cct(c0);
  }
gridRegionBegin:
  cogInfoPointCallCells++;
  wzAssert(grid->ccdef(c0));
  nl=wzCellTypePoints(t=grid->cct(c0)); cn0=grid->ccn(c0);
  for(i=0;i<nl;i++) ff[i]=grid->cnx(cn0[i]);
  sl=wzCellTypeSides(t);
  wzAssert(wzCellTypeCODIM(t)==wzIsRegion);
  i0 = 0;
gridRegionContinue:
  switch(wzCellTypeGDIM(t)){
  case 1:
    for(i=i0;i<sl;i++){
      if((vv[i]=sig[i]*(xx-ff[i][0]))<0) {goto gridRegionNext;}
    } break;
  case 2:
    for(i=i0;i<sl;i++){
      fj = ff[wzCellTypeEPoint1(t,i)];
      x1 = fj[0] - xx; y1 = fj[1] - yy;
      fj = ff[wzCellTypeEPoint2(t,i)];
      x2 = fj[0] - xx; y2 = fj[1] - yy;
      vv[i] = x1 * y2 - x2 * y1;
      if((vv[i]<-epsilon2)||((vv[i]<0)&&(cv0++<cvmax))) goto gridRegionNext;
    } break;
  case 3:
    for(i=i0;i<sl;i++){
      fj = ff[wzCellTypeSPoint1(t,i,0)];
      x1 = fj[0] - xx; y1 = fj[1] - yy; z1 = fj[2] - zz;
      fj = ff[wzCellTypeSPoint1(t,i,1)];
      x2 = fj[0] - xx; y2 = fj[1] - yy; z2 = fj[2] - zz;
      fj = ff[wzCellTypeSPoint1(t,i,2)];
      x3 = fj[0] - xx; y3 = fj[1] - yy; z3 = fj[2] - zz;
      vv[i] = -(x1*(y2*z3-y3*z2) + x2*(y3*z1-y1*z3) + x3*(y1*z2-y2*z1));
      if((vv[i]<-epsilon3)||((vv[i]<0)&&(cv0++<cvmax))) goto gridRegionNext;
    } break;
  }
  goto gridRegionEnd;
gridRegionNext:
  c1=grid->ccs(c0)[i];
  if(grid->ccund(c1)){
    i0 = i+1; goto gridRegionContinue;
  }else if(wzCellTypeCODIM(th=grid->cct(c1))==wzIsFace){
    if(c0==(ch=grid->ccright(c1,th)))	c1=grid->ccleft(c1,th);
    else					c1=ch;
    if(grid->ccund(c1)){
      i0 = i+1; goto gridRegionContinue;
    }
  }
  c0=c1;
  goto gridRegionBegin;
gridRegionEnd:
  if(floatvalues){
    wzFloat vvs=0,fv;
    for(i=0;i<sl;i++) vvs  += vv[i];
    for(i=0;i<sl;i++) vv[i] /= vvs;
    cn0=grid->ccn(c0);
    for(wzIndex fn=0;fn<floatvalues;fn++){
      fv=0;
      for(i=0;i<sl;i++) fv += vv[i]*(grid->cnf(cn0[i]))[fn];
      p0.f(fn) = fv;
    }
  }
  p0.S = (wzRegion) grid->ccu(c0);
  p0.I[wzCellOfPoint] = c0;
  ((CogenGrid*)this)->lastCellUsed = c0;
  return cogRCRegionFound;
}

static int suppress_unreachable_warning = 1;
wzIndex	CogenGrid::Line (cogFlag1& f, const cogLine& s) const
{
  int c0,c1,ch,c2,cb,cr,nl,i,ii,io,in,sl,*cn0,*cs0,mark[IBG2CSMAX];
  wzCellType t,th;
  wzRegion u,ur;
  wzFace ub;
  double vv0,vv1,vv2,vv3,vvv,p1,p2;
  double xb,xe,xbo,xeo,x1,xx0,xx1,xx2,xx3;
  double yb,ye,ybo,yeo,y1,yy0,yy1,yy2,yy3;
  double zb,ze,zbo,zeo,z1,zz0,zz1,zz2,zz3;
  wzFloat*	ff[IBG2CNMAX];
  register wzFloat *fj,dd;
  if(suppress_unreachable_warning) return Cogeometry::Line(f,s);
  cogInfoLineCalls++;
  int reverse=0;  // is 1 only if no boundary has been found starting from s.c0
  c0 = s.c0.I[wzCellOfPoint]; 
  c2 = s.c1.I[wzCellOfPoint];
  xbo = xb = s.c0[0];
  ybo = yb = s.c0[1];
  zbo = zb = s.c0[2];
  xeo = xe = s.c1[0];
  yeo = ye = s.c1[1];
  zeo = ze = s.c1[2];
  xx0 = xb-xe; yy0 = yb-ye; zz0 = zb-ze;
  u = grid->ccu(c0);
begin:
  wzAssert(grid->ccdef(c0));
  t=grid->cct(c0);  
  nl=wzCellTypePoints(t); cn0=grid->ccn(c0);
  sl=wzCellTypeSides (t); cs0=grid->ccs(c0);
  wzAssert(wzCellTypeCODIM(t)==wzIsRegion);
  for(i=0;i<nl;i++) ff[i]=grid->cnx(cn0[i]);
  // initializing the mark field: impossible neighbours marked.
  for(i=0;i<sl;i++){	
    if((ch=cs0[i])>0){
      switch(wzCellTypeCODIM(th=grid->cct(ch))){
      case wzIsFace:
	if((c1=grid->ccleft(ch,th))<=0){
	  {mark[i]=1;break;}
	}else if(c1==c0){
	  if((c1=grid->ccright(ch,th))<=0)
	    {mark[i]=1;break;}
	}
	if(grid->ccu(c1)<=0)
	  {mark[i]=1;break;}
      case wzIsRegion:
      default:
	mark[i]=0; break;
      }
    }else{
      // ???     mark[i]=0;
      mark[i]=1;
    }
  }
  //  now we want to find the optimal neighbour side ii and go to
  //  side_found. It should not be a marked side, the end node should be
  //  on the correct side, and if there remains some choice after this
  //  it is the side which intersects the line (xb,xe):
  switch(wzCellTypeGDIM(t)){
  case 1:
    if((vv0=(xe - ff[0][0]))<0) {ii=0; goto side_found;}
    if((vv0=(ff[1][0] - xe))<0) {ii=1; goto side_found;}
    goto end_is_inside;
  case 2:
    for(i=0;i<sl;i++){
      ii=i; io= -1;
    begin_2:
      fj  = ff[wzCellTypeSPoint1(t,ii,0)];
      xx1 = fj[0] - xe; yy1 = fj[1] - ye;
      fj  = ff[wzCellTypeSPoint1(t,ii,1)];
      xx2 = fj[0] - xe; yy2 = fj[1] - ye;
      vv0 = xx1*yy2 - xx2*yy1;
      if(vv0>=0) {
	mark[ii]=1; 
	if(io>-1) {ii=io; io= -1; goto begin_2;}
	else continue;
      }
      vv1 = xx0*yy2 - xx2*yy0;
      if(vv1>0) {
	in=wzCellTypeSSide(t,ii,1);
	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begin_2;}
      }
      vv2 = xx1*yy0 - xx0*yy1;
      if(vv2>0) {
	in=wzCellTypeSSide(t,ii,0);
	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begin_2;}
      }
      goto side_found;
    } 
    goto end_is_inside;
  case 3:
    for(i=0;i<sl;i++){
      ii=i; io= -1;
    begin_3:
      fj  = ff[wzCellTypeSPoint1(t,ii,0)];
      xx1 = fj[0] - xe; yy1 = fj[1] - ye; zz1 = fj[2] - ze;
      fj  = ff[wzCellTypeSPoint1(t,ii,1)];
      xx2 = fj[0] - xe; yy2 = fj[1] - ye; zz2 = fj[2] - ze;
      fj  = ff[wzCellTypeSPoint1(t,ii,2)];
      xx3 = fj[0] - xe; yy3 = fj[1] - ye; zz3 = fj[2] - ze;
      vv0 = xx1*(yy2*zz3-yy3*zz2)
	+   xx2*(yy3*zz1-yy1*zz3)
	+   xx3*(yy1*zz2-yy2*zz1);
      if(vv0<=0) {
	mark[ii]=1; 
	if(io>-1) {ii=io; io= -1; goto begin_3;}
	else continue;
      }
      vv1 = xx0*(yy2*zz3-yy3*zz2)
	+   xx2*(yy3*zz0-yy0*zz3)
	+   xx3*(yy0*zz2-yy2*zz0);
      if(vv1<0) {
	in=wzCellTypeSSide(t,ii,1);
	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begin_3;}
      }
      vv2 = xx1*(yy0*zz3-yy3*zz0)
	+   xx0*(yy3*zz1-yy1*zz3)
	+   xx3*(yy1*zz0-yy0*zz1);
      if(vv2<0) {
	in=wzCellTypeSSide(t,ii,2);
	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begin_3;}
      }
      vv3 = xx1*(yy2*zz0-yy0*zz2)
	+   xx2*(yy0*zz1-yy1*zz0)
	+   xx0*(yy1*zz2-yy2*zz1);
      if(vv3<0) {
	in=wzCellTypeSSide(t,ii,0);
	if(mark[in]==0) {mark[ii]=2; io=ii; ii=in; goto begin_3;}
      }
      goto side_found;
    } 
    goto end_is_inside;

  }
side_found:
  cb=grid->ccs(c0)[ii]; t=grid->cct(cb);
  if(wzCellTypeCODIM(t)==wzIsFace) {
    if(reverse){
      if(c0==(cr=grid->ccleft(cb,t))) cr=grid->ccright(cb,t);
      else wzAssert(c0==grid->ccright(cb,t));
      if(grid->ccu(cr)!=grid->ccu(c2)) {wzAssert(0); c0=cr; goto begin;}
    }
    goto intersection_found;
  }
  wzAssert(wzCellTypeCODIM(t)==wzIsRegion);
  if(cb==c2) {
    cogInfoFaceNotFound++;
    return cogRCFaceNotFound;
  }
  if(u!=grid->ccu(cr)){
    cr = cb; cb = 0; wzAssert(0); goto intersection_found;
  }
  c0=cb; goto begin;
end_is_inside:
  // xe is inside of c0 or behind an outer boundary of c0.
  // That's ok if the regions of c0 and c2 are identical:
  if(grid->ccu(c0)==grid->ccu(c2)){
    cogInfoFaceNotFound++;
    return cogRCFaceNotFound;
  }
  // We have not found any intersection, but the ends are in different regions.
  // Let's try the other direction.
  if(!reverse){
    cogInfoReverse++;
    reverse = 1;
    c0 = s.c1.I[wzCellOfPoint];
    c2 = s.c0.I[wzCellOfPoint];
    xbo = xb = s.c1[0];
    ybo = yb = s.c1[1];
    zbo = zb = s.c1[2];
    xeo = xe = s.c0[0];
    yeo = ye = s.c0[1];
    zeo = ze = s.c0[2];
    xx0 = xb-xe; yy0 = yb-ye; zz0 = zb-ze;
    u = grid->ccu(c0);
    goto begin;
  }
  // Seems, the whole line is on the boundary.  We choose the starting
  // point as the boundary point, but create an artificial boundary.
  // wzAssert(0);	// not yet implemented
  cogInfoReverse2++;
  ub = 1; 
  u  = grid->ccu(c2);
  ur = grid->ccu(c0);
  vv0= s.diameter(); vvv = -Delta;
  goto return_intersection;
intersection_found:
  ub = grid->ccu(cb);
  ur = grid->ccu(cr);
  switch(wzCellTypeGDIM(t)){
  case 1:
    if(ii)	vvv = xb - ff[0][0];
    else	vvv = ff[1][0] - xb;
    if(vvv<0) vvv=0;
    break;
  case 2:
    vvv	= (xx1-xx0)*(yy2-yy0) - (xx2-xx0)*(yy1-yy0);
    if(vvv<0) vvv=0;
    break;
  case 3:
    yy1 -= yy0; yy2 -= yy0; yy3 -= yy0;
    zz1 -= zz0; zz2 -= zz0; zz3 -= zz0;
    vvv	= (xx1-xx0)*(yy2*zz3-yy3*zz2)
      +   (xx2-xx0)*(yy3*zz1-yy1*zz3)
      +   (xx3-xx0)*(yy1*zz2-yy2*zz1);
    if(vvv>0) vvv=0;
    break;
  }
return_intersection:
  p1=vv0/(vv0-vvv); p2=vvv/(vvv-vv0); // correct if reverse???
  f.p1[0] = x1 = p1*xbo + p2*xeo;
  f.p1[1] = y1 = p1*ybo + p2*yeo;
  f.p1[2] = z1 = p1*zbo + p2*zeo;
  f.p1.I[wzCellOfPoint]	= c0;
  dd = 0.5*Delta/s.diameter();
  f.p0[0] = x1 + dd*(xbo-xeo);
  f.p0[1] = y1 + dd*(ybo-yeo);
  f.p0[2] = z1 + dd*(zbo-zeo);
  f.po[0] = x1 - dd*(xbo-xeo);
  f.po[1] = y1 - dd*(ybo-yeo);
  f.po[2] = z1 - dd*(zbo-zeo);
  f.p1.S = wzFace(ub);
  f.p0.S = wzFace(u);
  f.po.S = wzFace(ur);
#ifdef wzFloatValues
 cn0=ccn(cb); sl=wzCellTypePoints(t);
 switch(wzCellTypeGDIM(t)){
  case 1: vv[0]=1.0;
	break;
  case 2: x = p1*xb + p2*xe; y = p1*yb + p2*ye;
	f0=gnx(gg,cn0[0]); f1=gnx(gg,cn0[1]);
	xx0=f1[0]-f0[0]; yy0=f1[1]-f0[1];
	vv[0]=(f1[0]-x)*xx0 + (f1[1]-y)*yy0;
	vv[1]=(x-f0[0])*xx0 + (y-f0[1])*yy0;
	vvs  = vv[0]+vv[1];
	vv[0] /= vvs; vv[1] /= vvs;
	break;
  case 3:  x = p1*xb + p2*xe; y = p1*yb + p2*ye; z = p1*zb + p2*ze;
	f0=gnx(gg,cn0[0]); f1=gnx(gg,cn0[1]); f2=gnx(gg,cn0[2]);
	xx0=(f1[1]-f0[1])*(f2[2]-f0[2]) - (f1[2]-f0[2])*(f2[1]-f0[1]);
	yy0=(f1[2]-f0[2])*(f2[0]-f0[0]) - (f1[0]-f0[0])*(f2[2]-f0[2]);
	zz0=(f1[0]-f0[0])*(f2[1]-f0[1]) - (f1[1]-f0[1])*(f2[0]-f0[0]);
	vv[0] = (f1[0]-x)*(yy0*(f2[2]-z) - zz0*(f2[1]-y))
	      + (f1[1]-y)*(zz0*(f2[0]-x) - xx0*(f2[2]-z))	
	      + (f1[2]-z)*(xx0*(f2[1]-y) - yy0*(f2[0]-x));	
	vv[1] = (f2[0]-x)*(yy0*(f0[2]-z) - zz0*(f0[1]-y))
	      + (f2[1]-y)*(zz0*(f0[0]-x) - xx0*(f0[2]-z))	
	      + (f2[2]-z)*(xx0*(f0[1]-y) - yy0*(f0[0]-x));	
	vv[2] = (f0[0]-x)*(yy0*(f1[2]-z) - zz0*(f1[1]-y))
	      + (f0[1]-y)*(zz0*(f1[0]-x) - xx0*(f1[2]-z))	
	      + (f0[2]-z)*(xx0*(f1[1]-y) - yy0*(f1[0]-x));	
	vvs   = vv[0]+vv[1]+vv[2];
	vv[0] /= vvs; vv[1] /= vvs; vv[2] /= vvs;
	break;
 }
 for(f=0;f<grdf->fanz;f++){
	fn=grdf->flist[f];
	fv=0;
	for(i=0;i<sl;i++) fv += vv[i]*(ibg2pF(ccn(cn0[i]))[fn]);
	ibg2pF(*nint)[fn]=fv;
 }
#endif
 return cogRCFaceFound;
}