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