cogenprofile.cxx

Delaunay三角形的网格剖分程序

#include "cogenprofile.hxx"
#include "cogwarnings.hxx"
#include <stdio.h>
#include <math.h>
extern int ibgErrorneousRegion;

void CogenProfile::setRelativeBox(wzIndex yb, wzIndex ye, 
				  wzIndex zb, wzIndex ze, 
				  wzIndex mb, wzIndex me, 
				  wzFloat xminrel,wzFloat xmaxrel)
{
  bz=zb;ez=ze,by=yb;ey=ye;bm=mb;em=me;
  if(em>=nx) em = nx-1;
  if(ey>=ny) ey = ny-1;
  if(ez>=nz) ez = nz-1;
  if(bm> em) bm = 0;
  if(by> ey) by = 0;
  if(bz> ez) bz = 0;
  zmin = zz[bz]; zmax=zz[ez];
  ymin = yy[by]; ymax=yy[ey];
  wzFloat dd = xmax-xmin;
  if(xmaxrel==1) xmaxrel=1+1.e-7;
  if(xminrel==0) xminrel= -1.e-7;
  xmax = xmin + dd*xmaxrel;
  xmin = xmin + dd*xminrel;
}

void	CogenProfile::endInitialization()
{
  wzIndex bb,fb,bc;
  x0 = 0;
  y0 = ymin = yy[by];
  z0 = zmin = zz[bz];
  dx = 0.01*(xmax-xmin);
  if(ey>by) dy = ((ymax=yy[ey])-y0)/(ey-by); else dy=1;
  if(ez>bz) dz = ((zmax=zz[ez])-z0)/(ez-bz); else dz=1;
  for(bb=bc=ey-by,fb=2;fb<ey-by;fb*=2){
    bc = bc/2;
    if(fb*bc<bb){
      bb = fb*bc;
      addPlaneY(yy[by+bb]);
    }
  }
  for(bb=bc=ez-bz,fb=2;fb<ez-bz;fb*=2){
    bc = bc/2;
    if(fb*bc<bb){
      bb = fb*bc;
      addPlaneZ(zz[bz+bb]);
    }
  }
  setBorder(xmin,xmax,ymin,ymax,zmin,zmax);
  dx = 1;
  ixold = bm;
  relativeMinimalDistance[0] = 1.e-2;
  relativeMinimalDistance[1] = 1.e-1;
  relativeMinimalDistance[2] = 1.e-1;
  CogenOctree::endInitialization();
}

wzgrid CogenProfile::generateGrid(cogeometry g,wzpoints plist) const
{
  // shifting nodes up to separate them in the Delaunay algorithm
  wzFloat dd = -sig*materialShiftFactor;
  wzRange& range = *(plist.operator->());
  wzArray<wzFloat> shift(range);
  wzPoint PP; wzIndex save=PP.floatValues();
  PP.setFloatValues(save+1);
  //  moveMaterialsUp(*plist,plist->Point,dd);
  moveMaterialsUp(range,plist,shift,dd);
  plist->computeBox();
  wzdelaunayGenerator gen = new wzProfileDelaunayGenerator
    (plist->boxDim,plist->spaceDim);
  gen->Delaunay(plist);
  // shifting nodes back; Note the movement of the grid list: 6 for 3D;
  //moveMaterialsDown(gen->grid()->points,gen->grid()->Point,shift,6);
  moveMaterialsDown(gen->grid()->points,gen->grid()->Point,shift,
		    2*gen->grid()->gdim(),
		    &gen->grid()->cgmin(0),&gen->grid()->cgmax(0));
  try{
    gen->defineRegions(g);
    gen->createFaces(g);
  }catch(...){;}
  PP.setFloatValues(save);
  return gen->grid();
}

void CogenProfile::moveMaterialsDown(wzRange& range, wzArray<wzPoint>& Point, 
				 wzArray<wzFloat>& shift, wzIndex beg,
				 wzFloat *xmin, wzFloat *xmax) const
{
  wzIndex i; wzPoint PP;
  wzIndex save=PP.floatValues()-1;
  *xmin = wzInfty; *xmax = - wzInfty;
  wzRangeLoop(range,i){
    if(i>beg){
      Point[i].x() = Point[i].f(save);
      /*
      Point[i].x() *= horizontalFactor;
      if(Point[i].x() != Point[i].f(save) + shift[i-beg]){
	wzOutput& pout = wzOutput::Default;
	pout("node <>: <>-><> \n")<<
	  i,Point[i].x(),Point[i].x()-(Point[i].f(save) + shift[i-beg]);
      }
      Point[i].x() -= shift[i-beg];Point[i].f(save)
      */
      if(Point[i].x()<*xmin) *xmin = Point[i].x();
      if(Point[i].x()>*xmax) *xmax = Point[i].x();
    }
  }
}

