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