00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00026
00027
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
00034
00035
00036
00037
00038 void output_vtk(MeshS *pM, OutputS *pOut)
00039 {
00040 int nl,nd;
00041
00042
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
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
00068
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
00077 int big_end = ath_big_endian();
00078 int i,j,nx1,nx2;
00079 Real dmin, dmax;
00080 Real **data2d=NULL;
00081 float *data;
00082 double x1, x2, x3, dx1, dx2, dx3;
00083
00084
00085 data2d = OutData2(pGrid,pOut,&nx1,&nx2);
00086 if (data2d == NULL) return;
00087
00088
00089
00090
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
00106 if((pfile = fopen(fname,"w")) == NULL){
00107 ath_error("[output_vtk]: Unable to open vtk file %s\n",fname);
00108 }
00109
00110
00111 minmax2(data2d,nx2,nx1,&dmin,&dmax);
00112 pOut->gmin = MIN(dmin,pOut->gmin);
00113 pOut->gmax = MAX(dmax,pOut->gmax);
00114
00115
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
00123
00124
00125 fprintf(pfile,"# vtk DataFile Version 2.0\n");
00126
00127
00128
00129 fprintf(pfile,"Really cool Athena data at time= %e, level= %i, domain= %i\n",
00130 pGrid->time,nl,nd);
00131
00132
00133
00134 fprintf(pfile,"BINARY\n");
00135
00136
00137
00138
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
00152
00153 fprintf(pfile,"CELL_DATA %d \n", nx1*nx2);
00154
00155
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
00168
00169 fclose(pfile);
00170 free(data);
00171 free_2d_array(data2d);
00172 return;
00173 }
00174
00175
00176
00177
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
00186 int big_end = ath_big_endian();
00187 int nx1,nx2,nx3,i,j,k;
00188 Real dmin, dmax;
00189 Real ***data3d=NULL;
00190 float *data;
00191 double x1, x2, x3;
00192
00193
00194 data3d = OutData3(pGrid,pOut,&nx1,&nx2,&nx3);
00195
00196
00197
00198
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
00214 if((pfile = fopen(fname,"w")) == NULL){
00215 ath_error("[output_vtk]: Unable to open vtk file %s\n",fname);
00216 }
00217
00218
00219 minmax3(data3d,nx3,nx2,nx1,&dmin,&dmax);
00220 pOut->gmin = MIN(dmin,pOut->gmin);
00221 pOut->gmax = MAX(dmax,pOut->gmax);
00222
00223
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
00234
00235
00236 fprintf(pfile,"# vtk DataFile Version 2.0\n");
00237
00238
00239
00240 fprintf(pfile,"Really cool Athena data at time= %e, level= %i, domain= %i\n",
00241 pGrid->time,nl,nd);
00242
00243
00244
00245 fprintf(pfile,"BINARY\n");
00246
00247
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
00255
00256 fprintf(pfile,"CELL_DATA %d \n", nx1*nx2*nx3);
00257
00258
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
00273
00274 fclose(pfile);
00275 free(data);
00276 free_3d_array(data3d);
00277 return;
00278 }