ibgface.cxx

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

#include "wzgrid.hxx"
#include "cogwarnings.hxx"

// if >0, it defines the wzFace-number of outer face cells;
// else, outer face cells should not be created;
int ibgFaceOutside=1;
// if >0, the face number for errorneous face cells;
int ibgErrorneousFace=0;
// if >0, the face number for all face cells;
int ibgSetDefaultFace=0;
int ibgNumberOfFaceCells = 0;

void wzDelaunayGenerator::createFaces(cogeometry g)
{
  int c0,cs,cf,m0,ms,s,i;
  ibg2CellType t,ts,tf;
  int cl0,*cnf,*csf,*css,*cn0;
  if(gdim()<1)	return;
  cl0=Grid.cells.last();
  ibgNumberOfFaceCells = 0;
  for(c0=1;c0<=cl0;c0++) if(ccdef(c0)){
    if((m0=ccu(c0))<=0) continue;
    t=cct(c0);
    for(s=0;s<wzCellTypeSides(t);s++){
      cs=ccs(c0)[s];ts=cct(cs);
      if(cs<=0)		{if (ibgFaceOutside) ms=0; else continue;}
      else if(ccout(cs)){if (ibgFaceOutside) ms=0; else continue;}
      else if(wzCellTypeCODIM(ts)!=wzIsRegion)	continue;
      else if(m0<=(ms=ccu(cs)))		continue;
      switch(wzCellTypeSSize(ts,s)){
      case 1:	tf=wzCellType1Node;		break;
      case 2:	tf=wzCellType2Line;		break;
      case 3:	tf=wzCellType3Triangle;		break;
      case 4:	tf=wzCellType3Rectangle;	break;
      }
      wzAssert(wzCellTypePoints(tf)==wzCellTypeSSize(t,s));
      cf = createCell(tf);
      ibgNumberOfFaceCells++;
      cct(cf)=tf;
      cnf=ccn(cf); csf=ccs(cf); cn0=ccn(c0);
      for(i=0;i<wzCellTypePoints(tf);i++){
	cnf[i] = cn0[wzCellTypeSPoint1(t,s,i)];
      }
      for(i=0;i<wzCellTypeSides(tf);i++) csf[i]=0;
      cnf[wzCellTypeSExtR(tf)]=c0;
      cnf[wzCellTypeSExtL(tf)]=cs;
      ccs(c0)[s]=cf;
      if(cs){
	css=ccs(cs);
	for(i=0;i<wzCellTypeSides(ts);i++) if(css[i]==c0) {css[i]=cf; break;}
      }
      if(ms){
	defineFaceOfCell(cf, g);
      }else{
	defineFaceOutside(cf, g);
      }
    }
  }
  Grid.firstFace= cl0+1;
  Grid.lastFace = Grid.cells.last();
}

void wzDelaunayGenerator::defineFaceOutside(int c, cogeometry g)
{
  wzPoint pp1,pp0,ppo; cogFlag1 f(pp1,pp0,ppo);
  wzIndex i,j;
  wzInteger *cn;
  ibg2CellType t;
  t = cct(c); cn = ccn(c);
  for(j=0;j<wzPointDim;j++) pp1[j]=0;
  for(i=0;i<wzCellTypePoints(t);i++){
    wzFloat* pp = cnx(cn[i]);
    for(j=0;j<wzPointDim;j++) pp1[j] += pp[j];
  }
  for(j=0;j<wzPointDim;j++) pp0[j] = ppo[j] = (pp1[j] /= wzCellTypePoints(t));
  g->Point(pp0); g->Point(ppo);
  g->BoundaryCondition(f);
  ccu(c) = f.p1.segment().face();
}

void wzDelaunayGenerator::defineFaceOfCell(int c, cogeometry g)
{
  int rc,n,nn,cs,f=0,m,i,j,k,nspoints;
  wzPoint xx,yy,*ry,h1,h0,ho;
  cogFlag1 ff(h1,h0,ho);
  ibg2CellType t = cct(c);
  int npoints = wzCellTypePoints(t);
  int *cn=ccn(c),*csn;
  for(n=0;n<npoints;n++){
    switch(cnt(nn=cn[n])){
    case wzIsRegion:
      f = -1; break;
    case wzIsFace:
      if(f){
	if(cnu(nn) != f) f = -1;
      }else{
	f = cnu(nn);
      }
    }
  }
  if(f>0){ccu(c)=f; return;}
  if(f<0){
    cogInfoBadFaceCell();
    if(ibgErrorneousFace){
      ccu(c) = ibgErrorneousFace; return;
    }
  }
  for(k=0;k<wzPointDim;k++)  xx.X[k] = 0;
  for(n=0;n<npoints;n++){
    nn=cn[n]; for(k=0;k<wzPointDim;k++) xx.X[k] += cnx(nn)[k];
  }
  for(k=0;k<wzPointDim;k++)  xx.X[k] /= npoints;
  g->setStartpoint(*cnd(nn)); g->Point(xx);
  m=xx.segment().index();
  for(i=0;i<2;i++){
    cs = cn[wzCellTypeSExt(t)+i];
    wzAssert(ccu(cs)>0);
    if(ccu(cs)==m) continue;
    csn = ccn(cs); nspoints = wzCellTypePoints(cct(cs));
    ry = &yy;
    for(j=0;j<nspoints;j++){
      nn = csn[j];
      for(k=0;k<wzPointDim;k++) yy.X[k] += cnx(nn)[k];
      if(cnt(nn)==wzIsRegion){
	if(cnu(nn) != m){
	  ry = cnd(nn);
	}
      }
    }
    for(k=0;k<wzPointDim;k++)  yy.X[k] /= nspoints;
    if(ry == &yy) g->Point(*ry);
    wzAssert(ry->segment().index() != m);
    rc = g->Line(ff,cogLine(xx,*ry));
    wzAssert(rc == cogRCFaceFound);
    ccu(c) = ff.p1.segment().index();
    return;
  }
  wzAssert(0);
}