wzdelaunay.cxx

Delaunay三角形的网格剖分程序

	  dscl[s][dsi] = -1;
	}else{
	  wzAssert(wzCellTypeSSize(t0,s)==2);
	  c1 = createCell(wzCellType2Triangle);
	  cni=ccn(c1); csi=ccs(c1);
	  cn0 = ccn(c0);
	  cs0 = ccs(c0);
	  cni[0]=n0; cnc(n0)=c1; cni++; cnj=cni;
	  csi[0]=cc; csi++;
	  for(i=0;i<2;i++){
	    cni[i]  = nt = cn0[wzCellTypeSPoint1(t0,s,i)];
	    cnc(nt) = c1; csi[i] = -1;
	  }
	  cx1=cnx(*cnj++); cx2=cnx(*cnj++);
	  x1=cx1[0]-x; y1=cx1[1]-y;
	  x2=cx2[0]-x; y2=cx2[1]-y;
	  l1=x1*x1 + y1*y1;
	  l2=x2*x2 + y2*y2;
	  det=x1*y2 - x2*y1;
	  if((det<ibgDelaunayEps2)){
	    destroyCell(c1);
	    // cc has to be deleted, but it's neighbours don't know this if they are
	    // deleted before (v-test is ok for cc). test deleted neighbours for this effect
	    if(cc>=0){
	      csh=ccs(cc);
	      for(sh=0;sh<wzCellTypeSides(t);sh++){
		if(!(dsh= -ccu(ch=csh[sh]))) continue;
		if(ch==c0) continue;
		if(dsh>dsi) continue;
		for(ih=0;ih<wzCellTypeSides(cct(ch));ih++){
		  if(ccs(ch)[ih]==cc){
		    cl=(dscl[ih])[dsh];
		    wzAssert(cl>0);
		    destroyCell(cl);
		    //		 for(nsi=0;nsi<nsl;nsi++){
		    //		   if(nsc[nsi]==cl) nsco[nsi]=0;
		    //		 }
		    (dscl[ih])[dsh] = -1;
		  }
		}
	      }
	      del=1; goto deltest;
	    }else{
	      wzAssert(0);	// do the same stuff for 3D too!!!
	      continue;
	    }
	  }
	  mx=(l1*y2 - l2*y1)/det;
	  my=(l2*x1 - l1*x2)/det;
	  ccv(c1).x1[0]=x;
	  ccv(c1).x1[1]=y;
	  ccv(c1).x2[0]=x+mx;
	  ccv(c1).x2[1]=y+my;
	  //	 cvert = ccv(c1);
	  // Eintragen von c1 in die nsc- und dsc-Listen:
	  nsi = nsl; nsl++;
	  nsc[nsi]=c1; nscs[nsi]=s; nsco[nsi]=c0;
	  ccu(c1)=0; (dscl[s])[dsi]=c1;
	}
      }
      dsi++;
    }
    // Herstellung der Verbindungen der neuen Zellen untereinander: 
    for(nsi=0;nsi<nsl;nsi++){
      c1 = nsc[nsi]; cs1= ccs(c1);	//* neue Zelle 
      c0 =nsco[nsi]; s0 = nscs[nsi];	//* alte Zelle 
      ds0= -ccu(c0);
      if((dscl[s0])[ds0]!=c1) continue;
      t0 =cct(c0); t1=cct(c1); cn0=ccn(c0); cs0=ccs(c0);
      // Eintragen des Zeigers von cc auf c1: 
      if((cc=cs1[0])>0){
	se=(ss=ccs(cc))+wzCellTypeSides(t);
	for(;ss<se;ss++) if(ss[0]==c0)
	  {ss[0]=c1; break;}
      }
      for(i0=0;i0<wzCellTypeSSize(t0,s0);i0++){
	if(cs1[wzCellTypeSSide(t1,0,i0)]>0) continue;
	nh=cn0[wzCellTypeSPoint1(t0,s0,i0)];
	ct=(dscl[wzCellTypeSSide(t0,s0,i0)])[ds0];
	if(ct<=0){
	  /* "Ubergang zur Nachbarzelle */
	  ch=c0;th=t0;sh=s0;ih=i0;csh=cs0;
	  do{
	    co =ch;
	    ch =csh[wzCellTypeSSide(th,sh,ih)];
	    dsh= -ccu(ch);	wzAssert(dsh>0);
	    th =cct(ch);
	    cnh=ccn(ch);
	    csh=ccs(ch);
	    for(sh=0;;sh++)	{
	      if(csh[sh]==co) break;
	      wzAssert(sh<4);}
	    for(ih=0;;ih++) {
	      if(cnh[wzCellTypeSPoint1(th,sh,ih)]==nh)
		break;
	      wzAssert(ih<2);}
	    ct=dscl[wzCellTypeSSide(th,sh,ih)][dsh];
	  }while(ct<=0);
	}
	cs1[wzCellTypeSSide(t1,0,i0)]=ct;
	cnj = ccn(ct); csj = ccs(ct); tt=cct(ct);
	for(j=0;j<4;j++) if(cnj[wzCellTypeSPoint1(tt,0,j)]==nh){
	  csj[wzCellTypeSSide(tt,0,j)]= c1; break;
	}
      }
    }
    /* Eintragen der zerst"orten Zellen in die free-Listen */
    for(dsi=1;dsi<=dsl;dsi++){
      cc=dsc[dsi]; destroyCell(cc);
    }
  }
  return wzSuccess;
}

