cogenoctree.cxx

来自「有限元学习研究用源代码(老外的),供科研人员参考」· CXX 代码 · 共 309 行

#include "ibgoctree.hxx"
#include "ibg2.hxx"
#include <stdlib.h>
#include "cogenoctree.hxx"
#include "ibggenerator.hxx"
#include "cogwarnings.hxx"

wzFloat  CogenOctree::epsilon = 1.e-7;

cogenOctree	CogenOctree::cogenoctree() const {return (CogenOctree *) this;}

CogenOctree::CogenOctree()
  :blist(sizeof(CogenOctreeBox))
  ,flist(sizeof(CogenOctreeFace))
  ,rlist(sizeof(CogenOctreeRefineItem))
  ,box(blist.base)
  ,face(flist.base)
  ,refine(rlist.base)
  ,boundaryLevel(4)
{
  setBox();
}

void CogenOctree::setBorder(
			    wzFloat xmin, wzFloat xmax,
			    wzFloat ymin, wzFloat ymax,
			    wzFloat zmin, wzFloat zmax)
{
  min[0]=xmin;  max[0]=xmax;
  min[1]=ymin;  max[1]=ymax;
  min[2]=zmin;  max[2]=zmax;
  for(int i=0;i<wzPointDim;i++){
    if(min[i] != -wzInfty) addPlane(i,min[i]);
    if(max[i] !=  wzInfty) addPlane(i,max[i]);
  }
}

void CogenOctree::addPlane(wzIndex dir, wzFloat x)
{
  if(x>=min[dir] && x<=max[dir]) wzTensorGrid::add(dir,x);
}

void CogenOctree::setFaceInBox(wzFace f,
			      wzFloat xmin,wzFloat xmax,
			      wzFloat ymin,wzFloat ymax,
			      wzFloat zmin,wzFloat zmax)
{
  wzFloat Min[3]={xmin,ymin,zmin};
  wzFloat Max[3]={xmax,ymax,zmax};
  setFaceInBox(f,Min,Max);
}

void CogenOctree::setFaceInBox(wzFace f, wzFloat *xmin, wzFloat *xmax)
{
  int i;
  for(i=0;i<3;i++){
    if(xmax[i]<min[i]){cogInfoEmptyFace();}
    if(xmin[i]>max[i]){cogInfoEmptyFace();}
  }
  wzIndex l = flist.create();
  face[l].face = f;
  for(i=0;i<3;i++){
    if(xmin[i]<min[i])	face[l].min[i] = min[i];
    else		face[l].min[i] = xmin[i];
    if(xmax[i]>max[i])	face[l].max[i] = max[i];
    else		face[l].max[i] = xmax[i];
  }
}

wzIndex CogenOctree::addBox(wzRegion r,
			    wzFloat xmin,wzFloat xmax,
			    wzFloat ymin,wzFloat ymax,
			    wzFloat zmin,wzFloat zmax)
{
  wzFloat Min[3]={xmin,ymin,zmin};
  wzFloat Max[3]={xmax,ymax,zmax};
  return addBox(r,Min,Max);
}

wzIndex CogenOctree::addBox(wzRegion r,wzFloat *xmin, wzFloat *xmax)
{
  int i;
  for(i=0;i<3;i++){
    if(xmax[i]<min[i]){cogInfoEmptyRegion();}
    if(xmin[i]>max[i]){cogInfoEmptyRegion();}
  }
  wzIndex l = blist.create();
  box[l].region = r;
  for(i=0;i<3;i++){
    if(xmin[i]<min[i])	box[l].min[i] = min[i];
    else		box[l].min[i] = xmin[i];
    if(xmax[i]>max[i])	box[l].max[i] = max[i];
    else		box[l].max[i] = xmax[i];
  }
  return l;
}

static const wzFloat dfac = 1+1.e-7;
static const wzFloat dadd = 1/(wzInfty*wzInfty);

wzIndex CogenOctree::refinementInBox(wzFloat dx, wzFloat dy, wzFloat dz, wzIndex b)
{
  wzIndex r = rlist.create(); wzFloat dd;
  refine[r].box = b;
  dd = dfac*dx*dx + dadd;  refine[r].dist[0] = 1/dd;
  dd = dfac*dy*dy + dadd;  refine[r].dist[1] = 1/dd;
  dd = dfac*dz*dz + dadd;  refine[r].dist[2] = 1/dd;
  return r;
}

void CogenOctree::setFaceOfBox(wzFace f, wzIndex r)
{
  CogenOctreeBox& reg = *(box+r);
  setFaceInBox(f,reg.min,reg.max);
}

void CogenOctree::setFaceOfBoxSide(wzFace f, wzIndex r, wzIndex dir)
{
  CogenOctreeBox& R = *(box+r);
  switch(dir){
  case 0: setFaceInBox(f,R.max[0],R.max[0],R.min[1],R.max[1],R.min[2],R.max[2]); break;
  case 1: setFaceInBox(f,R.min[0],R.max[0],R.max[1],R.max[1],R.min[2],R.max[2]); break;
  case 2: setFaceInBox(f,R.min[0],R.max[0],R.min[1],R.max[1],R.max[2],R.max[2]); break;
  case 3: setFaceInBox(f,R.min[0],R.min[0],R.min[1],R.max[1],R.min[2],R.max[2]); break;
  case 4: setFaceInBox(f,R.min[0],R.max[0],R.min[1],R.min[1],R.min[2],R.max[2]); break;
  case 5: setFaceInBox(f,R.min[0],R.max[0],R.min[1],R.max[1],R.min[2],R.min[2]); break;
  }
}

void CogenOctree::endInitialization()
{
  int i,l;
  wzFloat delta;
  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]);
    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);
  wzTensorGrid::getBox(min,max);
  if(max[0]>min[0]) delta = max[0]-min[0];
  if(max[1]>min[1] && max[1]<min[1]+delta) delta = max[1]-min[1];
  if(max[2]>min[2] && max[2]<min[2]+delta) delta = max[2]-min[2];
  setDelta(Delta*delta);
  wzRangeLoop(flist,l){
    CogenOctreeFace& ff = face[l];
    for(i=0;i<3;i++){
      ff.min[i] -= Delta; 
      ff.max[i] += 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;
    }
  }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;
      }
      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;
      }
      ff = d.face;
    next:
      continue;
    }
  }catch(...){
    ff = DefaultFace;
  }
  f.p1.segment() = ff;
  return cogRCConditionFound;
}

wzIndex CogenOctree::Line (cogFlag1& f, const cogLine& s) const
{
  wzIndex rc = Cogeometry::Line(f,s);
  if(rc != cogRCFaceFound) return rc;
  BoundaryCondition(f);
  return cogRCFaceFound;
}

wzpoints CogenOctree::generatePoints(cogeometry g,wzmetric r,wzchart c) const
{
  wzFloat rel = 1.e-3;
  wzpoints plist = new wzPointList;
  ibgoctree gg  = ibgOctreeCreate((wzTensorGrid *)this,g,r,c);
  gg->setMinimalDistance(rel*(max[0]-min[0]),rel*(max[1]-min[1]),rel*(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::create(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(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);
}