isob.c

一个开源的火灾动力模拟的系统

      index = triangles[nn++];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);
      glVertex3f(x,y,z);
      glVertex3f(x0,y0,z0);
    }
    glEnd();
  }
  if(surfacedata->plottype==2){
    glColor3fv(color);
    glBegin(GL_POINTS);
    for (n = 0; n < nvert; n++) {
      x = xvert[n]; y = yvert[n]; z = zvert[n];
      glVertex3f(x,y,z);
    }
    glEnd();
  }

}
#endif

#ifndef pp_DRAWISO

/* ------------------ isoheader ------------------------ */

void CCisoheader(char *isofile, 
                 char *isolonglabel, char *isoshortlabel, char *isounits,
                 float *levels, int *nlevels, int *error){
  int version=1,option=1;
  FILE *isostream=NULL;
  int len[3];


  isostream=fopen(isofile,"wb");

  len[0]=strlen(isolonglabel)+1;
  len[1]=strlen(isoshortlabel)+1;
  len[2]=strlen(isounits)+1;

  fwrite(&version,1,4,isostream);
  fwrite(len,4,3,isostream);
  fwrite(isolonglabel, 1,len[0],isostream);
  fwrite(isoshortlabel,1,len[1],isostream);
  fwrite(isounits,     1,len[2],isostream);
  fwrite(&option,4,1,isostream);
  fwrite(nlevels,4,1,isostream);
  fwrite(levels,4,*nlevels,isostream);
  fclose(isostream);

}


/* ------------------ isoheader ------------------------ */

void CCtisoheader(char *isofile, 
                 char *isolonglabel, char *isoshortlabel, char *isounits,
                 float *levels, int *nlevels, int *error){
  int version=2,option=1;
  FILE *isostream=NULL;
  int len[3];
  int one=1;


  isostream=fopen(isofile,"wb");

  len[0]=strlen(isolonglabel)+1;
  len[1]=strlen(isoshortlabel)+1;
  len[2]=strlen(isounits)+1;

  fwrite(&one,4,1,isostream);  
  fwrite(&version,4,1,isostream);
  fwrite(len,4,3,isostream);
  fwrite(isolonglabel, 1,len[0],isostream);
  fwrite(isoshortlabel,1,len[1],isostream);
  fwrite(isounits,     1,len[2],isostream);
  fwrite(&option,4,1,isostream);
  fwrite(nlevels,4,1,isostream);
  fwrite(levels,4,*nlevels,isostream);
  fclose(isostream);

}

/* ------------------ isoout ------------------------ */

void isoout(FILE *isostream,float t, int timeindex, isosurface *surface, int *error){
/*			  unsigned short *vertices, unsigned short *tvertices, int nvertices, 
//			  int *trilist, int ntrilist, 
//			  int *error){  */
	unsigned char czero=0,*trilist1=NULL;
	unsigned short szero=0,*trilist2=NULL;
	int i;
  unsigned short *vertices, *tvertices;
  int nvertices, *trilist, ntrilist;

  vertices=surface->vertices;
  tvertices=surface->tvertices;
  nvertices=surface->nvertices;
  trilist=surface->triangles;
  ntrilist=surface->ntriangles;
/*      surface.vertices,NULL,surface.nvertices,
//      surface.triangles,surface.ntriangles,error); */


	if(timeindex==0)fwrite(&t,4,1,isostream);
	fwrite(&nvertices,4,1,isostream);
	fwrite(&ntrilist,4,1,isostream);
  if(nvertices>0){
    fwrite(vertices,2,3*nvertices,isostream);
    if(tvertices!=NULL){
      fwrite(&surface->tmin,4,1,isostream);
      fwrite(&surface->tmax,4,1,isostream);
      fwrite(tvertices,2,nvertices,isostream);
    }
  }
  if(ntrilist==0)return;
	if(nvertices<256){
    if(NewMemory((void **)&trilist1,sizeof(unsigned char)*ntrilist)==0){
      for(i=0;i<ntrilist;i++){
        fwrite(&czero,1,1,isostream);
      }
    }
    else{
	    for(i=0;i<ntrilist;i++){trilist1[i] = (unsigned char)trilist[i];}
	    fwrite(trilist1,1,ntrilist,isostream);
      FREEMEMORY(trilist1);
    }
	}
	else if(nvertices>=256&&nvertices<65536){
    if(NewMemory((void **)&trilist2,sizeof(unsigned short)*ntrilist)==0){
      for(i=0;i<ntrilist;i++){
        fwrite(&szero,2,1,isostream);
      }
    }
    else{
	    for(i=0;i<ntrilist;i++){trilist2[i] = (unsigned short)trilist[i];}
	    fwrite(trilist2,2,ntrilist,isostream);
      FREEMEMORY(trilist2);
    }
	}
	else{
	  fwrite(trilist,4,ntrilist,isostream);
	}
}

