Actual source code: grvtk.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsc/private/dmdaimpl.h>
  2: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>

  6: static PetscErrorCode DMDAVTKWriteAll_VTS(DM da,PetscViewer viewer)
  7: {
  8: #if defined(PETSC_USE_REAL_SINGLE)
  9:   const char precision[] = "Float32";
 10: #elif defined(PETSC_USE_REAL_DOUBLE)
 11:   const char precision[] = "Float64";
 12: #else
 13:   const char precision[] = "UnknownPrecision";
 14: #endif
 15:   MPI_Comm                 comm;
 16:   Vec                      Coords;
 17:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*)viewer->data;
 18:   PetscViewerVTKObjectLink link;
 19:   FILE                     *fp;
 20:   PetscMPIInt              rank,size,tag;
 21:   DMDALocalInfo            info;
 22:   PetscInt                 dim,mx,my,mz,cdim,bs,boffset,maxnnodes,i,j,k,f,r;
 23:   PetscInt                 rloc[6],(*grloc)[6] = NULL;
 24:   PetscScalar              *array,*array2;
 25:   PetscReal                gmin[3],gmax[3];
 26:   PetscErrorCode           ierr;

 29:   PetscObjectGetComm((PetscObject)da,&comm);
 30: #if defined(PETSC_USE_COMPLEX)
 31:   SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported");
 32: #endif
 33:   MPI_Comm_size(comm,&size);
 34:   MPI_Comm_rank(comm,&rank);
 35:   DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);
 36:   DMDAGetLocalInfo(da,&info);
 37:   DMDAGetBoundingBox(da,gmin,gmax);
 38:   DMGetCoordinates(da,&Coords);
 39:   if (Coords) {
 40:     PetscInt csize;
 41:     VecGetSize(Coords,&csize);
 42:     if (csize % (mx*my*mz)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate vector size mismatch");
 43:     cdim = csize/(mx*my*mz);
 44:     if (cdim < dim || cdim > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate vector size mismatch");
 45:   } else {
 46:     cdim = dim;
 47:   }

 49:   PetscFOpen(comm,vtk->filename,"wb",&fp);
 50:   PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
 51: #if defined(PETSC_WORDS_BIGENDIAN)
 52:   PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
 53: #else
 54:   PetscFPrintf(comm,fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
 55: #endif
 56:   PetscFPrintf(comm,fp,"  <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);

 58:   if (!rank) {PetscMalloc1(size*6,&grloc);}
 59:   rloc[0] = info.xs;
 60:   rloc[1] = info.xm;
 61:   rloc[2] = info.ys;
 62:   rloc[3] = info.ym;
 63:   rloc[4] = info.zs;
 64:   rloc[5] = info.zm;
 65:   MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);

 67:   /* Write XML header */
 68:   maxnnodes = 0;                /* Used for the temporary array size on rank 0 */
 69:   boffset   = 0;                /* Offset into binary file */
 70:   for (r=0; r<size; r++) {
 71:     PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
 72:     if (!rank) {
 73:       xs     = grloc[r][0];
 74:       xm     = grloc[r][1];
 75:       ys     = grloc[r][2];
 76:       ym     = grloc[r][3];
 77:       zs     = grloc[r][4];
 78:       zm     = grloc[r][5];
 79:       nnodes = xm*ym*zm;
 80:     }
 81:     maxnnodes = PetscMax(maxnnodes,nnodes);
 82: #if 0
 83:     switch (dim) {
 84:     case 1:
 85:       PetscFPrintf(comm,fp,"    <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,0,0,0,0);
 86:       break;
 87:     case 2:
 88:       PetscFPrintf(comm,fp,"    <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm,ys+ym-1,xs,xs+xm-1,0,0);
 89:       break;
 90:     case 3:
 91:       PetscFPrintf(comm,fp,"    <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
 92:       break;
 93:     default: SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"No support for dimension %D",dim);
 94:     }
 95: #endif
 96:     PetscFPrintf(comm,fp,"    <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
 97:     PetscFPrintf(comm,fp,"      <Points>\n");
 98:     PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);
 99:     boffset += 3*nnodes*sizeof(PetscScalar) + sizeof(int);
100:     PetscFPrintf(comm,fp,"      </Points>\n");

102:     PetscFPrintf(comm,fp,"      <PointData Scalars=\"ScalarPointData\">\n");
103:     for (link=vtk->link; link; link=link->next) {
104:       Vec        X        = (Vec)link->vec;
105:       const char *vecname = "";
106:       if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
107:         PetscObjectGetName((PetscObject)X,&vecname);
108:       }
109:       for (i=0; i<bs; i++) {
110:         char       buf[256];
111:         const char *fieldname;
112:         DMDAGetFieldName(da,i,&fieldname);
113:         if (!fieldname) {
114:           PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);
115:           fieldname = buf;
116:         }
117:         PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
118:         boffset += nnodes*sizeof(PetscScalar) + sizeof(int);
119:       }
120:     }
121:     PetscFPrintf(comm,fp,"      </PointData>\n");
122:     PetscFPrintf(comm,fp,"    </Piece>\n");
123:   }
124:   PetscFPrintf(comm,fp,"  </StructuredGrid>\n");
125:   PetscFPrintf(comm,fp,"  <AppendedData encoding=\"raw\">\n");
126:   PetscFPrintf(comm,fp,"_");

128:   /* Now write the arrays. */
129:   tag  = ((PetscObject)viewer)->tag;
130:   PetscMalloc2(maxnnodes*PetscMax(3,bs),&array,maxnnodes*3,&array2);
131:   for (r=0; r<size; r++) {
132:     MPI_Status status;
133:     PetscInt   xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
134:     if (!rank) {
135:       xs     = grloc[r][0];
136:       xm     = grloc[r][1];
137:       ys     = grloc[r][2];
138:       ym     = grloc[r][3];
139:       zs     = grloc[r][4];
140:       zm     = grloc[r][5];
141:       nnodes = xm*ym*zm;
142:     } else if (r == rank) {
143:       nnodes = info.xm*info.ym*info.zm;
144:     }

146:     /* Write the coordinates */
147:     if (Coords) {
148:       const PetscScalar *coords;
149:       VecGetArrayRead(Coords,&coords);
150:       if (!rank) {
151:         if (r) {
152:           PetscMPIInt nn;
153:           MPI_Recv(array,nnodes*cdim,MPIU_SCALAR,r,tag,comm,&status);
154:           MPI_Get_count(&status,MPIU_SCALAR,&nn);
155:           if (nn != nnodes*cdim) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
156:         } else {
157:           PetscMemcpy(array,coords,nnodes*cdim*sizeof(PetscScalar));
158:         }
159:         /* Transpose coordinates to VTK (C-style) ordering */
160:         for (k=0; k<zm; k++) {
161:           for (j=0; j<ym; j++) {
162:             for (i=0; i<xm; i++) {
163:               PetscInt Iloc = i+xm*(j+ym*k);
164:               array2[Iloc*3+0] = array[Iloc*cdim + 0];
165:               array2[Iloc*3+1] = cdim > 1 ? array[Iloc*cdim + 1] : 0.0;
166:               array2[Iloc*3+2] = cdim > 2 ? array[Iloc*cdim + 2] : 0.0;
167:             }
168:           }
169:         }
170:       } else if (r == rank) {
171:         MPI_Send((void*)coords,nnodes*cdim,MPIU_SCALAR,0,tag,comm);
172:       }
173:       VecRestoreArrayRead(Coords,&coords);
174:     } else {       /* Fabricate some coordinates using grid index */
175:       for (k=0; k<zm; k++) {
176:         for (j=0; j<ym; j++) {
177:           for (i=0; i<xm; i++) {
178:             PetscInt Iloc = i+xm*(j+ym*k);
179:             array2[Iloc*3+0] = xs+i;
180:             array2[Iloc*3+1] = ys+j;
181:             array2[Iloc*3+2] = zs+k;
182:           }
183:         }
184:       }
185:     }
186:     PetscViewerVTKFWrite(viewer,fp,array2,nnodes*3,PETSC_SCALAR);

188:     /* Write each of the objects queued up for this file */
189:     for (link=vtk->link; link; link=link->next) {
190:       Vec               X = (Vec)link->vec;
191:       const PetscScalar *x;

193:       VecGetArrayRead(X,&x);
194:       if (!rank) {
195:         if (r) {
196:           PetscMPIInt nn;
197:           MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);
198:           MPI_Get_count(&status,MPIU_SCALAR,&nn);
199:           if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r);
200:         } else {
201:           PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));
202:         }
203:         for (f=0; f<bs; f++) {
204:           /* Extract and transpose the f'th field */
205:           for (k=0; k<zm; k++) {
206:             for (j=0; j<ym; j++) {
207:               for (i=0; i<xm; i++) {
208:                 PetscInt Iloc = i+xm*(j+ym*k);
209:                 array2[Iloc] = array[Iloc*bs + f];
210:               }
211:             }
212:           }
213:           PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);
214:         }
215:       } else if (r == rank) {
216:         MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);
217:       }
218:       VecRestoreArrayRead(X,&x);
219:     }
220:   }
221:   PetscFree2(array,array2);
222:   PetscFree(grloc);

