cogengfz.cxx

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

#include "cogengfz.hxx"
#include "cogwarnings.hxx"
#include <stdio.h>
#include <math.h>
extern int ibgErrorneousRegion;

void CogenGFZ::setIndexBox(wzIndex mb, wzIndex me, wzIndex yb,wzIndex ye, wzIndex zb,wzIndex ze)
{
  bz=zb;ez=ze,by=yb;ey=ye;bm=mb;em=me;
  if(em>=nx) em = nx-1;
  if(ey>=ny) ey = ny-1;
  if(ez>=nz) ez = nz-1;
  if(bm> em) bm = 0;
  if(by> ey) by = 0;
  if(bz> ez) bz = 0;
  zmin = zz[bz]; zmax=zz[ez];
  ymin = yy[by]; ymax=yy[ey];
}

void	CogenGFZ::endInitialization()
{
  wzIndex bb,fb,bc;
  x0 = 0;
  y0 = ymin = yy[by];
  z0 = zmin = zz[bz];
  dx = 0.01*(xmax-xmin);
  if(ey>by) dy = ((ymax=yy[ey])-y0)/(ey-by); else dy=1;
  if(ez>bz) dz = ((zmax=zz[ez])-z0)/(ez-bz); else dz=1;
  for(bb=bc=ey-by,fb=2;fb<ey-by;fb*=2){
    bc = bc/2;
    if(fb*bc<bb){
      bb = fb*bc;
      addPlaneY(yy[by+bb]);
    }
  }
  for(bb=bc=ez-bz,fb=2;fb<ez-bz;fb*=2){
    bc = bc/2;
    if(fb*bc<bb){
      bb = fb*bc;
      addPlaneZ(zz[bz+bb]);
    }
  }
  setBorder(xmin-dx,xmax+dx,ymin,ymax,zmin,zmax);
  dx = 1;
  ixold = bm;
  CogenOctree::endInitialization();
}

wzgrid CogenGFZ::generateGrid(cogeometry g,wzpoints plist) const
{
  plist->computeBox();
  wzdelaunayGenerator gen = new wzHighlyAnisotropicalDelaunayGenerator(plist->boxDim,plist->spaceDim);
  gen->Delaunay(plist,g);
  return gen->grid();
}

wzIndex	CogenGFZ::Point(wzPoint& p0) const
{
 cogFloat px = (p0.x()-x0)/dx;
 cogFloat ry = (p0.y()-y0)/dy; 
 cogFloat rz = (p0.z()-z0)/dz; 
 cogInteger ix;
 cogInteger iy = (cogInteger) ry, jy = iy+1; ry -= iy;
 cogInteger iz = (cogInteger) rz, jz = iz+1; rz -= iz;
 cogFloat sz = 1-rz, sy = 1-ry, u0;
 if(iy<0){
        iy = jy = 0;
 }else if(jy>=ny){
        iy = jy = ny-1;
 }
 if(iz<0){
        iz = jz = 0;
 }else if(jz>=nz){
        iz = jz = nz-1;
 }
 wzIndex nii = nx*iy+np*iz;
  wzIndex nij = nx*iy+np*jz;
  wzIndex nji = nx*jy+np*iz;
 wzIndex njj = nx*jy+np*jz;
 ix = ixold;
 u0 = sz*(sy*xx[ix+nii] + ry*xx[ix+nji]) + rz*(sy*xx[ix+nij] + ry*xx[ix+njj]);
 if(px>u0){
   for(ix=ixold-1;ix>=bm;ix--){
     u0 = sz*(sy*xx[ix+nii] + ry*xx[ix+nji]) + rz*(sy*xx[ix+nij] + ry*xx[ix+njj]);
     if(px<=u0){
       ((CogenGFZ*)this)->ixold = ix;
       p0.segment() = wzRegion(ix+2);
       return cogRCRegionFound;
     }
   }
   ((CogenGFZ*)this)->ixold = bm;
   p0.segment() = wzRegion(bm+1);
   return cogRCRegionFound;
 }else{
   for(ix=ixold+1;ix<=em;ix++){
     u0 = sz*(sy*xx[ix+nii] + ry*xx[ix+nji]) + rz*(sy*xx[ix+nij] + ry*xx[ix+njj]);
     if(px>u0){
       ((CogenGFZ*)this)->ixold = ix;
       p0.segment() = wzRegion(ix+1);
       return cogRCRegionFound;
     }
   }
   ((CogenGFZ*)this)->ixold = em;
   p0.segment() = wzRegion(em+2);
   return cogRCRegionFound;
 }
}

cogIndex CogenGFZ::BoundaryCondition(cogFlag1& f) const
{
  wzIndex s0 = f.p0.segment().region();
  wzIndex s1 = f.po.segment().region();
  if(s0>s1){
    f.p1.segment() = wzFace(s0);
  }else{
    f.p1.segment() = wzFace(s1);
  }
  return cogRCConditionFound;
}