void CogenProfile::moveMaterialsUp(wzRange& range, wzpoints plist, 
				   wzArray<wzFloat>& shift, wzFloat dd) const
{
  wzIndex i; //,Onx,Ony,Onz,Onb;
  wzPoint PP; wzIndex save=PP.floatValues()-1;
  plist->computeBox();
  wzArray<wzPoint>& Point = plist->Point;
  wzFloat minx=plist->min[0]+epsilon, maxx=plist->max[0]-epsilon;
  wzFloat miny=plist->min[1]+epsilon, maxy=plist->max[1]-epsilon;
  wzFloat minz=plist->min[2]+epsilon, maxz=plist->max[2]-epsilon;
  wzFloat d = dd*(plist->max[0]-plist->min[0]);
  wzFloat sh;
  wzRangeLoop(range,i){
    wzPoint& P = Point[i];
    wzSegment s = P.segment();
    if(s.isRegion()){
      sh = s.index()*d; goto add;
    /*
    }
    if(P.x() <= minx || P.x() >= maxx)  Onx = wzTrue; else Onx = wzFalse;
    if(P.y() <= miny || P.y() >= maxy)  Ony = wzTrue; else Ony = wzFalse;
    if(P.z() <= minz || P.z() >= maxz)  Onz = wzTrue; else Onz = wzFalse;
    if(Onx || Ony || Onz) Onb = wzTrue; else Onb = wzFalse;
    if(Onb){
      if(s.isFace()){
	sh = s.index()*d; goto add;
      }else{
	sh = (s.index() - 0.5)*d; goto add;
      }
      */
    }else{
      sh = (s.index() - 0.5)*d; goto add;
    }
  add:
    P.f(save) = P.x();
    P.x() += sh;
    Point[i].x() /= horizontalFactor;
    shift[i]      = sh;
  }
}

wzIndex	CogenProfile::Point(wzPoint& p0) const
{
  cogFloat px = (p0.x()-x0)/dx;
  cogFloat ry = (p0.y()-y0)/dy; 
  cogFloat rz = (p0.z()-z0)/dz; 
  cogInteger ix;
  cogInteger iy = (cogInteger) ry, jy = iy+1; ry -= iy;
  cogInteger iz = (cogInteger) rz, jz = iz+1; rz -= iz;
  cogFloat sz = 1-rz, sy = 1-ry, u0;
  if(iy<0){
    iy = jy = 0;
  }else if(jy>=ny){
    iy = jy = ny-1;
  }
  if(iz<0){
    iz = jz = 0;
  }else if(jz>=nz){
    iz = jz = nz-1;
  }
  
  wzIndex nii = nx*iy+np*iz;
  wzIndex nij = nx*iy+np*jz;
  wzIndex nji = nx*jy+np*iz;
  wzIndex njj = nx*jy+np*jz;
  ix = ixold;
  u0 = sz*(sy*xx[ix+nii]+ry*xx[ix+nji]) + rz*(sy*xx[ix+nij]+ry*xx[ix+njj]);
  if(sig*(px-u0)>0){
    for(ix=ixold-1;ix>=bm;ix--){
      u0 = sz*(sy*xx[ix+nii]+ry*xx[ix+nji]) + rz*(sy*xx[ix+nij]+ry*xx[ix+njj]);
      if(sig*(px-u0)<=0){
	((CogenProfile*)this)->ixold = ix;
	p0.segment() = wzRegion(ix+2);
	return cogRCRegionFound;
      }
    }
    ((CogenProfile*)this)->ixold = bm;
    p0.segment() = wzRegion(bm+1);
    return cogRCRegionFound;
  }else{
    for(ix=ixold+1;ix<=em;ix++){
      u0 = sz*(sy*xx[ix+nii] + ry*xx[ix+nji]) + rz*(sy*xx[ix+nij] + ry*xx[ix+njj]);
      if(sig*(px-u0)>0){
	((CogenProfile*)this)->ixold = ix;
	p0.segment() = wzRegion(ix+1);
	return cogRCRegionFound;
      }
    }
    ((CogenProfile*)this)->ixold = em;
    p0.segment() = wzRegion(em+2);
    return cogRCRegionFound;
  }
}

