PLUTO
Macros | Functions
write_vtk.c File Reference

Write data in VTK format. More...

#include "pluto.h"
Include dependency graph for write_vtk.c:

Go to the source code of this file.

Macros

#define RECTILINEAR_GRID   14
 
#define STRUCTURED_GRID   35
 
#define VTK_FORMAT   STRUCTURED_GRID
 
#define VTK_TIME_INFO   NO /* VisIt to display data results */
 
#define VTK_HEADER_WRITE_STRING(header)   fprintf (fvtk,header);
 
#define VTK_HEADER_WRITE_FLTARR(arr, nelem)   fwrite(arr, sizeof(float), nelem, fvtk);
 
#define VTK_HEADER_WRITE_DBLARR(arr, nelem)   fwrite(arr, sizeof(double), nelem, fvtk);
 

Functions

void WriteVTK_Header (FILE *fvtk, Grid *grid)
 
void WriteVTK_Vector (FILE *fvtk, Data_Arr V, double unit, char *var_name, Grid *grid)
 
void WriteVTK_Scalar (FILE *fvtk, double ***V, double unit, char *var_name, Grid *grid)
 

Detailed Description

Write data in VTK format.

Collection of basic functions to write VTK files using the simple legacy format. Files can be written in serial or parallel mode and consists of the basic five parts :

  1. File version and identifier
  2. Header consisting of a string
  3. File format
  4. Dataset structure: describes the gwometry and topology of the dataset.
  5. Dataset attributes. This section is used to write the actual binary data as a vector or scalar data.

The mesh topology and the variable datasets are written using single precision (4 bytes) binary format. VTK file are usually written following big endian order. Therefore, we swap endianity only if local architecture has little endian ordering.

The WriteVTK_Header() function provides the basic functionality for steps 1, 2, 3 and 4. Only processor 0 does the actual writing. For cartesian/cylindrical geometries the default grid topology is "RECTILINEAR_GRIDS" whereas for polar/spherical we employ "STRUCTURED_GRID" to provide a convenient mapping to a cartesian mesh.
Note for 2D datasets: in order to produce a strictly 2D dataset we always set the third coordinate (x3) to zero. For this reason, in 2D spherical cordinates we swap the role of the "y" and "z" coordinates.

The WriteVTK_Vector() is fully parallel and is used to write data with the vector attribute (by default these include velocity and magnetic fields).

The WriteVTK_Scalar() is fully parallel and is used to write data with the scalar attribute (by default these include density, pressure, tracer and user-defined variables).

Reference

http://www.vtk.org/VTK/img/file-formats.pdf

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Aug 17, 2015

Definition in file write_vtk.c.

Macro Definition Documentation

#define RECTILINEAR_GRID   14

Definition at line 53 of file write_vtk.c.

#define STRUCTURED_GRID   35

Definition at line 54 of file write_vtk.c.

#define VTK_FORMAT   STRUCTURED_GRID

Definition at line 60 of file write_vtk.c.

#define VTK_HEADER_WRITE_DBLARR (   arr,
  nelem 
)    fwrite(arr, sizeof(double), nelem, fvtk);

Definition at line 86 of file write_vtk.c.

#define VTK_HEADER_WRITE_FLTARR (   arr,
  nelem 
)    fwrite(arr, sizeof(float), nelem, fvtk);

Definition at line 84 of file write_vtk.c.

#define VTK_HEADER_WRITE_STRING (   header)    fprintf (fvtk,header);

Definition at line 82 of file write_vtk.c.

#define VTK_TIME_INFO   NO /* VisIt to display data results */

Definition at line 65 of file write_vtk.c.

Function Documentation

void WriteVTK_Header ( FILE *  fvtk,
Grid *  grid 
)

Write VTK header in parallel or serial mode. In parallel mode only processor 0 does the actual writing (see al_io.c/AL_Write_header).

