Actual source code: plexvtk.c
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 == 0) {
124: PetscInt *remoteVertices, *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: DMPlexReorderCell(dm, c, vertices);
144: corners[numCells++] = nC;
145: PetscFPrintf(comm, fp, "%D ", nC);
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: PetscFPrintf(comm, fp, "%D ", nC);
164: for (v = 0; v < nC; ++v) {
165: PetscFPrintf(comm, fp, " %D", vertices[v]);
166: }
167: PetscFPrintf(comm, fp, "\n");
168: }
169: }
170: if (size > 1) {PetscFree(remoteVertices);}
171: PetscFree(vertices);
172: } else {
173: PetscInt *localVertices, numSend = numCells+numCorners, k = 0;
175: PetscMalloc1(numSend, &localVertices);
176: for (c = cStart, numCells = 0; c < cEnd; ++c) {
177: PetscInt *closure = NULL;
178: PetscInt closureSize, value, nC = 0;
180: if (label) {
181: DMLabelGetValue(label, c, &value);
182: if (value != 1) continue;
183: }
184: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
185: for (v = 0; v < closureSize*2; v += 2) {
186: if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
187: const PetscInt gv = gvertex[closure[v] - vStart];
188: closure[nC++] = gv < 0 ? -(gv+1) : gv;
189: }
190: }
191: corners[numCells++] = nC;
192: localVertices[k++] = nC;
193: for (v = 0; v < nC; ++v, ++k) {
194: localVertices[k] = closure[v];
195: }
196: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
197: DMPlexReorderCell(dm, c, localVertices+k-nC);
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 == 0) {
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, 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: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
249: hasLabel = numLabelCells > 0 ? PETSC_TRUE : PETSC_FALSE;
250: for (c = cStart; c < cEnd; ++c) {
251: if (hasLabel) {
252: PetscInt value;
254: DMGetLabelValue(dm, "vtk", c, &value);
255: if (value != 1) continue;
256: }
257: ++numCells;
258: }
259: if (rank == 0) {
260: for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", rank);}
261: for (proc = 1; proc < size; ++proc) {
262: MPI_Status status;
264: MPI_Recv(&numCells, 1, MPIU_INT, proc, tag, comm, &status);
265: for (c = 0; c < numCells; ++c) {PetscFPrintf(comm, fp, "%d\n", proc);}
266: }
267: } else {
268: MPI_Send(&numCells, 1, MPIU_INT, 0, tag, comm);
269: }
270: return(0);
271: }
273: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
274: typedef double PetscVTKReal;
275: #elif defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16)
276: typedef float PetscVTKReal;
277: #else
278: typedef PetscReal PetscVTKReal;
279: #endif
281: static PetscErrorCode DMPlexVTKWriteSection_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec v, FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscInt imag)
282: {
283: MPI_Comm comm;
284: const MPI_Datatype mpiType = MPIU_SCALAR;
285: PetscScalar *array;
286: PetscInt numDof = 0, maxDof;
287: PetscInt numLabelCells, cellHeight, cStart, cEnd, numLabelVertices, vStart, vEnd, pStart, pEnd, p;
288: PetscMPIInt size, rank, proc, tag;
289: PetscBool hasLabel;
290: PetscErrorCode ierr;
293: PetscObjectGetComm((PetscObject)dm,&comm);
296: if (precision < 0) precision = 6;
297: PetscCommGetNewTag(comm, &tag);
298: MPI_Comm_size(comm, &size);
299: MPI_Comm_rank(comm, &rank);
300: PetscSectionGetChart(section, &pStart, &pEnd);
301: /* VTK only wants the values at cells or vertices */
302: DMPlexGetVTKCellHeight(dm, &cellHeight);
303: DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
304: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
305: pStart = PetscMax(PetscMin(cStart, vStart), pStart);
306: pEnd = PetscMin(PetscMax(cEnd, vEnd), pEnd);
307: DMGetStratumSize(dm, "vtk", 1, &numLabelCells);
308: DMGetStratumSize(dm, "vtk", 2, &numLabelVertices);
309: hasLabel = numLabelCells > 0 || numLabelVertices > 0 ? PETSC_TRUE : PETSC_FALSE;
310: for (p = pStart; p < pEnd; ++p) {
311: /* Reject points not either cells or vertices */
312: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
313: if (hasLabel) {
314: PetscInt value;
316: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
317: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
318: DMGetLabelValue(dm, "vtk", p, &value);
319: if (value != 1) continue;
320: }
321: }
322: PetscSectionGetDof(section, p, &numDof);
323: if (numDof) break;
324: }
325: MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, comm);
326: enforceDof = PetscMax(enforceDof, maxDof);
327: VecGetArray(v, &array);
328: if (rank == 0) {
329: PetscVTKReal dval;
330: PetscScalar val;
331: char formatString[8];
333: PetscSNPrintf(formatString, 8, "%%.%de", precision);
334: for (p = pStart; p < pEnd; ++p) {
335: /* Here we lose a way to filter points by keeping them out of the Numbering */
336: PetscInt dof, off, goff, d;
338: /* Reject points not either cells or vertices */
339: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
340: if (hasLabel) {
341: PetscInt value;
343: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
344: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
345: DMGetLabelValue(dm, "vtk", p, &value);
346: if (value != 1) continue;
347: }
348: }
349: PetscSectionGetDof(section, p, &dof);
350: PetscSectionGetOffset(section, p, &off);
351: PetscSectionGetOffset(globalSection, p, &goff);
352: if (dof && goff >= 0) {
353: for (d = 0; d < dof; d++) {
354: if (d > 0) {
355: PetscFPrintf(comm, fp, " ");
356: }
357: val = array[off+d];
358: dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
359: PetscFPrintf(comm, fp, formatString, dval);
360: }
361: for (d = dof; d < enforceDof; d++) {
362: PetscFPrintf(comm, fp, " 0.0");
363: }
364: PetscFPrintf(comm, fp, "\n");
365: }
366: }
367: for (proc = 1; proc < size; ++proc) {
368: PetscScalar *remoteValues;
369: PetscInt size = 0, d;
370: MPI_Status status;
372: MPI_Recv(&size, 1, MPIU_INT, proc, tag, comm, &status);
373: PetscMalloc1(size, &remoteValues);
374: MPI_Recv(remoteValues, size, mpiType, proc, tag, comm, &status);
375: for (p = 0; p < size/maxDof; ++p) {
376: for (d = 0; d < maxDof; ++d) {
377: if (d > 0) {
378: PetscFPrintf(comm, fp, " ");
379: }
380: val = remoteValues[p*maxDof+d];
381: dval = (PetscVTKReal) ((imag ? PetscImaginaryPart(val) : PetscRealPart(val)) * scale);
382: PetscFPrintf(comm, fp, formatString, dval);
383: }
384: for (d = maxDof; d < enforceDof; ++d) {
385: PetscFPrintf(comm, fp, " 0.0");
386: }
387: PetscFPrintf(comm, fp, "\n");
388: }
389: PetscFree(remoteValues);
390: }
391: } else {
392: PetscScalar *localValues;
393: PetscInt size, k = 0;
395: PetscSectionGetStorageSize(section, &size);
396: PetscMalloc1(size, &localValues);
397: for (p = pStart; p < pEnd; ++p) {
398: PetscInt dof, off, goff, d;
400: /* Reject points not either cells or vertices */
401: if (((p < cStart) || (p >= cEnd)) && ((p < vStart) || (p >= vEnd))) continue;
402: if (hasLabel) {
403: PetscInt value;
405: if (((p >= cStart) && (p < cEnd) && numLabelCells) ||
406: ((p >= vStart) && (p < vEnd) && numLabelVertices)) {
407: DMGetLabelValue(dm, "vtk", p, &value);
408: if (value != 1) continue;
409: }
410: }
411: PetscSectionGetDof(section, p, &dof);
412: PetscSectionGetOffset(section, p, &off);
413: PetscSectionGetOffset(globalSection, p, &goff);
414: if (goff >= 0) {
415: for (d = 0; d < dof; ++d) {
416: localValues[k++] = array[off+d];
417: }
418: }
419: }
420: MPI_Send(&k, 1, MPIU_INT, 0, tag, comm);
421: MPI_Send(localValues, k, mpiType, 0, tag, comm);
422: PetscFree(localValues);
423: }
424: VecRestoreArray(v, &array);
425: return(0);
426: }
428: static PetscErrorCode DMPlexVTKWriteField_ASCII(DM dm, PetscSection section, PetscSection globalSection, Vec field, const char name[], FILE *fp, PetscInt enforceDof, PetscInt precision, PetscReal scale, PetscBool nameComplex, PetscInt imag)
429: {
430: MPI_Comm comm;
431: PetscInt numDof = 0, maxDof;
432: PetscInt pStart, pEnd, p;
436: PetscObjectGetComm((PetscObject)dm,&comm);
437: PetscSectionGetChart(section, &pStart, &pEnd);
438: for (p = pStart; p < pEnd; ++p) {
439: PetscSectionGetDof(section, p, &numDof);
440: if (numDof) break;
441: }
442: numDof = PetscMax(numDof, enforceDof);
443: MPIU_Allreduce(&numDof, &maxDof, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
444: if (!name) name = "Unknown";
445: if (maxDof == 3) {
446: if (nameComplex) {
447: PetscFPrintf(comm, fp, "VECTORS %s.%s double\n", name, imag ? "Im" : "Re");
448: } else {
449: PetscFPrintf(comm, fp, "VECTORS %s double\n", name);
450: }
451: } else {
452: if (nameComplex) {
453: PetscFPrintf(comm, fp, "SCALARS %s.%s double %D\n", name, imag ? "Im" : "Re", maxDof);
454: } else {
455: PetscFPrintf(comm, fp, "SCALARS %s double %D\n", name, maxDof);
456: }
457: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
458: }
459: DMPlexVTKWriteSection_ASCII(dm, section, globalSection, field, fp, enforceDof, precision, scale, imag);
460: return(0);
461: }
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 totVertices, totCells = 0, loops_per_scalar, l;
474: PetscBool hasPoint = PETSC_FALSE, hasCell = PETSC_FALSE, writePartition = PETSC_FALSE, localized, writeComplex;
475: const char *dmname;
476: PetscErrorCode ierr;
479: #if defined(PETSC_USE_COMPLEX)
480: loops_per_scalar = 2;
481: writeComplex = PETSC_TRUE;
482: #else
483: loops_per_scalar = 1;
484: writeComplex = PETSC_FALSE;
485: #endif
486: DMGetCoordinatesLocalized(dm,&localized);
487: if (localized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"VTK output with localized coordinates not yet supported");
488: PetscObjectGetComm((PetscObject)dm,&comm);
489: PetscFOpen(comm, vtk->filename, "wb", &fp);
490: PetscObjectGetName((PetscObject)dm, &dmname);
491: PetscFPrintf(comm, fp, "# vtk DataFile Version 2.0\n");
492: PetscFPrintf(comm, fp, "%s\n", dmname);
493: PetscFPrintf(comm, fp, "ASCII\n");
494: PetscFPrintf(comm, fp, "DATASET UNSTRUCTURED_GRID\n");
495: /* Vertices */
496: DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);
497: DMGetCoordinateSection(dm, &coordSection);
498: PetscSectionCreateGlobalSection(coordSection, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalCoordSection);
499: DMGetCoordinatesLocal(dm, &coordinates);
500: PetscSectionGetPointLayout(PetscObjectComm((PetscObject)dm), globalCoordSection, &vLayout);
501: PetscLayoutGetSize(vLayout, &totVertices);
502: PetscFPrintf(comm, fp, "POINTS %D double\n", totVertices);
503: DMPlexVTKWriteSection_ASCII(dm, coordSection, globalCoordSection, coordinates, fp, 3, PETSC_DETERMINE, lengthScale, 0);
504: /* Cells */
505: DMPlexVTKWriteCells_ASCII(dm, fp, &totCells);
506: /* Vertex fields */
507: for (link = vtk->link; link; link = link->next) {
508: if ((link->ft == PETSC_VTK_POINT_FIELD) || (link->ft == PETSC_VTK_POINT_VECTOR_FIELD)) hasPoint = PETSC_TRUE;
509: if ((link->ft == PETSC_VTK_CELL_FIELD) || (link->ft == PETSC_VTK_CELL_VECTOR_FIELD)) hasCell = PETSC_TRUE;
510: }
511: if (hasPoint) {
512: PetscFPrintf(comm, fp, "POINT_DATA %D\n", totVertices);
513: for (link = vtk->link; link; link = link->next) {
514: Vec X = (Vec) link->vec;
515: PetscSection section = NULL, globalSection, newSection = NULL;
516: char namebuf[256];
517: const char *name;
518: PetscInt enforceDof = PETSC_DETERMINE;
520: if ((link->ft != PETSC_VTK_POINT_FIELD) && (link->ft != PETSC_VTK_POINT_VECTOR_FIELD)) continue;
521: if (link->ft == PETSC_VTK_POINT_VECTOR_FIELD) enforceDof = 3;
522: PetscObjectGetName(link->vec, &name);
523: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
524: if (!section) {
525: DM dmX;
527: VecGetDM(X, &dmX);
528: if (dmX) {
529: DMLabel subpointMap, subpointMapX;
530: PetscInt dim, dimX, pStart, pEnd, qStart, qEnd;
532: DMGetLocalSection(dmX, §ion);
533: /* Here is where we check whether dmX is a submesh of dm */
534: DMGetDimension(dm, &dim);
535: DMGetDimension(dmX, &dimX);
536: DMPlexGetChart(dm, &pStart, &pEnd);
537: DMPlexGetChart(dmX, &qStart, &qEnd);
538: DMPlexGetSubpointMap(dm, &subpointMap);
539: DMPlexGetSubpointMap(dmX, &subpointMapX);
540: if (((dim != dimX) || ((pEnd-pStart) < (qEnd-qStart))) && subpointMap && !subpointMapX) {
541: const PetscInt *ind = NULL;
542: IS subpointIS;
543: PetscInt n = 0, q;
545: PetscSectionGetChart(section, &qStart, &qEnd);
546: DMPlexGetSubpointIS(dm, &subpointIS);
547: if (subpointIS) {
548: ISGetLocalSize(subpointIS, &n);
549: ISGetIndices(subpointIS, &ind);
550: }
551: PetscSectionCreate(comm, &newSection);
552: PetscSectionSetChart(newSection, pStart, pEnd);
553: for (q = qStart; q < qEnd; ++q) {
554: PetscInt dof, off, p;
556: PetscSectionGetDof(section, q, &dof);
557: if (dof) {
558: PetscFindInt(q, n, ind, &p);
559: if (p >= pStart) {
560: PetscSectionSetDof(newSection, p, dof);
561: PetscSectionGetOffset(section, q, &off);
562: PetscSectionSetOffset(newSection, p, off);
563: }
564: }
565: }
566: if (subpointIS) {
567: ISRestoreIndices(subpointIS, &ind);
568: }
569: /* No need to setup section */
570: section = newSection;
571: }
572: }
573: }
574: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
575: if (link->field >= 0) {
576: const char *fieldname;
578: PetscSectionGetFieldName(section, link->field, &fieldname);
579: PetscSectionGetField(section, link->field, §ion);
580: if (fieldname) {
581: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
582: } else {
583: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);
584: }
585: } else {
586: PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
587: }
588: PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
589: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
590: for (l = 0; l < loops_per_scalar; l++) {
591: DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
592: }
593: PetscSectionDestroy(&globalSection);
594: if (newSection) {PetscSectionDestroy(&newSection);}
595: }
596: }
597: /* Cell Fields */
598: PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_view_partition", &writePartition, NULL);
599: if (hasCell || writePartition) {
600: PetscFPrintf(comm, fp, "CELL_DATA %D\n", totCells);
601: for (link = vtk->link; link; link = link->next) {
602: Vec X = (Vec) link->vec;
603: PetscSection section = NULL, globalSection;
604: const char *name = "";
605: char namebuf[256];
606: PetscInt enforceDof = PETSC_DETERMINE;
608: if ((link->ft != PETSC_VTK_CELL_FIELD) && (link->ft != PETSC_VTK_CELL_VECTOR_FIELD)) continue;
609: if (link->ft == PETSC_VTK_CELL_VECTOR_FIELD) enforceDof = 3;
610: PetscObjectGetName(link->vec, &name);
611: PetscObjectQuery(link->vec, "section", (PetscObject*) §ion);
612: if (!section) {
613: DM dmX;
615: VecGetDM(X, &dmX);
616: if (dmX) {
617: DMGetLocalSection(dmX, §ion);
618: }
619: }
620: if (!section) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Vector %s had no PetscSection composed with it and could not create one from VecGetDM()", name);
621: if (link->field >= 0) {
622: const char *fieldname;
624: PetscSectionGetFieldName(section, link->field, &fieldname);
625: PetscSectionGetField(section, link->field, §ion);
626: if (fieldname) {
627: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%s", name, fieldname);
628: } else {
629: PetscSNPrintf(namebuf, sizeof(namebuf), "%s%D", name, link->field);
630: }
631: } else {
632: PetscSNPrintf(namebuf, sizeof(namebuf), "%s", name);
633: }
634: PetscViewerVTKSanitizeName_Internal(namebuf, sizeof(namebuf));
635: PetscSectionCreateGlobalSection(section, dm->sf, PETSC_FALSE, PETSC_FALSE, &globalSection);
636: for (l = 0; l < loops_per_scalar; l++) {
637: DMPlexVTKWriteField_ASCII(dm, section, globalSection, X, namebuf, fp, enforceDof, PETSC_DETERMINE, 1.0, writeComplex, l);
638: }
639: PetscSectionDestroy(&globalSection);
640: }
641: if (writePartition) {
642: PetscFPrintf(comm, fp, "SCALARS partition int 1\n");
643: PetscFPrintf(comm, fp, "LOOKUP_TABLE default\n");
644: DMPlexVTKWritePartition_ASCII(dm, fp);
645: }
646: }
647: /* Cleanup */
648: PetscSectionDestroy(&globalCoordSection);
649: PetscLayoutDestroy(&vLayout);
650: PetscFClose(comm, fp);
651: return(0);
652: }
654: /*@C
655: DMPlexVTKWriteAll - Write a file containing all the fields that have been provided to the viewer
657: Collective
659: Input Parameters:
660: + odm - The DMPlex specifying the mesh, passed as a PetscObject
661: - viewer - viewer of type VTK
663: Level: developer
665: Note:
666: This function is a callback used by the VTK viewer to actually write the file.
667: The reason for this odd model is that the VTK file format does not provide any way to write one field at a time.
668: Instead, metadata for the entire file needs to be available up-front before you can start writing the file.
670: .seealso: PETSCVIEWERVTK
671: @*/
672: PetscErrorCode DMPlexVTKWriteAll(PetscObject odm, PetscViewer viewer)
673: {
674: DM dm = (DM) odm;
675: PetscBool isvtk;
681: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
682: if (!isvtk) SETERRQ1(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_INCOMP, "Cannot use viewer type %s", ((PetscObject)viewer)->type_name);
683: switch (viewer->format) {
684: case PETSC_VIEWER_ASCII_VTK_DEPRECATED:
685: DMPlexVTKWriteAll_ASCII(dm, viewer);
686: break;
687: case PETSC_VIEWER_VTK_VTU:
688: DMPlexVTKWriteAll_VTU(dm, viewer);
689: break;
690: default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No support for format '%s'", PetscViewerFormats[viewer->format]);
691: }
692: return(0);
693: }