cogIndex CogenProfile::Line(cogFlag1& f, const cogLine& s) const
{
  return Cogeometry::Line(f,s);
}

cogIndex CogenProfile::BoundaryCondition(cogFlag1& f) const
{
  wzIndex s0 = f.p0.segment().region();
  wzIndex s1 = f.po.segment().region();
  if(s0>s1){
    f.p1.segment() = wzFace(s0);
  }else{
    f.p1.segment() = wzFace(s1);
  }
  return cogRCConditionFound;
}

wzIndex CogenProfile::Triangle(cogFlag2& f, const cogFlag1& o,
			       const cogTriangle& s) const
{
  wzIndex rc = CogenOctree::Triangle(f,o,s);
  wzPoint p1;
  cogFlag1 flag(p1,f.p0,f.po);
  BoundaryCondition(flag);
  f.p2.segment() = wzEdge(p1.segment().index());
  return rc;
}

wzIndex CogenProfile::Tetrahedron(cogFlag3& f, const cogFlag2& o,
				  const cogTetrahedron& s) const
{
  wzIndex rc = CogenOctree::Tetrahedron(f,o,s);
  wzPoint p1;
  cogFlag1 flag(p1,f.p0,f.po);
  BoundaryCondition(flag);
  f.p3.segment() = wzVertex(p1.segment().index());
  return rc;
}

CogenProfile::CogenProfile(wzString filename)
{
  int rc;
  wzIndex cc,iz,iy,ix,count=0;
  float z,y,x;
  wzFloat x1,xb;
  materialShiftFactor = 0.0;
  horizontalFactor = 1.0;
  omitEdgeComputation();
  omitFacesInDirection(0,1,1);
  dmmin = 0.001;
  const wzIndex buflen=256;
  char buffer[buflen];
  FILE *file = fopen(filename,"r"); 
  if(!file){throw wzFileOpenError();}
  fgets(buffer,buflen,file);
  rc = sscanf(buffer,"%d %d %d %d",&ny,&nz,&nx,&cc);
  zz.available(nz);
  yy.available(ny);
  xx.available(nz*ny*nx);
  np = nx*ny;
  if(rc<4) {throw wzFileReadError();}
  if(cc<3) {throw wzFileReadError();}
  xmin = wzInfty; xmax = -wzInfty;
  for(ix=0;ix<nx;ix++){
    for(iz=0;iz<nz;iz++){
      for(iy=0;iy<ny;iy++){
	if(! fgets(buffer,buflen,file)){throw wzFileReadError();}
	while(buffer[0]=='#') fgets(buffer,buflen,file);
	rc = sscanf(buffer,"%f %f %f",&y,&z,&x);
	if(rc<3){throw wzFileReadError();}
	xx[ix+nx*iy+np*iz] = x;
	if(xmin>x) xmin=x;
	if(xmax<x) xmax=x;
	/*
	if(iy || ix){
	  wzAssert(z==zz[iz]);
	}else{
	  zz[iz] = z;
	}
	if(iz){
	  wzAssert(y==y0);
	}else{
	  y0 = y;
	  if(ix){
	    wzAssert(y0==yy[iy]);
	  }else{
	    yy[iy] = y0;
	  }
	}
	*/
	if(iz || ix){
	  wzAssert(y==yy[iy]);
	}else{
	  yy[iy] = y;
	}
	if(iy){
	  wzAssert(z==z0);
	}else{
	  z0 = z;
	  if(ix){
	    wzAssert(z0==zz[iz]);
	  }else{
	    zz[iz] = z0;
	  }
	}
      }
    }
  }
  if(nx<=1) goto end;
  // remove the minimal artificial distances 
  if(xx[0]>xx[1]){
    sig = 1;
    for(iz=0;iz<nz;iz++){
      for(iy=0;iy<ny;iy++){
	x0 = xb = xx[nx*iy+np*iz];
	for(ix=1;ix<nx;ix++){
	  x1 = xx[ix+nx*iy+np*iz];
	  // wzAssert(x1<x0-0.099);
	  if(x1>x0-dmmin){
	    xx[ix+nx*iy+np*iz] = xb;
	    x0 = x1;
	    count++;
	  }else{
	    x0 = xb = x1;
	  }
	}
      }
    }
  }else if(xx[0]<xx[1]){
    sig = -1;
    for(iz=0;iz<nz;iz++){
      for(iy=0;iy<ny;iy++){
	x0 = xb = xx[nx*iy+np*iz];
	for(ix=1;ix<nx;ix++){
	  x1 = xx[ix+nx*iy+np*iz];
	  //	  wzAssert(x1>x0+0.099);
	  if(x1<x0+dmmin){
	    xx[ix+nx*iy+np*iz] = xb;
	    x0 = x1;
	    count++;
	  }else{
	    x0 = xb = x1;
	  }
	}
      }
    }
  }
end:
  ixold = 0;
  bz = 0; ez = nz-1;
  by = 0; ey = ny-1;
  bm = 0; em = nx-1;
  zmin = zz[bz]; zmax=zz[ez];
  ymin = yy[by]; ymax=yy[ey];
}