Parameters
[in]fvtkpointer to file
[in]gridpointer to an array of Grid structures
Todo:
Write the grid using several processors.

Definition at line 91 of file write_vtk.c.

103 {
104  int i, j, k;
105  int nx1, nx2, nx3;
106  char header[1024];
107  float x1, x2, x3;
108  static float ***node_coord, *xnode, *ynode, *znode;
109 
110 /* ------------------------------------------------
111  get global dimensions
112  ------------------------------------------------ */
113 
114  nx1 = grid[IDIR].gend + 1 - grid[IDIR].nghost;
115  nx2 = grid[JDIR].gend + 1 - grid[JDIR].nghost;
116  nx3 = grid[KDIR].gend + 1 - grid[KDIR].nghost;
117 
118 /* -------------------------------------------------------------
119  Allocate memory and define node coordinates only once.
120  ------------------------------------------------------------- */
121 
122  if (node_coord == NULL){
123  node_coord = ARRAY_3D(nx2 + JOFFSET, nx1 + IOFFSET, 3, float);
124 
125  #if VTK_FORMAT == RECTILINEAR_GRID
126  xnode = ARRAY_1D(nx1 + IOFFSET, float);
127  ynode = ARRAY_1D(nx2 + JOFFSET, float);
128  znode = ARRAY_1D(nx3 + KOFFSET, float);
129 
130  for (i = 0; i < nx1 + IOFFSET; i++){
131  x1 = (float)(grid[IDIR].xl_glob[i+IBEG]);
132  if (IsLittleEndian()) SWAP_VAR(x1);
133  xnode[i] = x1;
134  }
135  for (j = 0; j < nx2 + JOFFSET; j++){
136  x2 = (float)(grid[JDIR].xl_glob[j+JBEG]);
137  if (IsLittleEndian()) SWAP_VAR(x2);
138  ynode[j] = x2;
139  }
140  for (k = 0; k < nx3 + KOFFSET; k++){
141  x3 = (float)(grid[KDIR].xl_glob[k+KBEG]);
142  if (IsLittleEndian()) SWAP_VAR(x3);
143  #if DIMENSIONS == 2
144  znode[k] = 0.0;
145  #else
146  znode[k] = x3;
147  #endif
148  }
149  #endif
150  }
151 
152 /* ----------------------------------------------------------
153  Part I, II, III:
154  Write file header on string "header"
155  ---------------------------------------------------------- */
156 
157  sprintf(header,"# vtk DataFile Version 2.0\n");
158  sprintf(header+strlen(header),"PLUTO %s VTK Data\n",PLUTO_VERSION);
159  sprintf(header+strlen(header),"BINARY\n");
160  #if VTK_FORMAT == RECTILINEAR_GRID
161  sprintf(header+strlen(header),"DATASET %s\n","RECTILINEAR_GRID");
162  #elif VTK_FORMAT == STRUCTURED_GRID
163  sprintf(header+strlen(header),"DATASET %s\n","STRUCTURED_GRID");
164  #endif
165 
166  VTK_HEADER_WRITE_STRING(header);
167 
168  /* -- generate time info -- */
169 
170  #if VTK_TIME_INFO == YES
171  sprintf (header,"FIELD FieldData 1\n");
172  sprintf (header+strlen(header),"TIME 1 1 double\n");
173  double tt=g_time;
174  if (IsLittleEndian()) SWAP_VAR(tt);
175  VTK_HEADER_WRITE_STRING(header);
176  VTK_HEADER_WRITE_DBLARR(&tt, 1);
177  VTK_HEADER_WRITE_STRING("\n");
178 
179  #endif /* VTK_TIME_INFO */
180 
181  sprintf(header,"DIMENSIONS %d %d %d\n",
182  nx1 + IOFFSET, nx2 + JOFFSET, nx3 + KOFFSET);
183  VTK_HEADER_WRITE_STRING(header);
184 
185 #if VTK_FORMAT == RECTILINEAR_GRID
186 
187  /* -- reset header string and keep going -- */
188 
189  sprintf(header,"X_COORDINATES %d float\n", nx1 + IOFFSET);
190  VTK_HEADER_WRITE_STRING(header);
191  VTK_HEADER_WRITE_FLTARR(xnode, nx1 + IOFFSET);
192 
193  sprintf(header,"\nY_COORDINATES %d float\n", nx2 + JOFFSET);
194  VTK_HEADER_WRITE_STRING(header);
195  VTK_HEADER_WRITE_FLTARR(ynode, nx2 + JOFFSET);
196 
197  sprintf(header,"\nZ_COORDINATES %d float\n", nx3 + KOFFSET);
198  VTK_HEADER_WRITE_STRING(header);
199  VTK_HEADER_WRITE_FLTARR(znode, nx3 + KOFFSET);
200 
201  sprintf (header,"\nCELL_DATA %d\n", nx1*nx2*nx3);
202  VTK_HEADER_WRITE_STRING (header);
203 
204 #elif VTK_FORMAT == STRUCTURED_GRID
205 
206  sprintf(header,"POINTS %d float\n", (nx1+IOFFSET)*(nx2+JOFFSET)*(nx3+KOFFSET));
207  VTK_HEADER_WRITE_STRING(header);
208 
209 /* ---------------------------------------------------------------
210  Part IV: (structured) grid information
211  --------------------------------------------------------------- */
212 
213  x1 = x2 = x3 = 0.0;
214  for (k = 0; k < nx3 + KOFFSET; k++){
215  for (j = 0; j < nx2 + JOFFSET; j++){
216  for (i = 0; i < nx1 + IOFFSET; i++){
217  D_EXPAND(x1 = grid[IDIR].xl_glob[IBEG + i]; ,
218  x2 = grid[JDIR].xl_glob[JBEG + j]; ,
219  x3 = grid[KDIR].xl_glob[KBEG + k];)
220 
221  #if (GEOMETRY == CARTESIAN) || (GEOMETRY == CYLINDRICAL)
222  node_coord[j][i][0] = x1;
223  node_coord[j][i][1] = x2;
224  node_coord[j][i][2] = x3;
225  #elif GEOMETRY == POLAR
226  node_coord[j][i][0] = x1*cos(x2);
227  node_coord[j][i][1] = x1*sin(x2);
228  node_coord[j][i][2] = x3;
229  #elif GEOMETRY == SPHERICAL
230  #if DIMENSIONS == 2
231  node_coord[j][i][0] = x1*sin(x2);
232  node_coord[j][i][1] = x1*cos(x2);
233  node_coord[j][i][2] = 0.0;
234  #elif DIMENSIONS == 3
235  node_coord[j][i][0] = x1*sin(x2)*cos(x3);
236  node_coord[j][i][1] = x1*sin(x2)*sin(x3);
237  node_coord[j][i][2] = x1*cos(x2);
238  #endif
239  #endif
240 
241  if (IsLittleEndian()){
242  SWAP_VAR(node_coord[j][i][0]);
243  SWAP_VAR(node_coord[j][i][1]);
244  SWAP_VAR(node_coord[j][i][2]);
245  }
246  }}
247  VTK_HEADER_WRITE_FLTARR(node_coord[0][0],3*(nx1+IOFFSET)*(nx2+JOFFSET));
248  }
249 
250  sprintf (header,"\nCELL_DATA %d\n", nx1*nx2*nx3);
251  VTK_HEADER_WRITE_STRING(header);
252 
253 #endif
254 }
VTK_HEADER_WRITE_STRING
#define VTK_HEADER_WRITE_STRING(header)
Definition: write_vtk.c:82
VTK_HEADER_WRITE_FLTARR
#define VTK_HEADER_WRITE_FLTARR(arr, nelem)
Definition: write_vtk.c:84
PLUTO_VERSION
#define PLUTO_VERSION
Definition: pluto.h:16
GRID::gend
int gend
Global end index for the global array.
Definition: structs.h:114
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
KDIR
#define KDIR
Definition: pluto.h:195
SWAP_VAR
#define SWAP_VAR(x)
Definition: macros.h:115
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
IsLittleEndian
int IsLittleEndian(void)
Definition: tools.c:40
VTK_HEADER_WRITE_DBLARR
#define VTK_HEADER_WRITE_DBLARR(arr, nelem)
Definition: write_vtk.c:86
GRID::nghost
int nghost
Number of ghost zones.
Definition: structs.h:104
k
int k
Definition: analysis.c:2
GRID::xl_glob
double * xl_glob
Cell left interface.
Definition: structs.h:82
GEOMETRY
#define GEOMETRY
Definition: definitions_01.h:4
CARTESIAN
#define CARTESIAN
Definition: pluto.h:33
CYLINDRICAL
#define CYLINDRICAL
Definition: pluto.h:34
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
i
int i
Definition: analysis.c:2
g_time
double g_time
The current integration time.
Definition: globals.h:117
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35

