ibgorefine.cxx

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

#include "ibgoctree.hxx"
#include "cogwarnings.hxx"
#define is_undefined(x) ((x)<=0)
#define is_defined(x) ((x)>0)
#define is_outside(x) ((x)==ibgOctreeData::outside)
#define is_nothing(x) ((x)<0)

wzFloat ibggOref=0.49;
wzFloat ibggOinv=1.0/ibggOref;
wzFloat ibggOref2=ibggOref*ibggOref;
//static int make_compiler_happy_1 = 1;
wzFloat ibggPref=2.01;
wzIndex ibggPrefL=1;

ibgiNode ibgOctree::refineBoundaryLine(ibgIndex l)
{
  ibgIndex bnu = lNodeDown[l], d=lDirection(l);
  return refine(bnBase[bnu],d);
}

void ibgOctree::nInitializeMetric(ibgiNode nn)
{
  ref->getMetric(nMetric[nn],nPoint[nn]);
  nMetricalRefinement(nn);
}

void ibgOctree::nMetricalRefinement(ibgiNode nn)
{
  wzIndex rf,d;
  ibgiNode no;
  rf = 0;
  for(d=0;d<Lines;d++){
    if(is_undefined(no=nNext[d][nn]))	continue;
    if(ref->refineBetween(nMetric[nn],nPoint[nn],nPoint[no])) 
      if(is_defined(refine(nn,d))){
	rf=1;
	nStackLinearRegularization.append(no);
      }
  }
  if(rf){
    nStackMetric.append(nn);
  }else{
    chart->getMetric(nMetric[nn],nPoint[nn]);
    nStackOrthogonalRegularization.append(nn);
  }
}

void ibgOctree::nOrthogonalRegularization(ibgiNode nn)
{
  wzIndex dd,dr,ro,rr,ru,s;
  ibgiNode no,nr;
  wzFloat dx,dx2,dy,dy2;
  for(ro=0;ro<Lines;ro++){
    if(is_undefined(no=nNext[ro][nn])) continue;
    dd  = Up(ro); ru=Other(ro);
    dx  = nX[dd][no]-nX[dd][nn];
    dx2 = chart->g_ii(nMetric[nn],nPoint[nn],dd)*dx*dx;
    for(s=0;s<LineSides;s++){
      rr=lsLine[dd][s];
      if(is_undefined(nr=nNext[rr][no]))	continue;
      if(is_undefined(nNext[rr][nn])){
	dr  = Up(rr);
	dy  = (nX[dr][nr]-nX[dr][no]);
	dy2 = chart->g_ii(nMetric[nn],nPoint[nn],dr)*dy*dy/dx2;
	if(is_undefined(nNext[ru][nr])){
	  if(dy2>4.001) {
	    if (is_defined(plumbForced(nr,ru,rr))) 
	      nStackLinearRegularization.append(nr);
	  }else if(dy2<0.2499){
	    plumbForced(nn,rr,ro);
	  }
	}else{
	  if(dy2<0.999){
	    plumbForced(nn,rr,ro);
	  }
	}
      }else if(is_undefined(nNext[ru][nr])){
	dr  = Up(rr);
	dy  = (nX[dr][nr]-nX[dr][no]);
	dy2 = chart->g_ii(nMetric[nn],nPoint[nn],dr)*dy*dy/dx2;
	if(dy2>1.001){
	  if (is_defined(plumbForced(nr,ru,rr))) 
	    nStackLinearRegularization.append(nr);
	}
      }
    }
  }
  nStackLinearRegularization.append(nn);
}

void ibgOctree::nLinearRegularization(ibgiNode nn)
{
  wzIndex dd;
  ibgiNode no,nu;
  wzFloat d1,d2,x;
  //  goto end;
  for(dd=0;dd<GridDim;dd++){
    if(is_undefined(no=nUp[dd][nn])) continue;
    if(is_undefined(nu=nDown[dd][nn])) continue;
    d1 = nX[dd][no] - (x=nX[dd][nn]);
    d2 = (x - nX[dd][nu])/d1;
    if(d2>2.01){
      if(is_defined(refine(nu,dd))){
	nStackLinearRegularization.append(nn);
	nStackLinearRegularization.append(nu);
      }
    }else if (d2<0.499){
      if(is_defined(refine(nn,dd))){
	nStackLinearRegularization.append(nn);
	nStackLinearRegularization.append(no);
      }
    }
  }
end:
  nStackFace.append(nn);
}

