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