Here is the call graph for this function:

Here is the caller graph for this function:

void WriteVTK_Scalar ( FILE *  fvtk,
double ***  V,
double  unit,
char *  var_name,
Grid *  grid 
)

Write VTK scalar field.

Parameters
[in]fvtkpointer to file (handle)
[in]Vpointer to 3D data array
[in]unitthe corresponding cgs unit (if specified, 1 otherwise)
[in]var_namethe variable name appearing in the VTK file
[in]gridpointer to an array of Grid structures

Definition at line 341 of file write_vtk.c.

352 {
353  int i,j,k;
354  char header[512];
355  float ***Vflt;
356 
357  sprintf (header,"\nSCALARS %s float\n", var_name);
358  sprintf (header,"%sLOOKUP_TABLE default\n",header);
359 
360  #ifdef PARALLEL
361  MPI_Barrier (MPI_COMM_WORLD);
362  AL_Write_header (header, strlen(header), MPI_CHAR, SZ_float);
363  #else
364  fprintf (fvtk, "%s",header);
365  #endif
366 
367  Vflt = Convert_dbl2flt(V, unit, IsLittleEndian());
368  WriteBinaryArray (Vflt[0][0], sizeof(float), SZ_float, fvtk, -1);
369 }
AL_Write_header
int AL_Write_header(void *vbuffer, int nelem, AL_Datatype type, int sz_ptr)
Definition: al_io.c:285
SZ_float
int SZ_float
Definition: globals.h:27
Convert_dbl2flt
float *** Convert_dbl2flt(double ***Vdbl, double unit, int swap_endian)
Definition: bin_io.c:216
j
int j
Definition: analysis.c:2
IsLittleEndian
int IsLittleEndian(void)
Definition: tools.c:40
k
int k
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
WriteBinaryArray
void WriteBinaryArray(void *V, size_t dsize, int sz, FILE *fl, int istag)
Definition: bin_io.c:98

