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