/*
void ibgOctree::refineOld(ibgiNode nn)
{
  if(make_compiler_happy_1) return;
  for(d=0;d<Lines;d++){
    if(is_defined(no=nNext[d][nn])){
      for(o=0;o<LineSides;o++){
	rv=lsLine[d][o];
	if(is_defined(nv=nNext[rv][no])){
	  if(nPoint[nv].segment() != nPoint[nn].segment()){
	    if(is_undefined(nNext[rv][nn])) plumb(nn,rv,d);
	  }
	}
	if(nPoint[no].segment() != nPoint[nn].segment()){
       	  if(is_undefined(nNext[rv][nn])) plumb(nn,rv,d);
	}
      }
    }
  }
  nStackFace.append(nn);
}
*/
ibgiCell ibgOctree3D::qref2(ibgiCell qq, unsigned ro)
{
 int nu,no,nn,np,nur,nor,ou,oo;
 int o,ru=Other(ro);
 if(is_undefined(qq)) return 0;
 nn = cNode[firstCorner][qq];		np = cNode[lastCorner][qq];
 nu = cNode[(lsCell[ru][3])][qq];	no = cNode[(int)lsCell[ro][3]][qq];
 for(o=0;o<LineSides;o++){
 	ou = nu; nu = cNode[lsCell[ru][o]][qq];
 	oo = no; no = cNode[lsCell[ro][o]][qq];
 	if(no==nNext[ro][nu]){
 		refine(nu,ro);
 	 	if(nn!=cNode[firstCorner][qq]) return 1;
 		if(np!=cNode[lastCorner][qq]) return 1;
 	}
 	if((nor=nNext[lsLine[ro][o]][oo])==no) continue;
 	if((nur=nNext[lsLine[ru][o]][ou])==nu) continue;
        if(nor==nothing) return 0;
	if(nur==nothing) return 0;
  	if(nor==nNext[ro][nur]){
 		refine(nur,ro);
 	 	if(nn!=cNode[firstCorner][qq]) return 1;
 		if(np!=cNode[lastCorner][qq]) return 1;
 	}
 }
 return splitCell(qq,ro);
}

ibgiCell ibgOctree3D::qref1(ibgiCell qq, unsigned ro)
{
 int nu,no,nur,nor,ou,oo;
 int o,ru=Other(ro);
 if(is_undefined(qq)) return 0;
 nu = cNode[(lsCell[ru][3])][qq];  no = cNode[(int)lsCell[ro][3]][qq];
 for(o=0;o<LineSides;o++){
 	ou = nu; nu = cNode[lsCell[ru][o]][qq];
 	oo = no; no = cNode[lsCell[ro][o]][qq];
        if(no!=nNext[ro][nu])	return qref2(qq,ro);
        if((nor=nNext[lsLine[ro][o]][oo])==no) continue;
        if((nur=nNext[lsLine[ru][o]][ou])==nu) continue;
        if(nor!=nNext[ro][nur])	return qref2(qq,ro);
 }
 return 0;
}


