ibgoutput.cxx

Delaunay三角形的网格剖分程序

#include "wzgrid.hxx"
#include "cogwarnings.hxx"
#include "wzview.hxx"
// #include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>

static int writeMaterial = 0;

void wzGrid::write(char* filename)
{
  int* cnum;
  int  ic,ncells,nfaces,u,t;
  cnum = (int*)calloc((cells.last()+1),sizeof(int));
  ncells = 0;
  nfaces = 0;
  wzRangeLoop(cells,ic){
    t   = cct(ic);
    u   = ccu(ic);   
    if(u<=0) continue;
    switch(wzCellTypeCODIM(t)){
    case wzIsRegion:
      if(writeMaterial){
	u = wzRegion(u).material();
	if(u<=0) continue;
      }
      cnum[ic] = ++ncells; break;
    case wzIsFace:
      if(writeMaterial){
	u = wzFace(u).type();
	if(u<=0) continue;
      }
      cnum[ic] = ++nfaces; break;
    default:
      continue;
    }
  };
  write(filename, cnum, ncells,nfaces);
}

void wzGrid::write(char* filename, wzSegmentDescription& d)
{
  d.use();
  writeMaterial = 1;
  write(filename);
  writeMaterial = 0;
}  

void wzGrid::write(char* filename, int* cnum, int ncells, int nfaces)
{
  FILE *file;
  int dim,sdim,d,ic,cc,t1,u,in,nn,n,*cs,*cn,*nnum,i,
    nnodes,
    cl,cr;
  time_t date,t;
  date=time(&t);
  cogInfoWriteGrid("\n<>D-Grid with <> points, <> cells, <> face cells ");
  file = fopen(filename,"w");
  fprintf(file,"SimplexGrid");
  fprintf(file," ");
  fprintf(file,"2.0\n");
  fprintf(file,"#created by wzgrid (c) I.Schmelzer\n");
  fprintf(file,"#%s",ctime(&date));
  fprintf(file,"DIMENSION\n%d\n",dim=gdim());

  cogInfoWriteGrid << dim;
  nnum = (int*)calloc((points.last()+1),sizeof(int));
  wzRangeLoop(cells,ic){
    if(cnum[ic]==0) continue;
    t   = cct(ic);
    cn = ccn(ic);
    for(n=0;n<wzCellTypePoints(t);n++){
      nnum[cn[n]] = 1;
    }
  };
  nnodes = 0;
  wzRangeLoop(points,in){
    if(nnum[in]) nnum[in] = ++nnodes;
  }
  cogInfoWriteGrid << nnodes,ncells,nfaces;
  sdim=dim;
  if(sdim==0) sdim=xdim();
  fprintf(file,"NODES\n%d %d\n",nnodes,sdim);
  wzRangeLoop(points,in){
    if(nnum[in]==0) continue;
    for(d=0;d<sdim;d++) fprintf(file,"%f ",(float)(cnx(in)[d]));
    fprintf(file,"\n");
  }

  if(ncells){
    fprintf(file,"CELLS\n%d\n",ncells);
    wzRangeLoop(cells,ic){
      t   = cct(ic);
      u   = ccu(ic);
      if(wzCellTypeCODIM(t) != wzIsRegion) continue;
      if(cnum[ic]==0) continue;
      if(writeMaterial) u = wzRegion(u).material();
      cn = ccn(ic);
      for(i=0;i<wzCellTypePoints(t);i++){
	nn = cn[i];
	nn = nnum[nn];
	fprintf(file,"%d ",nn);
      }
      fprintf(file,"%d ",u);
      cs = ccs(ic);
      for(i=0;i<wzCellTypeSides(t);i++)  {
	cc = cs[i];
	if(cc>0){
	  t1 = cct(cc);
	  if(wzCellTypeCODIM(t1)==wzIsRegion)	  {
	    fprintf(file,"%d ",cnum[cc]);
	  }else if(wzCellTypeCODIM(t1)==wzIsFace)	  {
	    fprintf(file,"-%d ",cnum[cc]);
	  }
	  else	wzAssert(0);
	}else{
	  fprintf(file,"-0 ");
	}
      }
      fprintf(file,"\n");
    }
  }
  if(nfaces){
    fprintf(file,"FACES\n%d\n",nfaces);
    wzRangeLoop(cells,ic){
      t   = cct(ic);
      u   = ccu(ic);
      if(wzCellTypeCODIM(t) != wzIsFace) continue;
      if(cnum[ic]==0) continue;
      if(writeMaterial) u = wzFace(u).type();
      cn = ccn(ic);
      for(i=0;i<wzCellTypePoints(t);i++){
	nn = cn[i];
	nn = nnum[nn];
	fprintf(file,"%d ",nn);
      }
      fprintf(file,"%d ",u);
      if(sdim==3){
	cl = wzGridCellLeft(*this,ic,t);
	cr = wzGridCellRight(*this,ic,t); 
      }else{
	cr = wzGridCellLeft(*this,ic,t);
	cl = wzGridCellRight(*this,ic,t); 
      }
      if(ncells){
	if(cnum[cl]) fprintf(file,"%d ",cnum[cl]);
	else fprintf(file,"-%d ",ccu(cl));
	if(cnum[cr]) fprintf(file,"%d ",cnum[cr]);
	else fprintf(file,"-%d ",ccu(cr));
      }else{
	fprintf(file,"-%d ",ccu(cl));
	fprintf(file,"-%d ",ccu(cr));
      }
      cs = ccs(ic);
      for(i=0;i<wzCellTypeSides(t);i++){
	cc = cs[i];
	if(cc==0)
	  fprintf(file,"0 ");
	else if(wzCellTypeCODIM(t1 = cct(cc))==wzIsFace)  {
	  fprintf(file,"%d ",cnum[cc]);
	}else if(wzCellTypeCODIM(t1)==wzIsEdge){
	  fprintf(file,"0 ");
	}
	else	wzAssert(0);
      }
      fprintf(file,"\n");
    }
  }
  wzIndex nfuncs;
  if(nfuncs=floatvalues.last()){
    fprintf(file,"NODEFUNCTION\nnoname\n%d\n",nfuncs);
    wzRangeLoop(points,in){
      if(nnum[in]==0) continue;
      for(d=0;d<nfuncs;d++) fprintf(file,"%f ",(float)(cnf(in)[d]));
      fprintf(file,"\n");
    }
  }
  fprintf(file,"END\n");
  fclose(file);
  free(nnum);
  free(cnum);
  cogInfoWriteGrid("written into file <>;\n") << filename;
}

