xref: /petsc/src/dm/impls/plex/plexvtu.c (revision f37c82fbca2ef2d9869799d3305350a5f272bd86)
1 #include <petsc-private/dmpleximpl.h>
2 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
3 
4 typedef struct {
5   PetscInt nvertices;
6   PetscInt ncells;
7   PetscInt nconn;               /* number of entries in cell->vertex connectivity array */
8 } PieceInfo;
9 
10 #if defined(PETSC_USE_REAL_SINGLE)
11 static const char precision[] = "Float32";
12 #elif defined(PETSC_USE_REAL_DOUBLE)
13 static const char precision[] = "Float64";
14 #else
15 static const char precision[] = "UnknownPrecision";
16 #endif
17 
18 #undef __FUNCT__
19 #define __FUNCT__ "TransferWrite"
20 static PetscErrorCode TransferWrite(PetscViewer viewer,FILE *fp,PetscMPIInt srank,PetscMPIInt root,const void *send,void *recv,PetscMPIInt count,PetscDataType datatype,PetscMPIInt tag)
21 {
22   PetscMPIInt    rank;
23   PetscErrorCode ierr;
24   MPI_Comm       comm;
25   MPI_Datatype   mpidatatype;
26 
27   PetscFunctionBegin;
28   ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
29   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
30   ierr = PetscDataTypeToMPIDataType(datatype,&mpidatatype);CHKERRQ(ierr);
31 
32   if (rank == srank && rank != root) {
33     ierr = MPI_Send((void*)send,count,mpidatatype,root,tag,comm);CHKERRQ(ierr);
34   } else if (rank == root) {
35     const void *buffer;
36     if (root == srank) {        /* self */
37       buffer = send;
38     } else {
39       MPI_Status  status;
40       PetscMPIInt nrecv;
41       ierr = MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);CHKERRQ(ierr);
42       ierr = MPI_Get_count(&status,mpidatatype,&nrecv);CHKERRQ(ierr);
43       if (count != nrecv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
44       buffer = recv;
45     }
46     ierr = PetscViewerVTKFWrite(viewer,fp,buffer,count,datatype);CHKERRQ(ierr);
47   }
48   PetscFunctionReturn(0);
49 }
50 
51 #undef __FUNCT__
52 #define __FUNCT__ "DMPlexGetVTKConnectivity"
53 static PetscErrorCode DMPlexGetVTKConnectivity(DM dm,PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes)
54 {
55   PetscErrorCode ierr;
56   PetscVTKInt    *conn,*offsets;
57   PetscVTKType   *types;
58   PetscInt       dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,cMax,numLabelCells,hasLabel,c,v,countcell,countconn;
59 
60   PetscFunctionBegin;
61   ierr = PetscMalloc3(piece->nconn,PetscVTKInt,&conn,piece->ncells,PetscVTKInt,&offsets,piece->ncells,PetscVTKType,&types);CHKERRQ(ierr);
62 
63   ierr = DMPlexGetDimension(dm,&dim);CHKERRQ(ierr);
64   ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
65   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
66   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
67   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
68   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
69   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
70   ierr     = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
71   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
72 
73   countcell = 0;
74   countconn = 0;
75   for (c = cStart; c < cEnd; ++c) {
76     PetscInt *closure = NULL;
77     PetscInt  closureSize,nverts,celltype,startoffset;
78 
79     if (hasLabel) {
80       PetscInt value;
81 
82       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
83       if (value != 1) continue;
84     }
85     startoffset = countconn;
86     ierr        = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
87     for (v = 0; v < closureSize*2; v += 2) {
88       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
89         conn[countconn++] = closure[v] - vStart;
90       }
91     }
92     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
93 
94     offsets[countcell] = countconn;
95 
96     nverts = countconn - startoffset;
97     ierr   = DMPlexVTKGetCellType(dm,dim,nverts,&celltype);CHKERRQ(ierr);
98 
99     types[countcell] = celltype;
100     countcell++;
101   }
102   if (countcell != piece->ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count");
103   if (countconn != piece->nconn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count");
104   *oconn    = conn;
105   *ooffsets = offsets;
106   *otypes   = types;
107   PetscFunctionReturn(0);
108 }
109 
110 #undef __FUNCT__
111 #define __FUNCT__ "DMPlexVTKWriteAll_VTU"
112 /*
113   Write all fields that have been provided to the viewer
114   Multi-block XML format with binary appended data.
115 */
116 PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)
117 {
118   MPI_Comm                 comm;
119   PetscViewer_VTK          *vtk = (PetscViewer_VTK*)viewer->data;
120   PetscViewerVTKObjectLink link;
121   FILE                     *fp;
122   PetscMPIInt              rank,size,tag;
123   PetscErrorCode           ierr;
124   PetscInt                 dim,cellHeight,cStart,cEnd,vStart,vEnd,cMax,numLabelCells,hasLabel,c,v,r,i;
125   PieceInfo                piece,*gpiece = NULL;
126   void                     *buffer = NULL;
127 
128   PetscFunctionBegin;
129   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
130 #if defined(PETSC_USE_COMPLEX)
131   SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Complex values not supported");
132 #endif
133   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
134   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
135   tag  = ((PetscObject)viewer)->tag;
136 
137   ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr);
138   ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr);
139 #if defined(PETSC_WORDS_BIGENDIAN)
140   ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");CHKERRQ(ierr);
141 #else
142   ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");CHKERRQ(ierr);
143 #endif
144   ierr = PetscFPrintf(comm,fp,"  <UnstructuredGrid>\n");CHKERRQ(ierr);
145 
146   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
147   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
148   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
149   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
150   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
151   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
152   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
153 
154   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
155   piece.nvertices = vEnd - vStart;
156   piece.ncells    = 0;
157   piece.nconn     = 0;
158   for (c = cStart; c < cEnd; ++c) {
159     PetscInt *closure = NULL;
160     PetscInt closureSize;
161 
162     if (hasLabel) {
163       PetscInt value;
164 
165       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
166       if (value != 1) continue;
167     }
168     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
169     for (v = 0; v < closureSize*2; v += 2) {
170       if ((closure[v] >= vStart) && (closure[v] < vEnd)) piece.nconn++;
171     }
172     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
173     piece.ncells++;
174   }
175   if (!rank) {ierr = PetscMalloc(size*sizeof(piece),&gpiece);CHKERRQ(ierr);}
176   ierr = MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);CHKERRQ(ierr);
177 
178   /*
179    * Write file header
180    */
181   if (!rank) {
182     PetscInt boffset = 0;
183 
184     for (r=0; r<size; r++) {
185       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"    <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);CHKERRQ(ierr);
186       /* Coordinate positions */
187       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      <Points>\n");CHKERRQ(ierr);
188       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr);
189       boffset += gpiece[r].nvertices*3*sizeof(PetscScalar) + sizeof(int);
190       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      </Points>\n");CHKERRQ(ierr);
191       /* Cell connectivity */
192       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      <Cells>\n");CHKERRQ(ierr);
193       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
194       boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int);
195       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"offsets\"      NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
196       boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int);
197       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"UInt8\" Name=\"types\"        NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
198       boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int);
199       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      </Cells>\n");CHKERRQ(ierr);
200 
201       /*
202        * Cell Data headers
203        */
204       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      <CellData>\n");CHKERRQ(ierr);
205       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
206       boffset += gpiece[r].ncells*sizeof(int) + sizeof(int);
207       /* all the vectors */
208       for (link=vtk->link; link; link=link->next) {
209         Vec        X = (Vec)link->vec;
210         PetscInt   bs,nfields,field;
211         const char *vecname = "";
212         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
213         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. */
214           ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr);
215         }
216         ierr = PetscSectionGetDof(dm->defaultSection,cStart,&bs);CHKERRQ(ierr);
217         ierr = PetscSectionGetNumFields(dm->defaultSection,&nfields);CHKERRQ(ierr);
218         for (field=0,i=0; field<(nfields?nfields:1); field++) {
219           PetscInt   fbs,j;
220           const char *fieldname = NULL;
221           char       buf[256];
222           if (nfields) {        /* We have user-defined fields/components */
223             ierr = PetscSectionGetFieldDof(dm->defaultSection,cStart,field,&fbs);CHKERRQ(ierr);
224             ierr = PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);CHKERRQ(ierr);
225           } else fbs = bs;      /* Say we have one field with 'bs' components */
226           if (!fieldname) {
227             ierr = PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);CHKERRQ(ierr);
228             fieldname = buf;
229           }
230           for (j=0; j<fbs; j++) {
231             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);CHKERRQ(ierr);
232             boffset += gpiece[r].ncells*sizeof(PetscScalar) + sizeof(int);
233             i++;
234           }
235         }
236         if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs);
237       }
238       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"      </CellData>\n");CHKERRQ(ierr);
239 
240       /*
241        * Point Data headers
242        */
243       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"      <PointData>\n");CHKERRQ(ierr);
244       for (link=vtk->link; link; link=link->next) {
245         Vec        X = (Vec)link->vec;
246         PetscInt   bs,nfields,field;
247         const char *vecname = "";
248         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
249         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. */
250           ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr);
251         }
252         ierr = PetscSectionGetDof(dm->defaultSection,vStart,&bs);CHKERRQ(ierr);
253         ierr = PetscSectionGetNumFields(dm->defaultSection,&nfields);CHKERRQ(ierr);
254         for (field=0,i=0; field<(nfields?nfields:1); field++) {
255           PetscInt   fbs,j;
256           const char *fieldname = NULL;
257           char       buf[256];
258           if (nfields) {        /* We have user-defined fields/components */
259             ierr = PetscSectionGetFieldDof(dm->defaultSection,vStart,field,&fbs);CHKERRQ(ierr);
260             ierr = PetscSectionGetFieldName(dm->defaultSection,field,&fieldname);CHKERRQ(ierr);
261           } else fbs = bs;      /* Say we have one field with 'bs' components */
262           if (!fieldname) {
263             ierr = PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);CHKERRQ(ierr);
264             fieldname = buf;
265           }
266           for (j=0; j<fbs; j++) {
267             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.%D\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,j,boffset);CHKERRQ(ierr);
268             boffset += gpiece[r].nvertices*sizeof(PetscScalar) + sizeof(int);
269           }
270         }
271       }
272       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"      </PointData>\n");CHKERRQ(ierr);
273 
274       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"    </Piece>\n");CHKERRQ(ierr);
275     }
276   }
277 
278   ierr = PetscFPrintf(comm,fp,"  </UnstructuredGrid>\n");CHKERRQ(ierr);
279   ierr = PetscFPrintf(comm,fp,"  <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr);
280   ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr);
281 
282   if (!rank) {
283     PetscInt maxsize = 0;
284     for (r=0; r<size; r++) {
285       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscScalar)));
286       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*sizeof(PetscScalar)));
287       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
288     }
289     ierr = PetscMalloc(maxsize,&buffer);CHKERRQ(ierr);
290   }
291   for (r=0; r<size; r++) {
292     if (r == rank) {
293       PetscInt nsend;
294       {                         /* Position */
295         const PetscScalar *x;
296         PetscScalar       *y = NULL;
297         Vec               coords;
298         nsend = piece.nvertices*3;
299         ierr  = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr);
300         ierr  = VecGetArrayRead(coords,&x);CHKERRQ(ierr);
301         if (dim != 3) {
302           ierr = PetscMalloc(piece.nvertices*3*sizeof(PetscScalar),&y);CHKERRQ(ierr);
303           for (i=0; i<piece.nvertices; i++) {
304             y[i*3+0] = x[i*dim+0];
305             y[i*3+1] = (dim > 1) ? x[i*dim+1] : 0;
306             y[i*3+2] = 0;
307           }
308         }
309         ierr = TransferWrite(viewer,fp,r,0,y ? y : x,buffer,nsend,PETSC_SCALAR,tag);CHKERRQ(ierr);
310         ierr = PetscFree(y);CHKERRQ(ierr);
311         ierr = VecRestoreArrayRead(coords,&x);CHKERRQ(ierr);
312       }
313       {                           /* Connectivity, offsets, types */
314         PetscVTKInt  *connectivity = NULL,*offsets;
315         PetscVTKType *types;
316         ierr = DMPlexGetVTKConnectivity(dm,&piece,&connectivity,&offsets,&types);CHKERRQ(ierr);
317         ierr = TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,PETSC_INT32,tag);CHKERRQ(ierr);
318         ierr = TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,PETSC_INT32,tag);CHKERRQ(ierr);
319         ierr = TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,PETSC_UINT8,tag);CHKERRQ(ierr);
320         ierr = PetscFree3(connectivity,offsets,types);CHKERRQ(ierr);
321       }
322       {                         /* Owners (cell data) */
323         PetscVTKInt *owners;
324         ierr = PetscMalloc(piece.ncells*sizeof(PetscVTKInt),&owners);CHKERRQ(ierr);
325         for (i=0; i<piece.ncells; i++) owners[i] = rank;
326         ierr = TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,PETSC_INT32,tag);CHKERRQ(ierr);
327         ierr = PetscFree(owners);CHKERRQ(ierr);
328       }
329       /* Cell data */
330       for (link=vtk->link; link; link=link->next) {
331         Vec               X = (Vec)link->vec;
332         const PetscScalar *x;
333         PetscScalar       *y;
334         PetscInt          bs;
335         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
336         ierr = PetscSectionGetDof(dm->defaultSection,cStart,&bs);CHKERRQ(ierr);
337         ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
338         ierr = PetscMalloc(piece.ncells*sizeof(PetscScalar),&y);CHKERRQ(ierr);
339         for (i=0; i<bs; i++) {
340           PetscInt cnt;
341           for (c=cStart,cnt=0; c<cEnd; c++) {
342             const PetscScalar *xpoint;
343             if (hasLabel) {     /* Ignore some cells */
344               PetscInt value;
345               ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
346               if (value != 1) continue;
347             }
348             ierr     = DMPlexPointLocalRead(dm,c,x,&xpoint);CHKERRQ(ierr);
349             y[cnt++] = xpoint[i];
350           }
351           if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
352           ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,PETSC_SCALAR,tag);CHKERRQ(ierr);
353         }
354         ierr = PetscFree(y);CHKERRQ(ierr);
355         ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
356       }
357 
358       for (link=vtk->link; link; link=link->next) {
359         Vec               X = (Vec)link->vec;
360         const PetscScalar *x;
361         PetscScalar       *y;
362         PetscInt          bs;
363         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
364         ierr = PetscSectionGetDof(dm->defaultSection,vStart,&bs);CHKERRQ(ierr);
365         ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
366         ierr = PetscMalloc(piece.nvertices*sizeof(PetscScalar),&y);CHKERRQ(ierr);
367         for (i=0; i<bs; i++) {
368           PetscInt cnt;
369           for (v=vStart,cnt=0; v<vEnd; v++) {
370             const PetscScalar *xpoint;
371             ierr     = DMPlexPointLocalRead(dm,v,x,&xpoint);CHKERRQ(ierr);
372             y[cnt++] = xpoint[i];
373           }
374           if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
375           ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,PETSC_SCALAR,tag);CHKERRQ(ierr);
376         }
377         ierr = PetscFree(y);CHKERRQ(ierr);
378         ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
379       }
380     } else if (!rank) {
381       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,PETSC_SCALAR,tag);CHKERRQ(ierr); /* positions */
382       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,PETSC_INT32,tag);CHKERRQ(ierr); /* connectivity */
383       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag);CHKERRQ(ierr); /* offsets */
384       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_UINT8,tag);CHKERRQ(ierr); /* types */
385       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_INT32,tag);CHKERRQ(ierr); /* owner rank (cells) */
386       /* all cell data */
387       for (link=vtk->link; link; link=link->next) {
388         PetscInt bs;
389         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
390         ierr = PetscSectionGetDof(dm->defaultSection,cStart,&bs);CHKERRQ(ierr);
391         for (i=0; i<bs; i++) {
392           ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,PETSC_SCALAR,tag);CHKERRQ(ierr);
393         }
394       }
395       /* all point data */
396       for (link=vtk->link; link; link=link->next) {
397         PetscInt bs;
398         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
399         ierr = PetscSectionGetDof(dm->defaultSection,vStart,&bs);CHKERRQ(ierr);
400         for (i=0; i<bs; i++) {
401           ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,PETSC_SCALAR,tag);CHKERRQ(ierr);
402         }
403       }
404     }
405   }
406   ierr = PetscFree(gpiece);CHKERRQ(ierr);
407   ierr = PetscFree(buffer);CHKERRQ(ierr);
408   ierr = PetscFPrintf(comm,fp,"\n  </AppendedData>\n");CHKERRQ(ierr);
409   ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr);
410   PetscFunctionReturn(0);
411 }
412