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