isob.c

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

    edge = edges[n];
    v1 = case2[edge2vertex[edge][0]];
    v2 = case2[edge2vertex[edge][1]];
    val1 = vals[v1]-level; val2 = vals[v2]-level;
    denom = val2 - val1;
    factor = 0.5;
    if(denom!=0.0)factor = -val1/denom;
    if(factor<0.5){closestnodes[n]=nodeindexes[v1];}
    else{closestnodes[n]=nodeindexes[v2];}
    if(factor>1.0){/*factor=1;*/outofbounds=1;}
    if(factor<0.0){/*factor=0.0;*/outofbounds=1;}
    xx = xxval[v1]*(1.0-factor) + xxval[v2]*factor;
    yy = yyval[v1]*(1.0-factor) + yyval[v2]*factor;
    zz = zzval[v1]*(1.0-factor) + zzval[v2]*factor;
    xvert[n] = xx;
    yvert[n] = yy;
    zvert[n] = zz;
    if(tvert!=NULL)tvert[n]= tvals[v1]*(1.0-factor) + tvals[v2]*factor;

  }
  if(outofbounds==1){
    printf("*** warning - computed isosurface vertices are out of bounds for :\n");
    printf("case number=%i level=%f\n",casenum,level);
    printf("values=");
    for(n=0;n<8;n++){printf("%f ",vals[n]);}
    printf("\n");
    printf("x=%f %f y=%f %f z=%f %f\n\n",x[0],x[1],y[0],y[1],z[1],z[2]);
  }

  /* copy coordinates to output array */

  *nvert = nedges;
  *ntriangles = npath;
  for(n=0;n<npath;n++){triangles[n] = path[n];}
}

/* ------------------ calcNormal2 ------------------------ */

void calcNormal2(unsigned short *v1, unsigned short *v2, unsigned short *v3, float *out, float *area){
  float u[3], v[3];
  int i;


  for(i=0;i<3;i++){
    u[i]=v2[i]-v1[i];
    v[i]=v3[i]-v1[i];
  }


  out[0] = u[1]*v[2] - u[2]*v[1];
  out[1] = u[2]*v[0] - u[0]*v[2];
  out[2] = u[0]*v[1] - u[1]*v[0];
  *area = sqrt(out[0]*out[0]+out[1]*out[1]+out[2]*out[2])/2.0;

  ReduceToUnit(out);

}

/* ------------------ calcNormal ------------------------ */

void calcNormal(float *xx, float *yy, float *zz, float *out){
  float u[3],v[3], p1[3], p2[3], p3[3];
  static const int x = 0;
  static const int y = 1;
  static const int z = 2;

  p1[x]=xx[0]; p1[y]=yy[0]; p1[z]=zz[0];
  p2[x]=xx[1]; p2[y]=yy[1]; p2[z]=zz[1];
  p3[x]=xx[2]; p3[y]=yy[2]; p3[z]=zz[2];


  u[x] = p2[x] - p1[x];
  u[y] = p2[y] - p1[y];
  u[z] = p2[z] - p1[z];

  v[x] = p3[x] - p1[x];
  v[y] = p3[y] - p1[y];
  v[z] = p3[z] - p1[z];

  out[x] = u[y]*v[z] - u[z]*v[y];
  out[y] = u[z]*v[x] - u[x]*v[z];
  out[z] = u[x]*v[y] - u[y]*v[x];

  ReduceToUnit(out);

}

/* ------------------ ReduceToUnit ------------------------ */

void ReduceToUnit(float *v){

  float length;

  length = (float)sqrt((v[0]*v[0]+v[1]*v[1]+v[2]*v[2]));

  if(length==0.0f)length=1.0f;

  v[0] /= length;
  v[1] /= length;
  v[2] /= length;
}

/* ------------------ GetIsoSurface ------------------------ */

