xref: /petsc/src/dm/impls/plex/plexvtu.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94)
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) || defined(PETSC_USE_REAL___FP16)
11 /* output in float if single or half precision in memory */
12 static const char precision[] = "Float32";
13 typedef float PetscVTUReal;
14 #define MPIU_VTUREAL MPI_FLOAT
15 #elif defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
16 /* output in double if double or quad precision in memory */
17 static const char precision[] = "Float64";
18 typedef double PetscVTUReal;
19 #define MPIU_VTUREAL MPI_DOUBLE
20 #else
21 static const char precision[] = "UnknownPrecision";
22 typedef PetscReal PetscVTUReal;
23 #define MPIU_VTUREAL MPIU_REAL
24 #endif
25 
26 static PetscErrorCode TransferWrite(PetscViewer viewer,FILE *fp,PetscMPIInt srank,PetscMPIInt root,const void *send,void *recv,PetscMPIInt count,MPI_Datatype mpidatatype,PetscMPIInt tag)
27 {
28   PetscMPIInt    rank;
29   PetscErrorCode ierr;
30   MPI_Comm       comm;
31 
32   PetscFunctionBegin;
33   ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
34   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
35 
36   if (rank == srank && rank != root) {
37     ierr = MPI_Send((void*)send,count,mpidatatype,root,tag,comm);CHKERRQ(ierr);
38   } else if (rank == root) {
39     const void *buffer;
40     if (root == srank) {        /* self */
41       buffer = send;
42     } else {
43       MPI_Status  status;
44       PetscMPIInt nrecv;
45       ierr = MPI_Recv(recv,count,mpidatatype,srank,tag,comm,&status);CHKERRQ(ierr);
46       ierr = MPI_Get_count(&status,mpidatatype,&nrecv);CHKERRQ(ierr);
47       if (count != nrecv) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Array size mismatch");
48       buffer = recv;
49     }
50     ierr = PetscViewerVTKFWrite(viewer,fp,buffer,count,mpidatatype);CHKERRQ(ierr);
51   }
52   PetscFunctionReturn(0);
53 }
54 
55 static PetscErrorCode DMPlexGetVTKConnectivity(DM dm, PetscBool localized, PieceInfo *piece,PetscVTKInt **oconn,PetscVTKInt **ooffsets,PetscVTKType **otypes)
56 {
57   PetscErrorCode ierr;
58   PetscSection   coordSection;
59   PetscVTKInt    *conn,*offsets;
60   PetscVTKType   *types;
61   PetscInt       dim,vStart,vEnd,cStart,cEnd,pStart,pEnd,cellHeight,numLabelCells,hasLabel,c,v,countcell,countconn;
62 
63   PetscFunctionBegin;
64   ierr = PetscMalloc3(piece->nconn,&conn,piece->ncells,&offsets,piece->ncells,&types);CHKERRQ(ierr);
65   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
66   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
67   ierr = DMPlexGetChart(dm,&pStart,&pEnd);CHKERRQ(ierr);
68   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
69   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
70   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
71   ierr     = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
72   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
73 
74   countcell = 0;
75   countconn = 0;
76   for (c = cStart; c < cEnd; ++c) {
77     PetscInt nverts,dof = 0,celltype,startoffset,nC=0;
78 
79     if (hasLabel) {
80       PetscInt value;
81 
82       ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
83       if (value != 1) continue;
84     }
85     startoffset = countconn;
86     if (localized) {
87       ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr);
88     }
89     if (!dof) {
90       PetscInt *closure = NULL;
91       PetscInt  closureSize;
92 
93       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
94       for (v = 0; v < closureSize*2; v += 2) {
95         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
96           if (!localized) conn[countconn++] = closure[v] - vStart;
97           else conn[countconn++] = startoffset + nC;
98           ++nC;
99         }
100       }
101       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
102     } else {
103       for (nC = 0; nC < dof/dim; nC++) conn[countconn++] = startoffset + nC;
104     }
105     ierr = DMPlexInvertCell(dim, nC, &conn[countconn-nC]);CHKERRQ(ierr);
106 
107     offsets[countcell] = countconn;
108 
109     nverts = countconn - startoffset;
110     ierr   = DMPlexVTKGetCellType_Internal(dm,dim,nverts,&celltype);CHKERRQ(ierr);
111 
112     types[countcell] = celltype;
113     countcell++;
114   }
115   if (countcell != piece->ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent cell count");
116   if (countconn != piece->nconn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent connectivity count");
117   *oconn    = conn;
118   *ooffsets = offsets;
119   *otypes   = types;
120   PetscFunctionReturn(0);
121 }
122 
123 /*
124   Write all fields that have been provided to the viewer
125   Multi-block XML format with binary appended data.
126 */
127 PetscErrorCode DMPlexVTKWriteAll_VTU(DM dm,PetscViewer viewer)
128 {
129   MPI_Comm                 comm;
130   PetscSection             coordSection;
131   PetscViewer_VTK          *vtk = (PetscViewer_VTK*)viewer->data;
132   PetscViewerVTKObjectLink link;
133   FILE                     *fp;
134   PetscMPIInt              rank,size,tag;
135   PetscErrorCode           ierr;
136   PetscInt                 dimEmbed,cellHeight,cStart,cEnd,vStart,vEnd,numLabelCells,hasLabel,c,v,r,i;
137   PetscBool                localized;
138   PieceInfo                piece,*gpiece = NULL;
139   void                     *buffer = NULL;
140   const char               *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
141   PetscInt                 loops_per_scalar;
142 
143   PetscFunctionBegin;
144   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
145 #if defined(PETSC_USE_COMPLEX)
146   loops_per_scalar = 2;
147 #else
148   loops_per_scalar = 1;
149 #endif
150   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
151   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
152   tag  = ((PetscObject)viewer)->tag;
153 
154   ierr = PetscFOpen(comm,vtk->filename,"wb",&fp);CHKERRQ(ierr);
155   ierr = PetscFPrintf(comm,fp,"<?xml version=\"1.0\"?>\n");CHKERRQ(ierr);
156   ierr = PetscFPrintf(comm,fp,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order);CHKERRQ(ierr);
157   ierr = PetscFPrintf(comm,fp,"  <UnstructuredGrid>\n");CHKERRQ(ierr);
158 
159   ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr);
160   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
161   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
162   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
163   ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
164   ierr = DMGetCoordinatesLocalized(dm, &localized);CHKERRQ(ierr);
165   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
166 
167   hasLabel        = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
168   piece.nvertices = 0;
169   piece.ncells    = 0;
170   piece.nconn     = 0;
171   if (!localized) piece.nvertices = vEnd - vStart;
172   for (c = cStart; c < cEnd; ++c) {
173     PetscInt dof = 0;
174     if (hasLabel) {
175       PetscInt value;
176 
177       ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
178       if (value != 1) continue;
179     }
180     if (localized) {
181       ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr);
182     }
183     if (!dof) {
184       PetscInt *closure = NULL;
185       PetscInt closureSize;
186 
187       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
188       for (v = 0; v < closureSize*2; v += 2) {
189         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
190           piece.nconn++;
191           if (localized) piece.nvertices++;
192         }
193       }
194       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
195     } else {
196       piece.nvertices += dof/dimEmbed;
197       piece.nconn     += dof/dimEmbed;
198     }
199     piece.ncells++;
200   }
201   if (!rank) {ierr = PetscMalloc1(size,&gpiece);CHKERRQ(ierr);}
202   ierr = MPI_Gather((PetscInt*)&piece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,(PetscInt*)gpiece,sizeof(piece)/sizeof(PetscInt),MPIU_INT,0,comm);CHKERRQ(ierr);
203 
204   /*
205    * Write file header
206    */
207   if (!rank) {
208     PetscInt boffset = 0;
209 
210     for (r=0; r<size; r++) {
211       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"    <Piece NumberOfPoints=\"%D\" NumberOfCells=\"%D\">\n",gpiece[r].nvertices,gpiece[r].ncells);CHKERRQ(ierr);
212       /* Coordinate positions */
213       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      <Points>\n");CHKERRQ(ierr);
214       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"%s\" Name=\"Position\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,boffset);CHKERRQ(ierr);
215       boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
216       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      </Points>\n");CHKERRQ(ierr);
217       /* Cell connectivity */
218       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      <Cells>\n");CHKERRQ(ierr);
219       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
220       boffset += gpiece[r].nconn*sizeof(PetscInt) + sizeof(int);
221       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"offsets\"      NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
222       boffset += gpiece[r].ncells*sizeof(PetscInt) + sizeof(int);
223       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"UInt8\" Name=\"types\"        NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
224       boffset += gpiece[r].ncells*sizeof(unsigned char) + sizeof(int);
225       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      </Cells>\n");CHKERRQ(ierr);
226 
227       /*
228        * Cell Data headers
229        */
230       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"      <CellData>\n");CHKERRQ(ierr);
231       ierr     = PetscFPrintf(PETSC_COMM_SELF,fp,"        <DataArray type=\"Int32\" Name=\"Rank\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",boffset);CHKERRQ(ierr);
232       boffset += gpiece[r].ncells*sizeof(int) + sizeof(int);
233       /* all the vectors */
234       for (link=vtk->link; link; link=link->next) {
235         Vec        X = (Vec)link->vec;
236         DM         dmX = NULL;
237         PetscInt   bs,nfields,field;
238         const char *vecname = "";
239         PetscSection section;
240         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
241         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. */
242           ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr);
243         }
244         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
245         if (!dmX) dmX = dm;
246         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
247         if (!section) {ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);}
248         ierr = PetscSectionGetDof(section,cStart,&bs);CHKERRQ(ierr);
249         ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr);
250         field = 0;
251         if (link->field >= 0) {
252           field = link->field;
253           nfields = field + 1;
254         }
255         for (i=0; field<(nfields?nfields:1); field++) {
256           PetscInt     fbs,j;
257           PetscFV      fv = NULL;
258           PetscObject  f;
259           PetscClassId fClass;
260           const char *fieldname = NULL;
261           char       buf[256];
262           PetscBool    vector;
263           if (nfields) {        /* We have user-defined fields/components */
264             ierr = PetscSectionGetFieldDof(section,cStart,field,&fbs);CHKERRQ(ierr);
265             ierr = PetscSectionGetFieldName(section,field,&fieldname);CHKERRQ(ierr);
266           } else fbs = bs;      /* Say we have one field with 'bs' components */
267           ierr = DMGetField(dmX,field,NULL,&f);CHKERRQ(ierr);
268           ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr);
269           if (fClass == PETSCFV_CLASSID) {
270             fv = (PetscFV) f;
271           }
272           if (nfields && !fieldname) {
273             ierr = PetscSNPrintf(buf,sizeof(buf),"CellField%D",field);CHKERRQ(ierr);
274             fieldname = buf;
275           }
276           vector = PETSC_FALSE;
277           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
278             vector = PETSC_TRUE;
279             if (fbs > 3) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Cell vector fields can have at most 3 components, %D given\n", fbs);
280             for (j = 0; j < fbs; j++) {
281               const char *compName = NULL;
282               if (fv) {
283                 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr);
284                 if (compName) break;
285               }
286             }
287             if (j < fbs) vector = PETSC_FALSE;
288           }
289           if (vector) {
290 #if defined(PETSC_USE_COMPLEX)
291             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr);
292             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
293             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr);
294             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
295 #else
296             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr);
297             boffset += gpiece[r].ncells*3*sizeof(PetscVTUReal) + sizeof(int);
298 #endif
299             i+=fbs;
300           } else {
301             for (j=0; j<fbs; j++) {
302               const char *compName = NULL;
303               char finalname[256];
304               if (fv) {
305                 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr);
306               }
307               if (compName) {
308                 ierr = PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);CHKERRQ(ierr);
309               }
310               else if (fbs > 1) {
311                 ierr = PetscSNPrintf(finalname,255,"%s%s.%D",vecname,fieldname,j);CHKERRQ(ierr);
312               } else {
313                 ierr = PetscSNPrintf(finalname,255,"%s%s",vecname,fieldname);CHKERRQ(ierr);
314               }
315 #if defined(PETSC_USE_COMPLEX)
316               ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);CHKERRQ(ierr);
317               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
318               ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);CHKERRQ(ierr);
319               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
320 #else
321               ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);CHKERRQ(ierr);
322               boffset += gpiece[r].ncells*sizeof(PetscVTUReal) + sizeof(int);
323 #endif
324               i++;
325             }
326           }
327         }
328         if (i != bs) SETERRQ2(comm,PETSC_ERR_PLIB,"Total number of field components %D != block size %D",i,bs);
329       }
330       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"      </CellData>\n");CHKERRQ(ierr);
331 
332       /*
333        * Point Data headers
334        */
335       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"      <PointData>\n");CHKERRQ(ierr);
336       for (link=vtk->link; link; link=link->next) {
337         Vec        X = (Vec)link->vec;
338         DM         dmX;
339         PetscInt   bs,nfields,field;
340         const char *vecname = "";
341         PetscSection section;
342         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
343         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. */
344           ierr = PetscObjectGetName((PetscObject)X,&vecname);CHKERRQ(ierr);
345         }
346         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
347         if (!dmX) dmX = dm;
348         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
349         if (!section) {ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);}
350         ierr = PetscSectionGetDof(section,vStart,&bs);CHKERRQ(ierr);
351         ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr);
352         field = 0;
353         if (link->field >= 0) {
354           field = link->field;
355           nfields = field + 1;
356         }
357         for (i=0; field<(nfields?nfields:1); field++) {
358           PetscInt   fbs,j;
359           const char *fieldname = NULL;
360           char       buf[256];
361           if (nfields) {        /* We have user-defined fields/components */
362             ierr = PetscSectionGetFieldDof(section,vStart,field,&fbs);CHKERRQ(ierr);
363             ierr = PetscSectionGetFieldName(section,field,&fieldname);CHKERRQ(ierr);
364           } else fbs = bs;      /* Say we have one field with 'bs' components */
365           if (nfields && !fieldname) {
366             ierr = PetscSNPrintf(buf,sizeof(buf),"PointField%D",field);CHKERRQ(ierr);
367             fieldname = buf;
368           }
369           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
370             if (fbs > 3) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_ARG_SIZ,"Point vector fields can have at most 3 components, %D given\n", fbs);
371 #if defined(PETSC_USE_COMPLEX)
372             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Re\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr);
373             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
374             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s.Im\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr);
375             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
376 #else
377             ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s%s\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%D\" />\n",precision,vecname,fieldname,boffset);CHKERRQ(ierr);
378             boffset += gpiece[r].nvertices*3*sizeof(PetscVTUReal) + sizeof(int);
379 #endif
380           } else {
381             for (j=0; j<fbs; j++) {
382               const char *compName = NULL;
383               char finalname[256];
384               ierr = PetscSectionGetComponentName(section,field,j,&compName);CHKERRQ(ierr);
385               ierr = PetscSNPrintf(finalname,255,"%s%s.%s",vecname,fieldname,compName);CHKERRQ(ierr);
386 #if defined(PETSC_USE_COMPLEX)
387               ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Re\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);CHKERRQ(ierr);
388               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
389               ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s.Im\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);CHKERRQ(ierr);
390               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
391 #else
392               ierr = PetscFPrintf(comm,fp,"        <DataArray type=\"%s\" Name=\"%s\" NumberOfComponents=\"1\" format=\"appended\" offset=\"%D\" />\n",precision,finalname,boffset);CHKERRQ(ierr);
393               boffset += gpiece[r].nvertices*sizeof(PetscVTUReal) + sizeof(int);
394 #endif
395             }
396           }
397         }
398       }
399       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"      </PointData>\n");CHKERRQ(ierr);
400       ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"    </Piece>\n");CHKERRQ(ierr);
401     }
402   }
403 
404   ierr = PetscFPrintf(comm,fp,"  </UnstructuredGrid>\n");CHKERRQ(ierr);
405   ierr = PetscFPrintf(comm,fp,"  <AppendedData encoding=\"raw\">\n");CHKERRQ(ierr);
406   ierr = PetscFPrintf(comm,fp,"_");CHKERRQ(ierr);
407 
408   if (!rank) {
409     PetscInt maxsize = 0;
410     for (r=0; r<size; r++) {
411       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nvertices*3*sizeof(PetscVTUReal)));
412       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].ncells*3*sizeof(PetscVTUReal)));
413       maxsize = PetscMax(maxsize, (PetscInt) (gpiece[r].nconn*sizeof(PetscVTKInt)));
414     }
415     ierr = PetscMalloc(maxsize,&buffer);CHKERRQ(ierr);
416   }
417   for (r=0; r<size; r++) {
418     if (r == rank) {
419       PetscInt nsend;
420       {                         /* Position */
421         const PetscScalar *x;
422         PetscVTUReal      *y = NULL;
423         Vec               coords;
424         PetscBool         copy;
425 
426         ierr  = DMGetCoordinatesLocal(dm,&coords);CHKERRQ(ierr);
427         ierr  = VecGetArrayRead(coords,&x);CHKERRQ(ierr);
428 #if defined(PETSC_USE_COMPLEX)
429         copy = PETSC_TRUE;
430 #else
431         copy = (PetscBool) (dimEmbed != 3 || localized || (sizeof(PetscReal) != sizeof(PetscVTUReal)));
432 #endif
433         if (copy) {
434           ierr = PetscMalloc1(piece.nvertices*3,&y);CHKERRQ(ierr);
435           if (localized) {
436             PetscInt cnt;
437             for (c=cStart,cnt=0; c<cEnd; c++) {
438               PetscInt off, dof;
439 
440               ierr = PetscSectionGetDof(coordSection, c, &dof);CHKERRQ(ierr);
441               if (!dof) {
442                 PetscInt *closure = NULL;
443                 PetscInt closureSize;
444 
445                 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
446                 for (v = 0; v < closureSize*2; v += 2) {
447                   if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
448                     ierr = PetscSectionGetOffset(coordSection, closure[v], &off);CHKERRQ(ierr);
449                     if (dimEmbed != 3) {
450                       y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]);
451                       y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off+1]) : 0.0);
452                       y[cnt*3+2] = (PetscVTUReal) 0.0;
453                     } else {
454                       y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off+0]);
455                       y[cnt*3+1] = (PetscVTUReal) PetscRealPart(x[off+1]);
456                       y[cnt*3+2] = (PetscVTUReal) PetscRealPart(x[off+2]);
457                     }
458                     cnt++;
459                   }
460                 }
461                 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
462               } else {
463                 ierr = PetscSectionGetOffset(coordSection, c, &off);CHKERRQ(ierr);
464                 if (dimEmbed != 3) {
465                   for (i=0; i<dof/dimEmbed; i++) {
466                     y[cnt*3+0] = (PetscVTUReal) PetscRealPart(x[off + i*dimEmbed + 0]);
467                     y[cnt*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[off + i*dimEmbed + 1]) : 0.0);
468                     y[cnt*3+2] = (PetscVTUReal) 0.0;
469                     cnt++;
470                   }
471                 } else {
472                   for (i=0; i<dof; i ++) {
473                     y[cnt*3+i] = (PetscVTUReal) PetscRealPart(x[off + i]);
474                   }
475                   cnt += dof/dimEmbed;
476                 }
477               }
478             }
479             if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
480           } else {
481             for (i=0; i<piece.nvertices; i++) {
482               y[i*3+0] = (PetscVTUReal) PetscRealPart(x[i*dimEmbed+0]);
483               y[i*3+1] = (PetscVTUReal) ((dimEmbed > 1) ? PetscRealPart(x[i*dimEmbed+1]) : 0.);
484               y[i*3+2] = (PetscVTUReal) ((dimEmbed > 2) ? PetscRealPart(x[i*dimEmbed+2]) : 0.);
485             }
486           }
487         }
488         nsend = piece.nvertices*3;
489         ierr  = TransferWrite(viewer,fp,r,0,copy ? (const void *) y : (const void *) x,buffer,nsend,MPIU_VTUREAL,tag);CHKERRQ(ierr);
490         ierr  = PetscFree(y);CHKERRQ(ierr);
491         ierr  = VecRestoreArrayRead(coords,&x);CHKERRQ(ierr);
492       }
493       {                           /* Connectivity, offsets, types */
494         PetscVTKInt  *connectivity = NULL, *offsets = NULL;
495         PetscVTKType *types = NULL;
496         ierr = DMPlexGetVTKConnectivity(dm,localized,&piece,&connectivity,&offsets,&types);CHKERRQ(ierr);
497         ierr = TransferWrite(viewer,fp,r,0,connectivity,buffer,piece.nconn,MPI_INT,tag);CHKERRQ(ierr);
498         ierr = TransferWrite(viewer,fp,r,0,offsets,buffer,piece.ncells,MPI_INT,tag);CHKERRQ(ierr);
499         ierr = TransferWrite(viewer,fp,r,0,types,buffer,piece.ncells,MPI_CHAR,tag);CHKERRQ(ierr);
500         ierr = PetscFree3(connectivity,offsets,types);CHKERRQ(ierr);
501       }
502       {                         /* Owners (cell data) */
503         PetscVTKInt *owners;
504         ierr = PetscMalloc1(piece.ncells,&owners);CHKERRQ(ierr);
505         for (i=0; i<piece.ncells; i++) owners[i] = rank;
506         ierr = TransferWrite(viewer,fp,r,0,owners,buffer,piece.ncells,MPI_INT,tag);CHKERRQ(ierr);
507         ierr = PetscFree(owners);CHKERRQ(ierr);
508       }
509       /* Cell data */
510       for (link=vtk->link; link; link=link->next) {
511         Vec               X = (Vec)link->vec;
512         DM                dmX;
513         const PetscScalar *x;
514         PetscVTUReal      *y;
515         PetscInt          bs, nfields, field;
516         PetscSection      section = NULL;
517 
518         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
519         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
520         if (!dmX) dmX = dm;
521         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
522         if (!section) {ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);}
523         ierr = PetscSectionGetDof(section,cStart,&bs);CHKERRQ(ierr);
524         ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr);
525         field = 0;
526         if (link->field >= 0) {
527           field = link->field;
528           nfields = field + 1;
529         }
530         ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
531         ierr = PetscMalloc1(piece.ncells*3,&y);CHKERRQ(ierr);
532         for (i=0; field<(nfields?nfields:1); field++) {
533           PetscInt     fbs,j;
534           PetscFV      fv = NULL;
535           PetscObject  f;
536           PetscClassId fClass;
537           PetscBool    vector;
538           if (nfields) {        /* We have user-defined fields/components */
539             ierr = PetscSectionGetFieldDof(section,cStart,field,&fbs);CHKERRQ(ierr);
540           } else fbs = bs;      /* Say we have one field with 'bs' components */
541           ierr = DMGetField(dmX,field,NULL,&f);CHKERRQ(ierr);
542           ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr);
543           if (fClass == PETSCFV_CLASSID) {
544             fv = (PetscFV) f;
545           }
546           vector = PETSC_FALSE;
547           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
548             vector = PETSC_TRUE;
549             for (j = 0; j < fbs; j++) {
550               const char *compName = NULL;
551               if (fv) {
552                 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr);
553                 if (compName) break;
554               }
555             }
556             if (j < fbs) vector = PETSC_FALSE;
557           }
558           if (vector) {
559             PetscInt cnt, l;
560             for (l = 0; l < loops_per_scalar; l++) {
561               for (c=cStart,cnt=0; c<cEnd; c++) {
562                 const PetscScalar *xpoint;
563                 PetscInt off, j;
564 
565                 if (hasLabel) {     /* Ignore some cells */
566                   PetscInt value;
567                   ierr = DMGetLabelValue(dmX, "vtk", c, &value);CHKERRQ(ierr);
568                   if (value != 1) continue;
569                 }
570                 if (nfields) {
571                   ierr = PetscSectionGetFieldOffset(section, c, field, &off);CHKERRQ(ierr);
572                 } else {
573                   ierr = PetscSectionGetOffset(section, c, &off);CHKERRQ(ierr);
574                 }
575                 xpoint = &x[off];
576                 for (j = 0; j < fbs; j++) {
577                   y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
578                 }
579                 for (; j < 3; j++) y[cnt++] = 0.;
580               }
581               if (cnt != piece.ncells*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
582               ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells*3,MPIU_VTUREAL,tag);CHKERRQ(ierr);
583             }
584           } else {
585             for (i=0; i<fbs; i++) {
586               PetscInt cnt, l;
587               for (l = 0; l < loops_per_scalar; l++) {
588                 for (c=cStart,cnt=0; c<cEnd; c++) {
589                   const PetscScalar *xpoint;
590                   PetscInt off;
591 
592                   if (hasLabel) {     /* Ignore some cells */
593                     PetscInt value;
594                     ierr = DMGetLabelValue(dmX, "vtk", c, &value);CHKERRQ(ierr);
595                     if (value != 1) continue;
596                   }
597                   if (nfields) {
598                     ierr = PetscSectionGetFieldOffset(section, c, field, &off);CHKERRQ(ierr);
599                   } else {
600                     ierr = PetscSectionGetOffset(section, c, &off);CHKERRQ(ierr);
601                   }
602                   xpoint   = &x[off];
603                   y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
604                 }
605                 if (cnt != piece.ncells) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
606                 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.ncells,MPIU_VTUREAL,tag);CHKERRQ(ierr);
607               }
608             }
609           }
610         }
611         ierr = PetscFree(y);CHKERRQ(ierr);
612         ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
613       }
614       /* point data */
615       for (link=vtk->link; link; link=link->next) {
616         Vec               X = (Vec)link->vec;
617         DM                dmX;
618         const PetscScalar *x;
619         PetscVTUReal      *y;
620         PetscInt          bs, nfields, field;
621         PetscSection      section = NULL;
622 
623         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
624         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
625         if (!dmX) dmX = dm;
626         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
627         if (!section) {ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);}
628         ierr = PetscSectionGetDof(section,vStart,&bs);CHKERRQ(ierr);
629         ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr);
630         field = 0;
631         if (link->field >= 0) {
632           field = link->field;
633           nfields = field + 1;
634         }
635         ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
636         ierr = PetscMalloc1(piece.nvertices*3,&y);CHKERRQ(ierr);
637         for (i=0; field<(nfields?nfields:1); field++) {
638           PetscInt   fbs,j;
639           if (nfields) {        /* We have user-defined fields/components */
640             ierr = PetscSectionGetFieldDof(section,vStart,field,&fbs);CHKERRQ(ierr);
641           } else fbs = bs;      /* Say we have one field with 'bs' components */
642           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
643             PetscInt cnt, l;
644             for (l = 0; l < loops_per_scalar; l++) {
645               if (!localized) {
646                 for (v=vStart,cnt=0; v<vEnd; v++) {
647                   PetscInt    off;
648                   const PetscScalar *xpoint;
649 
650                   if (nfields) {
651                     ierr     = PetscSectionGetFieldOffset(section,v,field,&off);CHKERRQ(ierr);
652                   } else {
653                     ierr     = PetscSectionGetOffset(section,v,&off);CHKERRQ(ierr);
654                   }
655                   xpoint   = &x[off];
656                   for (j = 0; j < fbs; j++) {
657                     y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
658                   }
659                   for (; j < 3; j++) y[cnt++] = 0.;
660                 }
661               } else {
662                 for (c=cStart,cnt=0; c<cEnd; c++) {
663                   PetscInt *closure = NULL;
664                   PetscInt  closureSize, off;
665 
666                   ierr = DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
667                   for (v = 0, off = 0; v < closureSize*2; v += 2) {
668                     if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
669                       PetscInt    voff;
670                       const PetscScalar *xpoint;
671 
672                       if (nfields) {
673                         ierr = PetscSectionGetFieldOffset(section,closure[v],field,&voff);CHKERRQ(ierr);
674                       } else {
675                         ierr = PetscSectionGetOffset(section,closure[v],&voff);CHKERRQ(ierr);
676                       }
677                       xpoint         = &x[voff];
678                       for (j = 0; j < fbs; j++) {
679                         y[cnt + off++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[j]) : PetscRealPart(xpoint[j]));
680                       }
681                       for (; j < 3; j++) y[cnt + off++] = 0.;
682                     }
683                   }
684                   cnt += off;
685                   ierr = DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
686                 }
687               }
688               if (cnt != piece.nvertices*3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
689               ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices*3,MPIU_VTUREAL,tag);CHKERRQ(ierr);
690             }
691           } else {
692             for (i=0; i<fbs; i++) {
693               PetscInt cnt, l;
694               for (l = 0; l < loops_per_scalar; l++) {
695                 if (!localized) {
696                   for (v=vStart,cnt=0; v<vEnd; v++) {
697                     PetscInt    off;
698                     const PetscScalar *xpoint;
699 
700                     if (nfields) {
701                       ierr     = PetscSectionGetFieldOffset(section,v,field,&off);CHKERRQ(ierr);
702                     } else {
703                       ierr     = PetscSectionGetOffset(section,v,&off);CHKERRQ(ierr);
704                     }
705                     xpoint   = &x[off];
706                     y[cnt++] = (PetscVTUReal) (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
707                   }
708                 } else {
709                   for (c=cStart,cnt=0; c<cEnd; c++) {
710                     PetscInt *closure = NULL;
711                     PetscInt  closureSize, off;
712 
713                     ierr = DMPlexGetTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
714                     for (v = 0, off = 0; v < closureSize*2; v += 2) {
715                       if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
716                         PetscInt    voff;
717                         const PetscScalar *xpoint;
718 
719                         if (nfields) {
720                           ierr           = PetscSectionGetFieldOffset(section,closure[v],field,&voff);CHKERRQ(ierr);
721                         } else {
722                           ierr           = PetscSectionGetOffset(section,closure[v],&voff);CHKERRQ(ierr);
723                         }
724                         xpoint         = &x[voff];
725                         y[cnt + off++] = (l ? PetscImaginaryPart(xpoint[i]) : PetscRealPart(xpoint[i]));
726                       }
727                     }
728                     cnt += off;
729                     ierr = DMPlexRestoreTransitiveClosure(dmX, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
730                   }
731                 }
732                 if (cnt != piece.nvertices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not match");
733                 ierr = TransferWrite(viewer,fp,r,0,y,buffer,piece.nvertices,MPIU_VTUREAL,tag);CHKERRQ(ierr);
734               }
735             }
736           }
737         }
738         ierr = PetscFree(y);CHKERRQ(ierr);
739         ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
740       }
741     } else if (!rank) {
742       PetscInt l;
743 
744       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag);CHKERRQ(ierr); /* positions */
745       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nconn,MPI_INT,tag);CHKERRQ(ierr); /* connectivity */
746       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag);CHKERRQ(ierr); /* offsets */
747       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_CHAR,tag);CHKERRQ(ierr); /* types */
748       ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPI_INT,tag);CHKERRQ(ierr); /* owner rank (cells) */
749       /* all cell data */
750       for (link=vtk->link; link; link=link->next) {
751         Vec               X = (Vec)link->vec;
752         PetscInt bs, nfields, field;
753         DM           dmX;
754         PetscSection section = NULL;
755 
756         if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
757         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
758         if (!dmX) dmX = dm;
759         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
760         if (!section) {ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);}
761         ierr = PetscSectionGetDof(section,cStart,&bs);CHKERRQ(ierr);
762         ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr);
763         field = 0;
764         if (link->field >= 0) {
765           field = link->field;
766           nfields = field + 1;
767         }
768         for (i=0; field<(nfields?nfields:1); field++) {
769           PetscInt     fbs,j;
770           PetscFV      fv = NULL;
771           PetscObject  f;
772           PetscClassId fClass;
773           PetscBool    vector;
774           if (nfields) {        /* We have user-defined fields/components */
775             ierr = PetscSectionGetFieldDof(section,cStart,field,&fbs);CHKERRQ(ierr);
776           } else fbs = bs;      /* Say we have one field with 'bs' components */
777           ierr = DMGetField(dmX,field,NULL,&f);CHKERRQ(ierr);
778           ierr = PetscObjectGetClassId(f,&fClass);CHKERRQ(ierr);
779           if (fClass == PETSCFV_CLASSID) {
780             fv = (PetscFV) f;
781           }
782           vector = PETSC_FALSE;
783           if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) {
784             vector = PETSC_TRUE;
785             for (j = 0; j < fbs; j++) {
786               const char *compName = NULL;
787               if (fv) {
788                 ierr = PetscFVGetComponentName(fv,j,&compName);CHKERRQ(ierr);
789                 if (compName) break;
790               }
791             }
792             if (j < fbs) vector = PETSC_FALSE;
793           }
794           if (vector) {
795             for (l = 0; l < loops_per_scalar; l++) {
796               ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells*3,MPIU_VTUREAL,tag);CHKERRQ(ierr);
797             }
798           } else {
799             for (i=0; i<fbs; i++) {
800               for (l = 0; l < loops_per_scalar; l++) {
801                 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].ncells,MPIU_VTUREAL,tag);CHKERRQ(ierr);
802               }
803             }
804           }
805         }
806       }
807       /* all point data */
808       for (link=vtk->link; link; link=link->next) {
809         Vec               X = (Vec)link->vec;
810         DM                dmX;
811         PetscInt bs, nfields, field;
812         PetscSection section = NULL;
813 
814         if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
815         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
816         if (!dmX) dmX = dm;
817         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
818         if (!section) {ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);}
819         ierr = PetscSectionGetDof(section,vStart,&bs);CHKERRQ(ierr);
820         ierr = PetscSectionGetNumFields(section,&nfields);CHKERRQ(ierr);
821         field = 0;
822         if (link->field >= 0) {
823           field = link->field;
824           nfields = field + 1;
825         }
826         for (i=0; field<(nfields?nfields:1); field++) {
827           PetscInt   fbs;
828           if (nfields) {        /* We have user-defined fields/components */
829             ierr = PetscSectionGetFieldDof(section,vStart,field,&fbs);CHKERRQ(ierr);
830           } else fbs = bs;      /* Say we have one field with 'bs' components */
831           if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) {
832             for (l = 0; l < loops_per_scalar; l++) {
833               ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices*3,MPIU_VTUREAL,tag);CHKERRQ(ierr);
834             }
835           } else {
836             for (i=0; i<fbs; i++) {
837               for (l = 0; l < loops_per_scalar; l++) {
838                 ierr = TransferWrite(viewer,fp,r,0,NULL,buffer,gpiece[r].nvertices,MPIU_VTUREAL,tag);CHKERRQ(ierr);
839               }
840             }
841           }
842         }
843       }
844     }
845   }
846   ierr = PetscFree(gpiece);CHKERRQ(ierr);
847   ierr = PetscFree(buffer);CHKERRQ(ierr);
848   ierr = PetscFPrintf(comm,fp,"\n  </AppendedData>\n");CHKERRQ(ierr);
849   ierr = PetscFPrintf(comm,fp,"</VTKFile>\n");CHKERRQ(ierr);
850   ierr = PetscFClose(comm,fp);CHKERRQ(ierr);
851   PetscFunctionReturn(0);
852 }
853