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