/* ------------------ isosurface2file ------------------------ */
void CCisosurface2file(char *isofile, float *t, float *data, int *iblank, 
						float *level, int *nlevels,
                   float *xplt, int *nx, 
                   float *yplt, int *ny, 
                   float *zplt, int *nz,
                   int *isooffset, int *reduce_triangles, int *error
                   ){
  isosurface surface;
  int i;
  FILE *isostream=NULL;

#ifdef _DEBUG
  printf("before surface creation:");
  PrintMemoryInfo();
#endif
  isostream=fopen(isofile,"ab");
  *error = 0;
  for(i=0;i<*nlevels;i++){
    InitIsosurface(&surface,level[i],NULL,i);
    surface.dataflag=0;
    if(GetIsosurface(&surface,data,NULL,iblank,level[i],xplt,*nx,yplt,*ny,zplt,*nz,*isooffset)!=0){
      *error=1;
      return;
    }
    if(GetNormalSurface(&surface)!=0){
      *error=1;
      return;
    }
    if(CompressIsosurface(&surface,*reduce_triangles,
      xplt[0],xplt[*nx-1],
      yplt[0],yplt[*ny-1],
      zplt[0],zplt[*nz-1]
      )!=0){
      *error=1;
      return;
    }

    isoout(isostream,*t,i,&surface,error);
/*      surface.vertices,NULL,surface.nvertices,
//      surface.triangles,surface.ntriangles,error); */

    freesurface(&surface);
  }
  fclose(isostream);
#ifdef _DEBUG
  printf("after surface creation:");
  PrintMemoryInfo();
#endif
}

void CCisosurfacet2file(char *isofile, float *t, float *data, int *data2flag, float *data2, int *iblank, 
						float *level, int *nlevels,
                   float *xplt, int *nx, 
                   float *yplt, int *ny, 
                   float *zplt, int *nz,
                   int *isooffset, int *reduce_triangles, int *error
                   ){
  isosurface surface;
  int i;
  FILE *isostream=NULL;
  int dataflag=0;
  float *tdata=NULL;

  if(*data2flag==1){
    dataflag=1;
    tdata=data2;
  }


#ifdef _DEBUG
  printf("before surface creation:");
  PrintMemoryInfo();
#endif
  isostream=fopen(isofile,"ab");
  *error = 0;
  for(i=0;i<*nlevels;i++){
    InitIsosurface(&surface,level[i],NULL,i);
    surface.dataflag=dataflag;
    if(GetIsosurface(&surface,data,tdata,iblank,level[i],xplt,*nx,yplt,*ny,zplt,*nz,*isooffset)!=0){
      *error=1;
      return;
    }
    if(GetNormalSurface(&surface)!=0){
      *error=1;
      return;
    }
    if(CompressIsosurface(&surface,*reduce_triangles,
      xplt[0],xplt[*nx-1],
      yplt[0],yplt[*ny-1],
      zplt[0],zplt[*nz-1]
      )!=0){
      *error=1;
      return;
    }

    isoout(isostream,*t,i,&surface,error);

    freesurface(&surface);
  }
  fclose(isostream);
#ifdef _DEBUG
  printf("after surface creation:");
  PrintMemoryInfo();
#endif
}
#endif

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#ifndef pp_noappend
#define CCsmoke3dheader smoke3dheader_
#define CCsmoke3dtofile smoke3dtofile_
#else
#define CCsmoke3dheader smoke3dheader
#define CCsmoke3dtofile smoke3dtofile
#endif

unsigned int rle(unsigned char *buffer_in, int nchars_in, unsigned char *buffer_out);

void CCsmoke3dheader(char *file,int *is1, int *is2, int *js1, int *js2, int *ks1, int *ks2){
  FILE *binfile,*textfile;
  int nxyz[8], VERSION=0;
  char textfilename[1024];

  nxyz[0]=1;
  nxyz[1]=VERSION;
  nxyz[2]=*is1;
  nxyz[3]=*is2;
  nxyz[4]=*js1;
  nxyz[5]=*js2;
  nxyz[6]=*ks1;
  nxyz[7]=*ks2;



  binfile=fopen(file,"wb");
  if(binfile==NULL)return;
  fwrite(nxyz,  4,1,binfile);
  fwrite(nxyz+1,4,7,binfile);
  fclose(binfile);

  strcpy(textfilename,file);
  strcat(textfilename,".sz");
  textfile=fopen(textfilename,"w");
  fprintf(textfile,"%i\n",VERSION);
  fclose(textfile);
}

