xref: /petsc/src/dm/impls/plex/plexvtk.c (revision c4762a1b19cd2af06abeed90e8f9d34fb975dd94) 
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 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
275 typedef double PetscVTKReal;
276 #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
277 typedef float PetscVTKReal;
278 #else
279 typedef PetscReal PetscVTKReal;
280 #endif
281 
282 static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
283 {
284   MPI_Comm           comm;
285   const MPI_Datatype mpiType = MPIU_SCALAR;
286   PetscScalar        *array;
287   PetscInt           numDof = 0, maxDof;
288   PetscInt           numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
289   PetscMPIInt        size, rank, proc, tag;
290   PetscBool          hasLabel;
291   PetscErrorCode     ierr;
292 
293   PetscFunctionBegin;
294   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
295   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
296   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
297   if (precision < 0) precision = 6;
298   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
299   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
300   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
301   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
302   /* VTK only wants the values at cells or vertices */
303   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
304   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
305   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
306   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
307   if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
308   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
309   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
310   ierr     = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
311   ierr     = DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
312   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
313   for (p = pStart; p < pEnd; ++p) {
314     /* Reject points not either cells or vertices */
315     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
316     if (hasLabel) {
317       PetscInt value;
318 
319       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
320           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
321         ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
322         if (value != 1) continue;
323       }
324     }
325     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
326     if (numDof) break;
327   }
328   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
329   enforceDof = PetscMax(enforceDof, maxDof);
330   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
331   if (!rank) {
332     PetscVTKReal dval;
333     PetscScalar  val;
334     char formatString[8];
335 
336     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
337     for (p = pStart; p < pEnd; ++p) {
338       /* Here we lose a way to filter points by keeping them out of the Numbering */
339       PetscInt dof, off, goff, d;
340 
341       /* Reject points not either cells or vertices */
342       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
343       if (hasLabel) {
344         PetscInt value;
345 
346         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
347             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
348           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
349           if (value != 1) continue;
350         }
351       }
352       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
353       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
354       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
355       if (dof && goff >= 0) {
356         for (d = 0; d < dof; d++) {
357           if (d > 0) {
358             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
359           }
360           val = array[off+d];
361           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
362           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
363         }
364         for (d = dof; d < enforceDof; d++) {
365           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
366         }
367         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
368       }
369     }
370     for (proc = 1; proc < size; ++proc) {
371       PetscScalar *remoteValues;
372       PetscInt    size = 0, d;
373       MPI_Status  status;
374 
375       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
376       ierr = PetscMalloc1(size, &remoteValues);CHKERRQ(ierr);
377       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRQ(ierr);
378       for (p = 0; p < size/maxDof; ++p) {
379         for (d = 0; d < maxDof; ++d) {
380           if (d > 0) {
381             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
382           }
383           val = remoteValues[p*maxDof+d];
384           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
385           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
386         }
387         for (d = maxDof; d < enforceDof; ++d) {
388           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
389         }
390         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
391       }
392       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
393     }
394   } else {
395     PetscScalar *localValues;
396     PetscInt    size, k = 0;
397 
398     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
399     ierr = PetscMalloc1(size, &localValues);CHKERRQ(ierr);
400     for (p = pStart; p < pEnd; ++p) {
401       PetscInt dof, off, goff, d;
402 
403       /* Reject points not either cells or vertices */
404       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
405       if (hasLabel) {
406         PetscInt value;
407 
408         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
409             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
410           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
411           if (value != 1) continue;
412         }
413       }
414       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
415       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
416       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
417       if (goff >= 0) {
418         for (d = 0; d < dof; ++d) {
419           localValues[k++] = array[off+d];
420         }
421       }
422     }
423     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
424     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRQ(ierr);
425     ierr = PetscFree(localValues);CHKERRQ(ierr);
426   }
427   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
428   PetscFunctionReturn(0);
429 }
430 
431 static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscBool nameComplex, PetscInt imag)
432 {
433   MPI_Comm       comm;
434   PetscInt       numDof = 0, maxDof;
435   PetscInt       pStart, pEnd, p;
436   PetscErrorCode ierr;
437 
438   PetscFunctionBegin;
439   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
440   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
441   for (p = pStart; p < pEnd; ++p) {
442     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
443     if (numDof) break;
444   }
445   numDof = PetscMax(numDof, enforceDof);
446   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
447   if (!name) name = "Unknown";
448   if (maxDof == 3) {
449     if (nameComplex) {
450       ierr = PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");CHKERRQ(ierr);
451     } else {
452       ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
453     }
454   } else {
455     if (nameComplex) {
456       ierr = PetscFPrintf(comm, fp, "SCALARS %s.%s double %D\n", name, imag ? "Im" : "Re", maxDof);CHKERRQ(ierr);
457     } else {
458       ierr = PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);CHKERRQ(ierr);
459     }
460     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
461   }
462   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);CHKERRQ(ierr);
463   PetscFunctionReturn(0);
464 }
465 
466 static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
467 {
468   MPI_Comm                 comm;
469   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
470   FILE                     *fp;
471   PetscViewerVTKObjectLink link;
472   PetscSection             coordSection, globalCoordSection;
473   PetscLayout              vLayout;
474   Vec                      coordinates;
475   PetscReal                lengthScale;
476   PetscInt                 vMax, totVertices, totCells = 0, loops_per_scalar, l;
477   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
478   PetscErrorCode           ierr;
479 
480   PetscFunctionBegin;
481 #if defined(PETSC_USE_COMPLEX)
482   loops_per_scalar = 2;
483   writeComplex = PETSC_TRUE;
484 #else
485   loops_per_scalar = 1;
486   writeComplex = PETSC_FALSE;
487 #endif
488   ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr);
489   if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
490   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
491   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
492   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
493   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
494   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
495   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
496   /* Vertices */
497   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
498   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
499   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
500   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
501   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
502   if (vMax >= 0) {
503     PetscInt pStart, pEnd, p, localSize = 0;
504 
505     ierr = PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);CHKERRQ(ierr);
506     pEnd = PetscMin(pEnd, vMax);
507     for (p = pStart; p < pEnd; ++p) {
508       PetscInt dof;
509 
510       ierr = PetscSectionGetDof(globalCoordSection, p, &dof);CHKERRQ(ierr);
511       if (dof > 0) ++localSize;
512     }
513     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);CHKERRQ(ierr);
514     ierr = PetscLayoutSetLocalSize(vLayout, localSize);CHKERRQ(ierr);
515     ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
516     ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
517   } else {
518     ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);CHKERRQ(ierr);
519   }
520   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
521   ierr = PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);CHKERRQ(ierr);
522   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);CHKERRQ(ierr);
523   /* Cells */
524   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
525   /* Vertex fields */
526   for (link = vtk->link; link; link = link->next) {
527     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
528     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
529   }
530   if (hasPoint) {
531     ierr = PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);CHKERRQ(ierr);
532     for (link = vtk->link; link; link = link->next) {
533       Vec          X = (Vec) link->vec;
534       PetscSection section = NULL, globalSection, newSection = NULL;
535       char         namebuf[256];
536       const char   *name;
537       PetscInt     enforceDof = PETSC_DETERMINE;
538 
539       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
540       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
541       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
542       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
543       if (!section) {
544         DM           dmX;
545 
546         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
547         if (dmX) {
548           DMLabel  subpointMap, subpointMapX;
549           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
550 
551           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
552           /* Here is where we check whether dmX is a submesh of dm */
553           ierr = DMGetDimension(dm,  &dim);CHKERRQ(ierr);
554           ierr = DMGetDimension(dmX, &dimX);CHKERRQ(ierr);
555           ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
556           ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
557           ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
558           ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
559           if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
560             const PetscInt *ind = NULL;
561             IS              subpointIS;
562             PetscInt        n = 0, q;
563 
564             ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
565             ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
566             if (subpointIS) {
567               ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
568               ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
569             }
570             ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
571             ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
572             for (q = qStart; q < qEnd; ++q) {
573               PetscInt dof, off, p;
574 
575               ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
576               if (dof) {
577                 ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
578                 if (p >= pStart) {
579                   ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
580                   ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
581                   ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
582                 }
583               }
584             }
585             if (subpointIS) {
586               ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
587               ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
588             }
589             /* No need to setup section */
590             section = newSection;
591           }
592         }
593       }
594       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
595       if (link->field >= 0) {
596         const char *fieldname;
597 
598         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
599         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
600         if (fieldname) {
601           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);CHKERRQ(ierr);
602         } else {
603           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);CHKERRQ(ierr);
604         }
605       } else {
606         ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);CHKERRQ(ierr);
607       }
608       ierr = PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));CHKERRQ(ierr);
609       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
610       for (l = 0; l < loops_per_scalar; l++) {
611         ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);CHKERRQ(ierr);
612       }
613       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
614       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
615     }
616   }
617   /* Cell Fields */
618   ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
619   if (hasCell || writePartition) {
620     ierr = PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);CHKERRQ(ierr);
621     for (link = vtk->link; link; link = link->next) {
622       Vec          X = (Vec) link->vec;
623       PetscSection section = NULL, globalSection;
624       const char   *name = "";
625       char         namebuf[256];
626       PetscInt     enforceDof = PETSC_DETERMINE;
627 
628       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
629       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
630       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
631       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
632       if (!section) {
633         DM           dmX;
634 
635         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
636         if (dmX) {
637           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
638         }
639       }
640       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
641       if (link->field >= 0) {
642         const char *fieldname;
643 
644         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
645         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
646         if (fieldname) {
647           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);CHKERRQ(ierr);
648         } else {
649           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);CHKERRQ(ierr);
650         }
651       } else {
652         ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);CHKERRQ(ierr);
653       }
654       ierr = PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));CHKERRQ(ierr);
655       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
656       for (l = 0; l < loops_per_scalar; l++) {
657         ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);CHKERRQ(ierr);
658       }
659       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
660     }
661     if (writePartition) {
662       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
663       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
664       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
665     }
666   }
667   /* Cleanup */
668   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
669   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
670   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
671   PetscFunctionReturn(0);
672 }
673 
674 /*@C
675   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
676 
677   Collective
678 
679   Input Arguments:
680 + odm - The DMPlex specifying the mesh, passed as a PetscObject
681 - viewer - viewer of type VTK
682 
683   Level: developer
684 
685   Note:
686   This function is a callback used by the VTK viewer to actually write the file.
687   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
688   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
689 
690 .seealso: PETSCVIEWERVTK
691 @*/
692 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
693 {
694   DM             dm = (DM) odm;
695   PetscBool      isvtk;
696   PetscErrorCode ierr;
697 
698   PetscFunctionBegin;
699   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
700   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
701   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
702   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
703   switch (viewer->format) {
704   case PETSC_VIEWER_ASCII_VTK:
705     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
706     break;
707   case PETSC_VIEWER_VTK_VTU:
708     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
709     break;
710   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
711   }
712   PetscFunctionReturn(0);
713 }
714