ibgiNode ibgOctree3D::plumb(ibgiNode nn, ibgdLine rr, ibgdLine ro)
{
 int no,nu,nor,nur,noru,norr,nurr;
 wzAssert(is_undefined(nNext[rr][nn]));
 no =nNext[ro][nn];		if(is_undefined(no))	return nothing;
 nu =nNext[Other(ro)][nn];	if(is_undefined(nu))	return nothing;
 nor=nNext[rr][no];		if(is_undefined(nor))	return nothing;
 nur=nNext[rr][nu];		if(is_undefined(nur))	return nothing;
 noru=nNext[Other(ro)][nor];
 if(is_undefined(noru)){
#ifdef GLEVEL
 	dd=posit(rr); if(nLevel[dd][no]>=nLevel[dd][nor]) return nothing;
#endif
	norr=nNext[rr][nor];
	if(is_undefined(norr))		return nothing;
	noru=nNext[Other(ro)][norr];
	if(is_undefined(noru))	return nothing;
	if(noru==nur)		return refine(noru,ro);
	if(noru==nNext[rr][nur]){
				return nothing;
	}
	nurr=nNext[rr][nur];
	if(is_undefined(nurr))		return nothing;
	if(noru==nNext[ro][nurr]){
			return nothing;
	}
	else 			return nothing;
 }
 if(noru==nur)			return refine(noru,ro);
 if(noru==nNext[rr][nur])		return refine(noru,ro);
 if(noru==nNext[ro][nur]){
   // 	nNext[rr][nn] = noru;
   // 	nNext[Other(rr)][noru] = nn;
				return noru;
 }
#ifdef GLEVEL
 dd=posit(rr); if(nLevel[dd][nu]>=nLevel[dd][nur]) return nothing;
#endif
 nurr=nNext[rr][nur];
 if(is_undefined(nurr))			return nothing;
 if(noru==nNext[ro][nurr]){
   // 	nNext[rr][nn] = noru;
   //	nNext[Other(rr)][noru] = nn;
			   	return noru;
 }
 else	 			return nothing;
}

ibgiNode ibgOctree3D::plumbForced(ibgiNode nn, ibgdLine rr, ibgdLine ro)
{int rc,q0;
 if(is_defined(rc=plumb(nn,rr,ro)))	return rc;
 q0=nCell[llCell[rr][ro]][nn];
 wzAssert(q0==nCell[llCell[ro][rr]][nn]);
 qref2(q0,llOrtho[ro][rr]);
 if(is_defined(rc=nNext[rr][nn]))	return rc;
 if(is_defined(rc=plumb(nn,rr,ro)))	return rc;
 if(is_defined(rc=plumb(nn,rr,llOrtho[ro][rr])))	return rc;
 wzFailed(); 
 return 0;
}

ibgiNode	ibgOctree3D::refine(ibgiNode nu, ibgdLineUp ro)
{
  ibgiCell qq,qo,o,d,dh,eo,e1,rr,ru;
  int no,nh,nur,nurr,nor,nr,lo;
  if(isUp(ro))	{d=ro;	  ru=down(ro);no=nNext[ro][nu];}
  else		{d=up(ro);ru=ro;ro=d; no=nu;nu=nNext[ru][no];}
  wzAssert(is_defined(no)); wzAssert(is_defined(nu));
  for(o=0;o<LineSides;o++){
    if(is_undefined(qq= nCell[lo=lsCell[ro][o]][nu]))continue;
    if     (qq==nCell[lo][no])
      {if(!qref2(qq,ro)) return nothing;}
    else if(qq==nCell[lsCell[ru][o]][nu])
      {if(!qref2(qq,ro)) return nothing;}
  }
  if((nX[d][no]-nX[d][nu])<=dmin[d]){
    cogInfoMaximalRefinementReached(); return nothing;
  }
  // now, internal node creation part has started:
  nh = nodes.create();
  nNext[ro][nu] = nh; nNext[ru][nh] = nu; nNext[ru][no] = nh; nNext[ro][nh] = no;
  for(dh=0;dh<XDim;dh++){
    if(d!=dh){
      nX[dh][nh] = nX[dh][nu];
      nUp[dh][nh] = nDown[dh][nh] = nothing;
#ifdef GLEVEL
      nLevel[dh][nh] = nLevel[dh][nu];
#endif
   }else{
     nX[dh][nh] = 0.5*(nX[dh][nu]+nX[dh][no]);
#ifdef GLEVEL
     nLevel[dh][nh] = nLevel[dh][nu]+1;
     if(nLevel[dh][nh]<=nLevel[dh][no]) nLevel[dh][nh] = nLevel[dh][no]+1;
#endif
   }
  }
  geo->setStartpoint(nPoint[nu]); geo->Point(nPoint[nh]);
  nOld[nh]=nu;
  nStackInitial.append(nh);
  for(o=0;o<LineSides;o++){
    eo = lsCell[ro][o];
    nCell[eo][nh] = (nCell[lsCell[ru][o]][nh] = nCell[eo][nu]);
  }
  if(nBoundary.allocated()!=0) {splitBoundaryLine(nu,d,nh);}

  for(o=0;o<LineSides;o++){
    rr = lsLine[ro][o];
    if(is_undefined(nur=nNext[rr][nu]))		{nNext[rr][nh]=  nur;continue;}
    if(is_undefined(nr=nNext[ro][nur])){
      if(nCell[lsCell[ro][o]][nh ] !=
	 nCell[lsCell[ro][o]][nur])		{nNext[rr][nh]= -nur;continue;}
      if(is_undefined(nurr=nNext[rr][nur]))	{nNext[rr][nh]= -nur;continue;}
      if(is_undefined(nr  =nNext[ro][nurr])) 	{nNext[rr][nh]= -nur;continue;}
    }
    if(is_undefined(nor=nNext[rr][no]))		{nNext[rr][nh]=  nor;continue;}
    if(nr==nor)			 		{nNext[rr][nh]= -nor;continue;}
    if(nr==nNext[rr][nor])			{nNext[rr][nh]= -nor;continue;}
    if((nNext[ro][nr]!=nor)&&(nNext[ro][nr]!=nNext[rr][nor]))
      {nNext[rr][nh]= -nor;continue;}
    eo=lsCell[ro][o];	e1=lsCell[ro][o+1];
    if(nCell[eo][nu]==nCell[e1][nu])		{nNext[rr][nh]= -nor;continue;}
    nNext[Other(rr)][nr] = nh; nNext[rr][nh] = nr;
  }
  qq=nCell[lsCell[ro][3]][nh]; 
  for(o=0;o<LineSides;o++){
    qo=qq;
    if(qo==(qq=nCell[lsCell[ro][o]][nh])) continue;
    if(splitCell(qq,ro))	qq=nCell[lsCell[ro][o]][nh];
  }
  // t();
  return nh;
}


