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