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