void 	wzProfileDelaunayGenerator::defineFaceOfCell(int c, cogeometry g)
{
  wzIndex cl,cr,ul,ur;
  wzIndex t = cct(c);
  cl = wzGridCellLeft(*this,c,t);
  cr = wzGridCellRight(*this,c,t);
  if(ccdef(cl)) ul = ccu(cl); else ul=0;
  if(ccdef(cr)) ur = ccu(cr); else ur=0;
  if(ul>ur) ccu(c) = ul;
  else      ccu(c) = ur;
}

void 	wzProfileDelaunayGenerator::defineFaceOutside(int c, cogeometry g)
{
  wzIndex cl,cr,ul,ur;
  wzIndex t = cct(c);
  cl = wzGridCellLeft(*this,c,t);
  cr = wzGridCellRight(*this,c,t);
  if(ccdef(cl)) ul = ccu(cl); else ul=0;
  if(ccdef(cr)) ur = ccu(cr); else ur=0;
  if(ul>ur) ccu(c) = ul;
  else      ccu(c) = ur;
}

//static int suppress_unreachable_warning = 1;
void wzProfileDelaunayGenerator::defineRegionOfCell(int c,
						    cogeometry g)
{
  int n,nn,m=0,m1;
  wzPoint xx;
  wzArray<wzPoint>&  P = Grid.Point;
  //  wzFloat l0;
  wzPoint pp,*p0,*p1;
  wzCellType t = cct(c);
  int npoints = wzCellTypePoints(t);
  int *cn=ccn(c);
  for(n=0;n<npoints;n++){
    if(cnt(nn=cn[n])==wzIsRegion){
      if((m1=cnu(nn)) <= 0){ccu(c) = 0; return;}
      if(m && (m1 != m)) goto differentRegions;
      m = m1;
    }
  }
  if(m){ccu(c) = m; return;}
  //  if(suppress_unreachable_warning) goto useInnerPoint;
  //  l0 = wzInfty;
  for(n=0;n<npoints;n++){
    p0 = &P[cn[n]];
    for(m=n+1;m<npoints;m++){
      p1 = &P[cn[m]];
      if((*p0)[0]!=(*p1)[0]){
	if((*p0)[1]!=(*p1)[1]) continue;
	if((*p0)[2]!=(*p1)[2]) continue;
      }
      //      if((*p0)[1]!=(*p1)[1]){
      //	if((*p0)[2]!=(*p1)[2]) continue;
      //      }
      goto useLineNode;
    }
  }
  goto useInnerPoint;
useLineNode:
  cogLine(*p0,*p1).barycentre(pp);  
  goto ende;
differentRegions:
  cogInfoBadRegionCell();
  if(ibgErrorneousRegion){
    ccu(c) = ibgErrorneousRegion; return;
  }
useInnerPoint:
  {
    wzPoint& P0 = P[cn[0]];  // separate because of a compiler problem
    if(npoints==4)
      cogTetrahedron(P0,P[cn[1]],P[cn[2]],P[cn[3]]).barycentre(pp);
    else
      cogTriangle(P0,P[cn[1]],P[cn[2]]).barycentre(pp);
  }
ende:
  g->Point(pp);
  ccu(c) = pp.segment().region();
  return;
}

void	CogenProfile::getMetric(wzMetricData& data,const cogFlag1& f) const
{
  //  int i;
  wzFloat gxx = G_ii[0]*16;
  wzFloat gyy = G_ii[1];
  wzFloat gzz = G_ii[2];
  ((wzDiagonalMetricData*)(&data))->g_ii[0] = gxx;
  ((wzDiagonalMetricData*)(&data))->g_ii[1] = gyy;
  ((wzDiagonalMetricData*)(&data))->g_ii[2] = gzz;
}