ibgiNode	ibgOctree1D::refine(ibgiNode nu, ibgdLineUp ro)
{int d,ru;
 int no,nh;
 if(isUp(ro))	{d=ro;		ru=ro+GridDim;	no=nNext[ro][nu];}
 else		{d=ro-GridDim;	ru=ro;ro=d;	no=nu;nu=nNext[ru][no];}
 if((nX[d][no]-nX[d][nu])<=dmin[d]){cogInfoMaximalRefinementReached(); return nothing;}

 nh = nodes.create();
 nNext[0][nh] = nNext[1][nh] = nothing;
 nNext[ro][nu] = nh; nNext[ru][nh] = nu; nNext[ru][no] = nh; nNext[ro][nh] = no;

 nX[1][nh] = nX[1][nu];
 nX[2][nh] = nX[2][nu];
 nX[0][nh] = 0.5*(nX[0][nu]+nX[0][no]);
#ifdef GLEVEL
 nLevel[0][nh] = nLevel[0][nu]+1;
 if(nLevel[0][nh]<=nLevel[0][no]) nLevel[0][nh] = nLevel[0][no]+1;
#endif

 geo->setStartpoint(nPoint[nu]); geo->Point(nPoint[nh]);
 nOld[nh]=nu;
 if(nBoundary.allocated()!=0) {splitBoundaryLine(nu,d,nh);}
 nStackInitial.append(nh);

 return nh;
}

