xref: /petsc/src/dm/impls/plex/plexvtk.c (revision e8f6f0f66d69b00190dad407ead438ea5df5d8e1) 
1 #define PETSCDM_DLL
2 #include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
3 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexVTKGetCellType"
7 PetscErrorCode DMPlexVTKGetCellType(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
8 {
9   PetscFunctionBegin;
10   *cellType = -1;
11   switch (dim) {
12   case 0:
13     switch (corners) {
14     case 1:
15       *cellType = 1; /* VTK_VERTEX */
16       break;
17     default:
18       break;
19     }
20     break;
21   case 1:
22     switch (corners) {
23     case 2:
24       *cellType = 3; /* VTK_LINE */
25       break;
26     case 3:
27       *cellType = 21; /* VTK_QUADRATIC_EDGE */
28       break;
29     default:
30       break;
31     }
32     break;
33   case 2:
34     switch (corners) {
35     case 3:
36       *cellType = 5; /* VTK_TRIANGLE */
37       break;
38     case 4:
39       *cellType = 9; /* VTK_QUAD */
40       break;
41     case 6:
42       *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
43       break;
44     case 9:
45       *cellType = 23; /* VTK_QUADRATIC_QUAD */
46       break;
47     default:
48       break;
49     }
50     break;
51   case 3:
52     switch (corners) {
53     case 4:
54       *cellType = 10; /* VTK_TETRA */
55       break;
56     case 8:
57       *cellType = 12; /* VTK_HEXAHEDRON */
58       break;
59     case 10:
60       *cellType = 24; /* VTK_QUADRATIC_TETRA */
61       break;
62     case 27:
63       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
64       break;
65     default:
66       break;
67     }
68   }
69   PetscFunctionReturn(0);
70 }
71 
72 #undef __FUNCT__
73 #define __FUNCT__ "DMPlexVTKWriteCells_ASCII"
74 PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
75 {
76   MPI_Comm       comm;
77   IS             globalVertexNumbers;
78   const PetscInt *gvertex;
79   PetscInt       dim;
80   PetscInt       numCorners = 0, totCorners = 0, maxCorners, *corners;
81   PetscInt       numCells   = 0, totCells   = 0, maxCells, cellHeight;
82   PetscInt       numLabelCells, cMax, cStart, cEnd, c, vStart, vEnd, v;
83   PetscMPIInt    numProcs, rank, proc, tag;
84   PetscBool      hasLabel;
85   PetscErrorCode ierr;
86 
87   PetscFunctionBegin;
88   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
89   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
90   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
91   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
92   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
93   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
94   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
95   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
96   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
97   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
98   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
99 
100   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
101   for (c = cStart; c < cEnd; ++c) {
102     PetscInt *closure = NULL;
103     PetscInt closureSize;
104 
105     if (hasLabel) {
106       PetscInt value;
107 
108       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
109       if (value != 1) continue;
110     }
111     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
112     for (v = 0; v < closureSize*2; v += 2) {
113       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
114     }
115     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
116     ++numCells;
117   }
118   maxCells = numCells;
119   ierr     = MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
120   ierr     = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
121   ierr     = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
122   ierr     = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
123   ierr     = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
124   ierr     = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
125   ierr     = PetscMalloc(maxCells * sizeof(PetscInt), &corners);CHKERRQ(ierr);
126   ierr     = PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);CHKERRQ(ierr);
127   if (!rank) {
128     PetscInt *remoteVertices;
129 
130     for (c = cStart, numCells = 0; c < cEnd; ++c) {
131       PetscInt *closure = NULL;
132       PetscInt closureSize, nC = 0, tmp;
133 
134       if (hasLabel) {
135         PetscInt value;
136 
137         ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
138         if (value != 1) continue;
139       }
140       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
141       for (v = 0; v < closureSize*2; v += 2) {
142         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
143           const PetscInt gv = gvertex[closure[v] - vStart];
144           closure[nC++] = gv < 0 ? -(gv+1) : gv;
145         }
146       }
147       corners[numCells++] = nC;
148       ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
149       tmp        = closure[0];
150       closure[0] = closure[1];
151       closure[1] = tmp;
152       for (v = 0; v < nC; ++v) {
153         ierr = PetscFPrintf(comm, fp, " %d", closure[v]);CHKERRQ(ierr);
154       }
155       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
156       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);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         ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
168         for (v = 0; v < nC; ++v, ++c) {
169           ierr = PetscFPrintf(comm, fp, " %d", remoteVertices[c]);CHKERRQ(ierr);
170         }
171         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
172       }
173     }
174     if (numProcs > 1) {ierr = PetscFree(remoteVertices);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, nC = 0;
182 
183       if (hasLabel) {
184         PetscInt value;
185 
186         ierr = DMPlexGetLabelValue(dm, "vtk", 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 = DMPlexGetStratumSize(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 = DMPlexGetLabelValue(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     = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
313   ierr     = DMPlexGetStratumSize(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 = DMPlexGetLabelValue(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 = MPI_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 = DMPlexGetLabelValue(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, d;
371       MPI_Status  status;
372 
373       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
374       ierr = PetscMalloc(size * sizeof(PetscScalar), &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 = PetscMalloc(size * sizeof(PetscScalar), &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 = DMPlexGetLabelValue(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 = MPI_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 = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
482   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, 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 = DMPlexGetDimension(dm,  &dim);CHKERRQ(ierr);
533         ierr = DMPlexGetDimension(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;
540           IS              subpointIS;
541           PetscInt        n, q;
542 
543           ierr = PetscPrintf(PETSC_COMM_SELF, "Making translation PetscSection\n");CHKERRQ(ierr);
544           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
545           ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
546           ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
547           ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
548           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
549           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
550           for (q = qStart; q < qEnd; ++q) {
551             PetscInt dof, off, p;
552 
553             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
554             if (dof) {
555               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
556               if (p >= pStart) {
557                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
558                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
559                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
560               }
561             }
562           }
563           ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
564           ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
565           /* No need to setup section */
566           section = newSection;
567         }
568       } else {
569         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
570         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
571       }
572       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
573       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
574       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
575       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
576       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
577     }
578   }
579   /* Cell Fields */
580   ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
581   if (hasCell || writePartition) {
582     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
583     for (link = vtk->link; link; link = link->next) {
584       Vec          X = (Vec) link->vec;
585       DM           dmX;
586       PetscSection section, globalSection;
587       const char   *name;
588       PetscInt     enforceDof = PETSC_DETERMINE;
589 
590       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
591       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
592       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
593       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
594       if (dmX) {
595         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
596       } else {
597         PetscContainer c;
598 
599         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
600         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
601         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
602       }
603       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
604       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
605       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
606       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
607     }
608     if (writePartition) {
609       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
610       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
611       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
612     }
613   }
614   /* Cleanup */
615   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
616   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
617   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
618   PetscFunctionReturn(0);
619 }
620 
621 #undef __FUNCT__
622 #define __FUNCT__ "DMPlexVTKWriteAll"
623 /*@C
624   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
625 
626   Collective
627 
628   Input Arguments:
629 + odm - The DMPlex specifying the mesh, passed as a PetscObject
630 - viewer - viewer of type VTK
631 
632   Level: developer
633 
634   Note:
635   This function is a callback used by the VTK viewer to actually write the file.
636   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
637   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
638 
639 .seealso: PETSCVIEWERVTK
640 @*/
641 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
642 {
643   DM             dm = (DM) odm;
644   PetscBool      isvtk;
645   PetscErrorCode ierr;
646 
647   PetscFunctionBegin;
648   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
649   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
650   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
651   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
652   switch (viewer->format) {
653   case PETSC_VIEWER_ASCII_VTK:
654     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
655     break;
656   case PETSC_VIEWER_VTK_VTU:
657     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
658     break;
659   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
660   }
661   PetscFunctionReturn(0);
662 }
663