void	CogenProfile::getMetric(wzMetricData& data, const wzPoint& p0) const
{
  cogFloat px = (p0.x()-x0)/dx;
  cogFloat ry = (p0.y()-y0)/dy; 
  cogFloat rz = (p0.z()-z0)/dz; 
  cogInteger nearedge;
  cogInteger iy = (cogInteger) ry, jy = iy+1; ry -= iy;
  cogInteger iz = (cogInteger) rz, jz = iz+1; rz -= iz;
  cogFloat sz = 1-rz, sy = 1-ry;
  if(iy<0){
    iy = jy = 0;
  }else if(jy>=ny){
    iy = jy = ny-1;
  }
  if(iz<0){
    iz = jz = 0;
  }else if(jz>=nz){
    iz = jz = nz-1;
  }
  wzIndex nii = nx*iy+np*iz;
  wzIndex nij = nx*iy+np*jz;
  wzIndex nji = nx*jy+np*iz;
  wzIndex njj = nx*jy+np*jz;
  wzInteger m1  = ((wzIndex) p0.segment().region()) - 1, m0 = m1-1;
  wzFloat  yf[2],zf[2],yg[2],zg[2];
  if(m0>=bm){
    yf[0] = sy*xx[m0+nii] + ry*xx[m0+nji];
    yf[1] = sy*xx[m0+nij] + ry*xx[m0+njj];
    zf[0] = sz*xx[m0+nii] + rz*xx[m0+nij];
    zf[1] = sz*xx[m0+nji] + rz*xx[m0+njj];
  }else{
    yf[0] = yf[1] = zf[0] = zf[1] = -wzInfty;
  }
  //  wzIndex nearedge = 0;
  if(m1<=em){
    yg[0] = sy*xx[m1+nii] + ry*xx[m1+nji];
    yg[1] = sy*xx[m1+nij] + ry*xx[m1+njj];
    zg[0] = sz*xx[m1+nii] + rz*xx[m1+nij];
    zg[1] = sz*xx[m1+nji] + rz*xx[m1+njj];
    if(m0>=bm){
      if(xx[m0+nii]==xx[m1+nii]) nearedge = 1;
      if(xx[m0+nij]==xx[m1+nij]) nearedge = 1;
      if(xx[m0+nji]==xx[m1+nji]) nearedge = 1;
      if(xx[m0+njj]==xx[m1+njj]) nearedge = 1;
    }
  }else{
    yg[0] = yg[1] = zg[0] = zg[1] = wzInfty;
  }
  wzFloat deltam = sz*(yf[0]-yg[0]) + rz*(yf[1]-yg[1]);	if(deltam<0) deltam = -deltam;
  wzFloat deltay = (yg[0]+yf[0]) - (yg[1]+yf[1]);  	if(deltay<0) deltay = -deltay;
  wzFloat deltaz = (zg[0]+zf[0]) - (zg[1]+zf[1]);  	if(deltaz<0) deltaz = -deltaz;
  if(deltam<dmmin){
    wzAssert(nearedge);
    deltam = dmmin;
  }
  wzFloat gxx,gyy,gzz;
  //  wzFloat drx = 1/deltam; drx *= drx;
  wzFloat dry = deltay/(deltam*dy); dry *= dry;
  wzFloat drz = deltaz/(deltam*dz); drz *= drz;
  //  if(G_ii[0]<drx)  gxx = drx; else gxx = G_ii[0];
  if(G_ii[1]<dry)  gyy = dry; else gyy = G_ii[1];
  if(G_ii[2]<drz)  gzz = drz; else gzz = G_ii[2];
  //  if(G_ii[0]>16*drx)  gxx = 16*G_ii[0];
  gxx = G_ii[0];
  wzFloat ifac=8024.0;
  if(ifac*G_ii[1]<dry)  gyy = ifac*G_ii[1];
  if(ifac*G_ii[2]<drz)  gzz = ifac*G_ii[2];
  //  wzAssert(256*G_ii[2]>drz);
  ((wzDiagonalMetricData*)(&data))->g_ii[0] = gxx;
  ((wzDiagonalMetricData*)(&data))->g_ii[1] = gyy;
  ((wzDiagonalMetricData*)(&data))->g_ii[2] = gzz;
}