Here is the call graph for this function:

Here is the caller graph for this function:

void WriteVTK_Vector ( FILE *  fvtk,
Data_Arr  V,
double  unit,
char *  var_name,
Grid *  grid 
)

Write VTK vector field data. This is enabled only when VTK_VECTOR_DUMP is set to YES. For generality purposes, vectors are written always with 3 components, even when there're only 2 being used.

The following Maple script has been used to find vector components from cyl/sph to cartesian:

1 restart;
2 with(linalg);
3 Acyl := matrix (3,3,[ cos(phi), sin(phi), 0,
4  -sin(phi), cos(phi), 0,
5  0 , 0 , 1]);
6 Asph := matrix (3,3,[ sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta),
7  cos(theta)*cos(phi), cos(theta)*sin(phi), -sin(theta),
8  -sin(phi) , cos(phi) , 0]);
9 Bcyl := simplify(inverse(Acyl));
10 Bsph := simplify(inverse(Asph));
Parameters
[in]fvtkpointer to file
[in]Va 4D array [nv][k][j][i] containing the vector components (nv) defined at cell centers (k,j,i). The index nv = 0 marks the vector first component.
[in]unitthe corresponding cgs unit (if specified, 1 otherwise)
[in]var_namethe variable name appearing in the VTK file
[in]gridpointer to an array of Grid structures

