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