#define buflen 200

wzGrid::wzGrid(char *filename)
  :points(sizeof(wzPoint))
  ,Point(points.base)
  ,PointCell(points)
  ,cells(3*sizeof(wzInteger))
  ,CellSegment(cells.base)
  ,CellType(cells.base,1,sizeof(wzInteger))
  ,CellData(cells.base,2,sizeof(wzInteger))
  ,data(sizeof(wzInteger))
  ,Data(data.base)
  ,floatvalues(wzPointDimFloat)
{
  wzCurrentData::grid = this;
  wzCellTypeInit();
  FILE *file;
  file = fopen(filename,"r");
  sxConstructor();
  firstRegion=firstFace=firstEdge=firstNode=1;
  lastRegion=lastFace=lastEdge=lastNode=0;
  int *ccn,*ccs,cb,c,cl,cr,u,n,rc,i,j,k;
  int pnum,rpoints,fpoints,fsides,rcells,fcells,rdata,fdata,cdata;
  int num[9],gdim,xdim;
  wzCellType ftyp,rtyp;
  float x,y,z;
  wzFloat xmin = wzInfty, xmax = -wzInfty;
  wzFloat ymin = wzInfty, ymax = -wzInfty;
  wzFloat zmin = wzInfty, zmax = -wzInfty;
  wzPoint *point;
  char buffer[buflen];
  fgets(buffer,buflen,file);
  if(strncmp(buffer,"SimplexGrid 2.0",15)) goto readerror;
  fgets(buffer,buflen,file);
  while(buffer[0]=='#') fgets(buffer,buflen,file);
  if(strncmp(buffer,"DIMENSION",9)) goto readerror;
  fgets(buffer,buflen,file);
  rc = sscanf(buffer,"%d",&gdim);	// the grid dimension 
  if(rc != 1) goto readerror;
  fgets(buffer,buflen,file);
  if(strncmp(buffer,"NODES",5)) goto readerror;
  fgets(buffer,buflen,file);
  rc = sscanf(buffer,"%d %d",&pnum,&xdim);	// the number of points 
  if(rc <= 1) goto readerror;
  if(xdim>wzPointDim) xdim = wzPointDim;
  gridDim  	= gdim;
  spaceDim	= xdim;
  n = points.create(pnum);
  wzRangeLoop(points,n){
    fgets(buffer,buflen,file);
    switch(xdim){
    case 1:
      rc = sscanf(buffer,"%f",&x); y=0; z=0;
      if(rc != 1) goto readerror;
      break;
    case 2:
      rc = sscanf(buffer,"%f %f",&x,&y); z=0;
      if(rc != 2) goto readerror;
      break;
    case 3:
      rc = sscanf(buffer,"%f %f %f",&x,&y,&z);
      if(rc != 3) goto readerror;
      break;
    }
    point = &Point[n];
    if(x<xmin)	xmin = x;
    if(y<ymin)	ymin = y;
    if(z<zmin)	zmin = z;
    if(x>xmax)	xmax = x;
    if(y>ymax)	ymax = y;
    if(z>zmax)	zmax = z;
    point->x()=x;      point->y()=y;       point->z()=z; 
    point->segment() = wzRegion();
  }
  cgmin(0) = xmin;
  cgmin(1) = ymin;
  cgmin(2) = zmin;
  cgmax(0) = xmax;
  cgmax(1) = ymax;
  cgmax(2) = zmax;

  fgets(buffer,buflen,file);
  if(strncmp(buffer,"CELLS",5)) goto readerror;
  fgets(buffer,buflen,file);
  rc = sscanf(buffer,"%d",&rcells);	// the number of region cells 
  if(rc != 1) goto readerror;
  switch(gdim){
  case 1:
    rtyp = wzCellType1Line;	rpoints = 2;
    ftyp = wzCellType1Node;	fpoints = 1;  fsides = 0;
    break;
  case 2:
    rtyp = wzCellType2Triangle;	rpoints = 3;
    ftyp = wzCellType2Line;	fpoints = fsides = 2;
    break;
  case 3:
    rtyp = wzCellType3Tetrahedron;rpoints = 4;
    ftyp = wzCellType3Triangle;	fpoints = fsides = 3;
  }		     
  rdata = 2*rpoints;
  fdata = fpoints+fsides+2;
  c     = cells.create(rcells);
  cdata = data.create(rdata*rcells);
  ccn   = &Data[cdata]; 
  ccs   = ccn+rpoints;
  CellSegment[0] = 0;
  wzRangeLoop(cells,c){
    fgets(buffer,buflen,file);
    rc = sscanf(buffer,"%d %d %d %d %d %d %d %d %d ",
		&num[0],&num[1],&num[2],&num[3],
		&num[4],&num[5],&num[6],&num[7],&num[8]);
    CellType[c] = rtyp;
    CellData[c] = cdata;
    for(i=0;i<rpoints;i++){
      ccn[i] = num[i];	// at first the points 
    }
    CellSegment[c] = num[i++];	// now the segment 
    for(j=0;j<rpoints;i++,j++){		// now the neighbour cells: 
      if(num[i]>0)		ccs[j] = num[i];
      else if(num[i]<0)		ccs[j] = rcells - num[i];
      else 			ccs[j] = 0;
    }
    cdata += rdata;
    ccn   += rdata;
    ccs   += rdata;
  }

  fgets(buffer,buflen,file);
  if(strncmp(buffer,"FACES",5)) goto readerror;
  fgets(buffer,buflen,file);
  rc = sscanf(buffer,"%d",&fcells);	// the number of boundary face cells 
  if(rc != 1) goto readerror;
  c = cells.append(fcells);
  cdata = data.append(fdata*fcells);
  ccn   = &Data[cdata]; 
  ccs   = ccn+fpoints;
  for(cb=1;cb<=fcells;c++,cb++){
    for(i=0;i<9;i++) num[i]=0;
    fgets(buffer,buflen,file);
    rc = sscanf(buffer,"%d %d %d %d %d %d %d %d %d ",
		&num[0],&num[1],&num[2],&num[3],
		&num[4],&num[5],&num[6],&num[7],&num[8]);
    CellType[c] = ftyp;
    CellData[c] = cdata;
    for(i=0;i<fpoints;i++){
      ccn[i] = num[i];
    }
    u  = num[i++];		// segment  number 
    cr = num[i++];		// right neighbour (cell or segment number)
    cl = num[i++];		// left  neighbour (cell or segment number) 
    for(j=0;j<fsides;j++){
      if(k=num[i++]>0)
	ccs[j] = rcells+k;
      else
	ccs[j] = 0;
    }
    ccs[j++] = cl;
    ccs[j++] = cr;
    CellSegment[c] = u;
    cdata += fdata;
    ccn   += fdata;
    ccs   += fdata;
  }
  fgets(buffer,buflen,file);
  if(!strncmp(buffer,"END",3)) goto end;
  if(!strncmp(buffer,"NODEFUNCTION",12)){
    float u[6]; int nfuncs,in,d;
    fgets(buffer,buflen,file); // ignore the name;
    fgets(buffer,buflen,file);
    rc = sscanf(buffer,"%d",&nfuncs);
    if(rc != 1) goto readerror;
    if(nfuncs<1) goto readerror;
    if(nfuncs>6 || nfuncs>wzPointDimFloat){
      //      cogInfoMaximalFunctionValueNumber();
      nfuncs = wzPointDimFloat; 
      if(nfuncs>6) nfuncs = 6;
    }
    rc = floatvalues.create(nfuncs);
    wzAssert(rc==1);
    wzRangeLoop(points,in){
      fgets(buffer,buflen,file);
      rc = sscanf(buffer,"%f %f %f %f %f %f",&u[0],&u[1],&u[2],&u[3],&u[4],&u[5]);
      if(rc < nfuncs) goto readerror;
      for(d=0;d<nfuncs;d++){
	cnf(in)[d] = u[d];
      }
    }
  }
readerror:
end:
  fclose(file);
}