isob.c

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

  nnewvertices=0;
  for(i=1;i<nvertices+1;i++){
    iim1=ordered_closestnodes[i-1];
    if(i!=nvertices)ii=ordered_closestnodes[i];
    inverse_vertexmap[iim1] = nnewvertices;
    nn++;
    vx = vertices[3*iim1  ]; sumx += vx;
    vy = vertices[3*iim1+1]; sumy += vy;
    vz = vertices[3*iim1+2]; sumz += vz;
    if(surface->dataflag==1)sumt += tvertices[iim1];

    flag=0;
    if(i!=nvertices&&closestnodes[ii]!=closestnodes[iim1])flag=1;
    if(i==nvertices||flag==1){
      vertexmap[nnewvertices++] = iim1;
      vertices[3*iim1  ]=sumx/nn;
      vertices[3*iim1+1]=sumy/nn;
      vertices[3*iim1+2]=sumz/nn;
      if(surface->dataflag==1)tvertices[iim1]=sumt/nn;

      nn=0; 
      sumx = 0; sumy = 0;sumz = 0; sumt = 0;
    }
  }

  if(NewMemory((void **)&newvertices,3*nnewvertices*sizeof(unsigned short))==0){
    return 1;
  }
  if(surface->dataflag==1){
    if(NewMemory((void **)&newtvertices,nnewvertices*sizeof(unsigned short))==0){
      return 1;
    }
  }
  for(i=0;i<nnewvertices;i++){
    ii = vertexmap[i];
    newvertices[3*i  ] = vertices[3*ii  ];
    newvertices[3*i+1] = vertices[3*ii+1];
    newvertices[3*i+2] = vertices[3*ii+2];
    if(surface->dataflag==1)newtvertices[i]=tvertices[ii];
  }
  FREEMEMORY(vertices);
  surface->vertices = newvertices;
  surface->nvertices = nnewvertices;
  if(surface->dataflag==1){
    FREEMEMORY(tvertices);
    surface->tvertices=newtvertices;
  }
  triangles = surface->triangles;
  for(i=0;i<ntriangles;i++){
    triangles[i] = inverse_vertexmap[triangles[i]];
  }
  FREEMEMORY(ordered_closestnodes);
  FREEMEMORY(vertexmap);
  FREEMEMORY(inverse_vertexmap);

  return 0;

}

/* ------------------ UpdateIsosurface ------------------------ */

int UpdateIsosurface(isosurface *surface, 
                      float *xvert, float *yvert, float *zvert, float *tvert, 
                      int *closestnodes, int nvert, int *triangles, int ntriangles){
  int n,ns, noldvert, *is;
  float *xs=NULL, *ys=NULL, *zs=NULL, *ts=NULL;
  int *cn=NULL;

  if(ResizeSurface(surface,nvert,ntriangles,0)!=0)return 1;

  /* copy vertex data */

  if(nvert>0){
    noldvert = surface->nvertices;
    xs = surface->xvert + noldvert;
    ys = surface->yvert + noldvert;
    zs = surface->zvert + noldvert;
    if(tvert!=NULL)ts = surface->tvert + noldvert;
    cn = surface->closestnodes + noldvert;

    for(n=0;n<nvert;n++){
      xs[n] = xvert[n];
      ys[n] = yvert[n];
      zs[n] = zvert[n];
      if(tvert!=NULL)ts[n]=tvert[n];
      cn[n] = closestnodes[n];
    }
    surface->nvertices = noldvert + nvert;
  }

  /* copy triangle path data */

  if(ntriangles>0){
    ns = surface->ntriangles;
    is = surface->triangles + ns;
    for(n=0;n<ntriangles;n++){is[n] = triangles[n]+noldvert;}
    surface->ntriangles = ns + ntriangles;
  }
  return 0;
}

/* ------------------ ResizeSurface ------------------------ */

