xref: /petsc/src/dm/impls/plex/plexvtk.c (revision e630c359805a7909577b2c5515bf56b45bba1d70) 
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 PetscErrorCode DMPlexVTKGetCellType_Internal(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
6 {
7   PetscFunctionBegin;
8   *cellType = -1;
9   switch (dim) {
10   case 0:
11     switch (corners) {
12     case 1:
13       *cellType = 1; /* VTK_VERTEX */
14       break;
15     default:
16       break;
17     }
18     break;
19   case 1:
20     switch (corners) {
21     case 2:
22       *cellType = 3; /* VTK_LINE */
23       break;
24     case 3:
25       *cellType = 21; /* VTK_QUADRATIC_EDGE */
26       break;
27     default:
28       break;
29     }
30     break;
31   case 2:
32     switch (corners) {
33     case 3:
34       *cellType = 5; /* VTK_TRIANGLE */
35       break;
36     case 4:
37       *cellType = 9; /* VTK_QUAD */
38       break;
39     case 6:
40       *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
41       break;
42     case 9:
43       *cellType = 23; /* VTK_QUADRATIC_QUAD */
44       break;
45     default:
46       break;
47     }
48     break;
49   case 3:
50     switch (corners) {
51     case 4:
52       *cellType = 10; /* VTK_TETRA */
53       break;
54     case 6:
55       *cellType = 13; /* VTK_WEDGE */
56       break;
57     case 8:
58       *cellType = 12; /* VTK_HEXAHEDRON */
59       break;
60     case 10:
61       *cellType = 24; /* VTK_QUADRATIC_TETRA */
62       break;
63     case 27:
64       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
65       break;
66     default:
67       break;
68     }
69   }
70   PetscFunctionReturn(0);
71 }
72 
73 static PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
74 {
75   MPI_Comm       comm;
76   DMLabel        label;
77   IS             globalVertexNumbers = NULL;
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, maxLabelCells, cStart, cEnd, c, vStart, vEnd, v;
83   PetscMPIInt    size, rank, proc, tag;
84   PetscErrorCode ierr;
85 
86   PetscFunctionBegin;
87   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
88   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
89   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
90   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
91   ierr = DMGetDimension(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 = DMGetLabel(dm, "vtk", &label);CHKERRQ(ierr);
96   ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
97   ierr = MPIU_Allreduce(&numLabelCells, &maxLabelCells, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
98   if (!maxLabelCells) label = NULL;
99   for (c = cStart; c < cEnd; ++c) {
100     PetscInt *closure = NULL;
101     PetscInt closureSize, value;
102 
103     if (label) {
104       ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
105       if (value != 1) continue;
106     }
107     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
108     for (v = 0; v < closureSize*2; v += 2) {
109       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
110     }
111     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
112     ++numCells;
113   }
114   maxCells = numCells;
115   ierr     = MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
116   ierr     = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
117   ierr     = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
118   ierr     = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
119   ierr     = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
120   ierr     = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
121   ierr     = PetscMalloc1(maxCells, &corners);CHKERRQ(ierr);
122   ierr     = PetscFPrintf(comm, fp, "CELLS %D %D\n", totCells, totCorners+totCells);CHKERRQ(ierr);
123   if (!rank) {
124     PetscInt *remoteVertices;
125     int      *vertices;
126 
127     ierr = PetscMalloc1(maxCorners, &vertices);CHKERRQ(ierr);
128     for (c = cStart, numCells = 0; c < cEnd; ++c) {
129       PetscInt *closure = NULL;
130       PetscInt closureSize, value, nC = 0;
131 
132       if (label) {
133         ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
134         if (value != 1) continue;
135       }
136       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
137       for (v = 0; v < closureSize*2; v += 2) {
138         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
139           const PetscInt gv = gvertex[closure[v] - vStart];
140           vertices[nC++] = gv < 0 ? -(gv+1) : gv;
141         }
142       }
143       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
144       corners[numCells++] = nC;
145       ierr = PetscFPrintf(comm, fp, "%D ", nC);CHKERRQ(ierr);
146       ierr = DMPlexInvertCell(dim, nC, vertices);CHKERRQ(ierr);
147       for (v = 0; v < nC; ++v) {
148         ierr = PetscFPrintf(comm, fp, " %d", vertices[v]);CHKERRQ(ierr);
149       }
150       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
151     }
152     if (size > 1) {ierr = PetscMalloc1(maxCorners+maxCells, &remoteVertices);CHKERRQ(ierr);}
153     for (proc = 1; proc < size; ++proc) {
154       MPI_Status status;
155 
156       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
157       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
158       for (c = 0; c < numCorners;) {
159         PetscInt nC = remoteVertices[c++];
160 
161         for (v = 0; v < nC; ++v, ++c) {
162           vertices[v] = remoteVertices[c];
163         }
164         ierr = DMPlexInvertCell(dim, nC, vertices);CHKERRQ(ierr);
165         ierr = PetscFPrintf(comm, fp, "%D ", nC);CHKERRQ(ierr);
166         for (v = 0; v < nC; ++v) {
167           ierr = PetscFPrintf(comm, fp, " %d", vertices[v]);CHKERRQ(ierr);
168         }
169         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
170       }
171     }
172     if (size > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
173     ierr = PetscFree(vertices);CHKERRQ(ierr);
174   } else {
175     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
176 
177     ierr = PetscMalloc1(numSend, &localVertices);CHKERRQ(ierr);
178     for (c = cStart, numCells = 0; c < cEnd; ++c) {
179       PetscInt *closure = NULL;
180       PetscInt closureSize, value, nC = 0;
181 
182       if (label) {
183         ierr = DMLabelGetValue(label, 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_Internal(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
212       ierr = PetscFPrintf(comm, fp, "%D\n", cellType);CHKERRQ(ierr);
213     }
214     for (proc = 1; proc < size; ++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_Internal(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 static PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
234 {
235   MPI_Comm       comm;
236   PetscInt       numCells = 0, cellHeight;
237   PetscInt       numLabelCells, cStart, cEnd, c;
238   PetscMPIInt    size, rank, proc, tag;
239   PetscBool      hasLabel;
240   PetscErrorCode ierr;
241 
242   PetscFunctionBegin;
243   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
244   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
245   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
246   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
247   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
248   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
249   ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
250   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
251   for (c = cStart; c < cEnd; ++c) {
252     if (hasLabel) {
253       PetscInt value;
254 
255       ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
256       if (value != 1) continue;
257     }
258     ++numCells;
259   }
260   if (!rank) {
261     for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", rank);CHKERRQ(ierr);}
262     for (proc = 1; proc < size; ++proc) {
263       MPI_Status status;
264 
265       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
266       for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", proc);CHKERRQ(ierr);}
267     }
268   } else {
269     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
270   }
271   PetscFunctionReturn(0);
272 }
273 
274 static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
275 {
276   MPI_Comm           comm;
277   const MPI_Datatype mpiType = MPIU_SCALAR;
278   PetscScalar        *array;
279   PetscInt           numDof = 0, maxDof;
280   PetscInt           numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
281   PetscMPIInt        size, rank, proc, tag;
282   PetscBool          hasLabel;
283   PetscErrorCode     ierr;
284 
285   PetscFunctionBegin;
286   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
287   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
288   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
289   if (precision < 0) precision = 6;
290   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
291   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
292   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
293   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
294   /* VTK only wants the values at cells or vertices */
295   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
296   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
297   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
298   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
299   if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
300   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
301   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
302   ierr     = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
303   ierr     = DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
304   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
305   for (p = pStart; p < pEnd; ++p) {
306     /* Reject points not either cells or vertices */
307     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
308     if (hasLabel) {
309       PetscInt value;
310 
311       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
312           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
313         ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
314         if (value != 1) continue;
315       }
316     }
317     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
318     if (numDof) break;
319   }
320   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
321   enforceDof = PetscMax(enforceDof, maxDof);
322   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
323   if (!rank) {
324 #if defined(PETSC_USE_REAL___FLOAT128)
325     double dval;
326 #endif
327     char formatString[8];
328 
329     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
330     for (p = pStart; p < pEnd; ++p) {
331       /* Here we lose a way to filter points by keeping them out of the Numbering */
332       PetscInt dof, off, goff, d;
333 
334       /* Reject points not either cells or vertices */
335       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
336       if (hasLabel) {
337         PetscInt value;
338 
339         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
340             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
341           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
342           if (value != 1) continue;
343         }
344       }
345       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
346       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
347       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
348       if (dof && goff >= 0) {
349         for (d = 0; d < dof; d++) {
350           if (d > 0) {
351             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
352           }
353 #if defined(PETSC_USE_REAL___FLOAT128)
354           dval = (double)PetscRealPart(array[off+d])*scale;
355           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
356 #else
357           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);CHKERRQ(ierr);
358 #endif
359         }
360         for (d = dof; d < enforceDof; d++) {
361           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
362         }
363         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
364       }
365     }
366     for (proc = 1; proc < size; ++proc) {
367       PetscScalar *remoteValues;
368       PetscInt    size = 0, d;
369       MPI_Status  status;
370 
371       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
372       ierr = PetscMalloc1(size, &remoteValues);CHKERRQ(ierr);
373       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRQ(ierr);
374       for (p = 0; p < size/maxDof; ++p) {
375         for (d = 0; d < maxDof; ++d) {
376           if (d > 0) {
377             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
378           }
379 #if defined(PETSC_USE_REAL___FLOAT128)
380           dval = PetscRealPart(remoteValues[p*maxDof+d])*scale;
381           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
382 #else
383           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);CHKERRQ(ierr);
384 #endif
385         }
386         for (d = maxDof; d < enforceDof; ++d) {
387           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
388         }
389         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
390       }
391       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
392     }
393   } else {
394     PetscScalar *localValues;
395     PetscInt    size, k = 0;
396 
397     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
398     ierr = PetscMalloc1(size, &localValues);CHKERRQ(ierr);
399     for (p = pStart; p < pEnd; ++p) {
400       PetscInt dof, off, goff, d;
401 
402       /* Reject points not either cells or vertices */
403       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
404       if (hasLabel) {
405         PetscInt value;
406 
407         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
408             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
409           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
410           if (value != 1) continue;
411         }
412       }
413       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
414       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
415       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
416       if (goff >= 0) {
417         for (d = 0; d < dof; ++d) {
418           localValues[k++] = array[off+d];
419         }
420       }
421     }
422     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
423     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRQ(ierr);
424     ierr = PetscFree(localValues);CHKERRQ(ierr);
425   }
426   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
427   PetscFunctionReturn(0);
428 }
429 
430 static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
431 {
432   MPI_Comm       comm;
433   PetscInt       numDof = 0, maxDof;
434   PetscInt       pStart, pEnd, p;
435   PetscErrorCode ierr;
436 
437   PetscFunctionBegin;
438   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
439   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
440   for (p = pStart; p < pEnd; ++p) {
441     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
442     if (numDof) break;
443   }
444   numDof = PetscMax(numDof, enforceDof);
445   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
446   if (!name) name = "Unknown";
447   if (maxDof == 3) {
448     ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
449   } else {
450     ierr = PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);CHKERRQ(ierr);
451     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
452   }
453   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);CHKERRQ(ierr);
454   PetscFunctionReturn(0);
455 }
456 
457 static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
458 {
459   MPI_Comm                 comm;
460   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
461   FILE                     *fp;
462   PetscViewerVTKObjectLink link;
463   PetscSection             coordSection, globalCoordSection;
464   PetscLayout              vLayout;
465   Vec                      coordinates;
466   PetscReal                lengthScale;
467   PetscInt                 vMax, totVertices, totCells = 0;
468   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized;
469   PetscErrorCode           ierr;
470 
471   PetscFunctionBegin;
472   ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr);
473   if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
474   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
475   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
476   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
477   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
478   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
479   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
480   /* Vertices */
481   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
482   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
483   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
484   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
485   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
486   if (vMax >= 0) {
487     PetscInt pStart, pEnd, p, localSize = 0;
488 
489     ierr = PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);CHKERRQ(ierr);
490     pEnd = PetscMin(pEnd, vMax);
491     for (p = pStart; p < pEnd; ++p) {
492       PetscInt dof;
493 
494       ierr = PetscSectionGetDof(globalCoordSection, p, &dof);CHKERRQ(ierr);
495       if (dof > 0) ++localSize;
496     }
497     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);CHKERRQ(ierr);
498     ierr = PetscLayoutSetLocalSize(vLayout, localSize);CHKERRQ(ierr);
499     ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
500     ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
501   } else {
502     ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);CHKERRQ(ierr);
503   }
504   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
505   ierr = PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);CHKERRQ(ierr);
506   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);CHKERRQ(ierr);
507   /* Cells */
508   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
509   /* Vertex fields */
510   for (link = vtk->link; link; link = link->next) {
511     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
512     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
513   }
514   if (hasPoint) {
515     ierr = PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);CHKERRQ(ierr);
516     for (link = vtk->link; link; link = link->next) {
517       Vec          X = (Vec) link->vec;
518       PetscSection section = NULL, globalSection, newSection = NULL;
519       char         namebuf[256];
520       const char   *name;
521       PetscInt     enforceDof = PETSC_DETERMINE;
522 
523       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
524       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
525       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
526       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
527       if (!section) {
528         DM           dmX;
529 
530         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
531         if (dmX) {
532           DMLabel  subpointMap, subpointMapX;
533           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
534 
535           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
536           /* Here is where we check whether dmX is a submesh of dm */
537           ierr = DMGetDimension(dm,  &dim);CHKERRQ(ierr);
538           ierr = DMGetDimension(dmX, &dimX);CHKERRQ(ierr);
539           ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
540           ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
541           ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
542           ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
543           if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
544             const PetscInt *ind = NULL;
545             IS              subpointIS;
546             PetscInt        n = 0, q;
547 
548             ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
549             ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
550             if (subpointIS) {
551               ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
552               ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
553             }
554             ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
555             ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
556             for (q = qStart; q < qEnd; ++q) {
557               PetscInt dof, off, p;
558 
559               ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
560               if (dof) {
561                 ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
562                 if (p >= pStart) {
563                   ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
564                   ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
565                   ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
566                 }
567               }
568             }
569             if (subpointIS) {
570               ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
571               ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
572             }
573             /* No need to setup section */
574             section = newSection;
575           }
576         }
577       }
578       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
579       if (link->field >= 0) {
580         const char *fieldname;
581 
582         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
583         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
584         if (fieldname) {
585           ierr = PetscSNPrintf(namebuf, 255, "%s%s", name, fieldname);CHKERRQ(ierr);
586         } else {
587           ierr = PetscSNPrintf(namebuf, 255, "%s%D", name, link->field);CHKERRQ(ierr);
588         }
589         name = &namebuf[0];
590       }
591       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
592       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
593       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
594       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
595     }
596   }
597   /* Cell Fields */
598   ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
599   if (hasCell || writePartition) {
600     ierr = PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);CHKERRQ(ierr);
601     for (link = vtk->link; link; link = link->next) {
602       Vec          X = (Vec) link->vec;
603       PetscSection section = NULL, globalSection;
604       const char   *name;
605       char         namebuf[256];
606       PetscInt     enforceDof = PETSC_DETERMINE;
607 
608       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
609       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
610       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
611       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
612       if (!section) {
613         DM           dmX;
614 
615         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
616         if (dmX) {
617           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
618         }
619       }
620       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
621       if (link->field >= 0) {
622         const char *fieldname;
623 
624         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
625         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
626         if (fieldname) {
627           ierr = PetscSNPrintf(namebuf, 255, "%s%s", name, fieldname);CHKERRQ(ierr);
628         } else {
629           ierr = PetscSNPrintf(namebuf, 255, "%s%D", name, link->field);CHKERRQ(ierr);
630         }
631         name = &namebuf[0];
632       }
633       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
634       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
635       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
636     }
637     if (writePartition) {
638       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
639       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
640       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
641     }
642   }
643   /* Cleanup */
644   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
645   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
646   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
647   PetscFunctionReturn(0);
648 }
649 
650 /*@C
651   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
652 
653   Collective
654 
655   Input Arguments:
656 + odm - The DMPlex specifying the mesh, passed as a PetscObject
657 - viewer - viewer of type VTK
658 
659   Level: developer
660 
661   Note:
662   This function is a callback used by the VTK viewer to actually write the file.
663   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
664   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
665 
666 .seealso: PETSCVIEWERVTK
667 @*/
668 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
669 {
670   DM             dm = (DM) odm;
671   PetscBool      isvtk;
672   PetscErrorCode ierr;
673 
674   PetscFunctionBegin;
675   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
676   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
677   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
678   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
679   switch (viewer->format) {
680   case PETSC_VIEWER_ASCII_VTK:
681     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
682     break;
683   case PETSC_VIEWER_VTK_VTU:
684     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
685     break;
686   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
687   }
688   PetscFunctionReturn(0);
689 }
690