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 
DMPlexVTKGetCellType_Internal(DM dm,PetscInt dim,PetscInt corners,PetscInt * cellType)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 
DMPlexVTKWriteCells_ASCII(DM dm,FILE * fp,PetscInt * totalCells)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, *vertices;
125 
126     ierr = PetscMalloc1(maxCorners, &vertices);CHKERRQ(ierr);
127     for (c = cStart, numCells = 0; c < cEnd; ++c) {
128       PetscInt *closure = NULL;
129       PetscInt closureSize, value, nC = 0;
130 
131       if (label) {
132         ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
133         if (value != 1) continue;
134       }
135       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
136       for (v = 0; v < closureSize*2; v += 2) {
137         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
138           const PetscInt gv = gvertex[closure[v] - vStart];
139           vertices[nC++] = gv < 0 ? -(gv+1) : gv;
140         }
141       }
142       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
143       ierr = DMPlexReorderCell(dm, c, vertices);CHKERRQ(ierr);
144       corners[numCells++] = nC;
145       ierr = PetscFPrintf(comm, fp, "%D ", nC);CHKERRQ(ierr);
146       for (v = 0; v < nC; ++v) {
147         ierr = PetscFPrintf(comm, fp, " %D", vertices[v]);CHKERRQ(ierr);
148       }
149       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
150     }
151     if (size > 1) {ierr = PetscMalloc1(maxCorners+maxCells, &remoteVertices);CHKERRQ(ierr);}
152     for (proc = 1; proc < size; ++proc) {
153       MPI_Status status;
154 
155       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
156       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
157       for (c = 0; c < numCorners;) {
158         PetscInt nC = remoteVertices[c++];
159 
160         for (v = 0; v < nC; ++v, ++c) {
161           vertices[v] = remoteVertices[c];
162         }
163         ierr = PetscFPrintf(comm, fp, "%D ", nC);CHKERRQ(ierr);
164         for (v = 0; v < nC; ++v) {
165           ierr = PetscFPrintf(comm, fp, " %D", vertices[v]);CHKERRQ(ierr);
166         }
167         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
168       }
169     }
170     if (size > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
171     ierr = PetscFree(vertices);CHKERRQ(ierr);
172   } else {
173     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
174 
175     ierr = PetscMalloc1(numSend, &localVertices);CHKERRQ(ierr);
176     for (c = cStart, numCells = 0; c < cEnd; ++c) {
177       PetscInt *closure = NULL;
178       PetscInt closureSize, value, nC = 0;
179 
180       if (label) {
181         ierr = DMLabelGetValue(label, c, &value);CHKERRQ(ierr);
182         if (value != 1) continue;
183       }
184       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
185       for (v = 0; v < closureSize*2; v += 2) {
186         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
187           const PetscInt gv = gvertex[closure[v] - vStart];
188           closure[nC++] = gv < 0 ? -(gv+1) : gv;
189         }
190       }
191       corners[numCells++] = nC;
192       localVertices[k++]  = nC;
193       for (v = 0; v < nC; ++v, ++k) {
194         localVertices[k] = closure[v];
195       }
196       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
197       ierr = DMPlexReorderCell(dm, c, localVertices+k-nC);CHKERRQ(ierr);
198     }
199     if (k != numSend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %D should be %D", k, numSend);
200     ierr = MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
201     ierr = MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
202     ierr = PetscFree(localVertices);CHKERRQ(ierr);
203   }
204   ierr = ISRestoreIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
205   ierr = PetscFPrintf(comm, fp, "CELL_TYPES %D\n", totCells);CHKERRQ(ierr);
206   if (!rank) {
207     PetscInt cellType;
208 
209     for (c = 0; c < numCells; ++c) {
210       ierr = DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
211       ierr = PetscFPrintf(comm, fp, "%D\n", cellType);CHKERRQ(ierr);
212     }
213     for (proc = 1; proc < size; ++proc) {
214       MPI_Status status;
215 
216       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
217       ierr = MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
218       for (c = 0; c < numCells; ++c) {
219         ierr = DMPlexVTKGetCellType_Internal(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
220         ierr = PetscFPrintf(comm, fp, "%D\n", cellType);CHKERRQ(ierr);
221       }
222     }
223   } else {
224     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
225     ierr = MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
226   }
227   ierr        = PetscFree(corners);CHKERRQ(ierr);
228   *totalCells = totCells;
229   PetscFunctionReturn(0);
230 }
231 
DMPlexVTKWritePartition_ASCII(DM dm,FILE * fp)232 static PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
233 {
234   MPI_Comm       comm;
235   PetscInt       numCells = 0, cellHeight;
236   PetscInt       numLabelCells, cStart, cEnd, c;
237   PetscMPIInt    size, rank, proc, tag;
238   PetscBool      hasLabel;
239   PetscErrorCode ierr;
240 
241   PetscFunctionBegin;
242   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
243   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
244   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
245   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
246   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
247   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
248   ierr = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
249   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
250   for (c = cStart; c < cEnd; ++c) {
251     if (hasLabel) {
252       PetscInt value;
253 
254       ierr = DMGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
255       if (value != 1) continue;
256     }
257     ++numCells;
258   }
259   if (!rank) {
260     for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", rank);CHKERRQ(ierr);}
261     for (proc = 1; proc < size; ++proc) {
262       MPI_Status status;
263 
264       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
265       for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", proc);CHKERRQ(ierr);}
266     }
267   } else {
268     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
269   }
270   PetscFunctionReturn(0);
271 }
272 
273 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
274 typedef double PetscVTKReal;
275 #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
276 typedef float PetscVTKReal;
277 #else
278 typedef PetscReal PetscVTKReal;
279 #endif
280 
DMPlexVTKWriteSection_ASCII(DM dm,PetscSection section,PetscSection globalSection,Vec v,FILE * fp,PetscInt enforceDof,PetscInt precision,PetscReal scale,PetscInt imag)281 static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
282 {
283   MPI_Comm           comm;
284   const MPI_Datatype mpiType = MPIU_SCALAR;
285   PetscScalar        *array;
286   PetscInt           numDof = 0, maxDof;
287   PetscInt           numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vStart, vEnd, pStart, pEnd, p;
288   PetscMPIInt        size, rank, proc, tag;
289   PetscBool          hasLabel;
290   PetscErrorCode     ierr;
291 
292   PetscFunctionBegin;
293   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
294   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
295   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
296   if (precision < 0) precision = 6;
297   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
298   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
299   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
300   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
301   /* VTK only wants the values at cells or vertices */
302   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
303   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
304   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
305   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
306   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
307   ierr     = DMGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
308   ierr     = DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
309   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
310   for (p = pStart; p < pEnd; ++p) {
311     /* Reject points not either cells or vertices */
312     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
313     if (hasLabel) {
314       PetscInt value;
315 
316       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
317           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
318         ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
319         if (value != 1) continue;
320       }
321     }
322     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
323     if (numDof) break;
324   }
325   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
326   enforceDof = PetscMax(enforceDof, maxDof);
327   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
328   if (!rank) {
329     PetscVTKReal dval;
330     PetscScalar  val;
331     char formatString[8];
332 
333     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
334     for (p = pStart; p < pEnd; ++p) {
335       /* Here we lose a way to filter points by keeping them out of the Numbering */
336       PetscInt dof, off, goff, d;
337 
338       /* Reject points not either cells or vertices */
339       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
340       if (hasLabel) {
341         PetscInt value;
342 
343         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
344             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
345           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
346           if (value != 1) continue;
347         }
348       }
349       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
350       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
351       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
352       if (dof && goff >= 0) {
353         for (d = 0; d < dof; d++) {
354           if (d > 0) {
355             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
356           }
357           val = array[off+d];
358           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
359           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
360         }
361         for (d = dof; d < enforceDof; d++) {
362           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
363         }
364         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
365       }
366     }
367     for (proc = 1; proc < size; ++proc) {
368       PetscScalar *remoteValues;
369       PetscInt    size = 0, d;
370       MPI_Status  status;
371 
372       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
373       ierr = PetscMalloc1(size, &remoteValues);CHKERRQ(ierr);
374       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRQ(ierr);
375       for (p = 0; p < size/maxDof; ++p) {
376         for (d = 0; d < maxDof; ++d) {
377           if (d > 0) {
378             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
379           }
380           val = remoteValues[p*maxDof+d];
381           dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
382           ierr = PetscFPrintf(comm, fp, formatString, dval);CHKERRQ(ierr);
383         }
384         for (d = maxDof; d < enforceDof; ++d) {
385           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
386         }
387         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
388       }
389       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
390     }
391   } else {
392     PetscScalar *localValues;
393     PetscInt    size, k = 0;
394 
395     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
396     ierr = PetscMalloc1(size, &localValues);CHKERRQ(ierr);
397     for (p = pStart; p < pEnd; ++p) {
398       PetscInt dof, off, goff, d;
399 
400       /* Reject points not either cells or vertices */
401       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
402       if (hasLabel) {
403         PetscInt value;
404 
405         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
406             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
407           ierr = DMGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
408           if (value != 1) continue;
409         }
410       }
411       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
412       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
413       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
414       if (goff >= 0) {
415         for (d = 0; d < dof; ++d) {
416           localValues[k++] = array[off+d];
417         }
418       }
419     }
420     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
421     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRQ(ierr);
422     ierr = PetscFree(localValues);CHKERRQ(ierr);
423   }
424   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
425   PetscFunctionReturn(0);
426 }
427 
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)428 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)
429 {
430   MPI_Comm       comm;
431   PetscInt       numDof = 0, maxDof;
432   PetscInt       pStart, pEnd, p;
433   PetscErrorCode ierr;
434 
435   PetscFunctionBegin;
436   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
437   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
438   for (p = pStart; p < pEnd; ++p) {
439     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
440     if (numDof) break;
441   }
442   numDof = PetscMax(numDof, enforceDof);
443   ierr = MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
444   if (!name) name = "Unknown";
445   if (maxDof == 3) {
446     if (nameComplex) {
447       ierr = PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");CHKERRQ(ierr);
448     } else {
449       ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
450     }
451   } else {
452     if (nameComplex) {
453       ierr = PetscFPrintf(comm, fp, "SCALARS %s.%s double %D\n", name, imag ? "Im" : "Re", maxDof);CHKERRQ(ierr);
454     } else {
455       ierr = PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);CHKERRQ(ierr);
456     }
457     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
458   }
459   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);CHKERRQ(ierr);
460   PetscFunctionReturn(0);
461 }
462 
DMPlexVTKWriteAll_ASCII(DM dm,PetscViewer viewer)463 static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
464 {
465   MPI_Comm                 comm;
466   PetscViewer_VTK          *vtk = (PetscViewer_VTK*) viewer->data;
467   FILE                     *fp;
468   PetscViewerVTKObjectLink link;
469   PetscSection             coordSection, globalCoordSection;
470   PetscLayout              vLayout;
471   Vec                      coordinates;
472   PetscReal                lengthScale;
473   PetscInt                 totVertices, totCells = 0, loops_per_scalar, l;
474   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
475   PetscErrorCode           ierr;
476 
477   PetscFunctionBegin;
478 #if defined(PETSC_USE_COMPLEX)
479   loops_per_scalar = 2;
480   writeComplex = PETSC_TRUE;
481 #else
482   loops_per_scalar = 1;
483   writeComplex = PETSC_FALSE;
484 #endif
485   ierr = DMGetCoordinatesLocalized(dm,&localized);CHKERRQ(ierr);
486   if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
487   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
488   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
489   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
490   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
491   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
492   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
493   /* Vertices */
494   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
495   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
496   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
497   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
498   ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);CHKERRQ(ierr);
499   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
500   ierr = PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);CHKERRQ(ierr);
501   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);CHKERRQ(ierr);
502   /* Cells */
503   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
504   /* Vertex fields */
505   for (link = vtk->link; link; link = link->next) {
506     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
507     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
508   }
509   if (hasPoint) {
510     ierr = PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);CHKERRQ(ierr);
511     for (link = vtk->link; link; link = link->next) {
512       Vec          X = (Vec) link->vec;
513       PetscSection section = NULL, globalSection, newSection = NULL;
514       char         namebuf[256];
515       const char   *name;
516       PetscInt     enforceDof = PETSC_DETERMINE;
517 
518       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
519       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
520       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
521       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
522       if (!section) {
523         DM           dmX;
524 
525         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
526         if (dmX) {
527           DMLabel  subpointMap, subpointMapX;
528           PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
529 
530           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
531           /* Here is where we check whether dmX is a submesh of dm */
532           ierr = DMGetDimension(dm,  &dim);CHKERRQ(ierr);
533           ierr = DMGetDimension(dmX, &dimX);CHKERRQ(ierr);
534           ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
535           ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
536           ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
537           ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
538           if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
539             const PetscInt *ind = NULL;
540             IS              subpointIS;
541             PetscInt        n = 0, q;
542 
543             ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
544             ierr = DMPlexGetSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
545             if (subpointIS) {
546               ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
547               ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
548             }
549             ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
550             ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
551             for (q = qStart; q < qEnd; ++q) {
552               PetscInt dof, off, p;
553 
554               ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
555               if (dof) {
556                 ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
557                 if (p >= pStart) {
558                   ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
559                   ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
560                   ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
561                 }
562               }
563             }
564             if (subpointIS) {
565               ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
566             }
567             /* No need to setup section */
568             section = newSection;
569           }
570         }
571       }
572       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);
573       if (link->field >= 0) {
574         const char *fieldname;
575 
576         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
577         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
578         if (fieldname) {
579           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);CHKERRQ(ierr);
580         } else {
581           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);CHKERRQ(ierr);
582         }
583       } else {
584         ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);CHKERRQ(ierr);
585       }
586       ierr = PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));CHKERRQ(ierr);
587       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
588       for (l = 0; l < loops_per_scalar; l++) {
589         ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);CHKERRQ(ierr);
590       }
591       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
592       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
593     }
594   }
595   /* Cell Fields */
596   ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
597   if (hasCell || writePartition) {
598     ierr = PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);CHKERRQ(ierr);
599     for (link = vtk->link; link; link = link->next) {
600       Vec          X = (Vec) link->vec;
601       PetscSection section = NULL, globalSection;
602       const char   *name = "";
603       char         namebuf[256];
604       PetscInt     enforceDof = PETSC_DETERMINE;
605 
606       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
607       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
608       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
609       ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
610       if (!section) {
611         DM           dmX;
612 
613         ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
614         if (dmX) {
615           ierr = DMGetLocalSection(dmX, &section);CHKERRQ(ierr);
616         }
617       }
618       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);
619       if (link->field >= 0) {
620         const char *fieldname;
621 
622         ierr = PetscSectionGetFieldName(section, link->field, &fieldname);CHKERRQ(ierr);
623         ierr = PetscSectionGetField(section, link->field, &section);CHKERRQ(ierr);
624         if (fieldname) {
625           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);CHKERRQ(ierr);
626         } else {
627           ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);CHKERRQ(ierr);
628         }
629       } else {
630         ierr = PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);CHKERRQ(ierr);
631       }
632       ierr = PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));CHKERRQ(ierr);
633       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
634       for (l = 0; l < loops_per_scalar; l++) {
635         ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);CHKERRQ(ierr);
636       }
637       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
638     }
639     if (writePartition) {
640       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
641       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
642       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
643     }
644   }
645   /* Cleanup */
646   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
647   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
648   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
649   PetscFunctionReturn(0);
650 }
651 
652 /*@C
653   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
654 
655   Collective
656 
657   Input Arguments:
658 + odm - The DMPlex specifying the mesh, passed as a PetscObject
659 - viewer - viewer of type VTK
660 
661   Level: developer
662 
663   Note:
664   This function is a callback used by the VTK viewer to actually write the file.
665   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
666   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
667 
668 .seealso: PETSCVIEWERVTK
669 @*/
DMPlexVTKWriteAll(PetscObject odm,PetscViewer viewer)670 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
671 {
672   DM             dm = (DM) odm;
673   PetscBool      isvtk;
674   PetscErrorCode ierr;
675 
676   PetscFunctionBegin;
677   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
678   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
679   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
680   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
681   switch (viewer->format) {
682   case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
683     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
684     break;
685   case PETSC_VIEWER_VTK_VTU:
686     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
687     break;
688   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
689   }
690   PetscFunctionReturn(0);
691 }
692