CogenGFZ::CogenGFZ(wzString filename)
{
  int rc;
  wzIndex cc,iz,iy,ix,count=0;
  float z,y,x;
  wzFloat x1,xb;
  dmmin = 40;
  const wzIndex buflen=256;
  char buffer[buflen];
  FILE *file = fopen(filename,"r"); 
  if(!file){throw wzFileOpenError();}
  fgets(buffer,buflen,file);
  rc = sscanf(buffer,"%d %d %d %d",&nz,&ny,&nx,&cc);
  zz.available(nz);
  yy.available(ny);
  xx.available(nz*ny*nx);
  np = nx*ny;
  fgets(buffer,buflen,file);
  if(rc<4) {throw wzFileReadError();}
  if(cc<3) {throw wzFileReadError();}
  xmin = wzInfty; xmax = -wzInfty;
  for(ix=0;ix<nx;ix++){
    for(iy=0;iy<ny;iy++){
      for(iz=0;iz<nz;iz++){
	if(! fgets(buffer,buflen,file)){throw wzFileReadError();}
	rc = sscanf(buffer,"%f %f %f",&z,&y,&x);
	if(rc<3){throw wzFileReadError();}
	xx[ix+nx*iy+np*iz] = x;
	if(xmin>x) xmin=x;
	if(xmax<x) xmax=x;
	if(iy || ix){
	  wzAssert(z==zz[iz]);
	}else{
	  zz[iz] = z;
	}
	if(iz){
	  wzAssert(y==y0);
	}else{
	  y0 = y;
	  if(ix){
	    wzAssert(y0==yy[iy]);
	  }else{
	    yy[iy] = y0;
	  }
	}
      }
    }
  }
  // remove the 0.1 minimal artificial distances 
  for(iz=0;iz<nz;iz++){
    for(iy=0;iy<ny;iy++){
      x0 = xb = xx[nx*iy+np*iz];
      for(ix=1;ix<nx;ix++){
	x1 = xx[ix+nx*iy+np*iz];
	wzAssert(x1<x0-0.099);
	if(x1>x0-dmmin){
	  xx[ix+nx*iy+np*iz] = xb;
	  x0 = x1;
	  count++;
	}else{
	  x0 = xb = x1;
	}
      }
    }
  }
  ixold = 0;
  bz = 0; ez = nz-1;
  by = 0; ey = ny-1;
  bm = 0; em = nx-1;
  zmin = zz[bz]; zmax=zz[ez];
  ymin = yy[by]; ymax=yy[ey];
}

void 	wzHighlyAnisotropicalDelaunayGenerator::defineFaceOfCell(int c, cogeometry g)
{
  wzIndex cl,cr,ul,ur;
  wzIndex t = cct(c);
  cl = wzGridCellLeft(*this,c,t);
  cr = wzGridCellRight(*this,c,t);
  if(ccdef(cl)) ul = ccu(cl); else ul=0;
  if(ccdef(cr)) ur = ccu(cr); else ur=0;
  if(ul>ur) ccu(c) = ul;
  else      ccu(c) = ur;
}

void 	wzHighlyAnisotropicalDelaunayGenerator::defineFaceOutside(int c, cogeometry g)
{
  wzIndex cl,cr,ul,ur;
  wzIndex t = cct(c);
  cl = wzGridCellLeft(*this,c,t);
  cr = wzGridCellRight(*this,c,t);
  if(ccdef(cl)) ul = ccu(cl); else ul=0;
  if(ccdef(cr)) ur = ccu(cr); else ur=0;
  if(ul>ur) ccu(c) = ul;
  else      ccu(c) = ur;
}

static int suppress_unreachable_warning = 1;
void 	wzHighlyAnisotropicalDelaunayGenerator::defineRegionOfCell(int c, cogeometry g)
{
  int n,nn,m=0,m1;
  wzPoint xx;
  wzArray<wzPoint>&  P = Grid.Point;
  //  wzFloat l0;
  wzPoint pp,*p0,*p1;
  ibg2CellType t = cct(c);
  int npoints = wzCellTypePoints(t);
  int *cn=ccn(c);
  for(n=0;n<npoints;n++){
    if(cnt(nn=cn[n])==wzIsRegion){
      if((m1=cnu(nn)) <= 0){ccu(c) = 0; return;}
      if(m && (m1 != m)) goto differentRegions;
    }
  }
  if(m){ccu(c) = m; return;}
  if(suppress_unreachable_warning) goto useInnerPoint;
  //  l0 = wzInfty;
  for(n=0;n<npoints;n++){
    p0 = &P[cn[n]];
    for(m=n+1;m<npoints;m++){
      p1 = &P[cn[m]];
      if((*p0)[0]!=(*p1)[0]){
	if((*p0)[1]!=(*p1)[1]) continue;
	if((*p0)[2]!=(*p1)[2]) continue;
      }
      if((*p0)[1]!=(*p1)[1]){
	if((*p0)[2]!=(*p1)[2]) continue;
      }
      goto useLineNode;
    }
  }
  goto useInnerPoint;
useLineNode:
  cogLine(*p0,*p1).barycentre(pp);  
  goto ende;
differentRegions:
  cogInfoBadRegionCell();
  if(ibgErrorneousRegion){
    ccu(c) = ibgErrorneousRegion; return;
  }
useInnerPoint:
  {
    wzPoint& compiler_problem = P[cn[0]];
    cogTetrahedron(compiler_problem,P[cn[1]],P[cn[2]],P[cn[3]]).barycentre(pp);
  }
ende:
  g->Point(pp);
  ccu(c) = pp.segment().region();
  return;
}

