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