ibgoface.cxx

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

#include "ibgoctree.hxx"
#include "cogwarnings.hxx"

wzFloat doubleFaceFactor = 1.0/4;
wzFloat LongSideIntersectionMaximum = 100000;


wzFloat doubleFaceFactor2 = 0.99*doubleFaceFactor*doubleFaceFactor;
wzFloat LongSideIntersectionMaximum2 = LongSideIntersectionMaximum*LongSideIntersectionMaximum;

void ibgOctree::fInitializeMetric(wzIndex face)
{
  cogFlag1 f(fPoint[face],fPointDown[face],fPointUp[face]);
  ref->getMetric(fMetric[face],f);
  fMetricalRefinement(face);
}

void ibgOctree::fMetricalRefinement(wzIndex face)
{
  wzIndex l,no,nu;
  l  = fLine[face];
  no = bnBase[lNodeUp[l]];
  nu = bnBase[lNodeDown[l]];
  if(ref->refineBetween(fMetric[face],nPoint[no],nPoint[nu])){
    if(is_defined(refineBoundaryLine(l))){
      fStackMetric.append(face);
      nStackLinearRegularization.append(nu);
      nStackLinearRegularization.append(no);
      return;
    }
  }
  cogFlag1 f(fPoint[face],fPointDown[face],fPointUp[face]);
  chart->getMetric(fMetric[face],f);
  fStackDoubleIntersection.append(face);
}

void ibgOctree::fRegularization(wzIndex face)
{
  wzIndex l,d,dr,r,u,o,s,t;
  wzInteger no,nu,nro,nru,noo,nou,nh;
  l  = fLine[face]; d = lDirection(l); u = down(d);
  no = bnBase[lNodeUp[l]];
  nu = bnBase[lNodeDown[l]];
  for(s=0;s<LineSides;s++){
    r = lsLine[d][s]; dr = Up(r);
    if(is_nothing(nru=nNext[r][nu])){
      nh = plumbForced(nu,r,d); goto end;
    }else if(is_defined(nru)){
      if(is_nothing(nNext[d][nru])){
	nh=plumbForced(nru,d,r); goto end;
      }
      o = llOrtho[d][r];
      for(t=0;t<PlaneSides;t++){
	if(is_nothing(nou=nNext[o][nru])){
	  nh=plumbForced(nru,o,r); goto end;
	}else if(is_defined(nou)){
	  if(is_nothing(nNext[d][nou])){
	    nh=plumbForced(nou,d,o); goto end;
	  }
	}
	o = Other(o);
      }
    }
    if(is_nothing(nro=nNext[r][no])){
      nh=plumbForced(no,r,d); goto end;
    }else if(is_defined(nro)){
      if(is_nothing(nNext[u][nro])){
	nh=plumbForced(nro,u,r); goto end;
      }
      o = llOrtho[d][r];
      for(t=0;t<PlaneSides;t++){
	if(is_nothing(noo=nNext[o][nro])){
	  nh=plumbForced(nro,o,r); goto end;
	}else if(is_defined(noo)){
	  if(is_nothing(nNext[u][noo])){
	    nh=plumbForced(noo,u,o); goto end;
	  }
	}
	o = Other(o);
      }
    }
  }
  return;
end:
  if(is_defined(nh)){
    fStackRegularization.append(face);	
  }else{
    wzAssert(0);
  }
}

void ibgOctree::fDoubleIntersection(wzIndex face)
{
  wzIndex l,d,no,nu;
  wzFloat dx,l2;
  // at least two intersections of the line with faces.
  if(fUp[face]<0 || fDown[face]<0){
    l  = fLine[face]; d = lDirection(l);
    no = bnBase[lNodeUp[l]];
    nu = bnBase[lNodeDown[l]];
    dx = nPoint[no][d]-nPoint[nu][d];
    wzAssert(dx>0);
    l2 = ref->g_ii(fMetric[face],fPoint[face],d)*dx*dx;
    if(l2>doubleFaceFactor2){
      if(is_defined(refineBoundaryLine(l))){
	fStackDoubleIntersection.append(face);
	return;
      }
    }
  }
  fStackRegularization.append(face);
}

