cogenoctree.cxx

Delaunay三角形的网格剖分程序

	rr.min[i] = (1-epsilon)*rr.mino[i];
      else
	rr.min[i] = (1+epsilon)*rr.mino[i];
      if(rr.maxo[i]>=0)
	rr.max[i] = (1+epsilon)*rr.maxo[i];
      else
	rr.max[i] = (1-epsilon)*rr.maxo[i];
    }

// dangerous for large anisotropies - renormalization not related with
// the correct two sides, but with the largest side.
    rr.d = rr.dorig + epsilon*delta;  
  }
  if(Delta<10*epsilon) Delta = 10*epsilon; 
  if(Delta<1.e-10) Delta=1.e-10;
  if(Delta>1.e-4) Delta=1.e-4;
  setDelta(Delta*delta);
  wzRangeLoop(flist,l){
    CogenOctreeFace& ff = face[l];
    for(i=0;i<3;i++){
      ff.min[i] -= Delta; 
      ff.max[i] += Delta;
    }
    ff.d += Delta;
  }
}

void CogenOctree::endInitialization()
{
  int l;
  if(checked) return;
  //  wzFloat eps = epsilon;
  wzRangeLoop(blist,l){
    CogenOctreeBox& rr = box[l];
    addPlanesOfPoint(rr.min[0],rr.min[1],rr.min[2]);
    addPlanesOfPoint(rr.max[0],rr.max[1],rr.max[2]);
    /*
    int i;
    for(i=0;i<3;i++){
      if(rr.min[i]>=0)
	rr.min[i] -= eps*(rr.max[i]-rr.min[i]);
      else
	rr.min[i] += eps*rr.min[i];
      if(rr.max[i]>=0)
	rr.max[i] += eps*(rr.max[i]-rr.min[i]);
      else
	rr.max[i] -= eps*rr.max[i];
    }
    */
  }
  wzRangeLoop(flist,l){
    CogenOctreeFace& ff = face[l];
    addPlanesOfPoint(ff.min[0],ff.min[1],ff.min[2]);
    addPlanesOfPoint(ff.max[0],ff.max[1],ff.max[2]);
  }
  wzTensorGrid::endInitialization(min,max);	// now "checked" is true;
  wzTensorGrid::getBox(min,max);
  setEpsilon(epsilon);
  setDelta  (1.e-6);
}

void 	CogenOctree::setDelta(cogFloat delta)
{
  wzAssert(delta>2*epsilon);
  Cogeometry::setDelta(delta);
}

wzBoolean CogenOctree::inBox(wzIndex b, const wzPoint& p) const
{
  CogenOctreeBox& d = *(box+b);
  int i;
  wzFloat Y[wzPointDim];
  try{
    y(Y,p);
    for(i=0;i<wzPointDim;i++){
      if(d.min[i]>Y[i]) return wzFalse;
      if(d.max[i]<Y[i]) return wzFalse;
    }
    if(d.type == cogWedgeType){
      wzFloat d1 = d.d;
      for(i=0;i<wzPointDim;i++){
	d1 += Y[i]*d.a[i];
      }
      if(d1<0) return wzFalse;
    }
  }catch(...){return wzFalse;}
  return wzTrue;
}

void CogenOctree::getMetric(wzMetricData& data, const wzPoint& p) const
{
  wzIndex it,i,b;
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  for(i=0;i<wzPointDim;i++) d.g_ii[i] = G_ii[i];
  wzRangeLoop(rlist,it){
    b = refine[it].box;
    for(i=0;i<wzPointDim;i++){
      if(inBox(b,p)){
	for(i=0;i<wzPointDim;i++){
	  if(d.g_ii[i]<refine[it].dist[i]) d.g_ii[i] = refine[it].dist[i];
	}
      }
    }
  }
}

