Athena: output_vtk.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file output_vtk.c
00004  *  \brief Function to write a single variable in VTK "legacy" format. 
00005  *
00006  * PURPOSE: Function to write a single variable in VTK "legacy" format.  With
00007  *   SMR, dumps are made for all levels and domains, unless nlevel and ndomain
00008  *   are specified in <output> block.
00009  *
00010  * CONTAINS PUBLIC FUNCTIONS: 
00011  * - output_vtk() - writes VTK file (single variable).
00012  *
00013  * PRIVATE FUNCTION PROTOTYPES:
00014  * - output_vtk_2d() - write vtk file for 2D data
00015  * - output_vtk_3d() - write vtk file for 3D data
00016  *============================================================================*/
00017 #include <stdio.h>
00018 #include <stdlib.h>
00019 #include <string.h>
00020 #include "defs.h"
00021 #include "athena.h"
00022 #include "prototypes.h"
00023 
00024 /*==============================================================================
00025  * PRIVATE FUNCTION PROTOTYPES:
00026  *   output_vtk_2d() - write vtk file for 2D data
00027  *   output_vtk_3d() - write vtk file for 3D data
00028  *============================================================================*/
00029 
00030 static void output_vtk_2d(MeshS *pM, OutputS *pOut, int nl, int nd);
00031 static void output_vtk_3d(MeshS *pM, OutputS *pOut, int nl, int nd);
00032 
00033 /*=========================== PUBLIC FUNCTIONS ===============================*/
00034 /*----------------------------------------------------------------------------*/
00035 /*! \fn void output_vtk(MeshS *pM, OutputS *pOut)
00036  *  \brief Writes VTK file (single variable). */
00037 
00038 void output_vtk(MeshS *pM, OutputS *pOut)
00039 {
00040   int nl,nd;
00041 
00042 /* Loop over all Domains in Mesh, and output Grid data */
00043 
00044   for (nl=0; nl<(pM->NLevels); nl++){
00045     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00046       if (pM->Domain[nl][nd].Grid != NULL){
00047 
00048 /* write files if domain and level match input, or are not specified (-1) */
00049       if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00050           (pOut->ndomain == -1 || pOut->ndomain == nd)){
00051 
00052         if (pOut->ndim == 3) {
00053           output_vtk_3d(pM, pOut, nl, nd);
00054         } else if (pOut->ndim == 2) {
00055           output_vtk_2d(pM, pOut, nl, nd);
00056         } else {
00057           ath_error("[output_vtk]: Only able to output 2D or 3D");
00058         }
00059       }}
00060     }
00061   }
00062 
00063   return;
00064 }
00065 
00066 /*----------------------------------------------------------------------------*/
00067 /*! \fn static void output_vtk_2d(MeshS *pM, OutputS *pOut, int nl, int nd)
00068  *  \brief Writes 2D data  */
00069 
00070 static void output_vtk_2d(MeshS *pM, OutputS *pOut, int nl, int nd)
00071 {
00072   GridS *pGrid=pM->Domain[nl][nd].Grid;
00073   FILE *pfile;
00074   char *fname,*plev=NULL,*pdom=NULL;
00075   char levstr[8],domstr[8];
00076 /* Upper and Lower bounds on i,j,k for data dump */
00077   int big_end = ath_big_endian();
00078   int i,j,nx1,nx2;
00079   Real dmin, dmax;
00080   Real **data2d=NULL; /* 2D array of data to be dumped */
00081   float *data;        /* data actually output has to be floats */
00082   double x1, x2, x3, dx1, dx2, dx3;
00083 
00084 /* Allocate memory for and compute 2D array of data */
00085   data2d = OutData2(pGrid,pOut,&nx1,&nx2);
00086   if (data2d == NULL) return; /* data not in range of Grid */
00087 
00088 /* construct output filename.  pOut->id will either be name of variable,
00089  * if 'id=...' was included in <ouput> block, or 'outN' where N is number of
00090  * <output> block.  */
00091   if (nl>0) {
00092     plev = &levstr[0];
00093     sprintf(plev,"lev%d",nl);
00094   }
00095   if (nd>0) {
00096     pdom = &domstr[0];
00097     sprintf(pdom,"dom%d",nd);
00098   }
00099 
00100   if((fname = ath_fname(plev,pM->outfilename,plev,pdom,num_digit,
00101       pOut->num,pOut->id,"vtk")) == NULL){
00102     ath_error("[output_vtk]: Error constructing filename\n");
00103   }
00104 
00105 /* open output file */
00106   if((pfile = fopen(fname,"w")) == NULL){
00107     ath_error("[output_vtk]: Unable to open vtk file %s\n",fname);
00108   }
00109 
00110 /* Store the global min / max, for output at end of run */
00111   minmax2(data2d,nx2,nx1,&dmin,&dmax);
00112   pOut->gmin = MIN(dmin,pOut->gmin);
00113   pOut->gmax = MAX(dmax,pOut->gmax);
00114 
00115 /* Allocate memory for temporary array of floats */
00116 
00117   if((data = (float *)malloc(nx1*sizeof(float))) == NULL){
00118      ath_error("[output_vtk]: malloc failed for temporary array\n");
00119      return;
00120   }
00121 
00122 /* There are five basic parts to the VTK "legacy" file format.  */
00123 /*  1. Write file version and identifier */
00124 
00125   fprintf(pfile,"# vtk DataFile Version 2.0\n");
00126 
00127 /*  2. Header */
00128 
00129   fprintf(pfile,"Really cool Athena data at time= %e, level= %i, domain= %i\n",
00130     pGrid->time,nl,nd);
00131 
00132 /*  3. File format */
00133 
00134   fprintf(pfile,"BINARY\n");
00135 
00136 /*  4. Dataset structure */
00137 
00138 /* Set the Grid origin */
00139   x1 = pGrid->MinX[0];
00140   x2 = pGrid->MinX[1];
00141   x3 = pGrid->MinX[2];
00142   dx1 = (pOut->reduce_x1 == 1 ? pGrid->dx1 * pGrid->Nx[0] : pGrid->dx1);
00143   dx2 = (pOut->reduce_x2 == 1 ? pGrid->dx2 * pGrid->Nx[1] : pGrid->dx2);
00144   dx3 = (pOut->reduce_x3 == 1 ? pGrid->dx3 * pGrid->Nx[2] : pGrid->dx3);
00145 
00146   fprintf(pfile,"DATASET STRUCTURED_POINTS\n");
00147   fprintf(pfile,"DIMENSIONS %d %d %d\n",nx1+1,nx2+1,1);
00148   fprintf(pfile,"ORIGIN %e %e %e \n",x1,x2,x3);
00149   fprintf(pfile,"SPACING %e %e %e \n",dx1,dx2,dx3);
00150 
00151 /*  5. Data  */
00152 
00153   fprintf(pfile,"CELL_DATA %d \n", nx1*nx2);
00154 
00155 /* Write data */
00156 
00157   fprintf(pfile,"SCALARS %s float\n", pOut->id);
00158   fprintf(pfile,"LOOKUP_TABLE default\n");
00159   for (j=0; j<nx2; j++){
00160     for (i=0; i<nx1; i++){
00161       data[i] = (float)data2d[j][i];
00162     }
00163     if(!big_end) ath_bswap(data,sizeof(float),nx1);
00164     fwrite(data,sizeof(float),(size_t)nx1,pfile);
00165   }
00166 
00167 /* close file and free memory */
00168 
00169   fclose(pfile);
00170   free(data);
00171   free_2d_array(data2d);
00172   return;
00173 }
00174 
00175 /*----------------------------------------------------------------------------*/
00176 /*! \fn static void output_vtk_3d(MeshS *pM, OutputS *pOut, int nl, int nd)
00177  *  \brief Writes 3D data  */
00178 
00179 static void output_vtk_3d(MeshS *pM, OutputS *pOut, int nl, int nd)
00180 {
00181   GridS *pGrid=pM->Domain[nl][nd].Grid;
00182   FILE *pfile;
00183   char *fname,*plev=NULL,*pdom=NULL;
00184   char levstr[8],domstr[8];
00185 /* Upper and Lower bounds on i,j,k for data dump */
00186   int big_end = ath_big_endian();
00187   int nx1,nx2,nx3,i,j,k;
00188   Real dmin, dmax;
00189   Real ***data3d=NULL; /* 3D array of data to be dumped */
00190   float *data;         /* data actually output has to be floats */
00191   double x1, x2, x3;
00192 
00193 /* Allocate memory for and compute 3D array of data values */
00194   data3d = OutData3(pGrid,pOut,&nx1,&nx2,&nx3);
00195 
00196 /* construct output filename.  pOut->id will either be name of variable,
00197  * if 'id=...' was included in <ouput> block, or 'outN' where N is number of
00198  * <output> block.  */
00199   if (nl>0) {
00200     plev = &levstr[0];
00201     sprintf(plev,"lev%d",nl);
00202   }
00203   if (nd>0) {
00204     pdom = &domstr[0];
00205     sprintf(pdom,"dom%d",nd);
00206   }
00207 
00208   if((fname = ath_fname(plev,pM->outfilename,plev,pdom,num_digit,
00209       pOut->num,pOut->id,"vtk")) == NULL){
00210     ath_error("[output_vtk]: Error constructing filename\n");
00211   }
00212 
00213 /* open output file */
00214   if((pfile = fopen(fname,"w")) == NULL){
00215     ath_error("[output_vtk]: Unable to open vtk file %s\n",fname);
00216   }
00217 
00218 /* Store the global min / max, for output at end of run */
00219   minmax3(data3d,nx3,nx2,nx1,&dmin,&dmax);
00220   pOut->gmin = MIN(dmin,pOut->gmin);
00221   pOut->gmax = MAX(dmax,pOut->gmax);
00222 
00223 /* Allocate memory for temporary array of floats */
00224 
00225   x1 = pGrid->MinX[0];
00226   x2 = pGrid->MinX[1];
00227   x3 = pGrid->MinX[2];
00228   if((data = (float *)malloc(nx1*sizeof(float))) == NULL){
00229      ath_error("[output_vtk]: malloc failed for temporary array\n");
00230      return;
00231   }
00232 
00233 /* There are five basic parts to the VTK "legacy" file format.  */
00234 /*  1. Write file version and identifier */
00235 
00236   fprintf(pfile,"# vtk DataFile Version 2.0\n");
00237 
00238 /*  2. Header */
00239 
00240   fprintf(pfile,"Really cool Athena data at time= %e, level= %i, domain= %i\n",
00241     pGrid->time,nl,nd);
00242 
00243 /*  3. File format */
00244 
00245   fprintf(pfile,"BINARY\n");
00246 
00247 /*  4. Dataset structure */
00248 
00249   fprintf(pfile,"DATASET STRUCTURED_POINTS\n");
00250   fprintf(pfile,"DIMENSIONS %d %d %d\n",nx1+1,nx2+1,nx3+1);
00251   fprintf(pfile,"ORIGIN %e %e %e \n",x1,x2,x3);
00252   fprintf(pfile,"SPACING %e %e %e \n",pGrid->dx1,pGrid->dx2,pGrid->dx3);
00253 
00254 /*  5. Data  */
00255 
00256   fprintf(pfile,"CELL_DATA %d \n", nx1*nx2*nx3);
00257 
00258 /* Write data */
00259 
00260   fprintf(pfile,"SCALARS %s float\n", pOut->id);
00261   fprintf(pfile,"LOOKUP_TABLE default\n");
00262   for (k=0; k<nx3; k++) {
00263     for (j=0; j<nx2; j++) {
00264       for (i=0; i<nx1; i++) {
00265         data[i] = (float)data3d[k][j][i];
00266       }
00267       if(!big_end) ath_bswap(data,sizeof(float),nx1);
00268       fwrite(data,sizeof(float),(size_t)nx1,pfile);
00269     }
00270   }
00271 
00272 /* close file and free memory */
00273 
00274   fclose(pfile);
00275   free(data);
00276   free_3d_array(data3d);
00277   return;
00278 }