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