wzIndex CogenOctree::Point(wzPoint& p0) const
{
  wzFloat Y[wzPointDim];
  int i,l;
  wzIndex r = DefaultRegion;
  try{
    y(Y,p0);
    wzRangeLoop(blist,l){
      CogenOctreeBox& d = *(box+l);
      for(i=0;i<wzPointDim;i++){
	if(d.min[i]>Y[i]) goto next;
	if(d.max[i]<Y[i]) goto next;
      }
      if(d.type == cogWedgeType){
	wzFloat d1 = d.d;
	for(i=0;i<wzPointDim;i++){
	  d1 += Y[i]*d.a[i];
	}
	if(d1<0) goto next;
      }
      r = d.region;
    next:
      continue;
    }
    p0.segment()=wzRegion(r);
  }catch(...){
    p0.segment()=DefaultRegion;
  }
  return cogRCRegionFound;
}

cogIndex CogenOctree::BoundaryCondition(cogFlag1& f) const
{
  wzFloat Y[wzPointDim];
  int i,l;
  wzFace ff = DefaultFace;
  try{
    y(Y,f.p1);
    wzRangeLoop(flist,l){
      CogenOctreeFace& d = *(face+l);
      for(i=0;i<wzPointDim;i++){
	if(d.min[i]>Y[i]) goto next;
	if(d.max[i]<Y[i]) goto next;
      }
      if(d.type == cogWedgeType){
	wzFloat d1 = d.d;
	for(i=0;i<wzPointDim;i++){
	  d1 += Y[i]*d.a[i];
	}
	if(d1<0) goto next;
      }
      ff = d.face;
    next:
      continue;
    }
  }catch(...){
    ff = DefaultFace;
  }
  f.p1.segment() = ff;
  return cogRCConditionFound;
}

wzIndex CogenOctree::Line (cogFlag1& f, const cogLine& s) const
{
  if(SimpleBinaryIntersectionForLine){
    wzIndex rc = Cogeometry::Line(f,s);
    if(rc != cogRCFaceFound) return rc;
    BoundaryCondition(f);
    return cogRCFaceFound;
  }
  // wzPoint Y0,Y1;
  wzFloat p,q,qfound,D[wzPointDim],Y0[wzPointDim],Y1[wzPointDim],F1[wzPointDim];
  //  wzPoint F1;
  int i,j,l;
  wzIndex rc = Cogeometry::Line(f,s);
  try{
    y(Y0,s.c0);
    y(Y1,s.c1);
    for(i=0;i<wzPointDim;i++){D[i]=Y1[i]-Y0[i];}
    qfound = 1;
    wzRangeLoop(blist,l){
      const CogenOctreeBox& rr = box[l];
      for(i=0;i<wzPointDim;i++){
	if(D[i]!=0) q = (rr.mino[i]-Y0[i])/D[i]; else continue;
	if(q>=qfound || q<0) continue;
	p = (Y1[i]-rr.mino[i])/D[i];	// is 1 if mino == Y0
	for(j=0;j<wzPointDim;j++){F1[j]=p*Y0[j]+q*Y1[j];}
	F1[i]=rr.mino[i];
	if(!detectingFlag(f,F1,Y0,Y1,p,l)) continue;
	qfound = q;
      }
      for(i=0;i<wzPointDim;i++){
	if(D[i]!=0) q = (rr.maxo[i]-Y0[i])/D[i]; else continue;
	if(q>=qfound || q<0) continue;
	p = (Y1[i]-rr.maxo[i])/D[i];	// is 1 if maxo == Y0
	for(j=0;j<wzPointDim;j++){F1[j]=p*Y0[j]+q*Y1[j];}
	F1[i]=rr.maxo[i];
	if(!detectingFlag(f,F1,Y0,Y1,p,l)) continue;
	qfound = q;
      }
      if(rr.type == cogWedgeType){
	wzFloat d0 = rr.dorig;
	wzFloat d1 = rr.dorig;
	for(i=0;i<wzPointDim;i++){
	  d0 += Y0[i]*rr.a[i];
	  d1 += Y1[i]*rr.a[i];
	}
	if(d1!=d0) q = -d0/(d1-d0); else continue;
	if(q>=qfound || q<0) continue;
	p =  d1/(d1-d0);	// is 1 if Y0 is on the line (d0=0)
	for(j=0;j<wzPointDim;j++){F1[j]=p*Y0[j]+q*Y1[j];}
	if(!detectingFlag(f,F1,Y0,Y1,p,l)) continue;
	qfound = q;
      }    
    }
  }catch(...){
    wzAssert(0);
    qfound = 1;
  }
  if(qfound>=1){ // no exact intersection found: fallback to the default;
    //    wzAssert(0);
    wzIndex rc = Cogeometry::Line(f,s);
    if(rc != cogRCFaceFound) return rc;
  }
  BoundaryCondition(f);
  return cogRCFaceFound;
}

