ibgoinit.cxx

Delaunay三角形的网格剖分程序

#include "ibgoctree.hxx"
#include "wzview.hxx"

ibgIndex 		ibgOctreeData::maxFacesOnLine = 2;
const ibgInteger	ibgOctreeData::outside =  0;
const ibgInteger	ibgOctreeData::nothing = -999999999;
wzgrid		wzCurrentData::grid		= 0;
ibgoctree	wzCurrentData::octree		= 0;
wzpoints	wzCurrentData::pointList	= 0;
cogeometry	wzCurrentData::geometry		= 0;
cogenerator	wzCurrentData::generator	= 0;
wzmetric	wzCurrentData::metric		= 0;
wzCurrentData&	wzCurrentData::current		= * new wzCurrentData;
wzFloat		wzCurrentData::factor[3]	= {1,1,1};
wzFloat		wzCurrentData::shift[3]		= {0,0,0};

ibgOctreeData::ibgOctreeData(ibgIndex dim,
			     cogenOctree o,
			     wztensorgrid coarse,
			     cogeometry g, wzmetric r, wzchart c)
:ibgOctreeType(dim)
,nodes(sizeof(ibgPoint))
  //,cells(Corners*sizeof(ibgIndex))
,cells((Corners?Corners+1:0)*sizeof(ibgIndex))
,bnodes(2*sizeof(ibgIndex))
,lines (5*sizeof(ibgIndex))
,planes(8*sizeof(ibgIndex))
,bcells((3+ibgMaxEdgesOnCube)*sizeof(ibgIndex))
,faces(3*sizeof(ibgPoint))
,edges(4*sizeof(ibgPoint))
,vertices(5*sizeof(ibgPoint))
,nPoint(nodes.base)
,nx(nodes.base,wzOffsetOf(ibgPoint,X)),ny(nx,1),nz(ny,1)
,nOld(nodes)
,nBoundary(nodes)
,nStackInitial(nodes)
,nStackMetric(nStackInitial)
,nStackOrthogonalRegularization(nStackInitial)
,nStackLinearRegularization(nStackInitial)
,nStackFace(nStackInitial)
,nStackTest(nStackInitial)
,nMetric(nodes)
,bnBase(bnodes.base)
,bnType(bnBase,1)
,lNodeUp(lines.base)
,lNodeDown(lNodeUp,1)
,lFaceUp(lNodeUp,2)
,lFaceDown(lNodeUp,3)
,lType(lNodeUp,4)
,pEdge(planes.base)
,pType(pEdge,1)
,fPoint(faces.base)
,fPointUp(fPoint,1)
,fPointDown(fPoint,2)
,fUp(faces,4*sizeof(ibgInteger))
,fDown(fUp,1)
,fLine(fUp,2)
,fType(fUp,3)
,fStackInitial(faces)
,fStackMetric(fStackInitial)
,fStackDoubleIntersection(fStackInitial)
,fStackRegularization(fStackInitial)
,fMetric(faces)
,ePoint(edges.base)
,ePointDown(ePoint,1)
,ePointUp(ePoint,2)
,ePointFace(ePoint,3)
,eFace(edges,6*sizeof(ibgInteger))
,ePlane(eFace,1)
,eNext(eFace,2)
,eIdent(eFace,3)
,bcBase(bcells.base)
,bcType(bcBase,1)
,bcVertex(bcBase,2)
,vPoint(vertices.base)
,vPointEdge(vPoint,1)
,vPointFace(vPoint,2)
,vPointUp(vPoint,3)
,vPointDown(vPoint,4)
,vEdge(vertices,4*sizeof(ibgInteger))
,vCell(vEdge,1)
,vNext(vEdge,2)
,vIdent(vEdge,3)
,gen(o)
,geo(g)
,ref(r)
,chart(c)
{
 initialize(coarse);
}