224:   PetscFPrintf(comm,fp,"\n </AppendedData>\n");
225:   PetscFPrintf(comm,fp,"</VTKFile>\n");
226:   PetscFClose(comm,fp);
227:   return(0);
228: }


233: static PetscErrorCode DMDAVTKWriteAll_VTR(DM da,PetscViewer viewer)
234: {
235: #if defined(PETSC_USE_REAL_SINGLE)
236:   const char precision[] = "Float32";
237: #elif defined(PETSC_USE_REAL_DOUBLE)
238:   const char precision[] = "Float64";
239: #else
240:   const char precision[] = "UnknownPrecision";
241: #endif
242:   MPI_Comm                 comm;
243:   PetscViewer_VTK          *vtk = (PetscViewer_VTK*)viewer->data;
244:   PetscViewerVTKObjectLink link;
245:   FILE                     *fp;
246:   PetscMPIInt              rank,size,tag;
247:   DMDALocalInfo            info;
248:   PetscInt                 dim,mx,my,mz,bs,boffset,maxnnodes,i,j,k,f,r;
249:   PetscInt                 rloc[6],(*grloc)[6] = NULL;
250:   PetscScalar              *array,*array2;
251:   PetscReal                gmin[3],gmax[3];
252:   PetscErrorCode           ierr;

255:   PetscObjectGetComm((PetscObject)da,&comm);
256: #if defined(PETSC_USE_COMPLEX)
257:   SETERRQ(comm,PETSC_ERR_SUP,"Complex values not supported");
258: #endif
259:   MPI_Comm_size(comm,&size);
260:   MPI_Comm_rank(comm,&rank);
261:   DMDAGetInfo(da,&dim, &mx,&my,&mz, 0,0,0, &bs,0,0,0,0,0);
262:   DMDAGetLocalInfo(da,&info);
263:   DMDAGetBoundingBox(da,gmin,gmax);
264:   PetscFOpen(comm,vtk->filename,"wb",&fp);
265:   PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");
266: #if defined(PETSC_WORDS_BIGENDIAN)
267:   PetscFPrintf(comm,fp,"<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
268: #else
269:   PetscFPrintf(comm,fp,"<VTKFile type=\"RectilinearGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
270: #endif
271:   PetscFPrintf(comm,fp,"  <RectilinearGrid WholeExtent=\"%D %D %D %D %D %D\">\n",0,mx-1,0,my-1,0,mz-1);

273:   if (!rank) {PetscMalloc1(size*6,&grloc);}
274:   rloc[0] = info.xs;
275:   rloc[1] = info.xm;
276:   rloc[2] = info.ys;
277:   rloc[3] = info.ym;
278:   rloc[4] = info.zs;
279:   rloc[5] = info.zm;
280:   MPI_Gather(rloc,6,MPIU_INT,&grloc[0][0],6,MPIU_INT,0,comm);

282:   /* Write XML header */
283:   maxnnodes = 0;                /* Used for the temporary array size on rank 0 */
284:   boffset   = 0;                /* Offset into binary file */
285:   for (r=0; r<size; r++) {
286:     PetscInt xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
287:     if (!rank) {
288:       xs     = grloc[r][0];
289:       xm     = grloc[r][1];
290:       ys     = grloc[r][2];
291:       ym     = grloc[r][3];
292:       zs     = grloc[r][4];
293:       zm     = grloc[r][5];
294:       nnodes = xm*ym*zm;
295:     }
296:     maxnnodes = PetscMax(maxnnodes,nnodes);
297:     PetscFPrintf(comm,fp,"    <Piece Extent=\"%D %D %D %D %D %D\">\n",xs,xs+xm-1,ys,ys+ym-1,zs,zs+zm-1);
298:     PetscFPrintf(comm,fp,"      <Coordinates>\n");
299:     PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"Xcoord\"  format=\"appended\"  offset=\"%D\" />\n",precision,boffset);
300:     boffset += xm*sizeof(PetscScalar) + sizeof(int);
301:     PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"Ycoord\"  format=\"appended\"  offset=\"%D\" />\n",precision,boffset);
302:     boffset += ym*sizeof(PetscScalar) + sizeof(int);
303:     PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"Zcoord\"  format=\"appended\"  offset=\"%D\" />\n",precision,boffset);
304:     boffset += zm*sizeof(PetscScalar) + sizeof(int);
305:     PetscFPrintf(comm,fp,"      </Coordinates>\n");
306:     PetscFPrintf(comm,fp,"      <PointData Scalars=\"ScalarPointData\">\n");
307:     for (link=vtk->link; link; link=link->next) {
308:       Vec        X        = (Vec)link->vec;
309:       const char *vecname = "";
310:       if (((PetscObject)X)->name || link != vtk->link) { /* If the object is already named, use it. If it is past the first link, name it to disambiguate. */
311:         PetscObjectGetName((PetscObject)X,&vecname);
312:       }
313:       for (i=0; i<bs; i++) {
314:         char       buf[256];
315:         const char *fieldname;
316:         DMDAGetFieldName(da,i,&fieldname);
317:         if (!fieldname) {
318:           PetscSNPrintf(buf,sizeof(buf),"Unnamed%D",i);
319:           fieldname = buf;
320:         }
321:         PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);
322:         boffset += nnodes*sizeof(PetscScalar) + sizeof(int);
323:       }
324:     }
325:     PetscFPrintf(comm,fp,"      </PointData>\n");
326:     PetscFPrintf(comm,fp,"    </Piece>\n");
327:   }
328:   PetscFPrintf(comm,fp,"  </RectilinearGrid>\n");
329:   PetscFPrintf(comm,fp,"  <AppendedData encoding=\"raw\">\n");
330:   PetscFPrintf(comm,fp,"_");