Definition at line 259 of file write_vtk.c.

291 {
292  int i,j,k;
293  int vel_field, mag_field;
294  char header[512];
295  static Float_Vect ***vect3D;
296  double v[3], x1, x2, x3;
297 
298  if (vect3D == NULL){
299  vect3D = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, Float_Vect);
300  }
301 
302 /* --------------------------------------------------------
303  Write VTK vector fields
304  -------------------------------------------------------- */
305 
306  v[0] = v[1] = v[2] = 0.0;
307  x1 = x2 = x3 = 0.0;
308 
309  vel_field = (strcmp(var_name,"vx1") == 0);
310  mag_field = (strcmp(var_name,"bx1") == 0);
311  if (vel_field || mag_field) {
312  DOM_LOOP(k,j,i){
313  D_EXPAND(v[0] = V[0][k][j][i]; x1 = grid[IDIR].x[i]; ,
314  v[1] = V[1][k][j][i]; x2 = grid[JDIR].x[j]; ,
315  v[2] = V[2][k][j][i]; x3 = grid[KDIR].x[k];)
316 
317  VectorCartesianComponents(v, x1, x2, x3);
318  vect3D[k][j][i].v1 = (float)v[0]*unit;
319  vect3D[k][j][i].v2 = (float)v[1]*unit;
320  vect3D[k][j][i].v3 = (float)v[2]*unit;
321 
322  if (IsLittleEndian()){
323  SWAP_VAR(vect3D[k][j][i].v1);
324  SWAP_VAR(vect3D[k][j][i].v2);
325  SWAP_VAR(vect3D[k][j][i].v3);
326  }
327 
328  } /* endfor DOM_LOOP(k,j,i) */
329 
330  if (vel_field)
331  sprintf (header,"\nVECTORS %dD_Velocity_Field float\n", DIMENSIONS);
332  else
333  sprintf (header,"\nVECTORS %dD_Magnetic_Field float\n", DIMENSIONS);
334 
335  VTK_HEADER_WRITE_STRING(header);
336  WriteBinaryArray (vect3D[0][0], sizeof(Float_Vect), SZ_Float_Vect, fvtk, -1);
337  }
338 }
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
FLOAT_VECT::v1
float v1
Definition: structs.h:314
FLOAT_VECT::v2
float v2
Definition: structs.h:314
VTK_HEADER_WRITE_STRING
#define VTK_HEADER_WRITE_STRING(header)
Definition: write_vtk.c:82
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
FLOAT_VECT::v3
float v3
Definition: structs.h:314
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
KDIR
#define KDIR
Definition: pluto.h:195
SZ_Float_Vect
int SZ_Float_Vect
Definition: globals.h:29
SWAP_VAR
#define SWAP_VAR(x)
Definition: macros.h:115
IDIR
#define IDIR
Definition: pluto.h:193
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
IsLittleEndian
int IsLittleEndian(void)
Definition: tools.c:40
k
int k
Definition: analysis.c:2
FLOAT_VECT
Definition: structs.h:313
GRID::x
double * x
Definition: structs.h:80
VectorCartesianComponents
void VectorCartesianComponents(double *v, double x1, double x2, double x3)
Definition: math_misc.c:227
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
i
int i
Definition: analysis.c:2
JDIR
#define JDIR
Definition: pluto.h:194
WriteBinaryArray
void WriteBinaryArray(void *V, size_t dsize, int sz, FILE *fl, int istag)
Definition: bin_io.c:98
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the call graph for this function:

Here is the caller graph for this function: