xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 0adebc6ce4218bb5a986f93c65179d861410e0be) 
1 #define PETSCDM_DLL
2 #include <petsc-private/pleximpl.h>    /*I   "petscdmplex.h"   I*/
3 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexVTKGetCellType"
7 PetscErrorCode DMPlexVTKGetCellType(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
8 {
9   PetscFunctionBegin;
10   *cellType = -1;
11   switch(dim) {
12   case 0:
13     switch(corners) {
14     case 1:
15       *cellType = 1; /* VTK_VERTEX */
16       break;
17     default:
18       break;
19     }
20     break;
21   case 1:
22     switch(corners) {
23     case 2:
24       *cellType = 3; /* VTK_LINE */
25       break;
26     case 3:
27       *cellType = 21; /* VTK_QUADRATIC_EDGE */
28       break;
29     default:
30       break;
31     }
32     break;
33   case 2:
34     switch(corners) {
35     case 3:
36       *cellType = 5; /* VTK_TRIANGLE */
37       break;
38     case 4:
39       *cellType = 9; /* VTK_QUAD */
40       break;
41     case 6:
42       *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
43       break;
44     case 9:
45       *cellType = 23; /* VTK_QUADRATIC_QUAD */
46       break;
47     default:
48       break;
49     }
50     break;
51   case 3:
52     switch(corners) {
53     case 4:
54       *cellType = 10; /* VTK_TETRA */
55       break;
56     case 8:
57       *cellType = 12; /* VTK_HEXAHEDRON */
58       break;
59     case 10:
60       *cellType = 24; /* VTK_QUADRATIC_TETRA */
61       break;
62     case 27:
63       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
64       break;
65     default:
66       break;
67     }
68   }
69   PetscFunctionReturn(0);
70 }
71 
72 #undef __FUNCT__
73 #define __FUNCT__ "DMPlexVTKWriteCells_ASCII"
74 PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
75 {
76   MPI_Comm        comm = ((PetscObject) dm)->comm;
77   IS              globalVertexNumbers;
78   const PetscInt *gvertex;
79   PetscInt        dim;
80   PetscInt        numCorners = 0, totCorners = 0, maxCorners, *corners;
81   PetscInt        numCells   = 0, totCells   = 0, maxCells, cellHeight;
82   PetscInt        numLabelCells, cMax, cStart, cEnd, c, vStart, vEnd, v;
83   PetscMPIInt     numProcs, rank, proc, tag;
84   PetscBool       hasLabel;
85   PetscErrorCode  ierr;
86 
87   PetscFunctionBegin;
88   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
89   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
90   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
91   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
92   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
93   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
94   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
95   ierr = DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr);
96   if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
97   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
98   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
99   for(c = cStart; c < cEnd; ++c) {
100     PetscInt *closure = PETSC_NULL;
101     PetscInt  closureSize;
102 
103     if (hasLabel) {
104       PetscInt value;
105 
106       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
107       if (value != 1) continue;
108     }
109     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
110     for (v = 0; v < closureSize*2; v += 2) {
111       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
112     }
113     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
114     ++numCells;
115   }
116   maxCells = numCells;
117   ierr = MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
118   ierr = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
119   ierr = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
120   ierr = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
121   ierr = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
122   ierr = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
123   ierr = PetscMalloc(maxCells * sizeof(PetscInt), &corners);CHKERRQ(ierr);
124   ierr = PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);CHKERRQ(ierr);
125   if (!rank) {
126     PetscInt *remoteVertices;
127 
128     for(c = cStart, numCells = 0; c < cEnd; ++c) {
129       PetscInt *closure = PETSC_NULL;
130       PetscInt  closureSize, nC = 0;
131 
132       if (hasLabel) {
133         PetscInt value;
134 
135         ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
136         if (value != 1) continue;
137       }
138       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
139       for (v = 0; v < closureSize*2; v += 2) {
140         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
141           const PetscInt gv = gvertex[closure[v] - vStart];
142           closure[nC++] = gv < 0 ? -(gv+1) : gv;
143         }
144       }
145       corners[numCells++] = nC;
146       ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
147       for(v = 0; v < nC; ++v) {
148         ierr = PetscFPrintf(comm, fp, " %d", closure[v]);CHKERRQ(ierr);
149       }
150       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
151       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
152     }
153     if (numProcs > 1) {ierr = PetscMalloc((maxCorners+maxCells) * sizeof(PetscInt), &remoteVertices);CHKERRQ(ierr);}
154     for(proc = 1; proc < numProcs; ++proc) {
155       MPI_Status status;
156 
157       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
158       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
159       for(c = 0; c < numCorners;) {
160         PetscInt nC = remoteVertices[c++];
161 
162         ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
163         for(v = 0; v < nC; ++v, ++c) {
164           ierr = PetscFPrintf(comm, fp, " %d", remoteVertices[c]);CHKERRQ(ierr);
165         }
166         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
167       }
168     }
169     if (numProcs > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
170   } else {
171     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
172 
173     ierr = PetscMalloc(numSend * sizeof(PetscInt), &localVertices);CHKERRQ(ierr);
174     for (c = cStart, numCells = 0; c < cEnd; ++c) {
175       PetscInt *closure = PETSC_NULL;
176       PetscInt  closureSize, nC = 0;
177 
178       if (hasLabel) {
179         PetscInt value;
180 
181         ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
182         if (value != 1) continue;
183       }
184       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
185       for (v = 0; v < closureSize*2; v += 2) {
186         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
187           const PetscInt gv = gvertex[closure[v] - vStart];
188           closure[nC++] = gv < 0 ? -(gv+1) : gv;
189         }
190       }
191       corners[numCells++]  = nC;
192       localVertices[k++] = nC;
193       for(v = 0; v < nC; ++v, ++k) {
194         localVertices[k] = closure[v];
195       }
196       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
197     }
198     if (k != numSend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numSend);
199     ierr = MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
200     ierr = MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
201     ierr = PetscFree(localVertices);CHKERRQ(ierr);
202   }
203   ierr = ISRestoreIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
204   ierr = PetscFPrintf(comm, fp, "CELL_TYPES %d\n", totCells);CHKERRQ(ierr);
205   if (!rank) {
206     PetscInt cellType;
207 
208     for (c = 0; c < numCells; ++c) {
209       ierr = DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
210       ierr = PetscFPrintf(comm, fp, "%d\n", cellType);CHKERRQ(ierr);
211     }
212     for(proc = 1; proc < numProcs; ++proc) {
213       MPI_Status status;
214 
215       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
216       ierr = MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
217       for(c = 0; c < numCells; ++c) {
218         ierr = DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
219         ierr = PetscFPrintf(comm, fp, "%d\n", cellType);CHKERRQ(ierr);
220       }
221     }
222   } else {
223     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
224     ierr = MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
225   }
226   ierr = PetscFree(corners);CHKERRQ(ierr);
227   *totalCells = totCells;
228   PetscFunctionReturn(0);
229 }
230 
231 #undef __FUNCT__
232 #define __FUNCT__ "DMPlexVTKWriteSection_ASCII"
233 PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
234 {
235   MPI_Comm           comm    = ((PetscObject) dm)->comm;
236   const MPI_Datatype mpiType = MPIU_SCALAR;
237   PetscScalar       *array;
238   PetscInt           numDof = 0, maxDof;
239   PetscInt           numLabelCells, cellHeight, cMax, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
240   PetscMPIInt        numProcs, rank, proc, tag;
241   PetscBool          hasLabel;
242   PetscErrorCode     ierr;
243 
244   PetscFunctionBegin;
245   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
246   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
247   if (precision < 0) precision = 6;
248   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
249   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
250   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
251   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
252   /* VTK only wants the values at cells or vertices */
253   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
254   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
255   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
256   ierr = DMPlexGetHybridBounds(dm, &cMax, PETSC_NULL, PETSC_NULL, &vMax);CHKERRQ(ierr);
257   if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
258   if (vMax >= 0) {vEnd = PetscMin(vEnd, vMax);}
259   pStart = PetscMax(PetscMin(cStart, vStart), pStart);
260   pEnd   = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
261   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
262   ierr = DMPlexGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
263   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
264   for(p = pStart; p < pEnd; ++p) {
265     /* Reject points not either cells or vertices */
266     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
267     if (hasLabel) {
268       PetscInt value;
269 
270       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
271           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
272         ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
273         if (value != 1) continue;
274       }
275     }
276     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
277     if (numDof) {break;}
278   }
279   ierr = MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
280   enforceDof = PetscMax(enforceDof, maxDof);
281   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
282   if (!rank) {
283     char formatString[8];
284 
285     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
286     for (p = pStart; p < pEnd; ++p) {
287       /* Here we lose a way to filter points by keeping them out of the Numbering */
288       PetscInt dof, off, goff, d;
289 
290       /* Reject points not either cells or vertices */
291       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
292       if (hasLabel) {
293         PetscInt value;
294 
295         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
296             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
297           ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
298           if (value != 1) continue;
299         }
300       }
301       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
302       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
303       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
304       if (dof && goff >= 0) {
305         for (d = 0; d < dof; d++) {
306           if (d > 0) {
307             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
308           }
309           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);CHKERRQ(ierr);
310         }
311         for (d = dof; d < enforceDof; d++) {
312           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
313         }
314         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
315       }
316     }
317     for (proc = 1; proc < numProcs; ++proc) {
318       PetscScalar *remoteValues;
319       PetscInt     size, d;
320       MPI_Status   status;
321 
322       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
323       ierr = PetscMalloc(size * sizeof(PetscScalar), &remoteValues);CHKERRQ(ierr);
324       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRQ(ierr);
325       for (p = 0; p < size/maxDof; ++p) {
326         for (d = 0; d < maxDof; ++d) {
327           if (d > 0) {
328             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
329           }
330           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);CHKERRQ(ierr);
331         }
332         for (d = maxDof; d < enforceDof; ++d) {
333           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
334         }
335         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
336       }
337       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
338     }
339   } else {
340     PetscScalar *localValues;
341     PetscInt     size, k = 0;
342 
343     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
344     ierr = PetscMalloc(size * sizeof(PetscScalar), &localValues);CHKERRQ(ierr);
345     for (p = pStart; p < pEnd; ++p) {
346       PetscInt dof, off, goff, d;
347 
348       /* Reject points not either cells or vertices */
349       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
350       if (hasLabel) {
351         PetscInt value;
352 
353         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
354             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
355           ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
356           if (value != 1) continue;
357         }
358       }
359       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
360       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
361       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
362       if (goff >= 0) {
363         for (d = 0; d < dof; ++d) {
364           localValues[k++] = array[off+d];
365         }
366       }
367     }
368     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
369     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRQ(ierr);
370     ierr = PetscFree(localValues);CHKERRQ(ierr);
371   }
372   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
373   PetscFunctionReturn(0);
374 }
375 
376 #undef __FUNCT__
377 #define __FUNCT__ "DMPlexVTKWriteField_ASCII"
378 PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
379 {
380   MPI_Comm       comm = ((PetscObject) dm)->comm;
381   PetscInt       numDof = 0, maxDof;
382   PetscInt       pStart, pEnd, p;
383   PetscErrorCode ierr;
384 
385   PetscFunctionBegin;
386   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
387   for (p = pStart; p < pEnd; ++p) {
388     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
389     if (numDof) {break;}
390   }
391   numDof = PetscMax(numDof, enforceDof);
392   ierr = MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, ((PetscObject) dm)->comm);CHKERRQ(ierr);
393   if (!name) name = "Unknown";
394   if (maxDof == 3) {
395     ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
396   } else {
397     ierr = PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);CHKERRQ(ierr);
398     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
399   }
400   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);CHKERRQ(ierr);
401   PetscFunctionReturn(0);
402 }
403 
404 #undef __FUNCT__
405 #define __FUNCT__ "DMPlexVTKWriteAll_ASCII"
406 static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
407 {
408   MPI_Comm                 comm = ((PetscObject) dm)->comm;
409   PetscViewer_VTK         *vtk  = (PetscViewer_VTK *) viewer->data;
410   FILE                    *fp;
411   PetscViewerVTKObjectLink link;
412   PetscSection             coordSection, globalCoordSection;
413   PetscLayout              vLayout;
414   Vec                      coordinates;
415   PetscReal                lengthScale;
416   PetscInt                 vMax, totVertices, totCells;
417   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE;
418   PetscErrorCode           ierr;
419 
420   PetscFunctionBegin;
421   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
422   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
423   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
424   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
425   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
426   /* Vertices */
427   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
428   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
429   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
430   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
431   ierr = DMPlexGetHybridBounds(dm, PETSC_NULL, PETSC_NULL, PETSC_NULL, &vMax);CHKERRQ(ierr);
432   if (vMax >= 0) {
433     PetscInt pStart, pEnd, p, localSize = 0;
434 
435     ierr = PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);CHKERRQ(ierr);
436     pEnd = PetscMin(pEnd, vMax);
437     for (p = pStart; p < pEnd; ++p) {
438       PetscInt dof;
439 
440       ierr = PetscSectionGetDof(globalCoordSection, p, &dof);CHKERRQ(ierr);
441       if (dof > 0) {++localSize;}
442     }
443     ierr = PetscLayoutCreate(((PetscObject) dm)->comm, &vLayout);CHKERRQ(ierr);
444     ierr = PetscLayoutSetLocalSize(vLayout, localSize);CHKERRQ(ierr);
445     ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
446     ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
447   } else {
448     ierr = PetscSectionGetPointLayout(((PetscObject) dm)->comm, globalCoordSection, &vLayout);CHKERRQ(ierr);
449   }
450   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
451   ierr = PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);CHKERRQ(ierr);
452   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);CHKERRQ(ierr);
453   /* Cells */
454   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
455   /* Vertex fields */
456   for (link = vtk->link; link; link = link->next) {
457     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
458     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
459   }
460   if (hasPoint) {
461     ierr = PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);CHKERRQ(ierr);
462     for (link = vtk->link; link; link = link->next) {
463       Vec          X = (Vec) link->vec;
464       DM           dmX;
465       PetscSection section, globalSection;
466       const char  *name;
467       PetscInt     enforceDof = PETSC_DETERMINE;
468 
469       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
470       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
471       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
472       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
473       if (dmX) {
474         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
475       } else {
476         PetscContainer c;
477 
478         ierr = PetscObjectQuery(link->vec, "section", (PetscObject *) &c);CHKERRQ(ierr);
479         if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
480         ierr = PetscContainerGetPointer(c, (void **) &section);CHKERRQ(ierr);
481       }
482       if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
483       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
484       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
485       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
486     }
487   }
488   /* Cell Fields */
489   if (hasCell) {
490     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
491     for (link = vtk->link; link; link = link->next) {
492       Vec            X = (Vec) link->vec;
493       DM             dmX;
494       PetscSection   section, globalSection;
495       const char    *name;
496       PetscInt       enforceDof = PETSC_DETERMINE;
497 
498       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
499       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
500       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
501       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
502       if (dmX) {
503         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
504       } else {
505         PetscContainer c;
506 
507         ierr = PetscObjectQuery(link->vec, "section", (PetscObject *) &c);CHKERRQ(ierr);
508         if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
509         ierr = PetscContainerGetPointer(c, (void **) &section);CHKERRQ(ierr);
510       }
511       if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
512       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
513       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
514       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
515     }
516   }
517   /* Cleanup */
518   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
519   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
520   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
521   PetscFunctionReturn(0);
522 }
523 
524 #undef __FUNCT__
525 #define __FUNCT__ "DMPlexVTKWriteAll"
526 /*@C
527   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
528 
529   Collective
530 
531   Input Arguments:
532 + odm - The DMPlex specifying the mesh, passed as a PetscObject
533 - viewer - viewer of type VTK
534 
535   Level: developer
536 
537   Note:
538   This function is a callback used by the VTK viewer to actually write the file.
539   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
540   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
541 
542 .seealso: PETSCVIEWERVTK
543 @*/
544 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
545 {
546   DM              dm   = (DM) odm;
547   PetscBool       isvtk;
548   PetscErrorCode  ierr;
549 
550   PetscFunctionBegin;
551   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
552   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
553   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
554   if (!isvtk) SETERRQ1(((PetscObject) viewer)->comm, PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
555   switch (viewer->format) {
556   case PETSC_VIEWER_ASCII_VTK:
557     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
558     break;
559   case PETSC_VIEWER_VTK_VTU:
560     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
561     break;
562   default: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
563   }
564   PetscFunctionReturn(0);
565 }
566