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