xref: /petsc/src/dm/impls/plex/plexvtk.c (revision a30ec4eab21047868b4ea34d315622f953eb6ffa) 
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;
458   PetscErrorCode           ierr;
459 
460   PetscFunctionBegin;
461   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
462   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
463   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
464   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
465   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
466   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
467   /* Vertices */
468   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
469   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
470   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
471   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
472   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
473   if (vMax >= 0) {
474     PetscInt pStart, pEnd, p, localSize = 0;
475 
476     ierr = PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);CHKERRQ(ierr);
477     pEnd = PetscMin(pEnd, vMax);
478     for (p = pStart; p < pEnd; ++p) {
479       PetscInt dof;
480 
481       ierr = PetscSectionGetDof(globalCoordSection, p, &dof);CHKERRQ(ierr);
482       if (dof > 0) ++localSize;
483     }
484     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);CHKERRQ(ierr);
485     ierr = PetscLayoutSetLocalSize(vLayout, localSize);CHKERRQ(ierr);
486     ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
487     ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
488   } else {
489     ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);CHKERRQ(ierr);
490   }
491   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
492   ierr = PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);CHKERRQ(ierr);
493   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);CHKERRQ(ierr);
494   /* Cells */
495   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
496   /* Vertex fields */
497   for (link = vtk->link; link; link = link->next) {
498     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
499     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
500   }
501   if (hasPoint) {
502     ierr = PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);CHKERRQ(ierr);
503     for (link = vtk->link; link; link = link->next) {
504       Vec          X = (Vec) link->vec;
505       DM           dmX;
506       PetscSection section, globalSection, newSection = NULL;
507       const char   *name;
508       PetscInt     enforceDof = PETSC_DETERMINE;
509 
510       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
511       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
512       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
513       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
514       if (dmX) {
515         DMLabel  subpointMap, subpointMapX;
516         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
517 
518         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
519         /* Here is where we check whether dmX is a submesh of dm */
520         ierr = DMGetDimension(dm,  &dim);CHKERRQ(ierr);
521         ierr = DMGetDimension(dmX, &dimX);CHKERRQ(ierr);
522         ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
523         ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
524         ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
525         ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
526         if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
527           const PetscInt *ind = NULL;
528           IS              subpointIS;
529           PetscInt        n = 0, q;
530 
531           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
532           ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
533           if (subpointIS) {
534             ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
535             ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
536           }
537           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
538           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
539           for (q = qStart; q < qEnd; ++q) {
540             PetscInt dof, off, p;
541 
542             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
543             if (dof) {
544               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
545               if (p >= pStart) {
546                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
547                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
548                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
549               }
550             }
551           }
552           if (subpointIS) {
553             ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
554             ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
555           }
556           /* No need to setup section */
557           section = newSection;
558         }
559       } else {
560         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
561         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
562       }
563       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
564       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
565       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
566       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
567       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
568     }
569   }
570   /* Cell Fields */
571   ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
572   if (hasCell || writePartition) {
573     ierr = PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);CHKERRQ(ierr);
574     for (link = vtk->link; link; link = link->next) {
575       Vec          X = (Vec) link->vec;
576       DM           dmX;
577       PetscSection section, globalSection;
578       const char   *name;
579       PetscInt     enforceDof = PETSC_DETERMINE;
580 
581       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
582       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
583       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
584       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
585       if (dmX) {
586         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
587       } else {
588         PetscContainer c;
589 
590         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
591         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
592         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
593       }
594       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
595       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
596       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
597       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
598     }
599     if (writePartition) {
600       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
601       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
602       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
603     }
604   }
605   /* Cleanup */
606   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
607   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
608   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
609   PetscFunctionReturn(0);
610 }
611 
612 /*@C
613   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
614 
615   Collective
616 
617   Input Arguments:
618 + odm - The DMPlex specifying the mesh, passed as a PetscObject
619 - viewer - viewer of type VTK
620 
621   Level: developer
622 
623   Note:
624   This function is a callback used by the VTK viewer to actually write the file.
625   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
626   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
627 
628 .seealso: PETSCVIEWERVTK
629 @*/
630 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
631 {
632   DM             dm = (DM) odm;
633   PetscBool      isvtk;
634   PetscErrorCode ierr;
635 
636   PetscFunctionBegin;
637   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
638   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
639   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
640   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
641   switch (viewer->format) {
642   case PETSC_VIEWER_ASCII_VTK:
643     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
644     break;
645   case PETSC_VIEWER_VTK_VTU:
646     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
647     break;
648   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
649   }
650   PetscFunctionReturn(0);
651 }
652