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