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