ibgiNode	ibgOctree2D::refine(ibgiNode nu, ibgdLineUp ro)
{
  int o,u,d,rr,ru,no,nh,nur,nor,nr,nt;
  // ibgFloat dx,dy;
 if(isUp(ro))	{d=ro;		ru=ro+GridDim;	no=nNext[ro][nu];}
 else		{d=ro-GridDim;	ru=ro;ro=d;	no=nu;nu=nNext[ru][no];}
 if((nX[d][no]-nX[d][nu])<=dmin[d])
   {cogInfoMaximalRefinementReached(); return nothing;}
 o=lsLine[d][0];
 if(is_nothing(nNext[o][nu])) nt=plumb(nu,o,d);
 wzAssert(is_defined(nt));
 if(is_nothing(nNext[o][no])) nt=plumb(no,o,d);
 wzAssert(is_defined(nt));
 u=Other(o);
 if(is_nothing(nNext[u][nu])) nt=plumb(nu,u,d);
 wzAssert(is_defined(nt));
 if(is_nothing(nNext[u][no])) nt=plumb(no,u,d);
 wzAssert(is_defined(nt));

 nh = nodes.create();
 nNext[0][nh] = nNext[1][nh] = nNext[2][nh] = nNext[3][nh] = nothing;
 nNext[ro][nu] = nh; nNext[ru][nh] = nu; nNext[ru][no] = nh; nNext[ro][nh] = no;

 nX[o][nh] = nX[o][nu];
 nX[2][nh] = nX[2][nu];
 nX[d][nh] = 0.5*(nX[d][nu]+nX[d][no]);
#ifdef GLEVEL
 nLevel[o][nh] = nLevel[o][nu];
 nLevel[2][nh] = 0;
 nLevel[d][nh] = nLevel[d][nu]+1;
 if(nLevel[d][nh]<=nLevel[d][no]) nLevel[d][nh] = nLevel[d][no]+1;
#endif

 geo->setStartpoint(nPoint[nu]); geo->Point(nPoint[nh]);
 if(nBoundary.allocated()!=0) {splitBoundaryLine(nu,d,nh);}
 nOld[nh]=nu;
 nStackInitial.append(nh);

 for(o=0;o<LineSides;o++){
   rr=lsLine[ro][o];
   nur=nNext[rr][nu];
   if(nur==outside)	nNext[rr][nh] = outside;
   else{
     wzAssert(is_defined(nur));
     nr=nNext[ro][nur];
     if(is_nothing(nr)) {
       nur=nNext[rr][nur];	wzAssert(is_defined(nur));
       nr =nNext[ro][nur];	wzAssert(is_defined(nr));
     }
     nor=nNext[rr][no]; 	wzAssert(is_defined(nor));
     if(nr==nor)			nNext[rr][nh] = -nr;
     else if(nr==nNext[rr][nor])	nNext[rr][nh] = -nr;
     else {
       wzAssert((nNext[ro][nr]==nor)||(nNext[ro][nr]==nNext[rr][nor]));
       nNext[Other(rr)][nr] = nh;	nNext[rr][nh] = nr;
     }
   }
 }
 return nh;
}

ibgiNode ibgOctree2D::plumb(ibgiNode nn, ibgdLine rr, ibgdLine ro)
{
 int no,nu,nor,nur,noru,norr,nurr;
 wzAssert(is_undefined(nNext[rr][nn]));
 no =nNext[ro][nn];		if(is_undefined(no))	return nothing;
 nu =nNext[Other(ro)][nn];	if(is_undefined(nu))	return nothing;
 nor=nNext[rr][no];		if(is_undefined(nor))	return nothing;
 nur=nNext[rr][nu];		if(is_undefined(nur))	return nothing;
 noru=nNext[Other(ro)][nor];
 if(is_undefined(noru)){
#ifdef GLEVEL
 	dd=Up(rr); if(nLevel[dd][no]>=nLevel[dd][nor]) return nothing;
#endif
	norr=nNext[rr][nor];
	if(is_undefined(norr))		return nothing;
	noru=nNext[Other(ro)][norr];
	if(is_undefined(noru))	return nothing;
	if(noru==nur)		return refine(noru,ro);
	if(noru==nNext[rr][nur]){
		nNext[Other(ro)][nor] = nur;
		nNext[ro][nur]	  = nor;
				return refine(nur,ro);
	}
	nurr=nNext[rr][nur];
	if(is_undefined(nurr))		return nothing;
	if(noru==nNext[ro][nurr]){
		nNext[Other(ro)][nor] = nur;
		nNext[ro][nur]	  = nor;
				return refine(nur,ro);
	}
	else 			return nothing;
 }
 if(noru==nur)			return refine(noru,ro);
 if(noru==nNext[rr][nur])		return refine(noru,ro);
 if(noru==nNext[ro][nur]){
 	nNext[rr][nn] = noru;
 	nNext[Other(rr)][noru] = nn;
				return noru;
 }
#ifdef GLEVEL
 dd=Up(rr); if(nLevel[dd][nu]>=nLevel[dd][nur]) return nothing;
#endif
 nurr=nNext[rr][nur];
 if(is_undefined(nurr))		        return nothing;
 if(noru==nNext[ro][nurr]){
 	nNext[rr][nn] = noru;
	nNext[Other(rr)][noru] = nn;
				return noru;
 }
 else	 			return nothing;
}