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