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