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