int GetIsosurface(isosurface *surface, float *data, float *tdata, int *iblank_cell, float level,
                   float *xplt, int nx, float *yplt, int ny, float *zplt, int nz,
                   int isooffset
                   ){
#define ijk(i,j,k) ((i)+(j)*nx+(k)*nx*ny)
#define ijkcell(i,j,k) ((i)+(j)*ibar+(k)*ijbar)
#define ij(i,j) ((i)+(j)*nx)

  int ibar,ijbar;
  float xvert[12], yvert[12], zvert[12], tvert[12], *tvertptr=NULL;
  int triangles[18];
  int nvert, ntriangles;
  int i, j, k;
  float xxx[2], yyy[2], zzz[2], vals[8], tvals[8], *tvalsptr=NULL;
  int nodeindexes[8], closestnodes[18];
  float *xx, *yy, *zz;
  int ijkbase,ip1jk,ijkp1,ip1jkp1,ijp1k,ip1jp1k,ijp1kp1,ip1jp1kp1;
  int ijbase;
  int nxy;
  
  ibar = nx-1;
  ijbar = (nx-1)*(ny-1);
  xx = xxx;
  yy = yyy;
  zz = zzz;
  nxy = nx*ny;
  if(surface->defined==0)return 0;
  if(tdata!=NULL){
    tvalsptr=tvals;
    tvertptr=tvert;
  }
  for(i=0;i<nx-isooffset;){
    xx[0]=xplt[i];
    xx[1]=xplt[i+isooffset];
    for(j=0;j<ny-isooffset;){
      yy[0]=yplt[j];
      yy[1]=yplt[j+isooffset];
      ijbase = ij(i,j);
      for(k=0;k<nz-isooffset;){
        ijkbase = ijbase + k*nxy;
        ip1jk = ijkbase + isooffset;
        ijkp1 = ijkbase + isooffset*nxy;
        ip1jkp1 = ijkbase + isooffset*(1+nxy);
        ijp1k = ijkbase + isooffset*nx;
        ip1jp1k = ijkbase + isooffset*(1+nx);
        ijp1kp1 = ijkbase + isooffset*(nx+nxy);
        ip1jp1kp1 = ijkbase + isooffset*(1+nx+nxy);

        if(isooffset!=1||iblank_cell==NULL||iblank_cell[ijkcell(i,j,k)]!=0){
          zz[0]=zplt[k];
          zz[1]=zplt[k+isooffset];

          nodeindexes[0]=ijkbase;
          nodeindexes[1]=ijp1k;
          nodeindexes[2]=ip1jp1k;
          nodeindexes[3]=ip1jk;
          nodeindexes[4]=ijkp1;
          nodeindexes[5]=ijp1kp1;
          nodeindexes[6]=ip1jp1kp1;
          nodeindexes[7]=ip1jkp1;

          vals[0]=data[ijkbase];
          vals[1]=data[ijp1k];
          vals[2]=data[ip1jp1k];
          vals[3]=data[ip1jk];
          vals[4]=data[ijkp1];
          vals[5]=data[ijp1kp1];
          vals[6]=data[ip1jp1kp1];
          vals[7]=data[ip1jkp1];

          if(tdata!=NULL){
            tvals[0]=tdata[ijkbase];
            tvals[1]=tdata[ijp1k];
            tvals[2]=tdata[ip1jp1k];
            tvals[3]=tdata[ip1jk];
            tvals[4]=tdata[ijkp1];
            tvals[5]=tdata[ijp1kp1];
            tvals[6]=tdata[ip1jp1kp1];
            tvals[7]=tdata[ip1jkp1];
          }

          GetIsobox(xx, yy, zz, vals, tvalsptr, nodeindexes, level,
                    xvert, yvert, zvert, tvertptr, closestnodes, &nvert, triangles, &ntriangles);

          if(UpdateIsosurface(surface, xvert, yvert, zvert, tvertptr, 
                              closestnodes, nvert, triangles, ntriangles)!=0)return 1;

        }
        k+=isooffset;
      }
      j+=isooffset;
    }
    i+=isooffset;
  }
  surface->defined=1;
  return 0;
}

/* ------------------ compareisonodes ------------------------ */

int compareisonodes( const void *arg1, const void *arg2 ){
  int i, j;
  i=*(int *)arg1;
  j=*(int *)arg2;
  i *= 3;
  j *= 3;
  if(vertices[i]<vertices[j])return -1;
  if(vertices[i]>vertices[j])return 1;
  i++; j++;
  if(vertices[i]<vertices[j])return -1;
  if(vertices[i]>vertices[j])return 1;
  i++; j++;
  if(vertices[i]<vertices[j])return -1;
  if(vertices[i]>vertices[j])return 1;
  return 0;
}

/* ------------------ computerank ------------------------ */

