#include <stdio.h>#include <stdlib.h>#include "datadef.h"#include "visual.h"/*---------------------------------------------------------------------------*//* Writing U,V,P,PSI, and ZETA into "vecfile" for visualization              *//*---------------------------------------------------------------------------*/void OUTPUTVEC_bin(REAL **U,REAL **V,REAL **P,REAL **TEMP,		   REAL **PSI,REAL **ZETA,REAL **HEAT,int **FLAG,                   REAL xlength,REAL ylength,int imax,int jmax,                   char* vecfile){ int i,j; float temp; FILE *fp; fp = fopen(vecfile, "wb"); temp = xlength; fwrite(&temp, sizeof(float), 1, fp); temp = ylength; fwrite(&temp, sizeof(float), 1, fp); temp = imax; fwrite(&temp, sizeof(float), 1, fp); temp = jmax; fwrite(&temp, sizeof(float), 1, fp); for(j=1;j<=jmax;j+=1)  for(i=1;i<=imax;i+=1){   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) )       temp = (U[i][j]+U[i-1][j])/2.0;   else			 temp = 0.0;   fwrite(&temp, sizeof(float), 1, fp);  } for(j=1;j<=jmax;j+=1)  for(i=1;i<=imax;i+=1){   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) )       temp = (V[i][j]+V[i][j-1])/2.0;   else			 temp = 0.0;   fwrite(&temp, sizeof(float), 1, fp);  } for(j=1;j<=jmax;j+=1)  for(i=1;i<=imax;i+=1){   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) )       temp = P[i][j];   else			 temp = 0.0;   fwrite(&temp, sizeof(float), 1, fp);  } for(j=1;j<=jmax;j+=1)  for(i=1;i<=imax;i+=1){   if( (FLAG[i][j] & C_F) && (FLAG[i][j] < C_E) )         temp = TEMP[i][j];   else			 temp = -0.5;   temp = TEMP[i][j];   fwrite(&temp, sizeof(float), 1, fp);  } for(j=1;j<=jmax-1;j+=1)  for(i=1;i<=imax-1;i+=1){    temp = ZETA[i][j];    fwrite(&temp, sizeof(float), 1, fp);  } for(j=0;j<=jmax;j+=1)  for(i=0;i<=imax;i+=1){    temp = PSI[i][j];    fwrite(&temp, sizeof(float), 1, fp);  } for(j=0;j<=jmax;j+=1)  for(i=0;i<=imax;i+=1){    temp = HEAT[i][j];    fwrite(&temp, sizeof(float), 1, fp);  } fclose(fp);}/*-----------------------------------------------------*//* Computation of stream function and vorticity        *//*-----------------------------------------------------*/void COMPPSIZETA(REAL **U,REAL **V,REAL **PSI,REAL **ZETA,int **FLAG,                int imax,int jmax,REAL delx,REAL dely){ int i,j; /* Computation of the vorticity zeta at the upper right corner     */ /* of cell (i,j) (only if the corner is surrounded by fluid cells) */ /*-----------------------------------------------------------------*/ for (i=1;i<=imax-1;i++)    for (j=1;j<=jmax-1;j++)        if( ((FLAG[i][j] & C_F) && (FLAG[i][j] < C_E))  &&            ((FLAG[i+1][j] & C_F) && (FLAG[i+1][j] < C_E))  &&            ((FLAG[i][j+1] & C_F) && (FLAG[i][j+1] < C_E))  &&            ((FLAG[i+1][j+1] & C_F) && (FLAG[i+1][j+1] < C_E)) )          ZETA[i][j] = (U[i][j+1]-U[i][j])/dely - (V[i+1][j]-V[i][j])/delx;       else          ZETA[i][j] = 0.0; /* Computation of the stream function at the upper right corner    */ /* of cell (i,j) (only if bother lower cells are fluid cells)      */ /*-----------------------------------------------------------------*/ for (i=0;i<=imax;i++)   {    PSI[i][0] = 0.0;    for(j=1;j<=jmax;j++)        if( ((FLAG[i][j] & C_F) && (FLAG[i][j] < C_E))  ||            ((FLAG[i+1][j] & C_F) && (FLAG[i+1][j] < C_E)) )          PSI[i][j] = PSI[i][j-1] + U[i][j]*dely;       else          PSI[i][j] = PSI[i][j-1];  } }/*-----------------------------------------------------*//* Computation of the heat function                    *//*-----------------------------------------------------*/void COMP_HEAT(REAL **U,REAL **V,REAL **TEMP,REAL **HEAT,int **FLAG,               REAL Re,REAL Pr,int imax,int jmax,REAL delx,REAL dely){ int i,j; /* Computation at the upper right corner of cell (i,j) */ /*-----------------------------------------------------*/ for (i=0;i<=imax;i++)   {    HEAT[i][0] = 0.0;    for(j=1;j<=jmax;j++)       if( ((FLAG[i][j] & C_F) && (FLAG[i][j] < C_E))  ||           ((FLAG[i+1][j] & C_F) && (FLAG[i+1][j] < C_E)) )          HEAT[i][j] = HEAT[i][j-1] +                     dely*(U[i][j]*0.5*(1.0+TEMP[i+1][j]+TEMP[i][j])*Re*Pr-                          (TEMP[i+1][j]-TEMP[i][j])/delx );   }}/*-----------------------------------------------------------------*//* Allocate memory for a particle and set the coordinates to (x,y) *//*-----------------------------------------------------------------*/struct particle *PARTALLOC(REAL x, REAL y){  struct particle *part;  if((part=(struct particle *)malloc(sizeof(struct particle))) == NULL)    {	     printf("no memory\n");     exit(0);    }  part->x = x; part->y = y;  part->next = NULL;  return( part );}/*------------------------------------------------------------------*//* Allocate memory for a particleline and set coordinates where     *//* particles are injected.                                          *//*------------------------------------------------------------------*/struct particleline *SET_PARTICLES(int N,REAL pos1x,REAL pos1y,				   	 REAL pos2x,REAL pos2y){   int i;   REAL hx,hy;   struct particleline *Partlines;   if((Partlines=(struct particleline *)                 malloc((unsigned)(N) * sizeof(struct particleline))) == NULL)     {	      printf("no memory\n");      exit(0);     }   Partlines -= 1;     if (N>=2)     {      hx  = (pos2x-pos1x)/(N-1);      hy  = (pos2y-pos1y)/(N-1);      for(i=1; i<=N; i++){         Partlines[i].length = 0;         Partlines[i].Particles =              PARTALLOC(pos1x+hx*(i-1),pos1y+hy*(i-1));         Partlines[i].Particles->next =              PARTALLOC(pos1x+hx*(i-1),pos1y+hy*(i-1));         Partlines[i].length++;       }     }   return(Partlines);}/*End SET_PARTICLES*//*-------------------------------------------------------------*//* Moving particles                                            *//*-------------------------------------------------------------*/void ADVANCE_PARTICLES(int imax,int jmax,REAL delx,REAL dely,REAL delt,                       REAL **U,REAL **V,int **FLAG,                       int N,struct particleline *Partlines){  int i, j, k;  REAL x, y, x1, y1, x2, y2, u, v;   struct particle *part,*help;  for(k=1;k<=N;k++){     for(part=Partlines[k].Particles; part->next != NULL; part=part->next){        /* first element is only a dummy element           */        /*-------------------------------------------------*/        x = part->next->x; y = part->next->y;        /* Computation of new x-coordinates by discretizing dx/dt=u */        /*----------------------------------------------------------*/        i = (int)(x/delx)+1;  j = (int)((y+0.5*dely)/dely)+1;        x1 = (i-1)*delx;    y1 = ((j-1)-0.5)*dely;        x2 = i*delx;        y2 = (j-0.5)*dely;        /* bilinear interpolation */        /*------------------------*/        u = ((x2-x)*(y2-y)*U[i-1][j-1] +  	     (x-x1)*(y2-y)*U[i][j-1]   +	     (x2-x)*(y-y1)*U[i-1][j]   +	     (x-x1)*(y-y1)*U[i][j])/delx/dely;        /* Computation of new y-coordinates by discretizing dy/dt=v */        /*----------------------------------------------------------*/        i = (int)((x+0.5*delx)/delx)+1; j = (int)(y/dely)+1;        x1 = ((i-1)-0.5)*delx;    y1 = (j-1)*dely;        x2 = (i-0.5)*delx;        y2 = j*dely;        /* bilinear interpolation */        /*------------------------*/        v = ((x2-x)*(y2-y)*V[i-1][j-1] +	     (x-x1)*(y2-y)*V[i][j-1]   +	     (x2-x)*(y-y1)*V[i-1][j]   +	     (x-x1)*(y-y1)*V[i][j])/delx/dely;        x += delt*u;   y += delt*v;         /*-------------------------------------*/        /* determine new cell for the particle */        /*-------------------------------------*/        i = (int)(x/delx)+1;   j = (int)(y/dely)+1;        /* if particle left the fluid domain, delete it */        /*----------------------------------------------*/        if (x>=imax*delx || y>=jmax*dely || x<=0 || y<=0){          help = part->next->next;          free(part->next);          part->next = help;          Partlines[k].length--;	  if (help == NULL)	    break;         }        else{          /*-------------------------------------*/          /* special treatment if particle would */	  /* be in an inner obstacle cell        */