void CCsmoke3dtofile(char *file, float *time, float *dx, float *extcoef,int *type, float *xyz, int *nx, int *ny, int *nz){

  FILE *binfile,*textfile;
  unsigned char *buffer_in, *buffer_out;
  int nchars_in, nchars_out;
  int i;
  double xtype;
  int nchars[2];
  char textfilename[1024];
  int nxyz;
  double factor;

#define SOOT 1
#define FIRE 2
#define WATER 3

  nxyz=(*nx)*(*ny)*(*nz);
  if(nxyz<1)return;

  buffer_in=(unsigned char *)malloc(nxyz);
  if(buffer_in==NULL)return;

  buffer_out=(unsigned char *)malloc(nxyz);
  if(buffer_out==NULL){
    free(buffer_in);
    return;
  }

  binfile=fopen(file,"ab");
  if(binfile==NULL){
    free(buffer_in);
    free(buffer_out);
    return;
  }

  nchars_in=nxyz;
  switch (*type){
  case SOOT:
    xtype = *extcoef;
    factor = -xtype*(*dx)/1000.0;
    for(i=0;i<nxyz;i++){
      if(*xyz<0.0)*xyz=0.0;
      buffer_in[i]=254*(1.0-exp( factor*(*xyz++)) );
    }

    break;
  case FIRE:
    for(i=0;i<nxyz;i++){
      if(*xyz<0.0)*xyz=0.0;
      if(*xyz>1200.0)*xyz=1200.0;
      buffer_in[i]=254*(*xyz/1200.0);
      xyz++;
    }

    break;
  case WATER:
    factor=1.0/(0.1*0.5*(*dx));
    for(i=0;i<nxyz;i++){
      *xyz=*xyz-0.003;
      if(*xyz<0.0)*xyz=0.0;
      buffer_in[i]=254*(1.0-pow(0.5,*xyz*factor));
      xyz++;
    }

    break;
  }

  nchars_out=rle(buffer_in,nchars_in,buffer_out);

  nchars[0]=nchars_in;
  nchars[1]=nchars_out;

  fwrite(time,4,1,binfile);
  fwrite(nchars,4,2,binfile);
  if(nchars_out>0)fwrite(buffer_out,1,nchars_out,binfile);

  free(buffer_in);
  free(buffer_out);
  fclose(binfile);

  strcpy(textfilename,file);
  strcat(textfilename,".sz");
  textfile=fopen(textfilename,"a");
  fprintf(textfile,"%f %i %i\n",*time,nchars_in,nchars_out);
  fclose(textfile);
}

#define MARK 255

unsigned int rle(unsigned char *buffer_in, int nchars_in, unsigned char *buffer_out){
  unsigned char lastchar=MARK, cmark=MARK, thischar, *buffer_start;
  unsigned char *buffer_in_end;
  int nrepeats=1;

  buffer_start=buffer_out;
  buffer_in_end = buffer_in + nchars_in;

  while(buffer_in<buffer_in_end){
    thischar=*buffer_in;
    if(thischar==lastchar){
      nrepeats++;
    }
    else{
      nrepeats=1;
    }
    switch (nrepeats){
    case 1:
    case 2:
    case 3:
      *buffer_out=thischar;
      lastchar=thischar;
      break;
    case 4:
      buffer_out-=3;
      *buffer_out++=cmark;
      *buffer_out++=thischar;
    default:
      if(nrepeats!=4)buffer_out--;
      *buffer_out=nrepeats;
      if(nrepeats==254){
        nrepeats=1;
        lastchar=MARK;
      }
      break;
    }
    buffer_in++;
    buffer_out++;

  }
  return buffer_out-buffer_start;
}

unsigned int irle(unsigned char *buffer_in, int nchars_in, unsigned char *buffer_out){
  int nrepeats,nn;
  unsigned char thischar, *buffer_in_end;

  nn=0;
  buffer_in_end  = buffer_in  + nchars_in;

  while(buffer_in<buffer_in_end){
    if(*buffer_in==MARK){
      if(buffer_in+2>=buffer_in_end)break;
      buffer_in++;
      thischar=*buffer_in++;
      nrepeats=*buffer_in++;
      nn+=nrepeats;
      memset(buffer_out,thischar,nrepeats);
      buffer_out+=nrepeats;
    }
    else{
      *buffer_out++=*buffer_in++;
      nn++;
    }


  }
  return nn;
}