Actual source code: plexvtk.c
petsc-3.7.5 2017-01-01
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: DMLabel label;
78: IS globalVertexNumbers = NULL;
79: const PetscInt *gvertex;
80: PetscInt dim;
81: PetscInt numCorners = 0, totCorners = 0, maxCorners, *corners;
82: PetscInt numCells = 0, totCells = 0, maxCells, cellHeight;
83: PetscInt numLabelCells, maxLabelCells, cMax, cStart, cEnd, c, vStart, vEnd, v;
84: PetscMPIInt numProcs, rank, proc, tag;
88: PetscObjectGetComm((PetscObject)dm,&comm);
89: PetscCommGetNewTag(comm, &tag);
90: MPI_Comm_size(comm, &numProcs);
91: MPI_Comm_rank(comm, &rank);
92: DMGetDimension(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: DMGetLabel(dm, "vtk", &label);
99: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
100: MPIU_Allreduce(&numLabelCells, &maxLabelCells, 1, MPIU_INT, MPI_MAX, comm);
101: if (!maxLabelCells) label = NULL;
102: for (c = cStart; c < cEnd; ++c) {
103: PetscInt *closure = NULL;
104: PetscInt closureSize, value;
106: if (label) {
107: DMLabelGetValue(label, c, &value);
108: if (value != 1) continue;
109: }
110: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
111: for (v = 0; v < closureSize*2; v += 2) {
112: if ((closure[v] >= vStart) && (closure[v] < vEnd)) ++numCorners;
113: }
114: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
115: ++numCells;
116: }
117: maxCells = numCells;
118: MPI_Reduce(&numCells, &totCells, 1, MPIU_INT, MPI_SUM, 0, comm);
119: MPI_Reduce(&numCells, &maxCells, 1, MPIU_INT, MPI_MAX, 0, comm);
120: MPI_Reduce(&numCorners, &totCorners, 1, MPIU_INT, MPI_SUM, 0, comm);
121: MPI_Reduce(&numCorners, &maxCorners, 1, MPIU_INT, MPI_MAX, 0, comm);
122: DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
123: ISGetIndices(globalVertexNumbers, &gvertex);
124: PetscMalloc1(maxCells, &corners);
125: PetscFPrintf(comm, fp, "CELLS %d %d\n", totCells, totCorners+totCells);
126: if (!rank) {
127: PetscInt *remoteVertices;
128: int *vertices;
130: PetscMalloc1(maxCorners, &vertices);
131: for (c = cStart, numCells = 0; c < cEnd; ++c) {
132: PetscInt *closure = NULL;
133: PetscInt closureSize, value, nC = 0;
135: if (label) {
136: DMLabelGetValue(label, c, &value);
137: if (value != 1) continue;
138: }
139: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
140: for (v = 0; v < closureSize*2; v += 2) {
141: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
142: const PetscInt gv = gvertex[closure[v] - vStart];
143: vertices[nC++] = gv < 0 ? -(gv+1) : gv;
144: }
145: }
146: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
147: corners[numCells++] = nC;
148: PetscFPrintf(comm, fp, "%d ", nC);
149: DMPlexInvertCell(dim, nC, vertices);
150: for (v = 0; v < nC; ++v) {
151: PetscFPrintf(comm, fp, " %d", vertices[v]);
152: }
153: PetscFPrintf(comm, fp, "\n");
154: }
155: if (numProcs > 1) {PetscMalloc1(maxCorners+maxCells, &remoteVertices);}
156: for (proc = 1; proc < numProcs; ++proc) {
157: MPI_Status status;
159: MPI_Recv(&numCorners, 1, MPIU_INT, proc, tag, comm, &status);
160: MPI_Recv(remoteVertices, numCorners, MPIU_INT, proc, tag, comm, &status);
161: for (c = 0; c < numCorners;) {
162: PetscInt nC = remoteVertices[c++];
164: for (v = 0; v < nC; ++v, ++c) {
165: vertices[v] = remoteVertices[c];
166: }
167: DMPlexInvertCell(dim, nC, vertices);
168: PetscFPrintf(comm, fp, "%d ", nC);
169: for (v = 0; v < nC; ++v) {
170: PetscFPrintf(comm, fp, " %d", vertices[v]);
171: }
172: PetscFPrintf(comm, fp, "\n");
173: }
174: }
175: if (numProcs > 1) {PetscFree(remoteVertices);}
176: PetscFree(vertices);
177: } else {
178: PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
180: PetscMalloc1(numSend, &localVertices);
181: for (c = cStart, numCells = 0; c < cEnd; ++c) {
182: PetscInt *closure = NULL;
183: PetscInt closureSize, value, nC = 0;
185: if (label) {
186: DMLabelGetValue(label, 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 DMPlexVTKWritePartition_ASCII(DM dm, FILE *fp)
239: {
240: MPI_Comm comm;
241: PetscInt numCells = 0, cellHeight;
242: PetscInt numLabelCells, cMax, cStart, cEnd, c;
243: PetscMPIInt numProcs, rank, proc, tag;
244: PetscBool hasLabel;
248: PetscObjectGetComm((PetscObject)dm,&comm);
249: PetscCommGetNewTag(comm, &tag);
250: MPI_Comm_size(comm, &numProcs);
251: MPI_Comm_rank(comm, &rank);
252: DMPlexGetVTKCellHeight(dm, &cellHeight);
253: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
254: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
255: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
256: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
257: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
258: for (c = cStart; c < cEnd; ++c) {
259: if (hasLabel) {
260: PetscInt value;
262: DMGetLabelValue(dm, "vtk", c, &value);
263: if (value != 1) continue;
264: }
265: ++numCells;
266: }
267: if (!rank) {
268: for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", rank);}
269: for (proc = 1; proc < numProcs; ++proc) {
270: MPI_Status status;
272: MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
273: for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", proc);}
274: }
275: } else {
276: MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
277: }
278: return(0);
279: }
283: PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
284: {
285: MPI_Comm comm;
286: const MPI_Datatype mpiType = MPIU_SCALAR;
287: PetscScalar *array;
288: PetscInt numDof = 0, maxDof;
289: PetscInt numLabelCells, cellHeight, cMax, cStart, cEnd, numLabelVertices, vMax, vStart, vEnd, pStart, pEnd, p;
290: PetscMPIInt numProcs, rank, proc, tag;
291: PetscBool hasLabel;
292: PetscErrorCode ierr;
295: PetscObjectGetComm((PetscObject)dm,&comm);
298: if (precision < 0) precision = 6;
299: PetscCommGetNewTag(comm, &tag);
300: MPI_Comm_size(comm, &numProcs);
301: MPI_Comm_rank(comm, &rank);
302: PetscSectionGetChart(section, &pStart, &pEnd);
303: /* VTK only wants the values at cells or vertices */
304: DMPlexGetVTKCellHeight(dm, &cellHeight);
305: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
306: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
307: DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, &vMax);
308: if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
309: if (vMax >= 0) vEnd = PetscMin(vEnd, vMax);
310: pStart = PetscMax(PetscMin(cStart, vStart), pStart);
311: pEnd = PetscMin(PetscMax(cEnd, vEnd), pEnd);
312: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
313: DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);
314: hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
315: for (p = pStart; p < pEnd; ++p) {
316: /* Reject points not either cells or vertices */
317: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
318: if (hasLabel) {
319: PetscInt value;
321: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
322: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
323: DMGetLabelValue(dm, "vtk", p, &value);
324: if (value != 1) continue;
325: }
326: }
327: PetscSectionGetDof(section, p, &numDof);
328: if (numDof) break;
329: }
330: MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
331: enforceDof = PetscMax(enforceDof, maxDof);
332: VecGetArray(v, &array);
333: if (!rank) {
334: char formatString[8];
336: PetscSNPrintf(formatString, 8, "%%.%de", precision);
337: for (p = pStart; p < pEnd; ++p) {
338: /* Here we lose a way to filter points by keeping them out of the Numbering */
339: PetscInt dof, off, goff, d;
341: /* Reject points not either cells or vertices */
342: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
343: if (hasLabel) {
344: PetscInt value;
346: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
347: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
348: DMGetLabelValue(dm, "vtk", p, &value);
349: if (value != 1) continue;
350: }
351: }
352: PetscSectionGetDof(section, p, &dof);
353: PetscSectionGetOffset(section, p, &off);
354: PetscSectionGetOffset(globalSection, p, &goff);
355: if (dof && goff >= 0) {
356: for (d = 0; d < dof; d++) {
357: if (d > 0) {
358: PetscFPrintf(comm, fp, " ");
359: }
360: PetscFPrintf(comm, fp, formatString, PetscRealPart(array[off+d])*scale);
361: }
362: for (d = dof; d < enforceDof; d++) {
363: PetscFPrintf(comm, fp, " 0.0");
364: }
365: PetscFPrintf(comm, fp, "\n");
366: }
367: }
368: for (proc = 1; proc < numProcs; ++proc) {
369: PetscScalar *remoteValues;
370: PetscInt size = 0, d;
371: MPI_Status status;
373: MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
374: PetscMalloc1(size, &remoteValues);
375: MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
376: for (p = 0; p < size/maxDof; ++p) {
377: for (d = 0; d < maxDof; ++d) {
378: if (d > 0) {
379: PetscFPrintf(comm, fp, " ");
380: }
381: PetscFPrintf(comm, fp, formatString, PetscRealPart(remoteValues[p*maxDof+d])*scale);
382: }
383: for (d = maxDof; d < enforceDof; ++d) {
384: PetscFPrintf(comm, fp, " 0.0");
385: }
386: PetscFPrintf(comm, fp, "\n");
387: }
388: PetscFree(remoteValues);
389: }
390: } else {
391: PetscScalar *localValues;
392: PetscInt size, k = 0;
394: PetscSectionGetStorageSize(section, &size);
395: PetscMalloc1(size, &localValues);
396: for (p = pStart; p < pEnd; ++p) {
397: PetscInt dof, off, goff, d;
399: /* Reject points not either cells or vertices */
400: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
401: if (hasLabel) {
402: PetscInt value;
404: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
405: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
406: DMGetLabelValue(dm, "vtk", p, &value);
407: if (value != 1) continue;
408: }
409: }
410: PetscSectionGetDof(section, p, &dof);
411: PetscSectionGetOffset(section, p, &off);
412: PetscSectionGetOffset(globalSection, p, &goff);
413: if (goff >= 0) {
414: for (d = 0; d < dof; ++d) {
415: localValues[k++] = array[off+d];
416: }
417: }
418: }
419: MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
420: MPI_Send(localValues, k, mpiType, 0, tag, comm);
421: PetscFree(localValues);
422: }
423: VecRestoreArray(v, &array);
424: return(0);
425: }
429: PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale)
430: {
431: MPI_Comm comm;
432: PetscInt numDof = 0, maxDof;
433: PetscInt pStart, pEnd, p;
437: PetscObjectGetComm((PetscObject)dm,&comm);
438: PetscSectionGetChart(section, &pStart, &pEnd);
439: for (p = pStart; p < pEnd; ++p) {
440: PetscSectionGetDof(section, p, &numDof);
441: if (numDof) break;
442: }
443: numDof = PetscMax(numDof, enforceDof);
444: MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
445: if (!name) name = "Unknown";
446: if (maxDof == 3) {
447: PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
448: } else {
449: PetscFPrintf(comm, fp, "SCALARS %s double %d\n", name, maxDof);
450: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
451: }
452: DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale);
453: return(0);
454: }
458: static PetscErrorCode DMPlexVTKWriteAll_ASCII(DM dm, PetscViewer viewer)
459: {
460: MPI_Comm comm;
461: PetscViewer_VTK *vtk = (PetscViewer_VTK*) viewer->data;
462: FILE *fp;
463: PetscViewerVTKObjectLink link;
464: PetscSection coordSection, globalCoordSection;
465: PetscLayout vLayout;
466: Vec coordinates;
467: PetscReal lengthScale;
468: PetscInt vMax, totVertices, totCells;
469: PetscBool hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE;
470: PetscErrorCode ierr;
473: PetscObjectGetComm((PetscObject)dm,&comm);
474: PetscFOpen(comm, vtk->filename, "wb", &fp);
475: PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
476: PetscFPrintf(comm, fp, "Simplicial Mesh Example\n");
477: PetscFPrintf(comm, fp, "ASCII\n");
478: PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
479: /* Vertices */
480: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
481: DMGetCoordinateSection(dm, &coordSection);
482: PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
483: DMGetCoordinatesLocal(dm, &coordinates);
484: DMPlexGetHybridBounds(dm, NULL, NULL, NULL, &vMax);
485: if (vMax >= 0) {
486: PetscInt pStart, pEnd, p, localSize = 0;
488: PetscSectionGetChart(globalCoordSection, &pStart, &pEnd);
489: pEnd = PetscMin(pEnd, vMax);
490: for (p = pStart; p < pEnd; ++p) {
491: PetscInt dof;
493: PetscSectionGetDof(globalCoordSection, p, &dof);
494: if (dof > 0) ++localSize;
495: }
496: PetscLayoutCreate(PetscObjectComm((PetscObject)dm), &vLayout);
497: PetscLayoutSetLocalSize(vLayout, localSize);
498: PetscLayoutSetBlockSize(vLayout, 1);
499: PetscLayoutSetUp(vLayout);
500: } else {
501: PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);
502: }
503: PetscLayoutGetSize(vLayout, &totVertices);
504: PetscFPrintf(comm, fp, "POINTS %d double\n", totVertices);
505: DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale);
506: /* Cells */
507: DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
508: /* Vertex fields */
509: for (link = vtk->link; link; link = link->next) {
510: if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
511: if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
512: }
513: if (hasPoint) {
514: PetscFPrintf(comm, fp, "POINT_DATA %d\n", totVertices);
515: for (link = vtk->link; link; link = link->next) {
516: Vec X = (Vec) link->vec;
517: DM dmX;
518: PetscSection section, globalSection, newSection = NULL;
519: const char *name;
520: PetscInt enforceDof = PETSC_DETERMINE;
522: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
523: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
524: PetscObjectGetName(link->vec, &name);
525: VecGetDM(X, &dmX);
526: if (dmX) {
527: DMLabel subpointMap, subpointMapX;
528: PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
530: DMGetDefaultSection(dmX, §ion);
531: /* Here is where we check whether dmX is a submesh of dm */
532: DMGetDimension(dm, &dim);
533: DMGetDimension(dmX, &dimX);
534: DMPlexGetChart(dm, &pStart, &pEnd);
535: DMPlexGetChart(dmX, &qStart, &qEnd);
536: DMPlexGetSubpointMap(dm, &subpointMap);
537: DMPlexGetSubpointMap(dmX, &subpointMapX);
538: if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
539: const PetscInt *ind = NULL;
540: IS subpointIS;
541: PetscInt n = 0, q;
543: PetscSectionGetChart(section, &qStart, &qEnd);
544: DMPlexCreateSubpointIS(dm, &subpointIS);
545: if (subpointIS) {
546: ISGetLocalSize(subpointIS, &n);
547: ISGetIndices(subpointIS, &ind);
548: }
549: PetscSectionCreate(comm, &newSection);
550: PetscSectionSetChart(newSection, pStart, pEnd);
551: for (q = qStart; q < qEnd; ++q) {
552: PetscInt dof, off, p;
554: PetscSectionGetDof(section, q, &dof);
555: if (dof) {
556: PetscFindInt(q, n, ind, &p);
557: if (p >= pStart) {
558: PetscSectionSetDof(newSection, p, dof);
559: PetscSectionGetOffset(section, q, &off);
560: PetscSectionSetOffset(newSection, p, off);
561: }
562: }
563: }
564: if (subpointIS) {
565: ISRestoreIndices(subpointIS, &ind);
566: ISDestroy(&subpointIS);
567: }
568: /* No need to setup section */
569: section = newSection;
570: }
571: } else {
572: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
573: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
574: }
575: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
576: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
577: DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
578: PetscSectionDestroy(&globalSection);
579: if (newSection) {PetscSectionDestroy(&newSection);}
580: }
581: }
582: /* Cell Fields */
583: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);
584: if (hasCell || writePartition) {
585: PetscFPrintf(comm, fp, "CELL_DATA %d\n", totCells);
586: for (link = vtk->link; link; link = link->next) {
587: Vec X = (Vec) link->vec;
588: DM dmX;
589: PetscSection section, globalSection;
590: const char *name;
591: PetscInt enforceDof = PETSC_DETERMINE;
593: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
594: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
595: PetscObjectGetName(link->vec, &name);
596: VecGetDM(X, &dmX);
597: if (dmX) {
598: DMGetDefaultSection(dmX, §ion);
599: } else {
600: PetscContainer c;
602: PetscObjectQuery(link->vec, "section", (PetscObject*) &c);
603: if (!c) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
604: PetscContainerGetPointer(c, (void**) §ion);
605: }
606: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it", name);
607: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
608: DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, name, fp, enforceDof, PETSC_DETERMINE, 1.0);
609: PetscSectionDestroy(&globalSection);
610: }
611: if (writePartition) {
612: PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
613: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
614: DMPlexVTKWritePartition_ASCII(dm, fp);
615: }
616: }
617: /* Cleanup */
618: PetscSectionDestroy(&globalCoordSection);
619: PetscLayoutDestroy(&vLayout);
620: PetscFClose(comm, fp);
621: return(0);
622: }
626: /*@C
627: DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
629: Collective
631: Input Arguments:
632: + odm - The DMPlex specifying the mesh, passed as a PetscObject
633: - viewer - viewer of type VTK
635: Level: developer
637: Note:
638: This function is a callback used by the VTK viewer to actually write the file.
639: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
640: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
642: .seealso: PETSCVIEWERVTK
643: @*/
644: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
645: {
646: DM dm = (DM) odm;
647: PetscBool isvtk;
653: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
654: if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
655: switch (viewer->format) {
656: case PETSC_VIEWER_ASCII_VTK:
657: DMPlexVTKWriteAll_ASCII(dm, viewer);
658: break;
659: case PETSC_VIEWER_VTK_VTU:
660: DMPlexVTKWriteAll_VTU(dm, viewer);
661: break;
662: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
663: }
664: return(0);
665: }