cogenrefine.cxx

有限元学习研究用源代码(老外的),供科研人员参考

#include "cogenrefine.hxx"

static const wzFloat MinDelta = 1.e-16;
static const wzIndex nRefPoint = 0;
static const wzIndex nRefBox = 1;
//static const wzIndex nRefPoint2 = 2;
//static const wzIndex nRefPoint3 = 3;
static const wzIndex nRefRegion = 9;
static const wzIndex fRefFace = 10;
static const wzIndex fRefBoundary = 11;
static const wzIndex fRefBetween = 12;
  
cogenRefine	CogenRefine::cogenrefine() const {return (CogenRefine*) this;}

CogenRefine::CogenRefine(cogenerator gen)
  :CogenDecorator(gen)
  ,nList(sizeof(CogenRefineItem))
  ,nItem(nList.base)
  ,fList(sizeof(CogenRefineItem))
  ,fItem(fList.base)
  ,oldRefinement(gen->refinement())
  ,PointMode(nRefPoint)
{
}

void CogenRefine::setPointModeRectangular()
{PointMode = nRefBox;}
void CogenRefine::setPointModeSpherical()  
{PointMode = nRefPoint;}

// isotropical variants:

wzIndex CogenRefine::refinementIsotropicalInRegion(wzFloat d, wzRegion r)
{return refinementInRegion(d,d,d,r);}

wzIndex CogenRefine::refinementIsotropicalNearBoundary(wzFloat d, wzRegion r)
{return refinementNearBoundary(d,d,d,r);}

wzIndex CogenRefine::refinementIsotropicalBetween(wzFloat d, wzRegion r1, wzRegion r2)
{return refinementBetween(d,d,d,r1,r2);}

wzIndex CogenRefine::refinementIsotropicalNearFace(wzFloat d, wzFace f)
{return refinementNearFace(d,d,d,f);}

wzIndex CogenRefine::refinementIsotropicalNearPoint(wzFloat d, wzFloat n,
						    wzFloat xx, wzFloat yy, wzFloat zz)
{return refinementNearPoint(d,d,d,n,n,n,xx,yy,zz);}

wzIndex CogenRefine::refinementInRegion(wzFloat dx, wzFloat dy, wzFloat dz, wzRegion s)
{
  wzIndex i = nList.create();
  CogenRefineItem& r = nItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = nRefRegion;
  r.s    = s;
  return i;
}

wzIndex CogenRefine::refinementNearFace(wzFloat dx, wzFloat dy, wzFloat dz, wzFace s)
{
  wzIndex i = fList.create();
  CogenRefineItem& r = fItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = fRefFace;
  r.s    = s;
  return i;
}

wzIndex CogenRefine::refinementNearBoundary(wzFloat dx, wzFloat dy, wzFloat dz, wzRegion s)
{
  wzIndex i = fList.create();
  CogenRefineItem& r = fItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = fRefBoundary;
  r.s    = s;
  return i;
}

wzIndex CogenRefine::refinementBetween(wzFloat dx, wzFloat dy, wzFloat dz, 
				       wzRegion s1, wzRegion s2)
{
  wzIndex i = fList.create();
  CogenRefineItem& r = fItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  r.typ  = fRefBetween;
  r.s    = s1;
  r.so   = s2;
  return i;
}

wzIndex CogenRefine::refinementNearPoint(wzFloat dx, wzFloat dy, wzFloat dz,
					 wzFloat nx, wzFloat ny, wzFloat nz,
					 wzFloat xx, wzFloat yy, wzFloat zz)
{
  wzIndex i = nList.create();
  CogenRefineItem& r = nItem[i];
  if(dx<MinDelta) dx = MinDelta;
  if(dy<MinDelta) dy = MinDelta;
  if(dz<MinDelta) dz = MinDelta;
  r.dist[0] = 1.0/(dx*dx);
  r.dist[1] = 1.0/(dy*dy);
  r.dist[2] = 1.0/(dz*dz);
  if(nx<dx) nx=dx;
  if(ny<dy) ny=dy;
  if(nz<dz) nz=dz;
  r.near[0] = nx;
  r.near[1] = ny;
  r.near[2] = nz;
  r.x[0] = xx;
  r.x[1] = yy;
  r.x[2] = zz;
  r.typ  = PointMode;
  return i;
}

void CogenRefine::getMetric(wzMetricData& data, const cogFlag1& f) const
{
  wzIndex it,i;
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  wzDiagonalMetricSimple::getMetric(data,f);
  wzRangeLoop(fList,it){
    const CogenRefineItem& r = fItem[it];
    switch(r.typ){
    case fRefFace:
      if(r.s != f.p1.segment()) continue; 
      break;
    case fRefBetween:
      if(r.s == f.p0.segment()){
	if(r.so != f.po.segment()) continue; 
      }else if(r.so == f.p0.segment()){
	if(r.s != f.po.segment()) continue; 
      }else{
	continue;
      }
      break;
    case fRefBoundary:
      if(r.s != f.p0.segment()){
	if(r.s != f.po.segment()) continue; 
      }
      break;
    default:
      wzAssert(0); continue;
    }
    for(i=0;i<wzPointDim;i++){
      if(d.g_ii[i]<r.dist[i]) d.g_ii[i] = r.dist[i];
    }
  }
}

void CogenRefine::getMetric(wzMetricData& data, const wzPoint& p) const
{
  wzIndex it,i;
  wzFloat x,f,g;
  wzDiagonalMetricData& d = *((wzDiagonalMetricData*)(&data));
  oldRefinement->getMetric(data,p);
  for(i=0;i<wzPointDim;i++) 
    if(d.g_ii[i] < G_ii[i]) d.g_ii[i] = G_ii[i];
  wzRangeLoop(nList,it){
    const CogenRefineItem& r = nItem[it];
    switch(r.typ){
    case nRefRegion:
      if(r.s == p.segment()){
	for(i=0;i<wzPointDim;i++){
	  if(d.g_ii[i]<r.dist[i]) d.g_ii[i] = r.dist[i];
	}
      }
      break;
    case nRefBox:
      f = 1;
      for(i=0;i<wzPointDim;i++){
	x = (r.x[i] - p[i])/r.near[i];
	x *= x; x *= x; x *= x;
	f *= 1/(1+x*x);
      }
      for(i=0;i<wzPointDim;i++){
	g = f*r.dist[i];
	if(d.g_ii[i]<g) d.g_ii[i] = g;
      }
      break;
    case nRefPoint:
      f = 1;
      for(i=0;i<wzPointDim;i++){
	x = (r.x[i] - p[i])/r.near[i];
	f += x*x;
      }
      for(i=0;i<wzPointDim;i++){
	g = r.dist[i]/f;
	if(d.g_ii[i]<g) d.g_ii[i] = g;
      }
      break;
    }
  }
}