xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 6ac5842e34eedc6428162d8d42bedaaf46eae34c) 
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 #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         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
522         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
523       }
524       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
525       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
526       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
527       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
528       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
529     }
530   }
531   /* Cell Fields */
532   if (hasCell) {
533     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
534     for (link = vtk->link; link; link = link->next) {
535       Vec          X = (Vec) link->vec;
536       DM           dmX;
537       PetscSection section, globalSection;
538       const char   *name;
539       PetscInt     enforceDof = PETSC_DETERMINE;
540 
541       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
542       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
543       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
544       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
545       if (dmX) {
546         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
547       } else {
548         PetscContainer c;
549 
550         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
551         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
552         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
553       }
554       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
555       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
556       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
557       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
558     }
559   }
560   /* Cleanup */
561   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
562   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
563   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
564   PetscFunctionReturn(0);
565 }
566 
567 #undef __FUNCT__
568 #define __FUNCT__ "DMPlexVTKWriteAll"
569 /*@C
570   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
571 
572   Collective
573 
574   Input Arguments:
575 + odm - The DMPlex specifying the mesh, passed as a PetscObject
576 - viewer - viewer of type VTK
577 
578   Level: developer
579 
580   Note:
581   This function is a callback used by the VTK viewer to actually write the file.
582   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
583   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
584 
585 .seealso: PETSCVIEWERVTK
586 @*/
587 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
588 {
589   DM             dm = (DM) odm;
590   PetscBool      isvtk;
591   PetscErrorCode ierr;
592 
593   PetscFunctionBegin;
594   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
595   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
596   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
597   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
598   switch (viewer->format) {
599   case PETSC_VIEWER_ASCII_VTK:
600     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
601     break;
602   case PETSC_VIEWER_VTK_VTU:
603     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
604     break;
605   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
606   }
607   PetscFunctionReturn(0);
608 }
609