int computerank( const void *arg1, const void *arg2 ){
  int i, j;
  i=*(int *)arg1;
  j=*(int *)arg2;
  if(sortedlist[i]<sortedlist[j])return -1;
  if(sortedlist[i]>sortedlist[j])return 1;
  return 0;
}

int order_closestnodes( const void *arg1, const void *arg2 ){
  int i, j;
  int ii,jj;
  i=*(int *)arg1;
  j=*(int *)arg2;
  ii = closestnodes[i];
  jj = closestnodes[j];
  if(ii<jj)return -1;
  if(ii>jj)return 1;
  return 0;
}

/* ------------------ CompressIsosurface ------------------------ */

int CompressIsosurface(isosurface *surface, int reduce_triangles,
                        float xmin, float xmax, 
                        float ymin, float ymax, 
                        float zmin, float zmax){
  int i,j,nn;
  float *x=NULL, *y=NULL, *z=NULL, *t=NULL;
  int *map=NULL,*map2=NULL,nmap2,*vertexmap=NULL,*inverse_vertexmap=NULL;
  int nvertices;
  unsigned short *newvertices=NULL;
  int *triangles=NULL,*newtriangles=NULL,*newtriangles2=NULL,ntriangles;
  float xyzmaxdiff;
  int *cs=NULL, *ordered_closestnodes=NULL;
  int v1, v2, v3;
  int ii,iim1;
  unsigned short vx, vy, vz;
  int sumx, sumy, sumz, sumt;
  int nnewvertices;
  float tmin, tmax, tmaxmin;
  int flag;
  unsigned short *tvertices,*newtvertices;

  nvertices=surface->nvertices;
  if(nvertices==0)return 0;
  vertices=surface->vertices;

  x=surface->xvert; 
  y=surface->yvert; 
  z=surface->zvert; 
  if(surface->tvert!=NULL)t=surface->tvert;

  sortedlist=surface->sortedlist;
  rank=surface->rank;
  triangles=surface->triangles;
  ntriangles=surface->ntriangles;
  
  if(surface->dataflag==1){
    tvertices=surface->tvertices;
    FREEMEMORY(tvertices);
    if(NewMemory((void **)&tvertices,nvertices*sizeof(unsigned short))==0){
      FREEMEMORY(tvertices);
      return 1;
    }
    surface->tvertices=tvertices;
  }

  FREEMEMORY(vertices); FREEMEMORY(rank);
  if(NewMemory((void **)&vertices,3*nvertices*sizeof(unsigned short))==0||
     NewMemory((void **)&rank,nvertices*sizeof(int))==0||
     NewMemory((void **)&map,nvertices*sizeof(int))==0||
     NewMemory((void **)&map2,nvertices*sizeof(int))==0||
     NewMemory((void **)&sortedlist,nvertices*sizeof(int))==0){
    FREEMEMORY(vertices);
    FREEMEMORY(rank);
    FREEMEMORY(map);
    FREEMEMORY(map2);
    FREEMEMORY(sortedlist);
    return 1;
  }

  surface->vertices=vertices;
  surface->rank=rank;
  xyzmaxdiff = xmax - xmin;
  if(ymax-ymin>xyzmaxdiff)xyzmaxdiff=ymax-ymin;
  if(zmax-zmin>xyzmaxdiff)xyzmaxdiff=zmax-zmin;
  surface->xmin=xmin;
  surface->ymin=ymin;
  surface->zmin=zmin;
  surface->xyzmaxdiff=xyzmaxdiff;

  if(surface->dataflag==1){
    tmin=t[0];
    tmax=tmin;
    for(i=1;i<nvertices;i++){
      if(t[i]<tmin)tmin=t[i];
      if(t[i]>tmax)tmax=t[i];
    }
    surface->tmin=tmin;
    surface->tmax=tmax;
    tmaxmin=tmax-tmin;
    if(tmaxmin>0.0){
      for(i=0;i<nvertices;i++){
        tvertices[i]=65535*(t[i]-tmin)/tmaxmin;
      }
    }
    else{
      for(i=0;i<nvertices;i++){
        tvertices[i]=0;
      }
    }
  }

  for(i=0;i<nvertices;i++){
    vertices[3*i+0]=65535*(x[i]-xmin)/xyzmaxdiff;
    vertices[3*i+1]=65535*(y[i]-ymin)/xyzmaxdiff;
    vertices[3*i+2]=65535*(z[i]-zmin)/xyzmaxdiff;
    sortedlist[i]=i;
    map[i]=i;
  	rank[i]=i;
  }

  qsort((int *)sortedlist,(size_t)nvertices,sizeof(int),compareisonodes);
  qsort((int *)rank,(size_t)nvertices,sizeof(int),computerank);

  j=0;
  map[0]=0;
  map2[0]=0;
  nmap2=1;
  for(i=1;i<nvertices;i++){
    if(compareisonodes(sortedlist+i-1,sortedlist+i)!=0){
  	  j++;
  	  map2[j]=i;
  	  nmap2++;
    }
    map[i]=j;
  }

  nvertices=nmap2;

  if(NewMemory((void **)&newvertices,3*nvertices*sizeof(unsigned short))==0||
     NewMemory((void **)&cs,nvertices*sizeof(int))==0){
    FREEMEMORY(newvertices);
    FREEMEMORY(cs);
    return 1;
  }
  if(surface->dataflag==1){
    if(NewMemory((void **)&newtvertices,nvertices*sizeof(unsigned short))==0){
      FREEMEMORY(newtvertices);
      return 1;
    }
  }

  closestnodes=surface->closestnodes;

  for(i=0;i<nvertices;i++){
	  j=sortedlist[map2[i]];
	  newvertices[3*i+0]=vertices[3*j+0];
	  newvertices[3*i+1]=vertices[3*j+1];
	  newvertices[3*i+2]=vertices[3*j+2];
    cs[i] = closestnodes[j];
  }
  if(surface->dataflag==1){
    for(i=0;i<nvertices;i++){
	    j=sortedlist[map2[i]];
	    newtvertices[i]=tvertices[j];
    }
  }

  FREEMEMORY(closestnodes);
  surface->closestnodes = cs;

  if(NewMemory((void **)&newtriangles,ntriangles*sizeof(int))==0){
    return 1;
  }
  for(i=0;i<ntriangles;i++){newtriangles[i]=map[rank[triangles[i]]];}
  surface->triangles=newtriangles;
  surface->vertices=newvertices;
  surface->nvertices=nvertices;
  FREEMEMORY(vertices);
  if(surface->dataflag==1){
    surface->tvertices=newtvertices;
    FREEMEMORY(tvertices);
  }
  FREEMEMORY(triangles);
  FREEMEMORY(map);FREEMEMORY(map2);FREEMEMORY(sortedlist);

  if(reduce_triangles!=1)return 0;

  /* phase II compression, reduce the number of triangles */

  /* first, eliminate triangles whose nodes (2 or more) are closest to the same grid point */

  if(NewMemory((void **)&newtriangles2,ntriangles*sizeof(int))==0){
    return 1;
  }
  nn=0;
  for(i=0;i<ntriangles/3;i++){
    v1=newtriangles[3*i];
    v2=newtriangles[3*i+1];
    v3=newtriangles[3*i+2];
    if(cs[v1]!=cs[v2]&&cs[v1]!=cs[v3]&&cs[v2]!=cs[v3]){
      newtriangles2[nn++]=v1;
      newtriangles2[nn++]=v2;
      newtriangles2[nn++]=v3;
    }
  }
  FREEMEMORY(newtriangles);
  surface->triangles=newtriangles2;
  ntriangles=nn;
  surface->ntriangles=ntriangles;

  /* sort the closestnodes list */

  closestnodes=surface->closestnodes;
  if(NewMemory((void **)&ordered_closestnodes,nvertices*sizeof(int))==0){
    return 1;
  }
  for(i=0;i<nvertices;i++){ordered_closestnodes[i]=i;}
  qsort((int *)ordered_closestnodes,(size_t)nvertices,sizeof(int),order_closestnodes);
  if(NewMemory((void **)&vertexmap,nvertices*sizeof(int))==0||
     NewMemory((void **)&inverse_vertexmap,nvertices*sizeof(int))==0){
    FREEMEMORY(vertexmap);
    FREEMEMORY(inverse_vertexmap);
    return 1;
  }

  for(i=0;i<nvertices;i++){vertexmap[i]=i;}
  nn=0;
  sumx=0; sumy = 0; sumz = 0; sumt=0;

  /* average nodes */

  nn=0;
  vertices = surface->vertices;
  if(surface->dataflag==1)tvertices = surface->tvertices;