xref: /petsc/src/dm/impls/plex/plexvtk.c (revision c30958fdb06df61c40a76215b00e4d9018fd89ec) 
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, newSection = PETSC_NULL;
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         IS        subpointMap, subpointMapX;
475         PetscInt  dim, dimX, pStart, pEnd, qStart, qEnd;
476 
477         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
478         /* Here is where we check whether dmX is a submesh of dm */
479         ierr = DMPlexGetDimension(dm,  &dim);CHKERRQ(ierr);
480         ierr = DMPlexGetDimension(dmX, &dimX);CHKERRQ(ierr);
481         ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
482         ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
483         ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
484         ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
485         if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
486           const PetscInt *ind;
487           PetscInt        n, q;
488 
489           ierr = PetscPrintf(PETSC_COMM_SELF, "Making translation PetscSection\n");CHKERRQ(ierr);
490           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
491           ierr = ISGetLocalSize(subpointMap, &n);CHKERRQ(ierr);
492           ierr = ISGetIndices(subpointMap, &ind);CHKERRQ(ierr);
493           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
494           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
495           for(q = qStart; q < qEnd; ++q) {
496             PetscInt dof, off, p;
497 
498             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
499             if (dof) {
500               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
501               if (p >= pStart) {
502                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
503                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
504                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
505               }
506             }
507           }
508           ierr = ISRestoreIndices(subpointMap, &ind);CHKERRQ(ierr);
509           /* No need to setup section */
510           section = newSection;
511         }
512       } else {
513         PetscContainer c;
514 
515         ierr = PetscObjectQuery(link->vec, "section", (PetscObject *) &c);CHKERRQ(ierr);
516         if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
517         ierr = PetscContainerGetPointer(c, (void **) &section);CHKERRQ(ierr);
518       }
519       if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
520       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
521       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
522       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
523       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
524     }
525   }
526   /* Cell Fields */
527   if (hasCell) {
528     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
529     for (link = vtk->link; link; link = link->next) {
530       Vec            X = (Vec) link->vec;
531       DM             dmX;
532       PetscSection   section, globalSection;
533       const char    *name;
534       PetscInt       enforceDof = PETSC_DETERMINE;
535 
536       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
537       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
538       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
539       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
540       if (dmX) {
541         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
542       } else {
543         PetscContainer c;
544 
545         ierr = PetscObjectQuery(link->vec, "section", (PetscObject *) &c);CHKERRQ(ierr);
546         if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
547         ierr = PetscContainerGetPointer(c, (void **) &section);CHKERRQ(ierr);
548       }
549       if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
550       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
551       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
552       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
553     }
554   }
555   /* Cleanup */
556   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
557   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
558   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
559   PetscFunctionReturn(0);
560 }
561 
562 #undef __FUNCT__
563 #define __FUNCT__ "DMPlexVTKWriteAll"
564 /*@C
565   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
566 
567   Collective
568 
569   Input Arguments:
570 + odm - The DMPlex specifying the mesh, passed as a PetscObject
571 - viewer - viewer of type VTK
572 
573   Level: developer
574 
575   Note:
576   This function is a callback used by the VTK viewer to actually write the file.
577   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
578   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
579 
580 .seealso: PETSCVIEWERVTK
581 @*/
582 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
583 {
584   DM              dm   = (DM) odm;
585   PetscBool       isvtk;
586   PetscErrorCode  ierr;
587 
588   PetscFunctionBegin;
589   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
590   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
591   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
592   if (!isvtk) SETERRQ1(((PetscObject) viewer)->comm, PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
593   switch (viewer->format) {
594   case PETSC_VIEWER_ASCII_VTK:
595     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
596     break;
597   case PETSC_VIEWER_VTK_VTU:
598     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
599     break;
600   default: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
601   }
602   PetscFunctionReturn(0);
603 }
604