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