void ibgOctreeData::initialize(wztensorgrid coarse)
{unsigned i,pp0;
 unsigned ix,jx,kx,lx,mx,ppx,iy,jy,ky,ly,my,ppy,iz,jz,kz,lz,mz,ppz;
 int ii,i0,i1,i2,k0,k1,k2,kk,jj,o,r,fx,fy,fz;
 ibgFloat xx,dx,yy,dy,zz,dz;
 // initialization of the arrays of type wzArray<>
 nX[0].reinitialize(nx,0);
 nX[1].reinitialize(ny,0);
 nX[2].reinitialize(nz,0);
 fFlag[0].reinitialize(fPoint,0);
 fFlag[1].reinitialize(fPointUp,0);
 fFlag[2].reinitialize(fPointDown,0);
 eFlag[0].reinitialize(ePoint,0);
 eFlag[1].reinitialize(ePointFace,0);
 eFlag[2].reinitialize(ePointUp,0);
 eFlag[3].reinitialize(ePointDown,0);
 eVertex[0].reinitialize(eFace,4);
 eVertex[1].reinitialize(eFace,5);
 vFlag[0].reinitialize(vPoint,0);
 vFlag[1].reinitialize(vPointEdge,0);
 vFlag[2].reinitialize(vPointFace,0);
 vFlag[3].reinitialize(vPointUp,0);
 vFlag[4].reinitialize(vPointDown,0);
 fNext[0].reinitialize(fUp,0);
 fNext[1].reinitialize(fDown,0);
 lNode[0].reinitialize(lNodeUp,0);
 lNode[1].reinitialize(lNodeDown,0);
 pNode[0].reinitialize(pEdge,2);
 pNode[1].reinitialize(pEdge,3);
 pNode[2].reinitialize(pEdge,4);
 pNode[3].reinitialize(pEdge,5);
 pCell[0].reinitialize(pEdge,6);
 pCell[1].reinitialize(pEdge,7);
 nNext[0].reinitialize(nodes,2*GridDim);
 bnNext[0].reinitialize(nodes,2*GridDim);
 for(i=0;i<ibgMaxEdgesOnCube;i++){
   bcEdge[i].reinitialize(bcBase,3+i);
 }
 for(i=0;i<GridDim;i++){
   if(i>0) nNext[i].reinitialize(nNext[0],i);
   nNext[GridDim+i].reinitialize(nNext[0],GridDim+i);
   nUp[i].reinitialize(nNext[0],i);
   nDown[i].reinitialize(nNext[0],GridDim+i);
   if(i>0) bnNext[i].reinitialize(bnNext[0],i);
   bnNext[GridDim+i].reinitialize(bnNext[0],GridDim+i);
   bnUp[i].reinitialize(bnNext[0],i);
   bnDown[i].reinitialize(bnNext[0],GridDim+i);
 }
 bnPlane[0].reinitialize(bnodes,Planes);
 for(i=1;i<Planes;i++){
   bnPlane[i].reinitialize(bnPlane[0],i);
 }
 fEdge[0].reinitialize(faces,LineSides);
 lPlane[0].reinitialize(lines,LineSides);
 for(i=1;i<LineSides;i++){
   fEdge[i].reinitialize(fEdge[0],i);
   lPlane[i].reinitialize(lPlane[0],i);
 }
 mx=coarse->dimension(0);	fx=coarse->refinement(0)+1;
 my=coarse->dimension(1);	fy=coarse->refinement(1)+1;
 mz=coarse->dimension(2);  	fz=coarse->refinement(2)+1;
 lz=fz*(mz-1);	ppz=1;		//if(lz>100) CoarseGridTooFine();
 ly=fy*(my-1);	ppy=lz+1;	//if(ly>100) CoarseGridTooFine();
 lx=fx*(mx-1);	ppx=ppy*(ly+1);	//if(lx>100) CoarseGridTooFine();
 for(i=0;i<3;i++){
   min[i] = coarse->x[i][1];
   max[i] = coarse->x[i][coarse->dimension(i)];
 }
 i=1; 				pp0=ppx*(lx+1);
 nodes.available(pp0);
 iz=1; jz=0; kz=0; zz=coarse->x[2][iz]; if(mz>1) dz=(coarse->x[2][iz+1]-zz)/fz;
 iy=1; jy=0; ky=0; yy=coarse->x[1][iy]; if(my>1) dy=(coarse->x[1][iy+1]-yy)/fy;
 ix=1; jx=0; kx=0; xx=coarse->x[0][ix]; if(mx>1) dx=(coarse->x[0][ix+1]-xx)/fx;
 while(1){
	int i1 = nodes.create();
	nStackInitial.append(i);
	wzAssert(i1==i);
 	nx[i] = xx;	ny[i] = yy;	nz[i] = zz;
	geo->Point(nPoint[i]);
	nBoundary[i]=0;
#ifdef GLEVEL
	nLevel[0][i] =	nLevel[1][i] =	nLevel[2][i] = 0;
#endif
	switch(GridDim){
	case 3:
		if(kz<lz) nUp[2][i]	= i+ppz;
		else      nUp[2][i]	= outside;
		if(kz> 0) nDown[2][i]	= i-ppz;
		else      nDown[2][i]	= outside;
	case 2:
		if(ky<ly) nUp[1][i]	= i+ppy;
		else      nUp[1][i]	= outside;
		if(ky> 0) nDown[1][i]	= i-ppy;
		else      nDown[1][i]	= outside;
	case 1:
		if(kx<lx) nUp[0][i]	= i+ppx;
		else      nUp[0][i]	= outside;
		if(kx> 0) nDown[0][i]	= i-ppx;
		else      nDown[0][i]	= outside;
  	case 0: break;
	}
	if(kz==0){
		if(ky==0){
			if(kx==0)nOld[i] = 0;
			else	 nOld[i] = i-(ly+1)*(lz+1);
		}else		 nOld[i] = i-(lz+1);
	}else			 nOld[i] = i-1;
	i++;
	if(iz==mz){
		iz=1; jz=0; kz=0; zz=coarse->x[2][iz];
		if(mz>iz) dz=(coarse->x[2][iz+1]-zz)/fz;
		if(iy==my){
			iy=1; jy=0; ky=0; yy=coarse->x[1][iy];
			if(my>iy) dy=(coarse->x[1][iy+1]-yy)/fy;
			if(ix==mx) break;
			else if(jx==fx-1) {
				jx=0; kx++; ix++; xx=coarse->x[0][ix];
			        if(mx>ix) dx=(coarse->x[0][ix+1]-xx)/fx;
			}else{
				jx++; kx++; xx+=dx;
			}
		}else if(jy==fy-1) {
			jy=0; ky++; iy++; yy=coarse->x[1][iy];
			if(my>iy) dy=(coarse->x[1][iy+1]-yy)/fy;
		}else{
			jy++; ky++; yy+=dy;
		}
	}else if(jz==fz-1){
 		jz=0; iz++; kz++; zz=coarse->x[2][iz];
		if(mz>iz) dz=(coarse->x[2][iz+1]-zz)/fz;
	}else{
		jz++; kz++; zz+=dz;
	}
 }
 switch(GridDim){
 case 3:
   nCell[0].reinitialize(nodes,Corners);
   cNode[0].reinitialize(cells.base);
   for(o=1;o<Corners;o++){
     nCell[o].reinitialize(nCell[0],o);
     cNode[o].reinitialize(cNode[0],o);
   }
   cBoundary.reinitialize(cNode[0],Corners);
   cells.create(lx*ly*lz);
   for(ii=1,jj=1;ii<=pp0;ii++){
     for(r=0;r<Lines;r++){
       if(nNext[r][ii]==outside){
	 for(o=0;o<LineSides;o++) nCell[lsCell[r][o]][ii] = outside;
       }
     }
     if(is_undefined(i0=nNext[0][ii])) continue;
     if(is_undefined(i1=nNext[1][ii])) continue;
     if(is_undefined(i2=nNext[2][ii])) continue;
     k0=nNext[2][i1]; k1=nNext[0][i2]; k2=nNext[1][i0];
     kk=nNext[0][k0];
     cBoundary[jj] = 0;
     cNode[cdirnnn][jj] = ii; nCell[cdirppp][ii] = jj;
     cNode[cdirpnn][jj] = i0; nCell[cdirnpp][i0] = jj;
     cNode[cdirnpn][jj] = i1; nCell[cdirpnp][i1] = jj;
     cNode[cdirnnp][jj] = i2; nCell[cdirppn][i2] = jj;
     cNode[cdirnpp][jj] = k0; nCell[cdirpnn][k0] = jj;
     cNode[cdirpnp][jj] = k1; nCell[cdirnpn][k1] = jj;
     cNode[cdirppn][jj] = k2; nCell[cdirnnp][k2] = jj;
     cNode[cdirppp][jj] = kk; nCell[cdirnnn][kk] = jj;
     jj++;
   }
 default:break;
 }
}

