xref: /petsc/src/dm/impls/plex/plexvtk.c (revision 0dcf9c705d9e3a176e05c9fc6d3f909d0a86dcfd) 
1 #define PETSCDM_DLL
2 #include <petsc-private/dmpleximpl.h>    /*I   "petscdmplex.h"   I*/
3 #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexVTKGetCellType"
7 PetscErrorCode DMPlexVTKGetCellType(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
8 {
9   PetscFunctionBegin;
10   *cellType = -1;
11   switch (dim) {
12   case 0:
13     switch (corners) {
14     case 1:
15       *cellType = 1; /* VTK_VERTEX */
16       break;
17     default:
18       break;
19     }
20     break;
21   case 1:
22     switch (corners) {
23     case 2:
24       *cellType = 3; /* VTK_LINE */
25       break;
26     case 3:
27       *cellType = 21; /* VTK_QUADRATIC_EDGE */
28       break;
29     default:
30       break;
31     }
32     break;
33   case 2:
34     switch (corners) {
35     case 3:
36       *cellType = 5; /* VTK_TRIANGLE */
37       break;
38     case 4:
39       *cellType = 9; /* VTK_QUAD */
40       break;
41     case 6:
42       *cellType = 22; /* VTK_QUADRATIC_TRIANGLE */
43       break;
44     case 9:
45       *cellType = 23; /* VTK_QUADRATIC_QUAD */
46       break;
47     default:
48       break;
49     }
50     break;
51   case 3:
52     switch (corners) {
53     case 4:
54       *cellType = 10; /* VTK_TETRA */
55       break;
56     case 8:
57       *cellType = 12; /* VTK_HEXAHEDRON */
58       break;
59     case 10:
60       *cellType = 24; /* VTK_QUADRATIC_TETRA */
61       break;
62     case 27:
63       *cellType = 29; /* VTK_QUADRATIC_HEXAHEDRON */
64       break;
65     default:
66       break;
67     }
68   }
69   PetscFunctionReturn(0);
70 }
71 
72 #undef __FUNCT__
73 #define __FUNCT__ "DMPlexVTKWriteCells_ASCII"
74 PetscErrorCode DMPlexVTKWriteCells_ASCII(DM dm, FILE *fp, PetscInt *totalCells)
75 {
76   MPI_Comm       comm;
77   IS             globalVertexNumbers;
78   const PetscInt *gvertex;
79   PetscInt       dim;
80   PetscInt       numCorners = 0, totCorners = 0, maxCorners, *corners;
81   PetscInt       numCells   = 0, totCells   = 0, maxCells, cellHeight;
82   PetscInt       numLabelCells, cMax, cStart, cEnd, c, vStart, vEnd, v;
83   PetscMPIInt    numProcs, rank, proc, tag;
84   PetscBool      hasLabel;
85   PetscErrorCode ierr;
86 
87   PetscFunctionBegin;
88   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
89   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
90   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
91   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
92   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
93   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
94   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
95   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
96   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
97   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
98   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
99 
100   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
101   for (c = cStart; c < cEnd; ++c) {
102     PetscInt *closure = NULL;
103     PetscInt closureSize;
104 
105     if (hasLabel) {
106       PetscInt value;
107 
108       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
109       if (value != 1) continue;
110     }
111     ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
112     for (v = 0; v < closureSize*2; v += 2) {
113       if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
114     }
115     ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
116     ++numCells;
117   }
118   maxCells = numCells;
119   ierr     = MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
120   ierr     = MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
121   ierr     = MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);CHKERRQ(ierr);
122   ierr     = MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr);
123   ierr     = DMPlexGetVertexNumbering(dm, &globalVertexNumbers);CHKERRQ(ierr);
124   ierr     = ISGetIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
125   ierr     = PetscMalloc(maxCells * sizeof(PetscInt), &corners);CHKERRQ(ierr);
126   ierr     = PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);CHKERRQ(ierr);
127   if (!rank) {
128     PetscInt *remoteVertices;
129     int      *vertices;
130 
131     ierr = PetscMalloc(maxCorners * sizeof(int), &vertices);CHKERRQ(ierr);
132     for (c = cStart, numCells = 0; c < cEnd; ++c) {
133       PetscInt *closure = NULL;
134       PetscInt closureSize, nC = 0;
135 
136       if (hasLabel) {
137         PetscInt value;
138 
139         ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
140         if (value != 1) continue;
141       }
142       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
143       for (v = 0; v < closureSize*2; v += 2) {
144         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
145           const PetscInt gv = gvertex[closure[v] - vStart];
146           vertices[nC++] = gv < 0 ? -(gv+1) : gv;
147         }
148       }
149       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
150       corners[numCells++] = nC;
151       ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
152       ierr = DMPlexInvertCell(dim, nC, vertices);CHKERRQ(ierr);
153       for (v = 0; v < nC; ++v) {
154         ierr = PetscFPrintf(comm, fp, " %d", vertices[v]);CHKERRQ(ierr);
155       }
156       ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
157     }
158     if (numProcs > 1) {ierr = PetscMalloc((maxCorners+maxCells) * sizeof(PetscInt), &remoteVertices);CHKERRQ(ierr);}
159     for (proc = 1; proc < numProcs; ++proc) {
160       MPI_Status status;
161 
162       ierr = MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
163       ierr = MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
164       for (c = 0; c < numCorners;) {
165         PetscInt nC = remoteVertices[c++];
166 
167         for (v = 0; v < nC; ++v, ++c) {
168           vertices[v] = remoteVertices[c];
169         }
170         ierr = DMPlexInvertCell(dim, nC, vertices);CHKERRQ(ierr);
171         ierr = PetscFPrintf(comm, fp, "%d ", nC);CHKERRQ(ierr);
172         for (v = 0; v < nC; ++v) {
173           ierr = PetscFPrintf(comm, fp, " %d", vertices[v]);CHKERRQ(ierr);
174         }
175         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
176       }
177     }
178     if (numProcs > 1) {ierr = PetscFree(remoteVertices);CHKERRQ(ierr);}
179     ierr = PetscFree(vertices);CHKERRQ(ierr);
180   } else {
181     PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
182 
183     ierr = PetscMalloc(numSend * sizeof(PetscInt), &localVertices);CHKERRQ(ierr);
184     for (c = cStart, numCells = 0; c < cEnd; ++c) {
185       PetscInt *closure = NULL;
186       PetscInt closureSize, nC = 0;
187 
188       if (hasLabel) {
189         PetscInt value;
190 
191         ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
192         if (value != 1) continue;
193       }
194       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
195       for (v = 0; v < closureSize*2; v += 2) {
196         if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
197           const PetscInt gv = gvertex[closure[v] - vStart];
198           closure[nC++] = gv < 0 ? -(gv+1) : gv;
199         }
200       }
201       corners[numCells++] = nC;
202       localVertices[k++]  = nC;
203       for (v = 0; v < nC; ++v, ++k) {
204         localVertices[k] = closure[v];
205       }
206       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
207     }
208     if (k != numSend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numSend);
209     ierr = MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
210     ierr = MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
211     ierr = PetscFree(localVertices);CHKERRQ(ierr);
212   }
213   ierr = ISRestoreIndices(globalVertexNumbers, &gvertex);CHKERRQ(ierr);
214   ierr = PetscFPrintf(comm, fp, "CELL_TYPES %d\n", totCells);CHKERRQ(ierr);
215   if (!rank) {
216     PetscInt cellType;
217 
218     for (c = 0; c < numCells; ++c) {
219       ierr = DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
220       ierr = PetscFPrintf(comm, fp, "%d\n", cellType);CHKERRQ(ierr);
221     }
222     for (proc = 1; proc < numProcs; ++proc) {
223       MPI_Status status;
224 
225       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
226       ierr = MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
227       for (c = 0; c < numCells; ++c) {
228         ierr = DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);CHKERRQ(ierr);
229         ierr = PetscFPrintf(comm, fp, "%d\n", cellType);CHKERRQ(ierr);
230       }
231     }
232   } else {
233     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
234     ierr = MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
235   }
236   ierr        = PetscFree(corners);CHKERRQ(ierr);
237   *totalCells = totCells;
238   PetscFunctionReturn(0);
239 }
240 
241 #undef __FUNCT__
242 #define __FUNCT__ "DMPlexVTKWritePartition_ASCII"
243 PetscErrorCode DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
244 {
245   MPI_Comm       comm;
246   PetscInt       numCells = 0, cellHeight;
247   PetscInt       numLabelCells, cMax, cStart, cEnd, c;
248   PetscMPIInt    numProcs, rank, proc, tag;
249   PetscBool      hasLabel;
250   PetscErrorCode ierr;
251 
252   PetscFunctionBegin;
253   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
254   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
255   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
256   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
257   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
258   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
259   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);CHKERRQ(ierr);
260   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
261   ierr = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
262   hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
263   for (c = cStart; c < cEnd; ++c) {
264     if (hasLabel) {
265       PetscInt value;
266 
267       ierr = DMPlexGetLabelValue(dm, "vtk", c, &value);CHKERRQ(ierr);
268       if (value != 1) continue;
269     }
270     ++numCells;
271   }
272   if (!rank) {
273     for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", rank);CHKERRQ(ierr);}
274     for (proc = 1; proc < numProcs; ++proc) {
275       MPI_Status status;
276 
277       ierr = MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
278       for (c = 0; c < numCells; ++c) {ierr = PetscFPrintf(comm, fp, "%d\n", proc);CHKERRQ(ierr);}
279     }
280   } else {
281     ierr = MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
282   }
283   PetscFunctionReturn(0);
284 }
285 
286 #undef __FUNCT__
287 #define __FUNCT__ "DMPlexVTKWriteSection_ASCII"
288 PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
289 {
290   MPI_Comm           comm;
291   const MPI_Datatype mpiType = MPIU_SCALAR;
292   PetscScalar        *array;
293   PetscInt           numDof = 0, maxDof;
294   PetscInt           numLabelCells, cellHeight, cMax, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
295   PetscMPIInt        numProcs, rank, proc, tag;
296   PetscBool          hasLabel;
297   PetscErrorCode     ierr;
298 
299   PetscFunctionBegin;
300   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
301   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
302   PetscValidHeaderSpecific(v,VEC_CLASSID,4);
303   if (precision < 0) precision = 6;
304   ierr = PetscCommGetNewTag(comm, &tag);CHKERRQ(ierr);
305   ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
306   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
307   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
308   /* VTK only wants the values at cells or vertices */
309   ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr);
310   ierr = DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);CHKERRQ(ierr);
311   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
312   ierr = DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, &vMax);CHKERRQ(ierr);
313   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
314   if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
315   pStart   = PetscMax(PetscMin(cStart, vStart), pStart);
316   pEnd     = PetscMin(PetscMax(cEnd,   vEnd),   pEnd);
317   ierr     = DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);CHKERRQ(ierr);
318   ierr     = DMPlexGetStratumSize(dm, "vtk", 2, &numLabelVertices);CHKERRQ(ierr);
319   hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
320   for (p = pStart; p < pEnd; ++p) {
321     /* Reject points not either cells or vertices */
322     if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
323     if (hasLabel) {
324       PetscInt value;
325 
326       if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
327           ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
328         ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
329         if (value != 1) continue;
330       }
331     }
332     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
333     if (numDof) break;
334   }
335   ierr = MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr);
336   enforceDof = PetscMax(enforceDof, maxDof);
337   ierr = VecGetArray(v, &array);CHKERRQ(ierr);
338   if (!rank) {
339     char formatString[8];
340 
341     ierr = PetscSNPrintf(formatString, 8, "%%.%de", precision);CHKERRQ(ierr);
342     for (p = pStart; p < pEnd; ++p) {
343       /* Here we lose a way to filter points by keeping them out of the Numbering */
344       PetscInt dof, off, goff, d;
345 
346       /* Reject points not either cells or vertices */
347       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
348       if (hasLabel) {
349         PetscInt value;
350 
351         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
352             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
353           ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
354           if (value != 1) continue;
355         }
356       }
357       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
358       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
359       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
360       if (dof && goff >= 0) {
361         for (d = 0; d < dof; d++) {
362           if (d > 0) {
363             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
364           }
365           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);CHKERRQ(ierr);
366         }
367         for (d = dof; d < enforceDof; d++) {
368           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
369         }
370         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
371       }
372     }
373     for (proc = 1; proc < numProcs; ++proc) {
374       PetscScalar *remoteValues;
375       PetscInt    size, d;
376       MPI_Status  status;
377 
378       ierr = MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);CHKERRQ(ierr);
379       ierr = PetscMalloc(size * sizeof(PetscScalar), &remoteValues);CHKERRQ(ierr);
380       ierr = MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);CHKERRQ(ierr);
381       for (p = 0; p < size/maxDof; ++p) {
382         for (d = 0; d < maxDof; ++d) {
383           if (d > 0) {
384             ierr = PetscFPrintf(comm, fp, " ");CHKERRQ(ierr);
385           }
386           ierr = PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);CHKERRQ(ierr);
387         }
388         for (d = maxDof; d < enforceDof; ++d) {
389           ierr = PetscFPrintf(comm, fp, " 0.0");CHKERRQ(ierr);
390         }
391         ierr = PetscFPrintf(comm, fp, "\n");CHKERRQ(ierr);
392       }
393       ierr = PetscFree(remoteValues);CHKERRQ(ierr);
394     }
395   } else {
396     PetscScalar *localValues;
397     PetscInt    size, k = 0;
398 
399     ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr);
400     ierr = PetscMalloc(size * sizeof(PetscScalar), &localValues);CHKERRQ(ierr);
401     for (p = pStart; p < pEnd; ++p) {
402       PetscInt dof, off, goff, d;
403 
404       /* Reject points not either cells or vertices */
405       if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
406       if (hasLabel) {
407         PetscInt value;
408 
409         if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
410             ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
411           ierr = DMPlexGetLabelValue(dm, "vtk", p, &value);CHKERRQ(ierr);
412           if (value != 1) continue;
413         }
414       }
415       ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
416       ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
417       ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr);
418       if (goff >= 0) {
419         for (d = 0; d < dof; ++d) {
420           localValues[k++] = array[off+d];
421         }
422       }
423     }
424     ierr = MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);CHKERRQ(ierr);
425     ierr = MPI_Send(localValues, k, mpiType, 0, tag, comm);CHKERRQ(ierr);
426     ierr = PetscFree(localValues);CHKERRQ(ierr);
427   }
428   ierr = VecRestoreArray(v, &array);CHKERRQ(ierr);
429   PetscFunctionReturn(0);
430 }
431 
432 #undef __FUNCT__
433 #define __FUNCT__ "DMPlexVTKWriteField_ASCII"
434 PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
435 {
436   MPI_Comm       comm;
437   PetscInt       numDof = 0, maxDof;
438   PetscInt       pStart, pEnd, p;
439   PetscErrorCode ierr;
440 
441   PetscFunctionBegin;
442   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
443   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
444   for (p = pStart; p < pEnd; ++p) {
445     ierr = PetscSectionGetDof(section, p, &numDof);CHKERRQ(ierr);
446     if (numDof) break;
447   }
448   numDof = PetscMax(numDof, enforceDof);
449   ierr = MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
450   if (!name) name = "Unknown";
451   if (maxDof == 3) {
452     ierr = PetscFPrintf(comm, fp, "VECTORS %s double\n", name);CHKERRQ(ierr);
453   } else {
454     ierr = PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);CHKERRQ(ierr);
455     ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
456   }
457   ierr = DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);CHKERRQ(ierr);
458   PetscFunctionReturn(0);
459 }
460 
461 #undef __FUNCT__
462 #define __FUNCT__ "DMPlexVTKWriteAll_ASCII"
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                 vMax, totVertices, totCells;
474   PetscBool                hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE;
475   PetscErrorCode           ierr;
476 
477   PetscFunctionBegin;
478   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
479   ierr = PetscFOpen(comm, vtk->filename, "wb", &fp);CHKERRQ(ierr);
480   ierr = PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr);
481   ierr = PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");CHKERRQ(ierr);
482   ierr = PetscFPrintf(comm, fp, "ASCII\n");CHKERRQ(ierr);
483   ierr = PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");CHKERRQ(ierr);
484   /* Vertices */
485   ierr = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);CHKERRQ(ierr);
486   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
487   ierr = PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, &globalCoordSection);CHKERRQ(ierr);
488   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
489   ierr = DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);CHKERRQ(ierr);
490   if (vMax >= 0) {
491     PetscInt pStart, pEnd, p, localSize = 0;
492 
493     ierr = PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);CHKERRQ(ierr);
494     pEnd = PetscMin(pEnd, vMax);
495     for (p = pStart; p < pEnd; ++p) {
496       PetscInt dof;
497 
498       ierr = PetscSectionGetDof(globalCoordSection, p, &dof);CHKERRQ(ierr);
499       if (dof > 0) ++localSize;
500     }
501     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);CHKERRQ(ierr);
502     ierr = PetscLayoutSetLocalSize(vLayout, localSize);CHKERRQ(ierr);
503     ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
504     ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
505   } else {
506     ierr = PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);CHKERRQ(ierr);
507   }
508   ierr = PetscLayoutGetSize(vLayout, &totVertices);CHKERRQ(ierr);
509   ierr = PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);CHKERRQ(ierr);
510   ierr = DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);CHKERRQ(ierr);
511   /* Cells */
512   ierr = DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);CHKERRQ(ierr);
513   /* Vertex fields */
514   for (link = vtk->link; link; link = link->next) {
515     if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
516     if ((link->ft == PETSC_VTK_CELL_FIELD)  || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD))  hasCell  = PETSC_TRUE;
517   }
518   if (hasPoint) {
519     ierr = PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);CHKERRQ(ierr);
520     for (link = vtk->link; link; link = link->next) {
521       Vec          X = (Vec) link->vec;
522       DM           dmX;
523       PetscSection section, globalSection, newSection = NULL;
524       const char   *name;
525       PetscInt     enforceDof = PETSC_DETERMINE;
526 
527       if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
528       if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
529       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
530       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
531       if (dmX) {
532         DMLabel  subpointMap, subpointMapX;
533         PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
534 
535         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
536         /* Here is where we check whether dmX is a submesh of dm */
537         ierr = DMPlexGetDimension(dm,  &dim);CHKERRQ(ierr);
538         ierr = DMPlexGetDimension(dmX, &dimX);CHKERRQ(ierr);
539         ierr = DMPlexGetChart(dm,  &pStart, &pEnd);CHKERRQ(ierr);
540         ierr = DMPlexGetChart(dmX, &qStart, &qEnd);CHKERRQ(ierr);
541         ierr = DMPlexGetSubpointMap(dm,  &subpointMap);CHKERRQ(ierr);
542         ierr = DMPlexGetSubpointMap(dmX, &subpointMapX);CHKERRQ(ierr);
543         if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
544           const PetscInt *ind;
545           IS              subpointIS;
546           PetscInt        n, q;
547 
548           ierr = PetscPrintf(PETSC_COMM_SELF, "Making translation PetscSection\n");CHKERRQ(ierr);
549           ierr = PetscSectionGetChart(section, &qStart, &qEnd);CHKERRQ(ierr);
550           ierr = DMPlexCreateSubpointIS(dm, &subpointIS);CHKERRQ(ierr);
551           ierr = ISGetLocalSize(subpointIS, &n);CHKERRQ(ierr);
552           ierr = ISGetIndices(subpointIS, &ind);CHKERRQ(ierr);
553           ierr = PetscSectionCreate(comm, &newSection);CHKERRQ(ierr);
554           ierr = PetscSectionSetChart(newSection, pStart, pEnd);CHKERRQ(ierr);
555           for (q = qStart; q < qEnd; ++q) {
556             PetscInt dof, off, p;
557 
558             ierr = PetscSectionGetDof(section, q, &dof);CHKERRQ(ierr);
559             if (dof) {
560               ierr = PetscFindInt(q, n, ind, &p);CHKERRQ(ierr);
561               if (p >= pStart) {
562                 ierr = PetscSectionSetDof(newSection, p, dof);CHKERRQ(ierr);
563                 ierr = PetscSectionGetOffset(section, q, &off);CHKERRQ(ierr);
564                 ierr = PetscSectionSetOffset(newSection, p, off);CHKERRQ(ierr);
565               }
566             }
567           }
568           ierr = ISRestoreIndices(subpointIS, &ind);CHKERRQ(ierr);
569           ierr = ISDestroy(&subpointIS);CHKERRQ(ierr);
570           /* No need to setup section */
571           section = newSection;
572         }
573       } else {
574         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &section);CHKERRQ(ierr);
575         if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
576       }
577       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
578       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
579       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
580       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
581       if (newSection) {ierr = PetscSectionDestroy(&newSection);CHKERRQ(ierr);}
582     }
583   }
584   /* Cell Fields */
585   ierr = PetscOptionsGetBool(((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);CHKERRQ(ierr);
586   if (hasCell || writePartition) {
587     ierr = PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);CHKERRQ(ierr);
588     for (link = vtk->link; link; link = link->next) {
589       Vec          X = (Vec) link->vec;
590       DM           dmX;
591       PetscSection section, globalSection;
592       const char   *name;
593       PetscInt     enforceDof = PETSC_DETERMINE;
594 
595       if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
596       if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
597       ierr = PetscObjectGetName(link->vec, &name);CHKERRQ(ierr);
598       ierr = VecGetDM(X, &dmX);CHKERRQ(ierr);
599       if (dmX) {
600         ierr = DMGetDefaultSection(dmX, &section);CHKERRQ(ierr);
601       } else {
602         PetscContainer c;
603 
604         ierr = PetscObjectQuery(link->vec, "section", (PetscObject*) &c);CHKERRQ(ierr);
605         if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
606         ierr = PetscContainerGetPointer(c, (void**) &section);CHKERRQ(ierr);
607       }
608       if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
609       ierr = PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);CHKERRQ(ierr);
610       ierr = DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);CHKERRQ(ierr);
611       ierr = PetscSectionDestroy(&globalSection);CHKERRQ(ierr);
612     }
613     if (writePartition) {
614       ierr = PetscFPrintf(comm, fp, "SCALARS partition int 1\n");CHKERRQ(ierr);
615       ierr = PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");CHKERRQ(ierr);
616       ierr = DMPlexVTKWritePartition_ASCII(dm, fp);CHKERRQ(ierr);
617     }
618   }
619   /* Cleanup */
620   ierr = PetscSectionDestroy(&globalCoordSection);CHKERRQ(ierr);
621   ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
622   ierr = PetscFClose(comm, fp);CHKERRQ(ierr);
623   PetscFunctionReturn(0);
624 }
625 
626 #undef __FUNCT__
627 #define __FUNCT__ "DMPlexVTKWriteAll"
628 /*@C
629   DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
630 
631   Collective
632 
633   Input Arguments:
634 + odm - The DMPlex specifying the mesh, passed as a PetscObject
635 - viewer - viewer of type VTK
636 
637   Level: developer
638 
639   Note:
640   This function is a callback used by the VTK viewer to actually write the file.
641   The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
642   Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
643 
644 .seealso: PETSCVIEWERVTK
645 @*/
646 PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
647 {
648   DM             dm = (DM) odm;
649   PetscBool      isvtk;
650   PetscErrorCode ierr;
651 
652   PetscFunctionBegin;
653   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
654   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
655   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);CHKERRQ(ierr);
656   if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
657   switch (viewer->format) {
658   case PETSC_VIEWER_ASCII_VTK:
659     ierr = DMPlexVTKWriteAll_ASCII(dm, viewer);CHKERRQ(ierr);
660     break;
661   case PETSC_VIEWER_VTK_VTU:
662     ierr = DMPlexVTKWriteAll_VTU(dm, viewer);CHKERRQ(ierr);
663     break;
664   default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
665   }
666   PetscFunctionReturn(0);
667 }
668