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