xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 95c0884e6f7665b705eebf88174e89dc920c2fc0) 
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       DM           dmX;
519       PetscSection section, globalSection, newSection = NULL;
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 = VecGetDM(X, &dmX);CHKERRQ(ierr);
527       if (dmX) {
528         DMLabel  subpointMap, subpointMapX;
529         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
530 
531         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
532         /* Here is where we check whether dmX is a submesh of dm */
533         ierr = DMGetDimension(dm,  &dim);CHKERRQ(ierr);
534         ierr = DMGetDimension(dmX, &dimX);CHKERRQ(ierr);
535         ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
536         ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
537         ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
538         ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
539         if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
540           const PetscInt *ind = NULL;
541           IS              subpointIS;
542           PetscInt        n = 0, q;
543 
544           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
545           ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
546           if (subpointIS) {
547             ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
548             ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
549           }
550           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
551           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
552           for (q = qStart; q < qEnd; ++q) {
553             PetscInt dof, off, p;
554 
555             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
556             if (dof) {
557               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
558               if (p >= pStart) {
559                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
560                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
561                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
562               }
563             }
564           }
565           if (subpointIS) {
566             ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
567             ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
568           }
569           /* No need to setup section */
570           section = newSection;
571         }
572       } else {
573         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
574         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
575       }
576       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
577       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
578       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
579       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
580       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
581     }
582   }
583   /* Cell Fields */
584   ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
585   if (hasCell || writePartition) {
586     ierr = PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);CHKERRQ(ierr);
587     for (link = vtk->link; link; link = link->next) {
588       Vec          X = (Vec) link->vec;
589       DM           dmX;
590       PetscSection section, globalSection;
591       const char   *name;
592       PetscInt     enforceDof = PETSC_DETERMINE;
593 
594       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
595       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
596       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
597       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
598       if (dmX) {
599         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
600       } else {
601         PetscContainer c;
602 
603         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
604         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
605         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
606       }
607       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
608       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
609       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
610       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
611     }
612     if (writePartition) {
613       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
614       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
615       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
616     }
617   }
618   /* Cleanup */
619   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
620   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
621   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
622   PetscFunctionReturn(0);
623 }
624 
625 /*@C
626   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
627 
628   Collective
629 
630   Input Arguments:
631 + odm - The DMPlex specifying the mesh, passed as a PetscObject
632 - viewer - viewer of type VTK
633 
634   Level: developer
635 
636   Note:
637   This function is a callback used by the VTK viewer to actually write the file.
638   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
639   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
640 
641 .seealso: PETSCVIEWERVTK
642 @*/
643 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
644 {
645   DM             dm = (DM) odm;
646   PetscBool      isvtk;
647   PetscErrorCode ierr;
648 
649   PetscFunctionBegin;
650   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
651   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
652   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
653   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
654   switch (viewer->format) {
655   case PETSC_VIEWER_ASCII_VTK:
656     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
657     break;
658   case PETSC_VIEWER_VTK_VTU:
659     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
660     break;
661   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
662   }
663   PetscFunctionReturn(0);
664 }
665