int   wzDelaunayGenerator::Delaunay3(wzpoints list)
{
  int n,s,s0,sh,ms,i,i0,ih,j,del;
  int n0,nb,no,nold,nshift,nt,nh,c0,c1,cc,ct,ch,cl,co;
  int *cs0,*cs1,*csi,*csj,*csh,*cn0,*cni,*cnj,*cnh,*ss,*se;
  wzFloat *cx1,*cx2,*cx3;
  wzFloat x,y,z,x1,x2,x3,y1,y2,y3,z1,z2,z3;
  wzFloat l1,l2,l3,det,a1,a2,a3,mx,my,mz;
  wzCellType t,t0,t1,th,tt;
  wzPoint nd;
  ibgSphere v;
  int dsi,ds0,dsh,nsi;
  n0 = nb = Grid.points.append(list->length());
  nshift = nb - list->first(); n0--;
  nold = list->first()-1;
  wzRangeLoop(*list,n){
    wzAssert(n==nold+1);
    nold = n;
    n0++;
    Grid.Point[n0]=list->Point[n];
    // Suche nach einer Zelle, die n0 enth"alt: 
    no = list->Near[n] + nshift;
    wzAssert(list->Near[n]<n);
    if(no<nb) no = n0-1;
    c0=g3find(cnx(n0),cnc(no));
    // old code:
    x=cnx(n0)[0]; y=cnx(n0)[1]; z=cnx(n0)[2];
    nsl=0; dsl=0; dsi=1; dsc[++dsl]=c0; ccu(c0)= -dsl; 
    while(dsi<=dsl){
      c0  = dsc[dsi];
      cs0 = ccs(c0); ms=wzCellTypeSides((t0=cct(c0)));
      cn0 = ccn(c0);
      for(s=0;s<ms;s++){
	cc= cs0[s];
	if(cc<=0) {del=0; goto deltest;}
	if(ccu(cc)<0) {(dscl[s])[dsi]= -1; continue;}
	v = ccv(cc); t = cct(cc);
	if((     (v.x1[0]-x)*(v.x2[0]-x)
		 +    (v.x1[1]-y)*(v.x2[1]-y)
		 +    (v.x1[2]-z)*(v.x2[2]-z)) < ibgDelaunayEps2)	del=1;
	else							del=0;
      deltest:
	if(del){
	  dsc[++dsl] = cc; ccu(cc) = -dsl;
	  (dscl[s])[dsi] = -1;
	}else{
	  wzAssert(wzCellTypeSSize(t0,s)==3);
	  c1 = createCell(wzCellType3Tetrahedron);
	  cni=ccn(c1); csi=ccs(c1);
	  cn0 = ccn(c0);
	  cs0 = ccs(c0);
	  cni[0]=n0; cnc(n0)=c1; cni++; cnj=cni;
	  csi[0]=cc; csi++;
	  for(i=0;i<3;i++){
	    //	erfordert i+1 == wzCellTypeSPoint1(wzCellType3Tetrahedron,0,i)  
	    cni[i]  = nt = cn0[wzCellTypeSPoint1(t0,s,i)];
	    cnc(nt) = c1; csi[i] = -1;
	  }
	  // Berechnen und Eintragen der Daten der Umkugel: 
	  cx1=cnx(*cnj++); cx2=cnx(*cnj++); cx3=cnx(*cnj++);
	  x1=cx1[0]-x; y1=cx1[1]-y; z1=cx1[2]-z;
	  x2=cx2[0]-x; y2=cx2[1]-y; z2=cx2[2]-z;
	  x3=cx3[0]-x; y3=cx3[1]-y; z3=cx3[2]-z;
	  l1=x1*x1 + y1*y1 + z1*z1;
	  l2=x2*x2 + y2*y2 + z2*z2;
	  l3=x3*x3 + y3*y3 + z3*z3;
	  a1=y2*z3 - y3*z2;
	  a2=y3*z1 - y1*z3;
	  a3=y1*z2 - y2*z1;
	  det=x1*a1 + x2*a2 + x3*a3;
	  if((det>-ibgDelaunayEps3)){
	    destroyCell(c1);
	    // cc has to be deleted, but it's neighbours don't know this if they are
	    // deleted before (v-test is ok for cc). test deleted neighbours for this effect
	    csh=ccs(cc);
	    for(sh=0;sh<wzCellTypeSides(t);sh++){
	      if((ch=csh[sh])==0) continue;
	      if(!(dsh= -ccu(ch))) continue;
	      if(ch==c0) continue;
	      if(dsh>dsi) continue;
	      for(ih=0;ih<wzCellTypeSides(cct(ch));ih++){
		if(ccs(ch)[ih]==cc){
		  cl=(dscl[ih])[dsh];
		  wzAssert(cl>0);
		  destroyCell(cl);
		  (dscl[ih])[dsh] = -1;
		}
	      }
	    }
	    del=1;
	    goto deltest;
	  }
	  mx=(l1*a1+l2*a2+l3*a3)/det;
	  my=(l1*(z2*x3-z3*x2)+l2*(z3*x1-z1*x3)+l3*(z1*x2-z2*x1))/det;
	  mz=(l1*(x2*y3-x3*y2)+l2*(x3*y1-x1*y3)+l3*(x1*y2-x2*y1))/det;
	  ccv(c1).x1[0]=x;
	  ccv(c1).x1[1]=y;
	  ccv(c1).x1[2]=z;
	  ccv(c1).x2[0]=x+mx;
	  ccv(c1).x2[1]=y+my;
	  ccv(c1).x2[2]=z+mz;
	  // Eintragen von c1 in die nsc- und dsc-Listen: 
	  nsi = nsl; nsl++;
	  nsc[nsi]=c1; nscs[nsi]=s; nsco[nsi]=c0;
	  ccu(c1)=0; (dscl[s])[dsi]=c1;
	}
      }
      dsi++;
    }
    // Herstellung der Verbindungen der neuen Zellen untereinander: 
    for(nsi=0;nsi<nsl;nsi++){
      c1 = nsc[nsi]; cs1= ccs(c1);	/* neue Zelle */
      c0 =nsco[nsi]; s0 = nscs[nsi];	/* alte Zelle */
      ds0= -ccu(c0);
      if((dscl[s0])[ds0]!=c1) continue;
      t0 =cct(c0); t1=cct(c1); cn0=ccn(c0); cs0=ccs(c0);
      // Eintragen des Zeigers von cc auf c1: 
      if((cc=cs1[0])>0){
	se=(ss=ccs(cc))+wzCellTypeSides(t);
	for(;ss<se;ss++) if(ss[0]==c0)
	  {ss[0]=c1; break;}
      }
      for(i0=0;i0<wzCellTypeSSize(t0,s0);i0++){
	if(cs1[wzCellTypeSSide(t1,0,i0)]>0) continue;
	nh=cn0[wzCellTypeSPoint1(t0,s0,i0+1)];
	ct=(dscl[wzCellTypeSSide(t0,s0,i0)])[ds0];
	if(ct<=0){
	  /* "Ubergang zur Nachbarzelle */
	  ch=c0;th=t0;sh=s0;ih=i0;csh=cs0;
	  do{ co = ch;
	  ch = csh[wzCellTypeSSide(th,sh,ih)];
	  dsh= -ccu(ch);	wzAssert(dsh>0);
	  th = cct(ch);
	  cnh= ccn(ch);
	  csh= ccs(ch);
	  for(sh=0;;sh++)	{
	    if(csh[sh]==co) break;
	    wzAssert(sh<9);}
	  for(ih=0;;ih++) {
	    if(cnh[wzCellTypeSPoint1(th,sh,ih+1)]==nh)
	      break;
	    wzAssert(ih<5);}
	  ct=(dscl[wzCellTypeSSide(th,sh,ih)])[dsh];
	  }while(ct<=0);
	}
	cs1[wzCellTypeSSide(t1,0,i0)]=ct;
	cnj = ccn(ct); csj = ccs(ct); tt=cct(ct);
	for(j=0;j<4;j++) if(cnj[wzCellTypeSPoint1(tt,0,j)]==nh){
	  csj[wzCellTypeSSide(tt,0,j)]= c1; break;
	}
      }
    }
    for(dsi=1;dsi<=dsl;dsi++){
      cc=dsc[dsi]; destroyCell(cc);
    }
  }
  return wzSuccess;
}