void	CogenGFZ::getMetric(wzMetricData& data,const cogFlag1& f) const
{
  //  int i;
  wzFloat gxx = G_ii[0]*16;
  wzFloat gyy = G_ii[1];
  wzFloat gzz = G_ii[2];
  ((wzDiagonalMetricData*)(&data))->g_ii[0] = gxx;
  ((wzDiagonalMetricData*)(&data))->g_ii[1] = gyy;
  ((wzDiagonalMetricData*)(&data))->g_ii[2] = gzz;
}

void	CogenGFZ::getMetric(wzMetricData& data, const wzPoint& p0) const
{
  cogFloat px = (p0.x()-x0)/dx;
  cogFloat ry = (p0.y()-y0)/dy; 
  cogFloat rz = (p0.z()-z0)/dz; 
  cogInteger nearedge;
  cogInteger iy = (cogInteger) ry, jy = iy+1; ry -= iy;
  cogInteger iz = (cogInteger) rz, jz = iz+1; rz -= iz;
  cogFloat sz = 1-rz, sy = 1-ry;
  if(iy<0){
    iy = jy = 0;
  }else if(jy>=ny){
    iy = jy = ny-1;
  }
  if(iz<0){
    iz = jz = 0;
  }else if(jz>=nz){
    iz = jz = nz-1;
  }
  wzIndex nii = nx*iy+np*iz;
  wzIndex nij = nx*iy+np*jz;
  wzIndex nji = nx*jy+np*iz;
  wzIndex njj = nx*jy+np*jz;
  wzInteger m1  = ((wzIndex) p0.segment().region()) - 1, m0 = m1-1;
  wzFloat  yf[2],zf[2],yg[2],zg[2];
  if(m0>=bm){
    yf[0] = sy*xx[m0+nii] + ry*xx[m0+nji];
    yf[1] = sy*xx[m0+nij] + ry*xx[m0+njj];
    zf[0] = sz*xx[m0+nii] + rz*xx[m0+nij];
    zf[1] = sz*xx[m0+nji] + rz*xx[m0+njj];
  }else{
    yf[0] = yf[1] = zf[0] = zf[1] = -wzInfty;
  }
  //  wzIndex nearedge = 0;
  if(m1<=em){
    yg[0] = sy*xx[m1+nii] + ry*xx[m1+nji];
    yg[1] = sy*xx[m1+nij] + ry*xx[m1+njj];
    zg[0] = sz*xx[m1+nii] + rz*xx[m1+nij];
    zg[1] = sz*xx[m1+nji] + rz*xx[m1+njj];
    if(m0>=bm){
      if(xx[m0+nii]==xx[m1+nii]) nearedge = 1;
      if(xx[m0+nij]==xx[m1+nij]) nearedge = 1;
      if(xx[m0+nji]==xx[m1+nji]) nearedge = 1;
      if(xx[m0+njj]==xx[m1+njj]) nearedge = 1;
    }
  }else{
    yg[0] = yg[1] = zg[0] = zg[1] = wzInfty;
  }
  wzFloat deltam = sz*(yf[0]-yg[0]) + rz*(yf[1]-yg[1]);	if(deltam<0) deltam = -deltam;
  wzFloat deltay = (yg[0]+yf[0]) - (yg[1]+yf[1]);  	if(deltay<0) deltay = -deltay;
  wzFloat deltaz = (zg[0]+zf[0]) - (zg[1]+zf[1]);  	if(deltaz<0) deltaz = -deltaz;
  if(deltam<dmmin){
    wzAssert(nearedge);
    deltam = dmmin;
  }
  wzFloat gxx,gyy,gzz;
  //  wzFloat drx = 1/deltam; drx *= drx;
  wzFloat dry = deltay/(deltam*dy); dry *= dry;
  wzFloat drz = deltaz/(deltam*dz); drz *= drz;
  //  if(G_ii[0]<drx)  gxx = drx; else gxx = G_ii[0];
  if(G_ii[1]<dry)  gyy = dry; else gyy = G_ii[1];
  if(G_ii[2]<drz)  gzz = drz; else gzz = G_ii[2];
  //  if(G_ii[0]>16*drx)  gxx = 16*G_ii[0];
  gxx = G_ii[0];
  wzFloat ifac=8024.0;
  if(ifac*G_ii[1]<dry)  gyy = ifac*G_ii[1];
  if(ifac*G_ii[2]<drz)  gzz = ifac*G_ii[2];
  //  wzAssert(256*G_ii[2]>drz);
  ((wzDiagonalMetricData*)(&data))->g_ii[0] = gxx;
  ((wzDiagonalMetricData*)(&data))->g_ii[1] = gyy;
  ((wzDiagonalMetricData*)(&data))->g_ii[2] = gzz;
}