ibgOctree0D::ibgOctree0D(cogenOctree o,wztensorgrid gr,cogeometry g,wzmetric r,wzchart c)
:ibgOctree(0,o,gr,g,r,c)
{;}
ibgOctree1D::ibgOctree1D(cogenOctree o,wztensorgrid gr,cogeometry g,wzmetric r,wzchart c)
:ibgOctree(1,o,gr,g,r,c)
{;}
ibgOctree2D::ibgOctree2D(cogenOctree o,wztensorgrid gr,cogeometry g,wzmetric r,wzchart c)
:ibgOctree(2,o,gr,g,r,c)
{;}
ibgOctree3D::ibgOctree3D(cogenOctree o,wztensorgrid gr,cogeometry g,wzmetric r,wzchart c)
:ibgOctree(3,o,gr,g,r,c)
{;}

ibgoctree ibgOctreeCreate(cogenOctree o,wztensorgrid gr,
			  cogeometry g,wzmetric r,wzchart c)
{
  int i = gr->dimension();
  ibgoctree octree;
  wzCurrentData::geometry  = g;
  wzCurrentData::generator = o->generator();
  wzCurrentData::metric	   = r;
  switch(i){
  case 3:	octree = new ibgOctree3D(o,gr,g,r,c); break;
  case 2:	octree = new ibgOctree2D(o,gr,g,r,c); break;
  case 1:	octree = new ibgOctree1D(o,gr,g,r,c); break;
  case 0:	octree = new ibgOctree0D(o,gr,g,r,c); break;
  }
  wzCurrentData::octree = octree;
  return octree;
}

ibgOctree::ibgOctree(ibgIndex dim,
		     cogenOctree gen,
		     wztensorgrid coarse,
		     cogeometry g,wzmetric r,wzchart c)
  :ibgOctreeData(dim,gen,coarse,g,r,c)
{;}