xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 02281ff37adb6085fb26862da565c8baf3221123) 
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   DMLabel        label;
78   IS             globalVertexNumbers = NULL;
79   const PetscInt *gvertex;
80   PetscInt       dim;
81   PetscInt       numCorners = 0, totCorners = 0, maxCorners, *corners;
82   PetscInt       numCells   = 0, totCells   = 0, maxCells, cellHeight;
83   PetscInt       numLabelCells, cMax, cStart, cEnd, c, vStart, vEnd, v;
84   PetscMPIInt    numProcs, rank, proc, tag;
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 = DMPlexGetLabel(dm, "vtk", &label);CHKERRQ(ierr);
99   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
100   for (c = cStart; c < cEnd; ++c) {
101     PetscInt *closure = NULL;
102     PetscInt closureSize, value;
103 
104     if (label) {
105       ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
106       if (value != 1) continue;
107     }
108     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
109     for (v = 0; v < closureSize*2; v += 2) {
110       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
111     }
112     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
113     ++numCells;
114   }
115   maxCells = numCells;
116   ierr     = MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
117   ierr     = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
118   ierr     = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
119   ierr     = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
120   ierr     = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
121   ierr     = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
122   ierr     = PetscMalloc(maxCells * sizeof(PetscInt), &corners);CHKERRQ(ierr);
123   ierr     = PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);CHKERRQ(ierr);
124   if (!rank) {
125     PetscInt *remoteVertices;
126     int      *vertices;
127 
128     ierr = PetscMalloc(maxCorners * sizeof(int), &vertices);CHKERRQ(ierr);
129     for (c = cStart, numCells = 0; c < cEnd; ++c) {
130       PetscInt *closure = NULL;
131       PetscInt closureSize, value, nC = 0;
132 
133       if (label) {
134         ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
135         if (value != 1) continue;
136       }
137       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
138       for (v = 0; v < closureSize*2; v += 2) {
139         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
140           const PetscInt gv = gvertex[closure[v] - vStart];
141           vertices[nC++] = gv < 0 ? -(gv+1) : gv;
142         }
143       }
144       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
145       corners[numCells++] = nC;
146       ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
147       ierr = DMPlexInvertCell(dim, nC, vertices);CHKERRQ(ierr);
148       for (v = 0; v < nC; ++v) {
149         ierr = PetscFPrintf(comm, fp, " %d", vertices[v]);CHKERRQ(ierr);
150       }
151       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
152     }
153     if (numProcs > 1) {ierr = PetscMalloc((maxCorners+maxCells) * sizeof(PetscInt), &remoteVertices);CHKERRQ(ierr);}
154     for (proc = 1; proc < numProcs; ++proc) {
155       MPI_Status status;
156 
157       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
158       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
159       for (c = 0; c < numCorners;) {
160         PetscInt nC = remoteVertices[c++];
161 
162         for (v = 0; v < nC; ++v, ++c) {
163           vertices[v] = remoteVertices[c];
164         }
165         ierr = DMPlexInvertCell(dim, nC, vertices);CHKERRQ(ierr);
166         ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
167         for (v = 0; v < nC; ++v) {
168           ierr = PetscFPrintf(comm, fp, " %d", vertices[v]);CHKERRQ(ierr);
169         }
170         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
171       }
172     }
173     if (numProcs > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
174     ierr = PetscFree(vertices);CHKERRQ(ierr);
175   } else {
176     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
177 
178     ierr = PetscMalloc(numSend * sizeof(PetscInt), &localVertices);CHKERRQ(ierr);
179     for (c = cStart, numCells = 0; c < cEnd; ++c) {
180       PetscInt *closure = NULL;
181       PetscInt closureSize, value, nC = 0;
182 
183       if (label) {
184         ierr = DMLabelGetValue(label, 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 = 0, 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 = NULL;
538           IS              subpointIS;
539           PetscInt        n = 0, q;
540 
541           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
542           ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
543           if (subpointIS) {
544             ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
545             ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
546           }
547           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
548           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
549           for (q = qStart; q < qEnd; ++q) {
550             PetscInt dof, off, p;
551 
552             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
553             if (dof) {
554               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
555               if (p >= pStart) {
556                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
557                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
558                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
559               }
560             }
561           }
562           if (subpointIS) {
563             ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
564             ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
565           }
566           /* No need to setup section */
567           section = newSection;
568         }
569       } else {
570         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
571         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
572       }
573       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
574       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
575       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
576       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
577       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
578     }
579   }
580   /* Cell Fields */
581   ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
582   if (hasCell || writePartition) {
583     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
584     for (link = vtk->link; link; link = link->next) {
585       Vec          X = (Vec) link->vec;
586       DM           dmX;
587       PetscSection section, globalSection;
588       const char   *name;
589       PetscInt     enforceDof = PETSC_DETERMINE;
590 
591       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
592       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
593       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
594       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
595       if (dmX) {
596         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
597       } else {
598         PetscContainer c;
599 
600         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
601         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
602         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
603       }
604       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
605       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
606       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
607       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
608     }
609     if (writePartition) {
610       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
611       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
612       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
613     }
614   }
615   /* Cleanup */
616   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
617   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
618   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
619   PetscFunctionReturn(0);
620 }
621 
622 #undef __FUNCT__
623 #define __FUNCT__ "DMPlexVTKWriteAll"
624 /*@C
625   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
626 
627   Collective
628 
629   Input Arguments:
630 + odm - The DMPlex specifying the mesh, passed as a PetscObject
631 - viewer - viewer of type VTK
632 
633   Level: developer
634 
635   Note:
636   This function is a callback used by the VTK viewer to actually write the file.
637   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
638   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
639 
640 .seealso: PETSCVIEWERVTK
641 @*/
642 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
643 {
644   DM             dm = (DM) odm;
645   PetscBool      isvtk;
646   PetscErrorCode ierr;
647 
648   PetscFunctionBegin;
649   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
650   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
651   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
652   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
653   switch (viewer->format) {
654   case PETSC_VIEWER_ASCII_VTK:
655     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
656     break;
657   case PETSC_VIEWER_VTK_VTU:
658     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
659     break;
660   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
661   }
662   PetscFunctionReturn(0);
663 }
664