int ResizeSurface(isosurface *surfacedata, int incvert, int inctriangles, int incnorm){

  int maxnum, *itemp=NULL;
  float *temp=NULL;

  /* resize vertex data if necessary */

  maxnum = surfacedata->nvertices+incvert;
  if(maxnum>surfacedata->maxvertices){
    maxnum += INCPOINTS;
    surfacedata->maxvertices = maxnum;

    temp = surfacedata->xvert;
    if(temp==NULL){
      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    surfacedata->xvert = temp;

    temp = surfacedata->yvert;
    if(temp==NULL){
      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    surfacedata->yvert = temp;

    temp = surfacedata->zvert;
    if(temp==NULL){
      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    surfacedata->zvert = temp;

    if(surfacedata->dataflag==1){
      temp = surfacedata->tvert;
      if(temp==NULL){
        if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
      }
      else{
        if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
      }
      surfacedata->tvert = temp;
    }

    itemp = surfacedata->closestnodes;
    if(itemp==NULL){
      if(NewMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;
    }
    surfacedata->closestnodes = itemp;
  }

  /* resize norm data if necessary */

  maxnum = surfacedata->nnorm+incnorm;
  if(maxnum>surfacedata->maxnorm){
    maxnum += INCPOINTS;
    surfacedata->maxnorm = maxnum;

    temp = surfacedata->xnorm;
    if(temp==NULL){
      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    surfacedata->xnorm = temp;

    temp = surfacedata->ynorm;
    if(temp==NULL){
      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    surfacedata->ynorm = temp;

    temp = surfacedata->znorm;
    if(temp==NULL){
      if(NewMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&temp,maxnum*sizeof(float))==0)return 1;
    }
    surfacedata->znorm = temp;
  }

  /* resize triangles data if necessary */

  maxnum = surfacedata->ntriangles+inctriangles;
  if(maxnum>surfacedata->maxtriangles){
    maxnum += INCPOINTS;
    surfacedata->maxtriangles = maxnum;

    itemp = surfacedata->triangles;
    if(itemp==NULL){
      if(NewMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;
    }
    else{
      if(ResizeMemory((void **)&itemp,maxnum*sizeof(int))==0)return 1;
    }
    surfacedata->triangles = itemp;
  }
  return 0;

}

/* ------------------ InitIsosurface ------------------------ */

void InitIsosurface(isosurface *surfacedata, float level, float *color,int colorindex){
  surfacedata->level = level;
  surfacedata->color = color;
  surfacedata->colorindex=colorindex;
  surfacedata->defined=1;
  surfacedata->nvertices=0;
  surfacedata->maxvertices=-1;
  surfacedata->ntriangles=0;
  surfacedata->maxtriangles=-1;
  surfacedata->nnorm=0;
  surfacedata->maxnorm=-1;
  surfacedata->normtype=0;
  surfacedata->xvert=NULL;
  surfacedata->yvert=NULL;
  surfacedata->zvert=NULL;
  surfacedata->xnorm=NULL;
  surfacedata->ynorm=NULL;
  surfacedata->znorm=NULL;
  surfacedata->vertices=NULL;
  surfacedata->tvertices=NULL;
  surfacedata->rank=NULL;
  surfacedata->sortedlist=NULL;
  surfacedata->triangles=NULL;
  surfacedata->plottype = 0;
  surfacedata->closestnodes=NULL;
  surfacedata->tvert=NULL;
}

/* ------------------ freesurface ------------------------ */

void freesurface(isosurface *surfacedata){
  if(surfacedata->defined==0)return;
  FREEMEMORY(surfacedata->xvert);
  FREEMEMORY(surfacedata->yvert);
  FREEMEMORY(surfacedata->zvert);
  FREEMEMORY(surfacedata->xnorm);
  FREEMEMORY(surfacedata->ynorm);
  FREEMEMORY(surfacedata->znorm);
  FREEMEMORY(surfacedata->triangles);
  FREEMEMORY(surfacedata->vertices);
  FREEMEMORY(surfacedata->sortedlist);
  FREEMEMORY(surfacedata->rank);
  FREEMEMORY(surfacedata->closestnodes);
  FREEMEMORY(surfacedata->tvert);
  FREEMEMORY(surfacedata->tvertices);
  surfacedata->defined=0;
}

void SmoothIsoSurface(isosurface *surfacedata){
  int *triangles_i;
  unsigned short *vertices_i,*v1,*v2,*v3;
  int n,i1,i2,i3;
  int ntriangles_i,nvertices_i;
  float out[3],area;
  float *xyznorm;
  short *norm,*vertexnorm;

  triangles_i = surfacedata->triangles;
  vertices_i = surfacedata->vertices;
  ntriangles_i = surfacedata->ntriangles;
  nvertices_i = surfacedata->nvertices;

  if(nvertices_i==0||ntriangles_i==0)return;
  NewMemory((void **)&xyznorm,3*nvertices_i*sizeof(float));
  NewMemory((void **)&norm,ntriangles_i*sizeof(short));
  NewMemory((void **)&vertexnorm,3*nvertices_i*sizeof(short));
  surfacedata->norm=norm;
  surfacedata->vertexnorm=vertexnorm;
  for(n=0;n<3*nvertices_i;n++){xyznorm[n]=0.0;}
  for(n=0;n<ntriangles_i/3;n++){
    i1=3*triangles_i[3*n];
    i2=3*triangles_i[3*n+1];
    i3=3*triangles_i[3*n+2];
 /*   for(k=0;k<3;k++){
 //     vv1[k]=xyzmin[k]+xyzmaxdiff*v1[k]/65535.;
 //     vv2[k]=xyzmin[k]+xyzmaxdiff*v2[k]/65535.;
 //     vv3[k]=xyzmin[k]+xyzmaxdiff*v3[k]/65535.;
 //   } */
    v1=vertices_i+i1;
    v2=vertices_i+i2;
    v3=vertices_i+i3;
    calcNormal2(v1,v2,v3,out,&area);
    norm[3*n  ]=(short)(out[0]*32767);
    norm[3*n+1]=(short)(out[1]*32767);
    norm[3*n+2]=(short)(out[2]*32767);
	  xyznorm[i1  ] += out[0]*area;
    xyznorm[i1+1] += out[1]*area;
    xyznorm[i1+2] += out[2]*area;
    xyznorm[i2  ] += out[0]*area;
    xyznorm[i2+1] += out[1]*area;
    xyznorm[i2+2] += out[2]*area;
    xyznorm[i3  ] += out[0]*area;
    xyznorm[i3+1] += out[1]*area;
    xyznorm[i3+2] += out[2]*area;
  }
  for(n=0;n<nvertices_i;n++){
    ReduceToUnit(xyznorm+3*n);
    vertexnorm[3*n  ]=(short)(xyznorm[3*n  ]*32767);
    vertexnorm[3*n+1]=(short)(xyznorm[3*n+1]*32767);
    vertexnorm[3*n+2]=(short)(xyznorm[3*n+2]*32767);
  }
  FREEMEMORY(xyznorm);
}

/* ------------------ GetNormalSurface ------------------------ */

int GetNormalSurface(isosurface *surfacedata){

  float out[3];
  float vertx[3], verty[3], vertz[3];
  float *xnorm=NULL, *ynorm=NULL, *znorm=NULL;
  float *xvert=NULL, *yvert=NULL, *zvert=NULL;
  int *triangles=NULL;
  int ntriangles;
  int nn, n, index;
  
  ntriangles = (surfacedata->ntriangles)/3;
  if(ntriangles==0)return 0;

  xvert = surfacedata->xvert;
  yvert = surfacedata->yvert;
  zvert = surfacedata->zvert;

  triangles = surfacedata->triangles;
  FREEMEMORY(surfacedata->xnorm);
  FREEMEMORY(surfacedata->ynorm);
  FREEMEMORY(surfacedata->znorm);
  if(NewMemory((void **)&surfacedata->xnorm,ntriangles*sizeof(int))==0||
     NewMemory((void **)&surfacedata->ynorm,ntriangles*sizeof(int))==0||
     NewMemory((void **)&surfacedata->znorm,ntriangles*sizeof(int))==0){
    freesurface(surfacedata);
    return 1;
  }

  xnorm = surfacedata->xnorm;
  ynorm = surfacedata->ynorm;
  znorm = surfacedata->znorm;

  nn = 0;
  for(n=0;n<ntriangles;n++){
    index = triangles[nn++];
    vertx[0] = xvert[index]; verty[0] = yvert[index]; vertz[0] = zvert[index];
    index = triangles[nn++];
    vertx[1] = xvert[index]; verty[1] = yvert[index]; vertz[1] = zvert[index];
    index = triangles[nn++];
    vertx[2] = xvert[index]; verty[2] = yvert[index]; vertz[2] = zvert[index];
    calcNormal(vertx,verty,vertz,out);
    xnorm[n]=out[0];
    ynorm[n]=out[1];
    znorm[n]=out[2];
  }
  return 0;

}


/* ------------------ DrawIsosurface ------------------------ */

#ifdef pp_DRAWISO
void DrawIsosurface(isosurface *surfacedata){

  float shiny=10.;
  GLfloat *color=NULL;
  GLfloat specular[]={0.5,0.5,0.5,1.0};


  float *xvert, *yvert, *zvert, *xnorm, *ynorm, *znorm;
  float x0, y0, z0;
  float x, y, z;
  int *triangles, ntriangles;
  int index,n,nn,nvert;

  xvert = surfacedata->xvert;
  yvert = surfacedata->yvert;
  zvert = surfacedata->zvert;
  nvert = surfacedata->nvertices;
  xnorm = surfacedata->xnorm;
  ynorm = surfacedata->ynorm;
  znorm = surfacedata->znorm;
  triangles = surfacedata->triangles;
  ntriangles = surfacedata->ntriangles;

  nn = 0;
  color = surfacedata->color;
  if(surfacedata->plottype==0){

    glPushAttrib(GL_LIGHTING_BIT);

#ifdef ISO_DEBUG
    glLineWidth(1.0);
    glBegin(GL_LINES);
    nn=0;
    for(n=0;n<ntriangles/3;n++){
      index = triangles[nn];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);
      glVertex3f(x+xnorm[n],y+ynorm[n],z+znorm[n]);
      nn +=3;
    }
    glEnd(); 
#endif

    nn = 0;
    glEnable(GL_LIGHTING);
    glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,shiny);
    glBegin(GL_TRIANGLES);
    glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE,color);
    glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,specular);
    for (n = 0; n < ntriangles/3; n++) {
      glNormal3f(xnorm[n],ynorm[n],znorm[n]);

      index = triangles[nn++];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);

      index = triangles[nn++];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);

      index = triangles[nn++];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);
    }
    glEnd();
    glPopAttrib();

  }
  if(surfacedata->plottype==1){
    glColor3fv(color);
    glBegin(GL_LINES);
    for (n = 0; n < ntriangles/3; n++) {
      index = triangles[nn++];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);
      x0 = x; y0 = y; z0 = z;

      index = triangles[nn++];
      x = xvert[index]; y = yvert[index]; z = zvert[index];
      glVertex3f(x,y,z);
      glVertex3f(x,y,z);