xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 3c633725528e547aaaa9b672a746f5d686a276e1) 
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);CHKERRMPI(ierr);
90   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(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);CHKERRMPI(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);CHKERRMPI(ierr);
116   ierr     = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRMPI(ierr);
117   ierr     = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRMPI(ierr);
118   ierr     = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRMPI(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 == 0) {
124     PetscInt *remoteVertices, *vertices;
125 
126     ierr = PetscMalloc1(maxCorners, &vertices);CHKERRQ(ierr);
127     for (c = cStart, numCells = 0; c < cEnd; ++c) {
128       PetscInt *closure = NULL;
129       PetscInt closureSize, value, nC = 0;
130 
131       if (label) {
132         ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
133         if (value != 1) continue;
134       }
135       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
136       for (v = 0; v < closureSize*2; v += 2) {
137         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
138           const PetscInt gv = gvertex[closure[v] - vStart];
139           vertices[nC++] = gv < 0 ? -(gv+1) : gv;
140         }
141       }
142       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
143       ierr = DMPlexReorderCell(dm, c, vertices);CHKERRQ(ierr);
144       corners[numCells++] = nC;
145       ierr = PetscFPrintf(comm, fp, "%D ", nC);CHKERRQ(ierr);
146       for (v = 0; v < nC; ++v) {
147         ierr = PetscFPrintf(comm, fp, " %D", vertices[v]);CHKERRQ(ierr);
148       }
149       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
150     }
151     if (size > 1) {ierr = PetscMalloc1(maxCorners+maxCells, &remoteVertices);CHKERRQ(ierr);}
152     for (proc = 1; proc < size; ++proc) {
153       MPI_Status status;
154 
155       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRMPI(ierr);
156       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRMPI(ierr);
157       for (c = 0; c < numCorners;) {
158         PetscInt nC = remoteVertices[c++];
159 
160         for (v = 0; v < nC; ++v, ++c) {
161           vertices[v] = remoteVertices[c];
162         }
163         ierr = PetscFPrintf(comm, fp, "%D ", nC);CHKERRQ(ierr);
164         for (v = 0; v < nC; ++v) {
165           ierr = PetscFPrintf(comm, fp, " %D", vertices[v]);CHKERRQ(ierr);
166         }
167         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
168       }
169     }
170     if (size > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
171     ierr = PetscFree(vertices);CHKERRQ(ierr);
172   } else {
173     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
174 
175     ierr = PetscMalloc1(numSend, &localVertices);CHKERRQ(ierr);
176     for (c = cStart, numCells = 0; c < cEnd; ++c) {
177       PetscInt *closure = NULL;
178       PetscInt closureSize, value, nC = 0;
179 
180       if (label) {
181         ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
182         if (value != 1) continue;
183       }
184       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
185       for (v = 0; v < closureSize*2; v += 2) {
186         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
187           const PetscInt gv = gvertex[closure[v] - vStart];
188           closure[nC++] = gv < 0 ? -(gv+1) : gv;
189         }
190       }
191       corners[numCells++] = nC;
192       localVertices[k++]  = nC;
193       for (v = 0; v < nC; ++v, ++k) {
194         localVertices[k] = closure[v];
195       }
196       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
197       ierr = DMPlexReorderCell(dm, c, localVertices+k-nC);CHKERRQ(ierr);
198     }
199     PetscCheckFalse(k != numSend,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %D should be %D", k, numSend);
200     ierr = MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);CHKERRMPI(ierr);
201     ierr = MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);CHKERRMPI(ierr);
202     ierr = PetscFree(localVertices);CHKERRQ(ierr);
203   }
204   ierr = ISRestoreIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
205   ierr = PetscFPrintf(comm, fp, "CELL_TYPES %D\n", totCells);CHKERRQ(ierr);
206   if (rank == 0) {
207     PetscInt cellType;
208 
209     for (c = 0; c < numCells; ++c) {
210       ierr = DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
211       ierr = PetscFPrintf(comm, fp, "%D\n", cellType);CHKERRQ(ierr);
212     }
213     for (proc = 1; proc < size; ++proc) {
214       MPI_Status status;
215 
216       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRMPI(ierr);
217       ierr = MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);CHKERRMPI(ierr);
218       for (c = 0; c < numCells; ++c) {
219         ierr = DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
220         ierr = PetscFPrintf(comm, fp, "%D\n", cellType);CHKERRQ(ierr);
221       }
222     }
223   } else {
224     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRMPI(ierr);
225     ierr = MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);CHKERRMPI(ierr);
226   }
227   ierr        = PetscFree(corners);CHKERRQ(ierr);
228   *totalCells = totCells;
229   PetscFunctionReturn(0);
230 }
231 
232 static PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
233 {
234   MPI_Comm       comm;
235   PetscInt       numCells = 0, cellHeight;
236   PetscInt       numLabelCells, cStart, cEnd, c;
237   PetscMPIInt    size, rank, proc, tag;
238   PetscBool      hasLabel;
239   PetscErrorCode ierr;
240 
241   PetscFunctionBegin;
242   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
243   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
244   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
245   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
246   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
247   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
248   ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
249   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
250   for (c = cStart; c < cEnd; ++c) {
251     if (hasLabel) {
252       PetscInt value;
253 
254       ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
255       if (value != 1) continue;
256     }
257     ++numCells;
258   }
259   if (rank == 0) {
260     for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", rank);CHKERRQ(ierr);}
261     for (proc = 1; proc < size; ++proc) {
262       MPI_Status status;
263 
264       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRMPI(ierr);
265       for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", proc);CHKERRQ(ierr);}
266     }
267   } else {
268     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRMPI(ierr);
269   }
270   PetscFunctionReturn(0);
271 }
272 
273 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
274 typedef double PetscVTKReal;
275 #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
276 typedef float PetscVTKReal;
277 #else
278 typedef PetscReal PetscVTKReal;
279 #endif
280 
281 static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
282 {
283   MPI_Comm           comm;
284   const MPI_Datatype mpiType = MPIU_SCALAR;
285   PetscScalar        *array;
286   PetscInt           numDof = 0, maxDof;
287   PetscInt           numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vStart, vEnd, pStart, pEnd, p;
288   PetscMPIInt        size, rank, proc, tag;
289   PetscBool          hasLabel;
290   PetscErrorCode     ierr;
291 
292   PetscFunctionBegin;
293   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
294   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
295   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
296   if (precision < 0) precision = 6;
297   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
298   ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr);
299   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
300   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
301   /* VTK only wants the values at cells or vertices */
302   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
303   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
304   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
305   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
306   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
307   ierr     = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
308   ierr     = DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
309   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
310   for (p = pStart; p < pEnd; ++p) {
311     /* Reject points not either cells or vertices */
312     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
313     if (hasLabel) {
314       PetscInt value;
315 
316       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
317           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
318         ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
319         if (value != 1) continue;
320       }
321     }
322     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
323     if (numDof) break;
324   }
325   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRMPI(ierr);
326   enforceDof = PetscMax(enforceDof, maxDof);
327   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
328   if (rank == 0) {
329     PetscVTKReal dval;
330     PetscScalar  val;
331     char formatString[8];
332 
333     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
334     for (p = pStart; p < pEnd; ++p) {
335       /* Here we lose a way to filter points by keeping them out of the Numbering */
336       PetscInt dof, off, goff, d;
337 
338       /* Reject points not either cells or vertices */
339       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
340       if (hasLabel) {
341         PetscInt value;
342 
343         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
344             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
345           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
346           if (value != 1) continue;
347         }
348       }
349       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
350       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
351       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
352       if (dof && goff >= 0) {
353         for (d = 0; d < dof; d++) {
354           if (d > 0) {
355             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
356           }
357           val = array[off+d];
358           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
359           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
360         }
361         for (d = dof; d < enforceDof; d++) {
362           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
363         }
364         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
365       }
366     }
367     for (proc = 1; proc < size; ++proc) {
368       PetscScalar *remoteValues;
369       PetscInt    size = 0, d;
370       MPI_Status  status;
371 
372       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRMPI(ierr);
373       ierr = PetscMalloc1(size, &remoteValues);CHKERRQ(ierr);
374       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRMPI(ierr);
375       for (p = 0; p < size/maxDof; ++p) {
376         for (d = 0; d < maxDof; ++d) {
377           if (d > 0) {
378             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
379           }
380           val = remoteValues[p*maxDof+d];
381           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
382           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
383         }
384         for (d = maxDof; d < enforceDof; ++d) {
385           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
386         }
387         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
388       }
389       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
390     }
391   } else {
392     PetscScalar *localValues;
393     PetscInt    size, k = 0;
394 
395     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
396     ierr = PetscMalloc1(size, &localValues);CHKERRQ(ierr);
397     for (p = pStart; p < pEnd; ++p) {
398       PetscInt dof, off, goff, d;
399 
400       /* Reject points not either cells or vertices */
401       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
402       if (hasLabel) {
403         PetscInt value;
404 
405         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
406             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
407           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
408           if (value != 1) continue;
409         }
410       }
411       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
412       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
413       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
414       if (goff >= 0) {
415         for (d = 0; d < dof; ++d) {
416           localValues[k++] = array[off+d];
417         }
418       }
419     }
420     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRMPI(ierr);
421     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRMPI(ierr);
422     ierr = PetscFree(localValues);CHKERRQ(ierr);
423   }
424   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
425   PetscFunctionReturn(0);
426 }
427 
428 static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscBool nameComplex, PetscInt imag)
429 {
430   MPI_Comm       comm;
431   PetscInt       numDof = 0, maxDof;
432   PetscInt       pStart, pEnd, p;
433   PetscErrorCode ierr;
434 
435   PetscFunctionBegin;
436   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
437   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
438   for (p = pStart; p < pEnd; ++p) {
439     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
440     if (numDof) break;
441   }
442   numDof = PetscMax(numDof, enforceDof);
443   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr);
444   if (!name) name = "Unknown";
445   if (maxDof == 3) {
446     if (nameComplex) {
447       ierr = PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");CHKERRQ(ierr);
448     } else {
449       ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
450     }
451   } else {
452     if (nameComplex) {
453       ierr = PetscFPrintf(comm, fp, "SCALARS %s.%s double %D\n", name, imag ? "Im" : "Re", maxDof);CHKERRQ(ierr);
454     } else {
455       ierr = PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);CHKERRQ(ierr);
456     }
457     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
458   }
459   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);CHKERRQ(ierr);
460   PetscFunctionReturn(0);
461 }
462 
463 static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
464 {
465   MPI_Comm                 comm;
466   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
467   FILE                     *fp;
468   PetscViewerVTKObjectLink link;
469   PetscInt                 totVertices, totCells = 0, loops_per_scalar, l;
470   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
471   const char               *dmname;
472   PetscErrorCode           ierr;
473 
474   PetscFunctionBegin;
475 #if defined(PETSC_USE_COMPLEX)
476   loops_per_scalar = 2;
477   writeComplex = PETSC_TRUE;
478 #else
479   loops_per_scalar = 1;
480   writeComplex = PETSC_FALSE;
481 #endif
482   ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr);
483   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
484   PetscCheckFalse(localized,comm,PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
485   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
486   ierr = PetscObjectGetName((PetscObject)dm, &dmname);CHKERRQ(ierr);
487   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
488   ierr = PetscFPrintf(comm, fp, "%s\n", dmname);CHKERRQ(ierr);
489   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
490   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
491   /* Vertices */
492   {
493     PetscSection  section, coordSection, globalCoordSection;
494     Vec           coordinates;
495     PetscReal     lengthScale;
496     DMLabel       label;
497     IS            vStratumIS;
498     PetscLayout   vLayout;
499 
500     ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
501     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
502     ierr = DMPlexGetDepthLabel(dm, &label);CHKERRQ(ierr);
503     ierr = DMLabelGetStratumIS(label, 0, &vStratumIS);CHKERRQ(ierr);
504     ierr = DMGetCoordinateSection(dm, &section);CHKERRQ(ierr);                                 /* This section includes all points */
505     ierr = PetscSectionCreateSubmeshSection(section, vStratumIS, &coordSection);CHKERRQ(ierr); /* This one includes just vertices */
506     ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
507     ierr = PetscSectionGetPointLayout(comm, globalCoordSection, &vLayout);CHKERRQ(ierr);
508     ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
509     ierr = PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);CHKERRQ(ierr);
510     ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);CHKERRQ(ierr);
511     ierr = ISDestroy(&vStratumIS);CHKERRQ(ierr);
512     ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
513     ierr = PetscSectionDestroy(&coordSection);CHKERRQ(ierr);
514     ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
515   }
516   /* Cells */
517   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
518   /* Vertex fields */
519   for (link = vtk->link; link; link = link->next) {
520     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
521     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
522   }
523   if (hasPoint) {
524     ierr = PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);CHKERRQ(ierr);
525     for (link = vtk->link; link; link = link->next) {
526       Vec          X = (Vec) link->vec;
527       PetscSection section = NULL, globalSection, newSection = NULL;
528       char         namebuf[256];
529       const char   *name;
530       PetscInt     enforceDof = PETSC_DETERMINE;
531 
532       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
533       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
534       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
535       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
536       if (!section) {
537         DM           dmX;
538 
539         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
540         if (dmX) {
541           DMLabel  subpointMap, subpointMapX;
542           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
543 
544           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
545           /* Here is where we check whether dmX is a submesh of dm */
546           ierr = DMGetDimension(dm,  &dim);CHKERRQ(ierr);
547           ierr = DMGetDimension(dmX, &dimX);CHKERRQ(ierr);
548           ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
549           ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
550           ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
551           ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
552           if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
553             const PetscInt *ind = NULL;
554             IS              subpointIS;
555             PetscInt        n = 0, q;
556 
557             ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
558             ierr = DMPlexGetSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
559             if (subpointIS) {
560               ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
561               ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
562             }
563             ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
564             ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
565             for (q = qStart; q < qEnd; ++q) {
566               PetscInt dof, off, p;
567 
568               ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
569               if (dof) {
570                 ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
571                 if (p >= pStart) {
572                   ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
573                   ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
574                   ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
575                 }
576               }
577             }
578             if (subpointIS) {
579               ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
580             }
581             /* No need to setup section */
582             section = newSection;
583           }
584         }
585       }
586       PetscCheckFalse(!section,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
587       if (link->field >= 0) {
588         const char *fieldname;
589 
590         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
591         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
592         if (fieldname) {
593           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);CHKERRQ(ierr);
594         } else {
595           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);CHKERRQ(ierr);
596         }
597       } else {
598         ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);CHKERRQ(ierr);
599       }
600       ierr = PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));CHKERRQ(ierr);
601       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
602       for (l = 0; l < loops_per_scalar; l++) {
603         ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);CHKERRQ(ierr);
604       }
605       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
606       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
607     }
608   }
609   /* Cell Fields */
610   ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
611   if (hasCell || writePartition) {
612     ierr = PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);CHKERRQ(ierr);
613     for (link = vtk->link; link; link = link->next) {
614       Vec          X = (Vec) link->vec;
615       PetscSection section = NULL, globalSection;
616       const char   *name = "";
617       char         namebuf[256];
618       PetscInt     enforceDof = PETSC_DETERMINE;
619 
620       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
621       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
622       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
623       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
624       if (!section) {
625         DM           dmX;
626 
627         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
628         if (dmX) {
629           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
630         }
631       }
632       PetscCheckFalse(!section,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
633       if (link->field >= 0) {
634         const char *fieldname;
635 
636         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
637         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
638         if (fieldname) {
639           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);CHKERRQ(ierr);
640         } else {
641           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);CHKERRQ(ierr);
642         }
643       } else {
644         ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);CHKERRQ(ierr);
645       }
646       ierr = PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));CHKERRQ(ierr);
647       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
648       for (l = 0; l < loops_per_scalar; l++) {
649         ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);CHKERRQ(ierr);
650       }
651       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
652     }
653     if (writePartition) {
654       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
655       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
656       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
657     }
658   }
659   /* Cleanup */
660   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
661   PetscFunctionReturn(0);
662 }
663 
664 /*@C
665   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
666 
667   Collective
668 
669   Input Parameters:
670 + odm - The DMPlex specifying the mesh, passed as a PetscObject
671 - viewer - viewer of type VTK
672 
673   Level: developer
674 
675   Note:
676   This function is a callback used by the VTK viewer to actually write the file.
677   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
678   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
679 
680 .seealso: PETSCVIEWERVTK
681 @*/
682 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
683 {
684   DM             dm = (DM) odm;
685   PetscBool      isvtk;
686   PetscErrorCode ierr;
687 
688   PetscFunctionBegin;
689   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
690   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
691   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
692   PetscCheckFalse(!isvtk,PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
693   switch (viewer->format) {
694   case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
695     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
696     break;
697   case PETSC_VIEWER_VTK_VTU:
698     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
699     break;
700   default: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
701   }
702   PetscFunctionReturn(0);
703 }
704