wzBoolean CogenOctree::detectingFlag(cogFlag1 f, wzFloat* F1, 
				     wzFloat* Y0, 
				     wzFloat* Y1,
				     wzFloat p1, int box) const
{
  wzFloat eps = 1.e-5; // wrong - should be appropriately computed delta!!!
  wzFloat p0,q0,po,qo,q1,F0[wzPointDim],Fo[wzPointDim];
  int j;
  if(!inBox(box,f.p1)) return wzFalse;
  p0 = p1 + eps;  if (p0>1) {
    p0=1;p1=1-eps;q1=eps; for(j=0;j<wzPointDim;j++)F1[j]=p1*Y0[j]+q1*Y1[j];
  }
  po = p1 - eps;  if (po<0) {
    po=0;p1=eps;q1=1-eps; for(j=0;j<wzPointDim;j++)F1[j]=p1*Y0[j]+q1*Y1[j];
  }
  q0 = 1-p0;
  qo = 1-po;
  ((wzChart*)this)->x(f.p1,F1);
  for(j=0;j<wzPointDim;j++){
    F0[j]=p0*Y0[j]+q0*Y1[j];
    Fo[j]=po*Y0[j]+qo*Y1[j];
  }
  ((wzChart*)this)->x(f.p0,F0);	Point(f.p0);
  ((wzChart*)this)->x(f.po,Fo);	Point(f.po);
  if(f.p0.segment() == f.po.segment()){
    //    wzAssert(0);	// should happen only if something is hidden!
    return wzFalse;
  }
  return wzTrue;
}

void CogenOctree::omitFacesInDirection(wzBoolean x, wzBoolean y, wzBoolean z)
{
  suppressFaceIntersection[0] = x;
  suppressFaceIntersection[1] = y;
  suppressFaceIntersection[2] = z;
  suppressFaceIntersection[3] = x + y + z;
}

wzpoints CogenOctree::generatePoints(cogeometry g,wzmetric r,wzchart c) const
{
  //  wzFloat rel = 1.e-3;
  wzpoints plist = new wzPointList;
  ibgoctree gg = ibgOctreeCreate((CogenOctree *)this,(wzTensorGrid *)this,g,r,c);
  gg->setMinimalDistance(relativeMinimalDistance[0]*(max[0]-min[0]),
			 relativeMinimalDistance[1]*(max[1]-min[1]),
			 relativeMinimalDistance[2]*(max[2]-min[2]));
  //  if(boundaryLevel<2)
  //    gg->generateWithoutEdges();
  //  else
    gg->generate();
  mark(gg,g->Delta);
  ibgGetPoints(plist,gg);
  return plist;
}

wzFloat CogenOctree::chi(const wzPoint& p) const
{
  for(int i=0;i<wzPointDim;i++){
    if(p[i] < min[i]-epsilon)	return -1;
    if(p[i] > max[i]+epsilon) 	return -1;
  }
  return 1;
}

void CogenOctree::getMetric(wzMetricData& data, const cogFlag1& f) const
{
  wzDiagonalMetricSimple::getMetric(data,f);
}


// obsolete:

#include "ibggenerator.hxx"

wzgrid ibg3Generator::createGrid(cogeometry g,wztensorgrid box,wzmetric r,wzchart c)
{
  wzFloat rel = 1.e-3;
  wzpoints plist = new wzPointList; 
  wzFloat min[]={0,0,0};
  wzFloat max[]={0,0,0};
  box->endInitialization(min,max);
  ibgoctree gg  = ibgOctreeCreate(new CogenOctree(),box,g,r,c);
  gg->setMinimalDistance(rel*(max[0]-min[0]),rel*(max[1]-min[1]),rel*(max[2]-min[2]));
  gg->generate();
  //  mark(gg,g->Delta);
  ibgGetPoints(plist,gg);
  Cogenerator gen;
  return gen.generateGrid(g,plist);
}