ibgotests.cxx

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

#include "ibgoctree.hxx"

/* This file is a useful documentation. It tests most properties which
are supposed to be fulfilled in any valid state of the grid.  */

void ibgOctreeData::t()
{
  ibgiNode n,no,nb; ibgiCell q; int i,j,c,d,k;
  wzRangeLoop(nodes,n){
    for(i=0;i<LinesUp;i++){
      no = nUp[i][n];
      if(is_defined(no)){
	wzAssert(n==nDown[i][no]);
	for(d=0;d<GridDim;d++){
	  if(d==i)	wzAssert(nX[d][no] >nX[d][n]);
	  else		wzAssert(nX[d][no]==nX[d][n]);
	}
      }
    }
    for(i=0;i<Lines;i++){
      for(j=0;j<Lines;j++){
	if(llOrthogonal(i,j)){
	  k = llOrtho[i][j];
	  c = llCell[i][j];
	  q = nCell[c][n];
	  if(q==0){
	    wzAssert(
		     is_outside(nNext[i][n]) ||
		     is_outside(nNext[j][n]) ||
		     is_outside(nNext[k][n])
		     );
	    continue;
	  }
	  nb = cNode[cReverse[c]][q];
	  no = cNode[c][q];
	  for(d=0;d<LinesUp;d++){
	    if(nPoint[n][d]>nPoint[nb][d]){
	      wzAssert(nPoint[no][d]>nPoint[n][d]);
	    }else if(nPoint[n][d]<nPoint[nb][d]){
	      wzAssert(nPoint[no][d]<nPoint[n][d]);
	    }
	  }
	}
      }
    }
  }
  wzRangeLoop(cells,q){
    for(i=0;i<Corners;i++){
      n = cNode[i][q];
      j=cReverse[i];
      wzAssert(q==nCell[j][n]);
    }
  }
}


void ibgOctree::test()
{
  wzInteger n,nb,f,e,d,d1,c,c0,c1,i,j,k,l,p,pr,p1,p2,po,nu,no,nd,t,f1;
  //  return;
  for(i=0;i<Lines;i++){
    for(j=0;j<Lines;j++){
      p=llPlane[i][j];
      wzAssert(p==llPlane[j][i]);
      if(llOrthogonal(i,j)){
	if(GridDim>2){
	  k = llOrtho[i][j];
	  c = llCell[i][j];	// is the cell spanned by i,j,k:
	  wzAssert(c==llCell[j][k]);
	  wzAssert(c==llCell[k][i]);
	  p1 = pSide0[p];	p2 = pSide1[p];	po = pOrtho[p];
	  wzAssert(po==llOrtho[p1][p2]);
	  wzAssert(p==llPlane[p1][p2]);
	}
      }else{
	wzAssert(llParallel(i,j));
	wzAssert((i==j) || (i==Other(j)));
      }
    }
  }
  wzRangeLoop(nodes,n){
    if(nb=nBoundary[n]){
      wzAssert(bnBase[nb]==n);
    }
    for(d=0;d<LinesUp;d++){
      if(is_defined(no=nUp[d][n])){
	wzAssert(n==nDown[d][no]);
	for(d1=0;d1<LinesUp;d1++){
	  if(d1==d){
	    wzAssert(nX[d1][no]>nX[d1][n]);
	  }else{
	    wzAssert(nX[d1][no]==nX[d1][n]);
	  }
	}
      }
      if(is_defined(nu=nDown[d][n])){
	wzAssert(n==nUp[d][nu]);
	for(d1=0;d1<LinesUp;d1++){
	  if(d1==d){
	    wzAssert(nX[d1][nu]<nX[d1][n]);
	  }else{
	    wzAssert(nX[d1][nu]==nX[d1][n]);
	  }
	}
      }
    }
  }
  wzRangeLoop(cells,c){
    for(k=0;k<Corners;k++){
      n=cNode[k][c]; 
      wzAssert(nCell[cReverse[k]][n]==c);
    }
  }
  wzRangeLoop(bnodes,nb){
    n=bnBase[nb];
    wzAssert(nBoundary[n]==nb);
    for(i=0;i<Planes;i++){
      p = bnPlane[i][nb];
      if(p) wzAssert(pDirection(p)==pOrtho[i]);
    }
  }
  wzRangeLoop(lines,l){
    d = lDirection(l);
    nu = lNodeUp[l];
    nd = lNodeDown[l];
    wzAssert(bnUp[d][nd] == l);
    wzAssert(bnDown[d][nu] == l);
    wzAssert(nUp[d][bnBase[nd]] == bnBase[nu]);
  }
  wzRangeLoop(planes,p){
    t  = pDirection(p);
    if(GridDim>2){
      c0 = - pCell[0][p]; if(c0<0) c0 = bcBase[-c0];
      c1 = - pCell[1][p]; if(c1<0) c1 = bcBase[-c1];
    }
    for(i=0;i<4;i++){
      nb = pNode[i][p];
      wzAssert(bnPlane[pr=ptReverse[i][t]][nb]==p);
      if(GridDim>2){
	wzAssert(c0==nCell[llCell[pSide0[pr]][pSide1[pr]]][bnBase[nb]]);
	wzAssert(c1==nCell[llCell[pSide1[pr]][pSide0[pr]]][bnBase[nb]]);
      }
    }
  }
  wzRangeLoop(faces,f){
    e = fLine[f]; 
    d = lDirection(e);
    wzAssert(nPoint[bnBase[lNodeDown[e]]][d]<=fPoint[f][d]);
    wzAssert(nPoint[bnBase[lNodeUp[e]]][d]>=fPoint[f][d]);
    wzAssert(fPointDown[f][d]<fPoint[f][d]);
    wzAssert(fPointUp[f][d]>fPoint[f][d]);
    if((f1= -fDown[f])>0){
      wzAssert(fUp[f1]== -f);
      wzAssert(fLine[f1]==e);
      wzAssert(fPoint[f1][d]<fPoint[f][d]);
    }else{
      wzAssert(lFaceDown[e]==f);
    }
    if((f1=-fUp[f])>0){
      wzAssert(fDown[f1]== -f);
      wzAssert(fLine[f1]==e);
      wzAssert(fPoint[f1][d]>fPoint[f][d]);
    }else{
      wzAssert(lFaceUp[e]==f);
    }
  }
}