332:   /* Now write the arrays. */
333:   tag  = ((PetscObject)viewer)->tag;
334:   PetscMalloc2(PetscMax(maxnnodes*bs,info.xm+info.ym+info.zm),&array,PetscMax(maxnnodes*3,info.xm+info.ym+info.zm),&array2);
335:   for (r=0; r<size; r++) {
336:     MPI_Status status;
337:     PetscInt   xs=-1,xm=-1,ys=-1,ym=-1,zs=-1,zm=-1,nnodes = 0;
338:     if (!rank) {
339:       xs     = grloc[r][0];
340:       xm     = grloc[r][1];
341:       ys     = grloc[r][2];
342:       ym     = grloc[r][3];
343:       zs     = grloc[r][4];
344:       zm     = grloc[r][5];
345:       nnodes = xm*ym*zm;
346:     } else if (r == rank) {
347:       nnodes = info.xm*info.ym*info.zm;
348:     }
349:     {                           /* Write the coordinates */
350:       Vec Coords;
351:       DMGetCoordinates(da,&Coords);
352:       if (Coords) {
353:         const PetscScalar *coords;
354:         VecGetArrayRead(Coords,&coords);
355:         if (!rank) {
356:           if (r) {
357:             PetscMPIInt nn;
358:             MPI_Recv(array,xm+ym+zm,MPIU_SCALAR,r,tag,comm,&status);
359:             MPI_Get_count(&status,MPIU_SCALAR,&nn);
360:             if (nn != xm+ym+zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
361:           } else {
362:             /* x coordinates */
363:             for (j=0, k=0, i=0; i<xm; i++) {
364:               PetscInt Iloc = i+xm*(j+ym*k);
365:               array[i] = coords[Iloc*dim + 0];
366:             }
367:             /* y coordinates */
368:             for (i=0, k=0, j=0; j<ym; j++) {
369:               PetscInt Iloc = i+xm*(j+ym*k);
370:               array[j+xm] = dim > 1 ? coords[Iloc*dim + 1] : 0;
371:             }
372:             /* z coordinates */
373:             for (i=0, j=0, k=0; k<zm; k++) {
374:               PetscInt Iloc = i+xm*(j+ym*k);
375:               array[k+xm+ym] = dim > 2 ? coords[Iloc*dim + 2] : 0;
376:             }
377:           }
378:         } else if (r == rank) {
379:           xm = info.xm;
380:           ym = info.ym;
381:           zm = info.zm;
382:           for (j=0, k=0, i=0; i<xm; i++) {
383:             PetscInt Iloc = i+xm*(j+ym*k);
384:             array2[i] = coords[Iloc*dim + 0];
385:           }
386:           for (i=0, k=0, j=0; j<ym; j++) {
387:             PetscInt Iloc = i+xm*(j+ym*k);
388:             array2[j+xm] = dim > 1 ? coords[Iloc*dim + 1] : 0;
389:           }
390:           for (i=0, j=0, k=0; k<zm; k++) {
391:             PetscInt Iloc = i+xm*(j+ym*k);
392:             array2[k+xm+ym] = dim > 2 ? coords[Iloc*dim + 2] : 0;
393:           }
394:           MPI_Send((void*)array2,xm+ym+zm,MPIU_SCALAR,0,tag,comm);
395:         }
396:         VecRestoreArrayRead(Coords,&coords);
397:       } else {       /* Fabricate some coordinates using grid index */
398:         for (j=0, k=0, i=0; i<xm; i++) {
399:           array[i] = xs+i;
400:         }
401:         for (i=0, k=0, j=0; j<ym; j++) {
402:           array[j+xm] = ys+j;
403:         }
404:         for (i=0, j=0, k=0; k<zm; k++) {
405:           array[k+xm+ym] = zs+k;
406:         }
407:       }
408:       if (!rank) {
409:         PetscViewerVTKFWrite(viewer,fp,&(array[0]),xm,PETSC_SCALAR);
410:         PetscViewerVTKFWrite(viewer,fp,&(array[xm]),ym,PETSC_SCALAR);
411:         PetscViewerVTKFWrite(viewer,fp,&(array[xm+ym]),zm,PETSC_SCALAR);
412:       }
413:     }

415:     /* Write each of the objects queued up for this file */
416:     for (link=vtk->link; link; link=link->next) {
417:       Vec               X = (Vec)link->vec;
418:       const PetscScalar *x;

420:       VecGetArrayRead(X,&x);
421:       if (!rank) {
422:         if (r) {
423:           PetscMPIInt nn;
424:           MPI_Recv(array,nnodes*bs,MPIU_SCALAR,r,tag,comm,&status);
425:           MPI_Get_count(&status,MPIU_SCALAR,&nn);
426:           if (nn != nnodes*bs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch receiving from rank %D",r);
427:         } else {
428:           PetscMemcpy(array,x,nnodes*bs*sizeof(PetscScalar));
429:         }
430:         for (f=0; f<bs; f++) {
431:           /* Extract and transpose the f'th field */
432:           for (k=0; k<zm; k++) {
433:             for (j=0; j<ym; j++) {
434:               for (i=0; i<xm; i++) {
435:                 PetscInt Iloc = i+xm*(j+ym*k);
436:                 array2[Iloc] = array[Iloc*bs + f];
437:               }
438:             }
439:           }
440:           PetscViewerVTKFWrite(viewer,fp,array2,nnodes,PETSC_SCALAR);
441:         }
442:       } else if (r == rank) {
443:         MPI_Send((void*)x,nnodes*bs,MPIU_SCALAR,0,tag,comm);
444:       }
445:       VecRestoreArrayRead(X,&x);
446:     }
447:   }
448:   PetscFree2(array,array2);
449:   PetscFree(grloc);

451:   PetscFPrintf(comm,fp,"\n </AppendedData>\n");
452:   PetscFPrintf(comm,fp,"</VTKFile>\n");
453:   PetscFClose(comm,fp);
454:   return(0);
455: }

459: /*@C
460:    DMDAVTKWriteAll - Write a file containing all the fields that have been provided to the viewer

462:    Collective

464:    Input Arguments:
465:    odm - DM specifying the grid layout, passed as a PetscObject
466:    viewer - viewer of type VTK

468:    Level: developer

470:    Note:
471:    This function is a callback used by the VTK viewer to actually write the file.
472:    The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
473:    Instead, metadata for the entire file needs to be available up-front before you can start writing the file.

475: .seealso: PETSCVIEWERVTK
476: @*/
477: PetscErrorCode DMDAVTKWriteAll(PetscObject odm,PetscViewer viewer)
478: {
479:   DM             dm = (DM)odm;
480:   PetscBool      isvtk;

486:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);
487:   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_INCOMP,"Cannot use viewer type %s",((PetscObject)viewer)->type_name);
488:   switch (viewer->format) {
489:   case PETSC_VIEWER_VTK_VTS:
490:     DMDAVTKWriteAll_VTS(dm,viewer);
491:     break;
492:   case PETSC_VIEWER_VTK_VTR:
493:     DMDAVTKWriteAll_VTR(dm,viewer);
494:     break;
495:   default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for format '%s'",PetscViewerFormats[viewer->format]);
496:   }
497:   return(0);
498: }