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