void ibgOctree::fLongSideIntersection(wzIndex ff)
{
  wzIndex ll,ls,fs,bu,bo; //,bru,bro,bou,boo;
  wzInteger nu,no,nru;
  wzIndex d,u,r,dr,s;
  wzFloat dx,dx2,dy,dy2,df;
  wzSegment seg = fPoint[ff].segment();
  ll = fLine[ff];	 	
  d  = lDirection(ll);	u  = down(d);
  bo = lNodeUp[ll];	no = bnBase[bo];
  bu = lNodeDown[ll];	nu = bnBase[bu];
  dx = nPoint[no][d] - nPoint[nu][d];
  df = (fPoint[ff][d] - nPoint[nu][d])/dx;
  if(df<0.001 || df > 0.999) goto end;
  dx *= LongSideIntersectionFactor;
  dx2 = chart->g_ii(fMetric[ff],fPoint[ff],d)*dx*dx;
  for(s=0;s<LineSides;s++){
    r = lsLine[d][s]; dr = Up(r);
    if(is_defined(nru=nNext[r][nu])){
      dy = nPoint[nru][dr] - nPoint[nu][dr];
      dy2 = chart->g_ii(fMetric[ff],fPoint[ff],dr)*dy*dy;
      if(dy2<dx2){
	if(LongSideIntersectionMaximum2*dy2>dx2) goto OnLongSide;
      }
    }
    continue;
  }
  goto end;
OnLongSide:
  for(s=0;s<LineSides;s++){
    r = lsLine[d][s]; dr = Up(r);
    ls = bnNext[r][bu];
    if(ls && (fs=lFaceUp[ls])){
      if(fPoint[fs].segment()==seg) goto end;
    }
    ls = bnNext[r][bo];
    if(ls && (fs=lFaceUp[ls])){
      if(fPoint[fs].segment()==seg) goto end;
    }
  }
  if(is_defined(refineBoundaryLine(ll))){
    fStackDoubleIntersection.append(ff);
    return;
  }
end:
  fDoubleIntersection(ff);
  return;
}

ibgIndex ibgOctree::createSimpleBoundaryNode(ibgIndex node)
{
  wzIndex k,nb;
  wzAssert(nBoundary[node]==0);
  nb=nBoundary[node]=bnodes.create();
  bnBase[nb]=node;
  bnType[nb]=0;
  for(k=0;k<Lines;k++) bnNext[k][nb]=0;
  for(k=0;k<Planes;k++) bnPlane[k][nb]=0;
  return nb;
}

ibgIndex ibgOctree::createSimpleLine(ibgIndex bn0,ibgIndex bn1,ibgIndex d)
{
  wzIndex k,ll = lines.create();
  if(isUp(d)){
    lNodeUp[ll] = bn1; bnDown[d][bn1] = ll;
    lNodeDown[ll] = bn0; bnUp[d][bn0] = ll;
    wzAssert(nUp[d][bnBase[bn0]] == bnBase[bn1]);
  }else{
    d = up(d);
    lNodeUp[ll] = bn0; bnDown[d][bn0] = ll;
    lNodeDown[ll] = bn1; bnUp[d][bn1] = ll;
    wzAssert(nUp[d][bnBase[bn1]] == bnBase[bn0]);
  }
  lType[ll]=d;
  lFaceUp[ll]=lFaceDown[ll]=0;
  for(k=0;k<LineSides;k++) lPlane[k][ll]=0;
  return ll;
}

void ibgOctree::findFaces()
{
  wzIndex f,n,rc;
  freeFace = faces.create();
begin:
  rc = 0;
  while(1){
    if(n=nStackInitial.pop()){
      nInitializeMetric(n); continue;
    }
    if(n=nStackMetric.pop()){
      nMetricalRefinement(n); continue;
    }
    if(n=nStackOrthogonalRegularization.pop()){
      nOrthogonalRegularization(n); continue;
    }
    if(n=nStackLinearRegularization.pop()){
      nLinearRegularization(n); continue;
    }
    if(n=nStackFace.pop()){
      nFindFaces(n); continue;
    }
    if(f=fStackInitial.pop()){
      fInitializeMetric(f); continue;
    }
    if(f=fStackMetric.pop()){
      fMetricalRefinement(f); continue;
    }
    if(f=fStackDoubleIntersection.pop()){
      fLongSideIntersection(f); continue;
      //      fDoubleIntersection(f); continue;
    }
    if(f=fStackRegularization.pop()){
      fRegularization(f); continue;
    }
    break;
  }
  wzRangeLoop(nodes,n){
    if(nFindFaces(n)) rc=1;
  }
  if(rc) goto begin;
  faces.destroy(freeFace);
}

