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