xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 0298fd7132830bec7daee99a80be0eddb2b310a5) 
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, NULL, NULL, NULL);CHKERRQ(ierr);
96   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
97   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
98 
99   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
100   for (c = cStart; c < cEnd; ++c) {
101     PetscInt *closure = NULL;
102     PetscInt closureSize;
103 
104     if (hasLabel) {
105       PetscInt value;
106 
107       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
108       if (value != 1) continue;
109     }
110     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
111     for (v = 0; v < closureSize*2; v += 2) {
112       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
113     }
114     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
115     ++numCells;
116   }
117   maxCells = numCells;
118   ierr     = MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
119   ierr     = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
120   ierr     = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
121   ierr     = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
122   ierr     = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
123   ierr     = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
124   ierr     = PetscMalloc(maxCells * sizeof(PetscInt), &corners);CHKERRQ(ierr);
125   ierr     = PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);CHKERRQ(ierr);
126   if (!rank) {
127     PetscInt *remoteVertices;
128 
129     for (c = cStart, numCells = 0; c < cEnd; ++c) {
130       PetscInt *closure = NULL;
131       PetscInt closureSize, nC = 0;
132 
133       if (hasLabel) {
134         PetscInt value;
135 
136         ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
137         if (value != 1) continue;
138       }
139       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
140       for (v = 0; v < closureSize*2; v += 2) {
141         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
142           const PetscInt gv = gvertex[closure[v] - vStart];
143           closure[nC++] = gv < 0 ? -(gv+1) : gv;
144         }
145       }
146       corners[numCells++] = nC;
147       ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
148       for (v = 0; v < nC; ++v) {
149         ierr = PetscFPrintf(comm, fp, " %d", closure[v]);CHKERRQ(ierr);
150       }
151       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
152       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
153     }
154     if (numProcs > 1) {ierr = PetscMalloc((maxCorners+maxCells) * sizeof(PetscInt), &remoteVertices);CHKERRQ(ierr);}
155     for (proc = 1; proc < numProcs; ++proc) {
156       MPI_Status status;
157 
158       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
159       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
160       for (c = 0; c < numCorners;) {
161         PetscInt nC = remoteVertices[c++];
162 
163         ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
164         for (v = 0; v < nC; ++v, ++c) {
165           ierr = PetscFPrintf(comm, fp, " %d", remoteVertices[c]);CHKERRQ(ierr);
166         }
167         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
168       }
169     }
170     if (numProcs > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
171   } else {
172     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
173 
174     ierr = PetscMalloc(numSend * sizeof(PetscInt), &localVertices);CHKERRQ(ierr);
175     for (c = cStart, numCells = 0; c < cEnd; ++c) {
176       PetscInt *closure = NULL;
177       PetscInt closureSize, nC = 0;
178 
179       if (hasLabel) {
180         PetscInt value;
181 
182         ierr = DMPlexGetLabelValue(dm, "vtk", 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(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
211       ierr = PetscFPrintf(comm, fp, "%d\n", cellType);CHKERRQ(ierr);
212     }
213     for (proc = 1; proc < numProcs; ++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(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 #undef __FUNCT__
233 #define __FUNCT__ "DMPlexVTKWriteSection_ASCII"
234 PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
235 {
236   MPI_Comm           comm    = ((PetscObject) dm)->comm;
237   const MPI_Datatype mpiType = MPIU_SCALAR;
238   PetscScalar        *array;
239   PetscInt           numDof = 0, maxDof;
240   PetscInt           numLabelCells, cellHeight, cMax, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
241   PetscMPIInt        numProcs, rank, proc, tag;
242   PetscBool          hasLabel;
243   PetscErrorCode     ierr;
244 
245   PetscFunctionBegin;
246   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
247   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
248   if (precision < 0) precision = 6;
249   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
250   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
251   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
252   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
253   /* VTK only wants the values at cells or vertices */
254   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
255   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
256   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
257   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, &vMax);CHKERRQ(ierr);
258   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
259   if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
260   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
261   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
262   ierr     = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
263   ierr     = DMPlexGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
264   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
265   for (p = pStart; p < pEnd; ++p) {
266     /* Reject points not either cells or vertices */
267     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
268     if (hasLabel) {
269       PetscInt value;
270 
271       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
272           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
273         ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
274         if (value != 1) continue;
275       }
276     }
277     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
278     if (numDof) break;
279   }
280   ierr = MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
281   enforceDof = PetscMax(enforceDof, maxDof);
282   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
283   if (!rank) {
284     char formatString[8];
285 
286     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
287     for (p = pStart; p < pEnd; ++p) {
288       /* Here we lose a way to filter points by keeping them out of the Numbering */
289       PetscInt dof, off, goff, d;
290 
291       /* Reject points not either cells or vertices */
292       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
293       if (hasLabel) {
294         PetscInt value;
295 
296         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
297             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
298           ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
299           if (value != 1) continue;
300         }
301       }
302       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
303       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
304       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
305       if (dof && goff >= 0) {
306         for (d = 0; d < dof; d++) {
307           if (d > 0) {
308             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
309           }
310           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);CHKERRQ(ierr);
311         }
312         for (d = dof; d < enforceDof; d++) {
313           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
314         }
315         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
316       }
317     }
318     for (proc = 1; proc < numProcs; ++proc) {
319       PetscScalar *remoteValues;
320       PetscInt    size, d;
321       MPI_Status  status;
322 
323       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
324       ierr = PetscMalloc(size * sizeof(PetscScalar), &remoteValues);CHKERRQ(ierr);
325       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRQ(ierr);
326       for (p = 0; p < size/maxDof; ++p) {
327         for (d = 0; d < maxDof; ++d) {
328           if (d > 0) {
329             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
330           }
331           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);CHKERRQ(ierr);
332         }
333         for (d = maxDof; d < enforceDof; ++d) {
334           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
335         }
336         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
337       }
338       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
339     }
340   } else {
341     PetscScalar *localValues;
342     PetscInt    size, k = 0;
343 
344     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
345     ierr = PetscMalloc(size * sizeof(PetscScalar), &localValues);CHKERRQ(ierr);
346     for (p = pStart; p < pEnd; ++p) {
347       PetscInt dof, off, goff, d;
348 
349       /* Reject points not either cells or vertices */
350       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
351       if (hasLabel) {
352         PetscInt value;
353 
354         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
355             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
356           ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
357           if (value != 1) continue;
358         }
359       }
360       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
361       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
362       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
363       if (goff >= 0) {
364         for (d = 0; d < dof; ++d) {
365           localValues[k++] = array[off+d];
366         }
367       }
368     }
369     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
370     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRQ(ierr);
371     ierr = PetscFree(localValues);CHKERRQ(ierr);
372   }
373   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
374   PetscFunctionReturn(0);
375 }
376 
377 #undef __FUNCT__
378 #define __FUNCT__ "DMPlexVTKWriteField_ASCII"
379 PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
380 {
381   MPI_Comm       comm   = ((PetscObject) dm)->comm;
382   PetscInt       numDof = 0, maxDof;
383   PetscInt       pStart, pEnd, p;
384   PetscErrorCode ierr;
385 
386   PetscFunctionBegin;
387   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
388   for (p = pStart; p < pEnd; ++p) {
389     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
390     if (numDof) break;
391   }
392   numDof = PetscMax(numDof, enforceDof);
393   ierr = MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, ((PetscObject) dm)->comm);CHKERRQ(ierr);
394   if (!name) name = "Unknown";
395   if (maxDof == 3) {
396     ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
397   } else {
398     ierr = PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);CHKERRQ(ierr);
399     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
400   }
401   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);CHKERRQ(ierr);
402   PetscFunctionReturn(0);
403 }
404 
405 #undef __FUNCT__
406 #define __FUNCT__ "DMPlexVTKWriteAll_ASCII"
407 static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
408 {
409   MPI_Comm                 comm = ((PetscObject) dm)->comm;
410   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
411   FILE                     *fp;
412   PetscViewerVTKObjectLink link;
413   PetscSection             coordSection, globalCoordSection;
414   PetscLayout              vLayout;
415   Vec                      coordinates;
416   PetscReal                lengthScale;
417   PetscInt                 vMax, totVertices, totCells;
418   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE;
419   PetscErrorCode           ierr;
420 
421   PetscFunctionBegin;
422   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
423   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
424   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
425   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
426   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
427   /* Vertices */
428   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
429   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
430   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
431   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
432   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
433   if (vMax >= 0) {
434     PetscInt pStart, pEnd, p, localSize = 0;
435 
436     ierr = PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);CHKERRQ(ierr);
437     pEnd = PetscMin(pEnd, vMax);
438     for (p = pStart; p < pEnd; ++p) {
439       PetscInt dof;
440 
441       ierr = PetscSectionGetDof(globalCoordSection, p, &dof);CHKERRQ(ierr);
442       if (dof > 0) ++localSize;
443     }
444     ierr = PetscLayoutCreate(((PetscObject) dm)->comm, &vLayout);CHKERRQ(ierr);
445     ierr = PetscLayoutSetLocalSize(vLayout, localSize);CHKERRQ(ierr);
446     ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
447     ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
448   } else {
449     ierr = PetscSectionGetPointLayout(((PetscObject) dm)->comm, globalCoordSection, &vLayout);CHKERRQ(ierr);
450   }
451   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
452   ierr = PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);CHKERRQ(ierr);
453   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);CHKERRQ(ierr);
454   /* Cells */
455   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
456   /* Vertex fields */
457   for (link = vtk->link; link; link = link->next) {
458     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
459     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
460   }
461   if (hasPoint) {
462     ierr = PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);CHKERRQ(ierr);
463     for (link = vtk->link; link; link = link->next) {
464       Vec          X = (Vec) link->vec;
465       DM           dmX;
466       PetscSection section, globalSection, newSection = NULL;
467       const char   *name;
468       PetscInt     enforceDof = PETSC_DETERMINE;
469 
470       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
471       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
472       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
473       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
474       if (dmX) {
475         DMLabel  subpointMap, subpointMapX;
476         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
477 
478         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
479         /* Here is where we check whether dmX is a submesh of dm */
480         ierr = DMPlexGetDimension(dm,  &dim);CHKERRQ(ierr);
481         ierr = DMPlexGetDimension(dmX, &dimX);CHKERRQ(ierr);
482         ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
483         ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
484         ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
485         ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
486         if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
487           const PetscInt *ind;
488           IS              subpointIS;
489           PetscInt        n, q;
490 
491           ierr = PetscPrintf(PETSC_COMM_SELF, "Making translation PetscSection\n");CHKERRQ(ierr);
492           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
493           ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
494           ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
495           ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
496           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
497           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
498           for (q = qStart; q < qEnd; ++q) {
499             PetscInt dof, off, p;
500 
501             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
502             if (dof) {
503               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
504               if (p >= pStart) {
505                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
506                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
507                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
508               }
509             }
510           }
511           ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
512           ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
513           /* No need to setup section */
514           section = newSection;
515         }
516       } else {
517         PetscContainer c;
518 
519         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
520         if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
521         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
522       }
523       if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
524       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
525       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
526       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
527       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
528     }
529   }
530   /* Cell Fields */
531   if (hasCell) {
532     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
533     for (link = vtk->link; link; link = link->next) {
534       Vec          X = (Vec) link->vec;
535       DM           dmX;
536       PetscSection section, globalSection;
537       const char   *name;
538       PetscInt     enforceDof = PETSC_DETERMINE;
539 
540       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
541       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
542       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
543       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
544       if (dmX) {
545         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
546       } else {
547         PetscContainer c;
548 
549         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
550         if (!c) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
551         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
552       }
553       if (!section) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
554       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
555       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
556       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
557     }
558   }
559   /* Cleanup */
560   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
561   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
562   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
563   PetscFunctionReturn(0);
564 }
565 
566 #undef __FUNCT__
567 #define __FUNCT__ "DMPlexVTKWriteAll"
568 /*@C
569   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
570 
571   Collective
572 
573   Input Arguments:
574 + odm - The DMPlex specifying the mesh, passed as a PetscObject
575 - viewer - viewer of type VTK
576 
577   Level: developer
578 
579   Note:
580   This function is a callback used by the VTK viewer to actually write the file.
581   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
582   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
583 
584 .seealso: PETSCVIEWERVTK
585 @*/
586 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
587 {
588   DM             dm = (DM) odm;
589   PetscBool      isvtk;
590   PetscErrorCode ierr;
591 
592   PetscFunctionBegin;
593   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
594   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
595   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
596   if (!isvtk) SETERRQ1(((PetscObject) viewer)->comm, PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
597   switch (viewer->format) {
598   case PETSC_VIEWER_ASCII_VTK:
599     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
600     break;
601   case PETSC_VIEWER_VTK_VTU:
602     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
603     break;
604   default: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
605   }
606   PetscFunctionReturn(0);
607 }
608