wzIndex ibgOctree::lFindFaces(wzIndex ll)
{
  wzIndex bnu = lNodeDown[ll],nu=bnBase[bnu];
  wzIndex bno = lNodeUp[ll]  ,no=bnBase[bno];
  wzSegment su = nPoint[nu].segment();
  wzSegment so = nPoint[no].segment();
  wzIndex d   = lDirection(ll),j,k,rc=0;
  wzInteger fu  = lFaceDown[ll],fo,fb,fh;
  wzAssert(fu>=0);
  if(fu==0){ // nothing on the line:
    if(su==so) return rc;
    cogFlag1 f(fPoint[freeFace],fPointDown[freeFace],fPointUp[freeFace]);
    rc=geo->Line(f,cogLine(nPoint[nu],nPoint[no]));
    wzAssert(rc==cogRCFaceFound);
    fu = freeFace; freeFace = faces.create();
    fStackInitial.append(fu); rc=1;
    fLine[fu] = ll; lFaceDown[ll] = lFaceUp[ll] = fu;
    fUp[fu]   = bno; fDown[fu] = bnu;
    for(k=0;k<LineSides;k++) {fEdge[k][fu]=0;}
  }else if(su != fPointDown[fu].segment()){
    cogFlag1 f(fPoint[freeFace],fPointDown[freeFace],fPointUp[freeFace]);
    rc=geo->Line(f,cogLine(nPoint[nu],fPointDown[fu]));
    wzAssert(rc==cogRCFaceFound);
    fb = fu; fu = freeFace; freeFace = faces.create();
    fStackInitial.append(fu); rc=1;
    fLine[fu] = ll; lFaceDown[ll] = fu;
    fUp[fu]   = -fb; fDown[fb] = -fu; fDown[fu] = bnu;
    for(k=0;k<LineSides;k++) {fEdge[k][fu]=0;}
  }
  fo = lFaceUp[ll]; wzAssert(fo>0);
  if(so != fPointUp[fo].segment()){
    cogFlag1 f(fPoint[freeFace],fPointDown[freeFace],fPointUp[freeFace]);
    rc=geo->Line(f,cogLine(fPointUp[fo],nPoint[no]));
    wzAssert(rc==cogRCFaceFound);
    fb = fo; fo = freeFace; freeFace = faces.create();
    fStackInitial.append(fo); rc=1;
    fLine[fo] = ll; lFaceUp[ll] = fo;
    fUp[fb]   = -fo; fDown[fo] = -fb; fUp[fo] = bno;
    for(k=0;k<LineSides;k++) {fEdge[k][fo]=0;}
  }
  for(j=0;j<20;j++){
    if(fu==fo) return rc;
    fb = -fUp[fu]; wzAssert(fb>0); 
    if(fPointUp[fu].segment() != fPointDown[fb].segment()){
      cogFlag1 f(fPoint[freeFace],fPointDown[freeFace],fPointUp[freeFace]);
      rc=geo->Line(f,cogLine(fPointUp[fu],fPointDown[fb]));
      wzAssert(rc==cogRCFaceFound);
      fh = fu; fu = freeFace; freeFace = faces.create();
      fStackInitial.append(fu); rc=1;
      fLine[fu] = ll;
      fUp[fh]   = -fu; fDown[fu] = -fh; 
      fUp[fu]   = -fb; fDown[fb] = -fu;
      for(k=0;k<LineSides;k++) {fEdge[k][fu]=0;}
    }else{
      fu = fb;
    }
  }
  wzAssert(0);
  return rc;
}

wzIndex ibgOctree::nFindFaces(ibgiNode n)
{
  wzIndex d,rc=0;
  wzSegment s=nPoint[n].segment(),so;
  wzInteger nb=nBoundary[n],nob,no,ll;
  for(d=0;d<Lines;d++){
    if(is_undefined(no=nNext[d][n]))	continue;
    if(s==(so=nPoint[no].segment())) 	continue;
    if(!(nb))			nb  = createSimpleBoundaryNode(n);
    if(!(ll=bnNext[d][nb])){
      if(!(nob=nBoundary[no]))	nob = createSimpleBoundaryNode(no);
      ll  = createSimpleLine(nb,nob,d);
    }
    if(lFindFaces(ll)) rc=1;
  }
  return rc;
}