Actual source code: plexvtk.c
petsc-3.4.5 2014-06-29
1: #define PETSCDM_DLL
2: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
3: #include <../src/sys/classes/viewer/impls/vtk/vtkvimpl.h>
7: PetscErrorCode DMPlexVTKGetCellType(DM dm, PetscInt dim, PetscInt corners, PetscInt *cellType)
8: {
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: return(0);
70: }
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;
88: PetscObjectGetComm((PetscObject)dm,&comm);
89: PetscCommGetNewTag(comm, &tag);
90: MPI_Comm_size(comm, &numProcs);
91: MPI_Comm_rank(comm, &rank);
92: DMPlexGetDimension(dm, &dim);
93: DMPlexGetVTKCellHeight(dm, &cellHeight);
94: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
95: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
96: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
97: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
98: DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
100: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
101: for (c = cStart; c < cEnd; ++c) {
102: PetscInt *closure = NULL;
103: PetscInt closureSize;
105: if (hasLabel) {
106: PetscInt value;
108: DMPlexGetLabelValue(dm, "vtk", c, &value);
109: if (value != 1) continue;
110: }
111: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
112: for (v = 0; v < closureSize*2; v += 2) {
113: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
114: }
115: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
116: ++numCells;
117: }
118: maxCells = numCells;
119: MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);
120: MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);
121: MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);
122: MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);
123: DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
124: ISGetIndices(globalVertexNumbers, &gvertex);
125: PetscMalloc(maxCells * sizeof(PetscInt), &corners);
126: PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);
127: if (!rank) {
128: PetscInt *remoteVertices;
130: for (c = cStart, numCells = 0; c < cEnd; ++c) {
131: PetscInt *closure = NULL;
132: PetscInt closureSize, nC = 0, tmp;
134: if (hasLabel) {
135: PetscInt value;
137: DMPlexGetLabelValue(dm, "vtk", c, &value);
138: if (value != 1) continue;
139: }
140: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
141: for (v = 0; v < closureSize*2; v += 2) {
142: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
143: const PetscInt gv = gvertex[closure[v] - vStart];
144: closure[nC++] = gv < 0 ? -(gv+1) : gv;
145: }
146: }
147: corners[numCells++] = nC;
148: PetscFPrintf(comm, fp, "%d ", nC);
149: tmp = closure[0];
150: closure[0] = closure[1];
151: closure[1] = tmp;
152: for (v = 0; v < nC; ++v) {
153: PetscFPrintf(comm, fp, " %d", closure[v]);
154: }
155: PetscFPrintf(comm, fp, "\n");
156: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
157: }
158: if (numProcs > 1) {PetscMalloc((maxCorners+maxCells) * sizeof(PetscInt), &remoteVertices);}
159: for (proc = 1; proc < numProcs; ++proc) {
160: MPI_Status status;
162: MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);
163: MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);
164: for (c = 0; c < numCorners;) {
165: PetscInt nC = remoteVertices[c++];
167: PetscFPrintf(comm, fp, "%d ", nC);
168: for (v = 0; v < nC; ++v, ++c) {
169: PetscFPrintf(comm, fp, " %d", remoteVertices[c]);
170: }
171: PetscFPrintf(comm, fp, "\n");
172: }
173: }
174: if (numProcs > 1) {PetscFree(remoteVertices);}
175: } else {
176: PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
178: PetscMalloc(numSend * sizeof(PetscInt), &localVertices);
179: for (c = cStart, numCells = 0; c < cEnd; ++c) {
180: PetscInt *closure = NULL;
181: PetscInt closureSize, nC = 0;
183: if (hasLabel) {
184: PetscInt value;
186: DMPlexGetLabelValue(dm, "vtk", c, &value);
187: if (value != 1) continue;
188: }
189: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
190: for (v = 0; v < closureSize*2; v += 2) {
191: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
192: const PetscInt gv = gvertex[closure[v] - vStart];
193: closure[nC++] = gv < 0 ? -(gv+1) : gv;
194: }
195: }
196: corners[numCells++] = nC;
197: localVertices[k++] = nC;
198: for (v = 0; v < nC; ++v, ++k) {
199: localVertices[k] = closure[v];
200: }
201: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
202: }
203: if (k != numSend) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numSend);
204: MPI_Send(&numSend, 1, MPIU_INT, 0, tag, comm);
205: MPI_Send(localVertices, numSend, MPIU_INT, 0, tag, comm);
206: PetscFree(localVertices);
207: }
208: ISRestoreIndices(globalVertexNumbers, &gvertex);
209: PetscFPrintf(comm, fp, "CELL_TYPES %d\n", totCells);
210: if (!rank) {
211: PetscInt cellType;
213: for (c = 0; c < numCells; ++c) {
214: DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);
215: PetscFPrintf(comm, fp, "%d\n", cellType);
216: }
217: for (proc = 1; proc < numProcs; ++proc) {
218: MPI_Status status;
220: MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
221: MPI_Recv(corners, numCells, MPIU_INT, proc, tag, comm, &status);
222: for (c = 0; c < numCells; ++c) {
223: DMPlexVTKGetCellType(dm, dim, corners[c], &cellType);
224: PetscFPrintf(comm, fp, "%d\n", cellType);
225: }
226: }
227: } else {
228: MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
229: MPI_Send(corners, numCells, MPIU_INT, 0, tag, comm);
230: }
231: PetscFree(corners);
232: *totalCells = totCells;
233: return(0);
234: }
238: PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
239: {
240: MPI_Comm comm;
241: const MPI_Datatype mpiType = MPIU_SCALAR;
242: PetscScalar *array;
243: PetscInt numDof = 0, maxDof;
244: PetscInt numLabelCells, cellHeight, cMax, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
245: PetscMPIInt numProcs, rank, proc, tag;
246: PetscBool hasLabel;
247: PetscErrorCode ierr;
250: PetscObjectGetComm((PetscObject)dm,&comm);
253: if (precision < 0) precision = 6;
254: PetscCommGetNewTag(comm, &tag);
255: MPI_Comm_size(comm, &numProcs);
256: MPI_Comm_rank(comm, &rank);
257: PetscSectionGetChart(section, &pStart, &pEnd);
258: /* VTK only wants the values at cells or vertices */
259: DMPlexGetVTKCellHeight(dm, &cellHeight);
260: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
261: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
262: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, &vMax);
263: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
264: if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
265: pStart = PetscMax(PetscMin(cStart, vStart), pStart);
266: pEnd = PetscMin(PetscMax(cEnd, vEnd), pEnd);
267: DMPlexGetStratumSize(dm, "vtk", 1, &numLabelCells);
268: DMPlexGetStratumSize(dm, "vtk", 2, &numLabelVertices);
269: hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
270: for (p = pStart; p < pEnd; ++p) {
271: /* Reject points not either cells or vertices */
272: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
273: if (hasLabel) {
274: PetscInt value;
276: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
277: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
278: DMPlexGetLabelValue(dm, "vtk", p, &value);
279: if (value != 1) continue;
280: }
281: }
282: PetscSectionGetDof(section, p, &numDof);
283: if (numDof) break;
284: }
285: MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
286: enforceDof = PetscMax(enforceDof, maxDof);
287: VecGetArray(v, &array);
288: if (!rank) {
289: char formatString[8];
291: PetscSNPrintf(formatString, 8, "%%.%de", precision);
292: for (p = pStart; p < pEnd; ++p) {
293: /* Here we lose a way to filter points by keeping them out of the Numbering */
294: PetscInt dof, off, goff, d;
296: /* Reject points not either cells or vertices */
297: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
298: if (hasLabel) {
299: PetscInt value;
301: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
302: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
303: DMPlexGetLabelValue(dm, "vtk", p, &value);
304: if (value != 1) continue;
305: }
306: }
307: PetscSectionGetDof(section, p, &dof);
308: PetscSectionGetOffset(section, p, &off);
309: PetscSectionGetOffset(globalSection, p, &goff);
310: if (dof && goff >= 0) {
311: for (d = 0; d < dof; d++) {
312: if (d > 0) {
313: PetscFPrintf(comm, fp, " ");
314: }
315: PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);
316: }
317: for (d = dof; d < enforceDof; d++) {
318: PetscFPrintf(comm, fp, " 0.0");
319: }
320: PetscFPrintf(comm, fp, "\n");
321: }
322: }
323: for (proc = 1; proc < numProcs; ++proc) {
324: PetscScalar *remoteValues;
325: PetscInt size = 0, d;
326: MPI_Status status;
328: MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
329: PetscMalloc(size * sizeof(PetscScalar), &remoteValues);
330: MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
331: for (p = 0; p < size/maxDof; ++p) {
332: for (d = 0; d < maxDof; ++d) {
333: if (d > 0) {
334: PetscFPrintf(comm, fp, " ");
335: }
336: PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);
337: }
338: for (d = maxDof; d < enforceDof; ++d) {
339: PetscFPrintf(comm, fp, " 0.0");
340: }
341: PetscFPrintf(comm, fp, "\n");
342: }
343: PetscFree(remoteValues);
344: }
345: } else {
346: PetscScalar *localValues;
347: PetscInt size, k = 0;
349: PetscSectionGetStorageSize(section, &size);
350: PetscMalloc(size * sizeof(PetscScalar), &localValues);
351: for (p = pStart; p < pEnd; ++p) {
352: PetscInt dof, off, goff, d;
354: /* Reject points not either cells or vertices */
355: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
356: if (hasLabel) {
357: PetscInt value;
359: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
360: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
361: DMPlexGetLabelValue(dm, "vtk", p, &value);
362: if (value != 1) continue;
363: }
364: }
365: PetscSectionGetDof(section, p, &dof);
366: PetscSectionGetOffset(section, p, &off);
367: PetscSectionGetOffset(globalSection, p, &goff);
368: if (goff >= 0) {
369: for (d = 0; d < dof; ++d) {
370: localValues[k++] = array[off+d];
371: }
372: }
373: }
374: MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
375: MPI_Send(localValues, k, mpiType, 0, tag, comm);
376: PetscFree(localValues);
377: }
378: VecRestoreArray(v, &array);
379: return(0);
380: }
384: PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
385: {
386: MPI_Comm comm;
387: PetscInt numDof = 0, maxDof;
388: PetscInt pStart, pEnd, p;
392: PetscObjectGetComm((PetscObject)dm,&comm);
393: PetscSectionGetChart(section, &pStart, &pEnd);
394: for (p = pStart; p < pEnd; ++p) {
395: PetscSectionGetDof(section, p, &numDof);
396: if (numDof) break;
397: }
398: numDof = PetscMax(numDof, enforceDof);
399: MPI_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
400: if (!name) name = "Unknown";
401: if (maxDof == 3) {
402: PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
403: } else {
404: PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);
405: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
406: }
407: DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);
408: return(0);
409: }
413: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
414: {
415: MPI_Comm comm;
416: PetscViewer_VTK *vtk = (PetscViewer_VTK*) viewer->data;
417: FILE *fp;
418: PetscViewerVTKObjectLink link;
419: PetscSection coordSection, globalCoordSection;
420: PetscLayout vLayout;
421: Vec coordinates;
422: PetscReal lengthScale;
423: PetscInt vMax, totVertices, totCells;
424: PetscBool hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE;
425: PetscErrorCode ierr;
428: PetscObjectGetComm((PetscObject)dm,&comm);
429: PetscFOpen(comm, vtk->filename, "wb", &fp);
430: PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
431: PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");
432: PetscFPrintf(comm, fp, "ASCII\n");
433: PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
434: /* Vertices */
435: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
436: DMPlexGetCoordinateSection(dm, &coordSection);
437: PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, &globalCoordSection);
438: DMGetCoordinatesLocal(dm, &coordinates);
439: DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
440: if (vMax >= 0) {
441: PetscInt pStart, pEnd, p, localSize = 0;
443: PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);
444: pEnd = PetscMin(pEnd, vMax);
445: for (p = pStart; p < pEnd; ++p) {
446: PetscInt dof;
448: PetscSectionGetDof(globalCoordSection, p, &dof);
449: if (dof > 0) ++localSize;
450: }
451: PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);
452: PetscLayoutSetLocalSize(vLayout, localSize);
453: PetscLayoutSetBlockSize(vLayout, 1);
454: PetscLayoutSetUp(vLayout);
455: } else {
456: PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);
457: }
458: PetscLayoutGetSize(vLayout, &totVertices);
459: PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);
460: DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);
461: /* Cells */
462: DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
463: /* Vertex fields */
464: for (link = vtk->link; link; link = link->next) {
465: if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
466: if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
467: }
468: if (hasPoint) {
469: PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);
470: for (link = vtk->link; link; link = link->next) {
471: Vec X = (Vec) link->vec;
472: DM dmX;
473: PetscSection section, globalSection, newSection = NULL;
474: const char *name;
475: PetscInt enforceDof = PETSC_DETERMINE;
477: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
478: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
479: PetscObjectGetName(link->vec, &name);
480: VecGetDM(X, &dmX);
481: if (dmX) {
482: DMLabel subpointMap, subpointMapX;
483: PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
485: DMGetDefaultSection(dmX, §ion);
486: /* Here is where we check whether dmX is a submesh of dm */
487: DMPlexGetDimension(dm, &dim);
488: DMPlexGetDimension(dmX, &dimX);
489: DMPlexGetChart(dm, &pStart, &pEnd);
490: DMPlexGetChart(dmX, &qStart, &qEnd);
491: DMPlexGetSubpointMap(dm, &subpointMap);
492: DMPlexGetSubpointMap(dmX, &subpointMapX);
493: if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
494: const PetscInt *ind;
495: IS subpointIS;
496: PetscInt n, q;
498: PetscPrintf(PETSC_COMM_SELF, "Making translation PetscSection\n");
499: PetscSectionGetChart(section, &qStart, &qEnd);
500: DMPlexCreateSubpointIS(dm, &subpointIS);
501: ISGetLocalSize(subpointIS, &n);
502: ISGetIndices(subpointIS, &ind);
503: PetscSectionCreate(comm, &newSection);
504: PetscSectionSetChart(newSection, pStart, pEnd);
505: for (q = qStart; q < qEnd; ++q) {
506: PetscInt dof, off, p;
508: PetscSectionGetDof(section, q, &dof);
509: if (dof) {
510: PetscFindInt(q, n, ind, &p);
511: if (p >= pStart) {
512: PetscSectionSetDof(newSection, p, dof);
513: PetscSectionGetOffset(section, q, &off);
514: PetscSectionSetOffset(newSection, p, off);
515: }
516: }
517: }
518: ISRestoreIndices(subpointIS, &ind);
519: ISDestroy(&subpointIS);
520: /* No need to setup section */
521: section = newSection;
522: }
523: } else {
524: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
525: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
526: }
527: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
528: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);
529: DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
530: PetscSectionDestroy(&globalSection);
531: if (newSection) {PetscSectionDestroy(&newSection);}
532: }
533: }
534: /* Cell Fields */
535: if (hasCell) {
536: PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);
537: for (link = vtk->link; link; link = link->next) {
538: Vec X = (Vec) link->vec;
539: DM dmX;
540: PetscSection section, globalSection;
541: const char *name;
542: PetscInt enforceDof = PETSC_DETERMINE;
544: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
545: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
546: PetscObjectGetName(link->vec, &name);
547: VecGetDM(X, &dmX);
548: if (dmX) {
549: DMGetDefaultSection(dmX, §ion);
550: } else {
551: PetscContainer c;
553: PetscObjectQuery(link->vec, "section", (PetscObject*) &c);
554: if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
555: PetscContainerGetPointer(c, (void**) §ion);
556: }
557: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
558: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, &globalSection);
559: DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
560: PetscSectionDestroy(&globalSection);
561: }
562: }
563: /* Cleanup */
564: PetscSectionDestroy(&globalCoordSection);
565: PetscLayoutDestroy(&vLayout);
566: PetscFClose(comm, fp);
567: return(0);
568: }
572: /*@C
573: DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
575: Collective
577: Input Arguments:
578: + odm - The DMPlex specifying the mesh, passed as a PetscObject
579: - viewer - viewer of type VTK
581: Level: developer
583: Note:
584: This function is a callback used by the VTK viewer to actually write the file.
585: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
586: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
588: .seealso: PETSCVIEWERVTK
589: @*/
590: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
591: {
592: DM dm = (DM) odm;
593: PetscBool isvtk;
599: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
600: if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
601: switch (viewer->format) {
602: case PETSC_VIEWER_ASCII_VTK:
603: DMPlexVTKWriteAll_ASCII(dm, viewer);
604: break;
605: case PETSC_VIEWER_VTK_VTU:
606: DMPlexVTKWriteAll_VTU(dm, viewer);
607: break;
608: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
609: }
610: return(0);
611: }